text:
push!(LOAD_PATH, "/Users/thomasriddick/Documents/workspace/Dynamic_HD_Code/Dynamic_HD_Julia_Code/")
using ArgParse
using HDDriverModule: drive_hd_model,drive_hd_and_lake_model
using IOModule: load_river_parameters, load_lake_parameters
using IOModule: load_drainage_fields, load_runoff_fields
using IOModule: load_lake_evaporation_fields
using IOModule: load_lake_initial_values,load_river_initial_values
using GridModule: Grid, LatLonGrid
using FieldModule: Field,LatLonField,repeat,divide
using Profile
using ProfileView
using Serialization
using InteractiveUtils

function pass_arguments()
  settings = ArgParseSettings()
  @add_arg_table settings begin
  "--hd-para-file", "-p"
    help = "Filepath to HD parameters file"
    required = true
  "--hd-init-file", "-i"
    help = "Filepath to HD initialisation file (optional)"
  "--lake-para-file", "-l"
    help = "Filepath to lake model parameters file (optional)"
  "--lake-init-file", "-n"
    help = "Filepath to lake model initialisation file (optional)"
  "--drainage-file", "-d"
    help = "Filepath to file containing drainage values"
  "--runoff-file", "-r"
    help = "Filepath to file containing runoff values"
  "--lake-evaporation-file", "-e"
    help = "Filepath to file containing values for evaporation from lakes"
  "--timesteps", "-t"
    help = "Number of timesteps to run"
    arg_type = Int
    default=100
  end
  return parse_args(settings)
end

function main()
  args = pass_arguments()
  timesteps = args["timesteps"]
  grid = LatLonGrid(360,720,true)
  lake_grid = LatLonGrid(1080,2160,true)
  surface_model_grid = LatLonGrid(48,96,true)
  river_parameters = load_river_parameters(args["hd-para-file"],grid)
  local drainages::Array{Field{Float64},1}
  local runoffs::Array{Field{Float64},1}
  local lake_evaporations::Array{Field{Float64},1}
  if args["drainage-file"] != nothing
    drainages = load_drainage_fields(args["drainage-file"],grid,
                                     last_timestep=12)
    drainages = [divide(x,30.0) for x in drainages]
  else
    drainage::Field{Float64} = LatLonField{Float64}(river_parameters.grid,269747790.0*2.25/100.0)
    drainages = repeat(drainage,12*51)
  end
  if args["runoff-file"] != nothing
    runoffs = load_runoff_fields(args["runoff-file"],grid,
                                 last_timestep=12)
    runoffs = [divide(x,30.0) for x in runoffs]
  else
    runoffs = deepcopy(drainages)
  end
  if args["lake-evaporation-file"] != nothing
    lake_evaporations = load_lake_evaporation_fields(args["lake-evaporation-file"],grid,
                                                     last_timestep=12)
    lake_evaporations = [divide(x,30.0) for x in lake_evaporations]
  else
    lake_evaporation::Field{Float64} = LatLonField{Float64}(river_parameters.grid,0.0)
    lake_evaporations = repeat(lake_evaporation,12*51)
  end
  if args["lake-para-file"] != nothing
    lake_parameters = load_lake_parameters(args["lake-para-file"],lake_grid,grid,surface_model_grid)
    drainages_copy = deepcopy(drainages)
    runoffs_copy = deepcopy(runoffs)
    lake_evaporations_copy = deepcopy(lake_evaporations)
    if args["hd-init-file"] != nothing
      river_fields = load_river_initial_values(args["hd-init-file"],grid,river_parameters)
    end
    if args["lake-init-file"] != nothing
      initial_water_to_lake_centers::LatLonField{Float64},
      initial_spillover_to_rivers::LatLonField{Float64} =
        load_lake_initial_values(args["lake-init-file"],lake_grid,grid)
      if args["hd-init-file"] != nothing
        drive_hd_and_lake_model(river_parameters,river_fields,
                                lake_parameters,drainages,runoffs,lake_evaporations,
                                timesteps,true,initial_water_to_lake_centers,
                                initial_spillover_to_rivers;output_timestep=10,
                                print_timestep_results=false)
        # Profile.clear()
        # Profile.init(delay=0.01)
        # @time drive_hd_and_lake_model(river_parameters,river_fields,
        #                               lake_parameters,drainages_copy,runoffs_copy,
        #                               lake_evaporations_copy,timesteps,
        #                               true,initial_water_to_lake_centers,
        #                               initial_spillover_to_rivers;print_timestep_results=true)
        # Profile.print()
        # r = Profile.retrieve();
        # f = open("/Users/thomasriddick/Downloads/profile.bin", "w")
        # Serialization.serialize(f, r)
        # close(f)
      else
        drive_hd_and_lake_model(river_parameters,lake_parameters,drainages,runoffs,
                                lake_evaporations,timesteps,true,initial_water_to_lake_centers,
                                initial_spillover_to_rivers;output_timestep=10,
                                print_timestep_results=true)
      end
    else
      if args["hd-init-file"] != nothing
        drive_hd_and_lake_model(river_parameters,river_fields,
                                lake_parameters,drainages,runoffs,
                                lake_evaporations,timesteps;output_timestep=10,
                                print_timestep_results=true)
      else
        drive_hd_and_lake_model(river_parameters,lake_parameters,
                                drainages,runoffs,lake_evaporations,
                                timesteps;output_timestep=10,print_timestep_results=true)
        # @time drive_hd_and_lake_model(river_parameters,lake_parameters,
        #                               drainages_copy,runoffs_copy,
        #                               lake_evaporations_copy,timesteps;
        #                               print_timestep_results=true)
      end
    end

  else
    if args["hd-init-file"] != nothing
      river_fields = load_river_initial_values(args["hd-init-file"],grid,river_parameters)
      drive_hd_model(river_parameters,river_fields,
                     drainages,runoffs,lake_evaporations,
                     timesteps;output_timestep=10,print_timestep_results=false)
      @time drive_hd_model(river_parameters,river_fields,
                           drainages_copy,runoffs_copy,lake_evaporations_copy,
                           timesteps;output_timestep=10,print_timestep_results=false)
    else
      drainages_copy = deepcopy(drainages)
      runoffs_copy = deepcopy(runoffs)
      lake_evaporations_copy = deepcopy(lake_evaporations)
      drive_hd_model(river_parameters,drainages,runoffs,
                     lake_evaporations,timesteps;output_timestep=10,
                     print_timestep_results=true)
      Profile.clear()
      Profile.init(delay=0.0001)
      @time drive_hd_model(river_parameters,drainages_copy,
                           runoffs_copy,lake_evaporations_copy,
                           timesteps;output_timestep=10,print_timestep_results=false)
      Profile.print()
      r = Profile.retrieve();
      f = open("/Users/thomasriddick/Downloads/profile.bin", "w")
      Serialization.serialize(f, r)
      close(f)
    end
  end
end

# empty!(ARGS)
#push!(ARGS,"-p/Users/thomasriddick/Documents/data/HDdata/hdfiles/hdpara_file_from_current_model.nc")
# push!(ARGS,"-p/Users/thomasriddick/Documents/data/HDdata/hdfiles/generated/hd_file_prepare_basins_from_glac1D_20190925_225029_1400.nc")
# push!(ARGS,"-l/Users/thomasriddick/Documents/data/HDdata/lakeparafiles/lakeparas_prepare_basins_from_glac1D_20190925_225029_1400.nc")
#push!(ARGS,"-n/Users/thomasriddick/Documents/data/temp/transient_sim_1/results_for_1400/lake_model_start_1400.nc")
#push!(ARGS,"-i/Users/thomasriddick/Documents/data/temp/transient_sim_1/results_for_1400/hdstart_1400.nc")
# push!(ARGS,"-t6000")
push!(ARGS,"-p/Users/thomasriddick/Documents/data/temp/8510/hdpara_8510k.nc")
push!(ARGS,"-l/Users/thomasriddick/Documents/data/temp/8510/lakepara_8510k.nc")
push!(ARGS,"-n/Users/thomasriddick/Documents/data/temp/8510/lakestart_8510k.nc")
push!(ARGS,"-i/Users/thomasriddick/Documents/data/temp/8510/hdrestart_8510k.nc")
#push!(ARGS,"-n/Users/thomasriddick/Documents/data/temp/transient_sim_1/results_for_1400/lake_model_start_1400.nc")
#push!(ARGS,"-i/Users/thomasriddick/Documents/data/temp/transient_sim_1/results_for_1400/hdstart_1400.nc")
push!(ARGS,"-t1000")
main()

max_stars_repo_name: ThomasRiddick/DynamicHD

max_issues_repo_name: ThomasRiddick/DynamicHD

max_forks_repo_name: ThomasRiddick/DynamicHD


text:
using ToyFHE
using OffsetArrays
using Test
using Primes
using GaloisFields
using Random

N = 2^4
q₀ = nextprime(Int128(2)^40 + 1, 1; interval=2N)
ps = nextprime(q₀ + 2N, 1; interval=2N)

ℛ = let CT = ToyFHE.CRTEncoded{2, Tuple{GaloisField.((q₀,ps))...}}
    ζ₂n = GaloisFields.minimal_primitive_root(CT, 2N)
    ToyFHE.NegacyclicRing{CT, N}(ζ₂n)
end

scale = 2^60
Tscale = FixedRational{scale}

plain = CKKSEncoding{Tscale}(zero(ℛ))
plain .= OffsetArray(1:div(N,2), 0:div(N,2)-1)
plain[0] += 1im

re = convert(ToyFHE.NTT.RingElement, plain)
@test CKKSEncoding{Tscale}(ToyFHE.NTT.apply_galois_element(re, 3)) ≈ circshift(plain, -1)

# The same, but encrypted
ckks_params = CKKSParams(ℛ, 1, 3.2)
kp = keygen(ckks_params)

function ToyFHE.NTT.apply_galois_element(c::CipherText{Enc}, galois_element) where {Enc}
    CipherText{Enc}(c.params, map(re->ToyFHE.NTT.apply_galois_element(re, galois_element), c.cs))
end

rt = let c = encrypt(kp, plain)
    c′ = ToyFHE.NTT.apply_galois_element(c, 3)
    ek = ToyFHE.make_eval_key(Random.GLOBAL_RNG, ToyFHE.NTT.apply_galois_element(kp.priv.secret, 3)=>kp.priv)
    c′ = ToyFHE.keyswitch(ek, c′)
    decrypt(kp, c′)
end
@test rt ≈ circshift(plain, -1)

let gk = ToyFHE.keygen(GaloisKey, kp.priv; steps=1)
    @test decrypt(kp, ToyFHE.rotate(gk, encrypt(kp, plain))) ≈ circshift(plain, 1)
end

ext: jl


text:
using Genie, Genie.Router, Genie.WebChannels, Genie.Util, Genie.Configuration
using Revise
using Distributed, Logging

Genie.config.websockets_server = true

const WEBCHANNEL_NAME = "autoreload"
const GENIE_AUTORELOAD = true
const WATCHED_FOLDERS = ["public", pwd()]

function collect_watched_files(folders::Vector{String} = String[])
  result = String[]

  for f in folders
    try
      push!(result, Genie.Util.walk_dir(f, only_extensions = ["jl", "html", "md", "js", "css"])...)
    catch ex
      @error ex
    end
  end

  result
end

function watch()
  @info "Watching $WATCHED_FOLDERS"

  entr(collect_watched_files(WATCHED_FOLDERS)) do
    @info "Reloading!"

    try
      Genie.WebChannels.message(WEBCHANNEL_NAME, "autoreload:full")
    catch ex
      # @error ex
    end
  end
end

if Genie.Configuration.isdev()
  channel("/$WEBCHANNEL_NAME/subscribe") do
    WebChannels.subscribe(@params(:WS_CLIENT), WEBCHANNEL_NAME)
    @show "Subscription OK"
  end

  @spawn watch()
  @spawn WebChannels.unsubscribe_disconnected_clients(WEBCHANNEL_NAME)
end


text:
function setup_etc(backend::GPUBackend, state::SoftlabelSoftmaxLossLayerState, inputs::Vector{Blob})
  state.etc = make_blob(backend, eltype(inputs[1]), size(inputs[1]))
end
function shutdown_etc(backend::GPUBackend, state::SoftlabelSoftmaxLossLayerState)
  destroy(state.etc)
end
function forward(backend::GPUBackend, state::SoftlabelSoftmaxLossLayerState, inputs::Vector{Blob})
  pred = inputs[1]
  label = inputs[2]

  forward(backend, state.softmax, Blob[pred])
  copy!(state.etc, state.softmax.blobs[1])
  prob = state.etc

  dims = size(prob)
  data_type = eltype(prob)

  CuVec.log!(backend, prob)
  loss = -CuBLAS.dot(backend.cublas_ctx, data_type, length(prob), prob.ptr, 1, label.ptr, 1)
  state.loss = state.layer.weight * loss / (prod(dims) / dims[state.op_dim])
end

function backward(backend::GPUBackend, state::SoftlabelSoftmaxLossLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  diff = diffs[1]

  if isa(diff, CuTensorBlob)
    label = inputs[2]
    copy!(diff, state.softmax.blobs[1])
    data_type = eltype(diff)
    dims = size(diff)

    CuBLAS.axpy(backend.cublas_ctx, length(diff), -one(data_type), label.ptr, 1, diff.ptr, 1)
    CuBLAS.scal(backend.cublas_ctx, length(diff), convert(data_type, state.layer.weight * dims[state.op_dim]/prod(dims)), diff.ptr, 1)
  end
end


text:
using Weave
using Test

function pljtest(source, resfile, doctype)
  weave("documents/$source", out_path = "documents/plotsjl/$resfile", doctype=doctype)
  result =  read("documents/plotsjl/$resfile", String)
  ref =  read("documents/plotsjl/$resfile.ref", String)
  @test result == ref
  rm("documents/plotsjl/$resfile")
end

pljtest("plotsjl_test_gr.jmd", "plotsjl_test_gr.md", "pandoc")
pljtest("plotsjl_test_gr.jmd", "plotsjl_test_gr.tex", "tex")

# test cache with plots
isdir("documents/cache") && rm("documents/cache", recursive = true)
weave("documents/plotsjl_test_gr.jmd", cache=:all)
result =  read("documents/plotsjl_test_gr.html", String)
rm("documents/plotsjl_test_gr.html")
weave("documents/plotsjl_test_gr.jmd", cache=:all)
cached_result =  read("documents/plotsjl_test_gr.html", String)
@test result == cached_result


text:
# This file was generated, do not modify it. # hide
eltype(boston2.Chas)


text:
"""
Perform a MetropolisHastings step on a given Markovian state m, using a user
supplied proposal function q. Typically q will use x_t as the proposal mean and
generate x' using some pdf.
N.B. We never constrain variables to be positive/on some interval so this needs
to be done using transforms in the m.target_pdf function. e.g. exp(σ) instead of
σ.
"""
function MetropolisHastings!(m::MarkovChain, q; logs = true)
    # rename everything from `parameters` to x and x' to match notes/maths
    x = m.parameters

    x_prime = rand(q(x), 1)[:, 1] # draw proposal
    # evaluate g(x | x'). N.B. x' in first argument but conditioning
    # on x' still
    g_x = pdf(q(x_prime), x)
    g_x_prime = pdf(q(x), x_prime)

    # Calculate likelihood/probability whatever you want to call it at target
    # and proposal params.
    P_x_prime = m.target_pdf(x_prime, m.data)
    P_x = m.target_pdf(x, m.data)
    # In exp/log space to prevent over/underflow
    # IF LOGS = TRUE THEN m.target_pdf BETTER BE IN LOGS EDWARD
    if logs == true
        A = exp(P_x_prime - P_x + log(g_x) - log(g_x_prime))
    else
        A = (P_x_prime * g_x) / (P_x * g_x_prime)
    end
    # If accept, update params using x'. Otherwise, keep params == x.
    if rand() < A
        m.parameters[:] = x_prime
    end
    return m
end

"""
Call MetropolisHastings! N times with `burn` burn in draws on a given Markov
state with proposal function q.
"""
function MetropolisHastings(N::Int,
                            m::MarkovChain,
                            q;
                            logs = true,
                            burn = Int(floor(N/5)))
    MCMC_draws = Vector{Vector}(undef, N)
    MCMC_diagnostics = Matrix{Float64}(undef, N, 1)
    for i in 1:burn
        MetropolisHastings!(m, q, logs = logs)
    end
    old_params = m.parameters
    for i in 1:N
        new_params = MetropolisHastings!(m, q, logs = logs).parameters[:]
        MCMC_draws[i] = new_params
        MCMC_diagnostics[i, 1] = !isequal(old_params, new_params)
        old_params = new_params
    end
    return MCMC_draws, MCMC_diagnostics
end

"""
Perform a MetropolisHastings step on a given Markovian state m, using a user
supplied proposal function q. Typically q will use x_t as the proposal mean and
generate x' using some pdf.
N.B. We never constrain variables to be positive/on some interval so this needs
to be done using transforms in the m.target_pdf function. e.g. exp(σ) instead of
σ.
"""
function MetropolisHastings!(m::MarkovChainOverSample,
                             q,
                             f;
                             logs = true)

    S = m.S
    # rename everything from `parameters` to x and x' to match notes/maths
    x_single = m.single_sample_parameters
    x_over = m.oversample_parameters
    # No need to recompute this each time

    x_prime_single = rand(q(x_single), 1)[:, 1] # draw proposal
    x_prime_over = rand(f(x_prime_single), S)

    # N.B. g_x doesn't use over sampled param cond density.
    g_x = pdf(q(x_prime_single), x_single)
    g_x_prime = pdf(q(x_single), x_prime_single)

    P_x_prime = m.target_pdf(x_prime_single, x_prime_over, m.data)
    P_x = m.target_pdf(x_single, x_over, m.data)

    # In exp/log space to prevent over/underflow
    # IF LOGS = TRUE THEN m.target_pdf BETTER BE IN LOGS EDWARD
    if logs == true
        A = exp(P_x_prime - P_x + log(g_x) - log(g_x_prime))
    else
        A = (P_x_prime * g_x) / (P_x * g_x_prime)
    end

    if rand() < A
        m.single_sample_parameters[:] = x_prime_single
        m.oversample_parameters[:, :] = x_prime_over
    # If accept, update params using x'. Otherwise, keep params == x.
    end

    return m
end

function MetropolisHastings(N::Int,
                            m::MarkovChainOverSample,
                            q,
                            f;
                            logs = true,
                            burn = Int(floor(N/5)))
    MCMC_draws = Vector{Vector}(undef, N)
    MCMC_diagnostics = Matrix{Float64}(undef, N, 1)
    for i in 1:burn
        MetropolisHastings!(m, q, f, logs = logs)
    end
    old_params = m.single_sample_parameters
    for i in 1:N
        new_params = MetropolisHastings!(m, q, f, logs = logs).single_sample_parameters[:]
        MCMC_draws[i] = new_params
        MCMC_diagnostics[i, 1] = !isequal(old_params, new_params)
        old_params = new_params
    end
    return MCMC_draws, MCMC_diagnostics
end

max_stars_repo_name: EdJeeOnGitHub/Archibald

max_issues_repo_name: EdJeeOnGitHub/Archibald

max_forks_repo_name: EdJeeOnGitHub/Archibald


text:
# forward-mode evaluation of an expression tree
# tree is represented as Vector{NodeData} and adjacency matrix from adjmat()
# assumes values of expressions have already been computed
# partials_storage[k] is the partial derivative of nd[k].parent with respect to the value
# of node k. It's efficient to compute this at the same time as the value of the parent.
# Since we use it in reverse mode and in dual forward mode.
# Note that partials_storage makes a subtle assumption that we have a tree instead of
# a general DAG. If we have a DAG, then need to associate storage with each edge of the DAG.
# user_input_buffer and user_output_buffer are used as temporary storage
# when handling user-defined functions
function forward_eval{T}(storage::Vector{T},partials_storage::Vector{T},nd::Vector{NodeData},adj,const_values,parameter_values,x_values::Vector{T},subexpression_values,user_input_buffer=[],user_output_buffer=[])

    @assert length(storage) >= length(nd)
    @assert length(partials_storage) >= length(nd)

    # nd is already in order such that parents always appear before children
    # so a backwards pass through nd is a forward pass through the tree

    children_arr = rowvals(adj)

    for k in length(nd):-1:1
        # compute the value of node k
        @inbounds nod = nd[k]
        partials_storage[k] = zero(T)
        if nod.nodetype == VARIABLE
            @inbounds storage[k] = x_values[nod.index]
        elseif nod.nodetype == VALUE
            @inbounds storage[k] = const_values[nod.index]
        elseif nod.nodetype == SUBEXPRESSION
            @inbounds storage[k] = subexpression_values[nod.index]
        elseif nod.nodetype == PARAMETER
            @inbounds storage[k] = parameter_values[nod.index]
        elseif nod.nodetype == CALL
            op = nod.index
            @inbounds children_idx = nzrange(adj,k)
            #@show children_idx
            n_children = length(children_idx)
            if op == 1 # :+
                tmp_sum = zero(T)
                # sum over children
                for c_idx in children_idx
                    @inbounds ix = children_arr[c_idx]
                    @inbounds partials_storage[ix] = one(T)
                    @inbounds tmp_sum += storage[ix]
                end
                storage[k] = tmp_sum
            elseif op == 2 # :-
                child1 = first(children_idx)
                @assert n_children == 2
                @inbounds ix1 = children_arr[child1]
                @inbounds ix2 = children_arr[child1+1]
                @inbounds tmp_sub = storage[ix1]
                @inbounds tmp_sub -= storage[ix2]
                @inbounds partials_storage[ix1] = one(T)
                @inbounds partials_storage[ix2] = -one(T)
                storage[k] = tmp_sub
            elseif op == 3 # :*
                tmp_prod = one(T)
                for c_idx in children_idx
                    @inbounds tmp_prod *= storage[children_arr[c_idx]]
                end
                if tmp_prod == zero(T) # inefficient
                    for c_idx in children_idx
                        prod_others = one(T)
                        for c_idx2 in children_idx
                            (c_idx == c_idx2) && continue
                            prod_others *= storage[children_arr[c_idx2]]
                        end
                        partials_storage[children_arr[c_idx]] = prod_others
                    end
                else
                    for c_idx in children_idx
                        ix = children_arr[c_idx]
                        partials_storage[ix] = tmp_prod/storage[ix]
                    end
                end
                @inbounds storage[k] = tmp_prod
            elseif op == 4 # :^
                @assert n_children == 2
                idx1 = first(children_idx)
                idx2 = last(children_idx)
                @inbounds ix1 = children_arr[idx1]
                @inbounds ix2 = children_arr[idx2]
                @inbounds base = storage[ix1]
                @inbounds exponent = storage[ix2]
                if exponent == 2
                    @inbounds storage[k] = base*base
                    @inbounds partials_storage[ix1] = 2*base
                else
                    storage[k] = pow(base,exponent)
                    partials_storage[ix1] = exponent*pow(base,exponent-1)
                end
                partials_storage[ix2] = storage[k]*log(base)
            elseif op == 5 # :/
                @assert n_children == 2
                idx1 = first(children_idx)
                idx2 = last(children_idx)
                @inbounds ix1 = children_arr[idx1]
                @inbounds ix2 = children_arr[idx2]
                @inbounds numerator = storage[ix1]
                @inbounds denominator = storage[ix2]
                recip_denominator = 1/denominator
                @inbounds partials_storage[ix1] = recip_denominator
                partials_storage[ix2] = -numerator*recip_denominator*recip_denominator
                storage[k] = numerator*recip_denominator
            elseif op == 6 # ifelse
                @assert n_children == 3
                idx1 = first(children_idx)
                @inbounds condition = storage[children_arr[idx1]]
                @inbounds lhs = storage[children_arr[idx1+1]]
                @inbounds rhs = storage[children_arr[idx1+2]]
                @inbounds partials_storage[children_arr[idx1+1]] = condition == 1
                @inbounds partials_storage[children_arr[idx1+2]] = !(condition == 1)
                storage[k] = ifelse(condition == 1, lhs, rhs)
            elseif op >= USER_OPERATOR_ID_START
                evaluator = user_operator_map[op]
                f_input = Compat.view(user_input_buffer, 1:n_children)
                grad_output = Compat.view(user_output_buffer, 1:n_children)
                r = 1
                for c_idx in children_idx
                    ix = children_arr[c_idx]
                    f_input[r] = storage[ix]
                    grad_output[r] = 0.0
                    r += 1
                end
                fval = MathProgBase.eval_f(evaluator, f_input)::T
                MathProgBase.eval_grad_f(evaluator, grad_output, f_input)
                storage[k] = fval
                r = 1
                for c_idx in children_idx
                    ix = children_arr[c_idx]
                    partials_storage[ix] = grad_output[r]
                    r += 1
                end
            else
                error("Unsupported operation $(operators[op])")
            end
        elseif nod.nodetype == CALLUNIVAR # univariate function
            op = nod.index
            @inbounds child_idx = children_arr[adj.colptr[k]]
            #@assert child_idx == children_arr[first(nzrange(adj,k))]
            child_val = storage[child_idx]
            if op >= USER_UNIVAR_OPERATOR_ID_START
                f = user_univariate_operator_f[op]
                fprime = user_univariate_operator_fprime[op]
                fval = f(child_val)::T
                fprimeval = fprime(child_val)::T
            else
                fval, fprimeval = eval_univariate(op, child_val)
            end
            @inbounds partials_storage[child_idx] = fprimeval
            @inbounds storage[k] = fval
        elseif nod.nodetype == COMPARISON
            op = nod.index

            @inbounds children_idx = nzrange(adj,k)
            n_children = length(children_idx)
            result = true
            for r in 1:n_children-1
                partials_storage[children_arr[children_idx[r]]] = zero(T)
                cval_lhs = storage[children_arr[children_idx[r]]]
                cval_rhs = storage[children_arr[children_idx[r+1]]]
                if op == 1
                    result &= cval_lhs <= cval_rhs
                elseif op == 2
                    result &= cval_lhs == cval_rhs
                elseif op == 3
                    result &= cval_lhs >= cval_rhs
                elseif op == 4
                    result &= cval_lhs < cval_rhs
                elseif op == 5
                    result &= cval_lhs > cval_rhs
                end
            end
            partials_storage[children_arr[children_idx[n_children]]] = zero(T)
            storage[k] = result
        elseif nod.nodetype == LOGIC
            op = nod.index
            @inbounds children_idx = nzrange(adj,k)
            # boolean values are stored as floats
            partials_storage[children_arr[first(children_idx)]] = zero(T)
            partials_storage[children_arr[last(children_idx)]] = zero(T)
            cval_lhs = (storage[children_arr[first(children_idx)]] == 1)
            cval_rhs = (storage[children_arr[last(children_idx)]] == 1)
            if op == 1
                storage[k] = cval_lhs && cval_rhs
            elseif op == 2
                storage[k] = cval_lhs || cval_rhs
            end
        end

    end
    #@show storage

    return storage[1]

end

export forward_eval

# Evaluate directional derivatives of an expression tree.
# This is equivalent to evaluating the expression tree using DualNumbers or ForwardDiff,
# but instead we keep the storage for the epsilon components separate so that we don't
# need to recompute the real components.
# Computes partials_storage_ϵ as well
# We assume that forward_eval has already been called.
function forward_eval_ϵ{N,T}(storage::Vector{T},storage_ϵ::DenseVector{ForwardDiff.Partials{N,T}},partials_storage::Vector{T},partials_storage_ϵ::DenseVector{ForwardDiff.Partials{N,T}},nd::Vector{NodeData},adj,x_values_ϵ,subexpression_values_ϵ)

    @assert length(storage_ϵ) >= length(nd)
    @assert length(partials_storage_ϵ) >= length(nd)

    zero_ϵ = zero(ForwardDiff.Partials{N,T})

    # nd is already in order such that parents always appear before children
    # so a backwards pass through nd is a forward pass through the tree

    children_arr = rowvals(adj)

    for k in length(nd):-1:1
        # compute the value of node k
        @inbounds nod = nd[k]
        partials_storage_ϵ[k] = zero_ϵ
        if nod.nodetype == VARIABLE
            @inbounds storage_ϵ[k] = x_values_ϵ[nod.index]
        elseif nod.nodetype == VALUE
            @inbounds storage_ϵ[k] = zero_ϵ
        elseif nod.nodetype == SUBEXPRESSION
            @inbounds storage_ϵ[k] = subexpression_values_ϵ[nod.index]
        elseif nod.nodetype == PARAMETER
            @inbounds storage_ϵ[k] = zero_ϵ
        else
            ϵtmp = zero_ϵ
            @inbounds children_idx = nzrange(adj,k)
            for c_idx in children_idx
                ix = children_arr[c_idx]
                @inbounds partial = partials_storage[ix]
                if !isfinite(partials_storage[ix]) && storage_ϵ[ix] == zero_ϵ
                    continue
                end
                ϵtmp += storage_ϵ[ix]*partials_storage[ix]
            end
            storage_ϵ[k] = ϵtmp

            if nod.nodetype == CALL
                op = nod.index
                n_children = length(children_idx)
                if op == 3 # :*
                    # Lazy approach for now
                    anyzero = false
                    tmp_prod = one(ForwardDiff.Dual{N,T})
                    for c_idx in children_idx
                        ix = children_arr[c_idx]
                        sval = storage[ix]
                        gnum = ForwardDiff.Dual(sval,storage_ϵ[ix])
                        tmp_prod *= gnum
                        anyzero = ifelse(sval*sval == zero(T), true, anyzero)
                    end
                    # By a quirk of floating-point numbers, we can have
                    # anyzero == true && ForwardDiff.value(tmp_prod) != zero(T)
                    if anyzero # inefficient
                        for c_idx in children_idx
                            prod_others = one(ForwardDiff.Dual{N,T})
                            for c_idx2 in children_idx
                                (c_idx == c_idx2) && continue
                                ix = children_arr[c_idx2]
                                gnum = ForwardDiff.Dual(storage[ix],storage_ϵ[ix])
                                prod_others *= gnum
                            end
                            partials_storage_ϵ[children_arr[c_idx]] = ForwardDiff.partials(prod_others)
                        end
                    else
                        for c_idx in children_idx
                            ix = children_arr[c_idx]
                            prod_others = tmp_prod/ForwardDiff.Dual(storage[ix],storage_ϵ[ix])
                            partials_storage_ϵ[ix] = ForwardDiff.partials(prod_others)
                        end
                    end
                elseif op == 4 # :^
                    @assert n_children == 2
                    idx1 = first(children_idx)
                    idx2 = last(children_idx)
                    @inbounds ix1 = children_arr[idx1]
                    @inbounds ix2 = children_arr[idx2]
                    @inbounds base = storage[ix1]
                    @inbounds base_ϵ = storage_ϵ[ix1]
                    @inbounds exponent = storage[ix2]
                    @inbounds exponent_ϵ = storage_ϵ[ix2]
                    base_gnum = ForwardDiff.Dual(base,base_ϵ)
                    exponent_gnum = ForwardDiff.Dual(exponent,exponent_ϵ)
                    if exponent == 2
                        partials_storage_ϵ[ix1] = 2*base_ϵ
                    else
                        partials_storage_ϵ[ix1] = ForwardDiff.partials(exponent_gnum*pow(base_gnum,exponent_gnum-1))
                    end
                    result_gnum = ForwardDiff.Dual(storage[k],storage_ϵ[k])
                    partials_storage_ϵ[ix2] = ForwardDiff.partials(result_gnum*log(base_gnum))
                elseif op == 5 # :/
                    @assert n_children == 2
                    idx1 = first(children_idx)
                    idx2 = last(children_idx)
                    @inbounds ix1 = children_arr[idx1]
                    @inbounds ix2 = children_arr[idx2]
                    @inbounds numerator = storage[ix1]
                    @inbounds numerator_ϵ = storage_ϵ[ix1]
                    @inbounds denominator = storage[ix2]
                    @inbounds denominator_ϵ = storage_ϵ[ix2]
                    recip_denominator = 1/ForwardDiff.Dual(denominator,denominator_ϵ)
                    partials_storage_ϵ[ix1] = ForwardDiff.partials(recip_denominator)
                    partials_storage_ϵ[ix2] = ForwardDiff.partials(-ForwardDiff.Dual(numerator,numerator_ϵ)*recip_denominator*recip_denominator)
                elseif op >= USER_OPERATOR_ID_START
                    error("User-defined operators not supported for hessian computations")
                end
            elseif nod.nodetype == CALLUNIVAR # univariate function
                op = nod.index
                @inbounds child_idx = children_arr[adj.colptr[k]]
                child_val = storage[child_idx]
                if op >= USER_UNIVAR_OPERATOR_ID_START
                    fprimeprime = user_univariate_operator_fprimeprime[op](child_val)::T
                else
                    fprimeprime = eval_univariate_2nd_deriv(op, child_val,storage[k])
                end
                partials_storage_ϵ[child_idx] = fprimeprime*storage_ϵ[child_idx]
            end
        end

    end

    return storage_ϵ[1]

end

export forward_eval_ϵ

switchblock = Expr(:block)
for i = 1:length(univariate_operators)
    op = univariate_operators[i]
    deriv_expr = univariate_operator_deriv[i]
	ex = :(return $op(x), $deriv_expr::T)
    push!(switchblock.args,i,ex)
end
switchexpr = Expr(:macrocall, Expr(:.,:Lazy,quot(Symbol("@switch"))), :operator_id,switchblock)

@eval @inline function eval_univariate{T}(operator_id,x::T)
    $switchexpr
end

# TODO: optimize sin/cos/exp
switchblock = Expr(:block)
for i = 1:length(univariate_operators)
    op = univariate_operators[i]
    if op == :asec || op == :acsc || op == :asecd || op == :acscd || op == :acsch || op == :trigamma
        # there's an abs in the derivative that Calculus can't symbolically differentiate
        continue
    end
    if i in 1:3 # :+, :-, :abs
        deriv_expr = :(zero(T))
    elseif op == :sin || op == :cos
        deriv_expr = :(-fval)
    elseif op == :exp
        deriv_expr = :(fval)
    else
        deriv_expr = Calculus.differentiate(univariate_operator_deriv[i],:x)
    end
	ex = :(return $deriv_expr::T)
    push!(switchblock.args,i,ex)
end
switchexpr = Expr(:macrocall, Expr(:.,:Lazy,quot(Symbol("@switch"))), :operator_id,switchblock)

@eval @inline function eval_univariate_2nd_deriv{T}(operator_id,x::T,fval::T)
    $switchexpr
end