text:
import ForwardDiff

using Celeste: Model, SensitiveFloats, DeterministicVI
using Distributions
import DeterministicVI: BvnComponent, GalaxyCacheComponent
import DeterministicVI: eval_bvn_pdf!, get_bvn_derivs!, transform_bvn_derivs!
using DerivativeTestUtils

"""
This is the function of which get_bvn_derivs!() returns the derivatives.
It is only used for testing.
"""
function eval_bvn_log_density{NumType <: Number}(
        elbo_vars::DeterministicVI.ElboIntermediateVariables{NumType},
        bvn::BvnComponent{NumType}, x::Vector{Float64})

    DeterministicVI.eval_bvn_pdf!(elbo_vars.bvn_derivs, bvn, x)

    -0.5 * (
        (x[1] - bvn.the_mean[1]) * elbo_vars.bvn_derivs.py1[1] +
        (x[2] - bvn.the_mean[2]) * elbo_vars.bvn_derivs.py2[1] -
        log(bvn.precision[1, 1] * bvn.precision[2, 2] - bvn.precision[1, 2] ^ 2))
end

function test_process_active_pixels()
    blob, ea, bodies = gen_two_body_dataset()

    # Choose four pixels only to keep the test fast.
    active_pixels = Array(DeterministicVI.ActivePixel, 4)
    active_pixels[1] = DeterministicVI.ActivePixel(1, 1, 10, 11)
    active_pixels[2] = DeterministicVI.ActivePixel(1, 1, 11, 10)
    active_pixels[3] = DeterministicVI.ActivePixel(5, 1, 10, 11)
    active_pixels[4] = DeterministicVI.ActivePixel(5, 1, 11, 10)

    function tile_lik_wrapper_fun{NumType <: Number}(
            ea::ElboArgs{NumType}, calculate_derivs::Bool)

        elbo_vars = DeterministicVI.ElboIntermediateVariables(
            NumType, ea.S,
            length(ea.active_sources),
            calculate_derivs=calculate_derivs,
            calculate_hessian=calculate_derivs)
        DeterministicVI.process_active_pixels!(elbo_vars, ea, active_pixels)
        deepcopy(elbo_vars.elbo)
    end

    function tile_lik_value_wrapper{NumType <: Number}(x::Vector{NumType})
        ea_local = unwrap_vp_vector(x, ea)
        tile_lik_wrapper_fun(ea_local, false).v[1]
    end

    elbo = tile_lik_wrapper_fun(ea, true)

    x = wrap_vp_vector(ea, true)
    test_with_autodiff(tile_lik_value_wrapper, x, elbo)
end

function test_add_log_term()
    blob, ea, bodies = gen_two_body_dataset()

    # Test this pixel
    h, w = (10, 10)

    for b = 1:5
        println("Testing log term for band $b.")
        x_nbm = 70.
        tile = ea.images[b].tiles[1,1]
        tile_source_map = ea.tile_source_map[b][1,1]

        iota = median(blob[b].iota_vec)

        function add_log_term_wrapper_fun{NumType <: Number}(
                ea::ElboArgs{NumType}, calculate_derivs::Bool)

            star_mcs, gal_mcs =
                DeterministicVI.load_bvn_mixtures(ea, b, calculate_derivs=calculate_derivs)
            sbs = DeterministicVI.SourceBrightness{NumType}[
                DeterministicVI.SourceBrightness(ea.vp[s], calculate_derivs=calculate_derivs)
                for s in 1:ea.S]

            elbo_vars_loc = DeterministicVI.ElboIntermediateVariables(NumType, ea.S, ea.S)
            elbo_vars_loc.calculate_derivs = calculate_derivs
            DeterministicVI.populate_fsm_vecs!(
                elbo_vars_loc, ea, tile_source_map, tile, h, w, gal_mcs, star_mcs)
            DeterministicVI.combine_pixel_sources!(
                elbo_vars_loc, ea, tile_source_map, tile, sbs)

            DeterministicVI.add_elbo_log_term!(elbo_vars_loc, x_nbm, iota)

            deepcopy(elbo_vars_loc.elbo)
        end

        function ad_wrapper_fun{NumType <: Number}(x::Vector{NumType})
            ea_local = unwrap_vp_vector(x, ea)
            add_log_term_wrapper_fun(ea_local, false).v[1]
        end

        x = wrap_vp_vector(ea, true)
        elbo = add_log_term_wrapper_fun(ea, true)
        test_with_autodiff(ad_wrapper_fun, x, elbo)
    end
end

function test_combine_pixel_sources()
    blob, ea, bodies = gen_two_body_dataset()

    S = length(ea.active_sources)
    P = length(CanonicalParams)
    h = 10
    w = 10

    for test_var = [false, true], b=1:5
        test_var_string = test_var ? "E_G" : "var_G"
        println("Testing $(test_var_string), band $b")

        tile = ea.images[b].tiles[1,1]; # Note: only one tile in this simulated dataset.
        tile_source_map = ea.tile_source_map[b][1,1]

        function e_g_wrapper_fun{NumType <: Number}(
                ea::ElboArgs{NumType}; calculate_derivs=true)

            star_mcs, gal_mcs =
                DeterministicVI.load_bvn_mixtures(ea, b, calculate_derivs=calculate_derivs)
            sbs = DeterministicVI.SourceBrightness{NumType}[
                DeterministicVI.SourceBrightness(ea.vp[s], calculate_derivs=calculate_derivs)
                for s in 1:ea.S]

            elbo_vars_loc = DeterministicVI.ElboIntermediateVariables(NumType, ea.S, ea.S)
            elbo_vars_loc.calculate_derivs = calculate_derivs
            DeterministicVI.populate_fsm_vecs!(
                elbo_vars_loc, ea, tile_source_map, tile, h, w, gal_mcs, star_mcs)
            DeterministicVI.combine_pixel_sources!(
                elbo_vars_loc, ea, tile_source_map, tile, sbs)
            deepcopy(elbo_vars_loc)
        end

        function wrapper_fun{NumType <: Number}(x::Vector{NumType})
            ea_local = unwrap_vp_vector(x, ea)
            elbo_vars_local = e_g_wrapper_fun(ea_local, calculate_derivs=false)
            test_var ? elbo_vars_local.var_G.v[1] : elbo_vars_local.E_G.v[1]
        end

        x = wrap_vp_vector(ea, true)
        elbo_vars = e_g_wrapper_fun(ea)
        sf = test_var ? deepcopy(elbo_vars.var_G) : deepcopy(elbo_vars.E_G)

        test_with_autodiff(wrapper_fun, x, sf)
    end
end

function test_e_g_s_functions()
    blob, ea, bodies = gen_two_body_dataset()

    # S = length(ea.active_sources)
    P = length(CanonicalParams)
    h = 10
    w = 10
    s = 1

    for test_var = [false, true], b=1:5
        test_var_string = test_var ? "E_G" : "var_G"
        println("Testing $(test_var_string), band $b")

        tile = ea.images[b].tiles[1,1]; # Note: only one tile in this simulated dataset.
        tile_source_map = ea.tile_source_map[b][1,1]

        function e_g_wrapper_fun{NumType <: Number}(
                ea::ElboArgs{NumType}; calculate_derivs=true)

            star_mcs, gal_mcs =
                DeterministicVI.load_bvn_mixtures(ea, b, calculate_derivs=calculate_derivs)
            sbs = DeterministicVI.SourceBrightness{NumType}[
                DeterministicVI.SourceBrightness(ea.vp[s], calculate_derivs=calculate_derivs)
                for s in 1:ea.S]

            elbo_vars_loc = DeterministicVI.ElboIntermediateVariables(
                NumType, ea.S, length(ea.active_sources))
            elbo_vars_loc.calculate_derivs = calculate_derivs
            DeterministicVI.populate_fsm_vecs!(
                elbo_vars_loc, ea, tile_source_map, tile, h, w, gal_mcs, star_mcs)
            DeterministicVI.accumulate_source_brightness!(elbo_vars_loc, ea, sbs, s, b)
            deepcopy(elbo_vars_loc)
        end

        function wrapper_fun{NumType <: Number}(x::Vector{NumType})
            @assert length(x) == P
            ea_local = forward_diff_model_params(NumType, ea)
            ea_local.vp[s] = x
            elbo_vars_local = e_g_wrapper_fun(ea_local, calculate_derivs=false)
            test_var ? elbo_vars_local.var_G_s.v[1] : elbo_vars_local.E_G_s.v[1]
        end

        x = ea.vp[s]

        elbo_vars = e_g_wrapper_fun(ea)

        # Sanity check the variance value.
        @test_approx_eq(elbo_vars.var_G_s.v,
                                        elbo_vars.E_G2_s.v[1] - (elbo_vars.E_G_s.v[1] ^ 2))

        sf = test_var ? deepcopy(elbo_vars.var_G_s) : deepcopy(elbo_vars.E_G_s)

        test_with_autodiff(wrapper_fun, x, sf)
    end
end

function test_fs1m_derivatives()
    # TODO: test with a real and asymmetric wcs jacobian.
    blob, ea, three_bodies = gen_three_body_dataset()
    omitted_ids = Int[]
    kept_ids = setdiff(1:length(ids), omitted_ids)

    s = 1
    b = 1

    patch = ea.patches[s, b]
    u = ea.vp[s][ids.u]
    u_pix = Model.linear_world_to_pix(
        patch.wcs_jacobian, patch.center, patch.pixel_center, u)
    x = ceil(u_pix + [1.0, 2.0])

    elbo_vars = DeterministicVI.ElboIntermediateVariables(Float64, 1, 1)

    ###########################
    # Galaxies

    # Pick out a single galaxy component for testing.
    # The index is (psf, galaxy, gal type, source)
    for psf_k=1:psf_K, type_i = 1:2, gal_j in 1:[8,6][type_i]
        gcc_ind = (psf_k, gal_j, type_i, s)
        function f_wrap_gal{T <: Number}(par::Vector{T})
            # This uses ea, x, wcs_jacobian, and gcc_ind from the enclosing namespace.
            ea_fd = forward_diff_model_params(T, ea)
            elbo_vars_fd = DeterministicVI.ElboIntermediateVariables(T, 1, 1)

            # Make sure par is as long as the galaxy parameters.
            @assert length(par) == length(shape_standard_alignment[2])
            for p1 in 1:length(par)
                    p0 = shape_standard_alignment[2][p1]
                    ea_fd.vp[s][p0] = par[p1]
            end
            star_mcs, gal_mcs =
                DeterministicVI.load_bvn_mixtures(ea_fd, b, calculate_derivs=false)

            # Raw:
            gcc = gal_mcs[gcc_ind...]
            eval_bvn_pdf!(elbo_vars_fd.bvn_derivs, gcc.bmc, x)
            elbo_vars_fd.bvn_derivs.f_pre[1] * gcc.e_dev_i
        end

        function ea_to_par_gal(ea::ElboArgs{Float64})
            par = zeros(length(shape_standard_alignment[2]))
            for p1 in 1:length(par)
                    p0 = shape_standard_alignment[2][p1]
                    par[p1] = ea.vp[s][p0]
            end
            par
        end

        par_gal = ea_to_par_gal(ea)

        star_mcs, gal_mcs = DeterministicVI.load_bvn_mixtures(ea, b)
        clear!(elbo_vars.fs1m_vec[s])
        DeterministicVI.accum_galaxy_pos!(
            elbo_vars, s, gal_mcs[gcc_ind...], x, patch.wcs_jacobian, true)
        fs1m = deepcopy(elbo_vars.fs1m_vec[s])

        # Two sanity checks.
        gcc = gal_mcs[gcc_ind...]
        clear!(elbo_vars.fs1m_vec[s])
        v = eval_bvn_log_density(elbo_vars, gcc.bmc, x)
        gc = galaxy_prototypes[gcc_ind[3]][gcc_ind[2]]
        pc = ea.patches[s, b].psf[gcc_ind[1]]

        @test_approx_eq(
            pc.alphaBar * gc.etaBar * gcc.e_dev_i * exp(v) / (2 * pi),
            fs1m.v)

        test_with_autodiff(f_wrap_gal, par_gal, fs1m)
    end
end

function test_fs0m_derivatives()
    # TODO: test with a real and asymmetric wcs jacobian.
    blob, ea, three_bodies = gen_three_body_dataset()
    omitted_ids = Int[]
    kept_ids = setdiff(1:length(ids), omitted_ids)

    s = 1
    b = 1

    patch = ea.patches[s, b]
    u = ea.vp[s][ids.u]
    u_pix = Model.linear_world_to_pix(
        patch.wcs_jacobian, patch.center, patch.pixel_center, u)
    x = ceil(u_pix + [1.0, 2.0])

    elbo_vars = DeterministicVI.ElboIntermediateVariables(Float64, 1, 1)

    ###########################
    # Stars

    # Pick out a single star component for testing.
    # The index is psf, source
    for psf_k=1:psf_K
        bmc_ind = (psf_k, s)
        function f_wrap_star{T <: Number}(par::Vector{T})
            # This uses ea, x, wcs_jacobian, and gcc_ind from the enclosing namespace.
            ea_fd = forward_diff_model_params(T, ea)

            # Make sure par is as long as the galaxy parameters.
            @assert length(par) == length(ids.u)
            for p1 in 1:2
                    p0 = ids.u[p1]
                    ea_fd.vp[s][p0] = par[p1]
            end
            star_mcs, gal_mcs = DeterministicVI.load_bvn_mixtures(ea_fd, b)
            elbo_vars_fd = DeterministicVI.ElboIntermediateVariables(T, 1, 1)
            DeterministicVI.accum_star_pos!(
                elbo_vars_fd, s, star_mcs[bmc_ind...], x, patch.wcs_jacobian, true)
            elbo_vars_fd.fs0m_vec[s].v[1]
        end

        function ea_to_par_star(ea::ElboArgs{Float64})
            par = zeros(2)
            for p1 in 1:length(par)
                    par[p1] = ea.vp[s][ids.u[p1]]
            end
            par
        end

        par_star = ea_to_par_star(ea)

        clear!(elbo_vars.fs0m_vec[s])
        star_mcs, gal_mcs = DeterministicVI.load_bvn_mixtures(ea, b)
        DeterministicVI.accum_star_pos!(
            elbo_vars, s, star_mcs[bmc_ind...], x, patch.wcs_jacobian, true)
        fs0m = deepcopy(elbo_vars.fs0m_vec[s])

        test_with_autodiff(f_wrap_star, par_star, fs0m)
    end
end

function test_bvn_derivatives()
    # Test log(bvn prob) / d(mean, sigma)

    x = Float64[2.0, 3.0]

    e_angle, e_axis, e_scale = (1.1, 0.02, 4.8)
    sigma = DeterministicVI.get_bvn_cov(e_axis, e_angle, e_scale)

    offset = Float64[0.5, 0.25]

    # Note that get_bvn_derivs doesn't use the weight, so set it to something
    # strange to check that it doesn't matter.
    weight = 0.724

    bvn = BvnComponent{Float64}(offset, sigma, weight)
    elbo_vars = DeterministicVI.ElboIntermediateVariables(Float64, 1, 1)
    DeterministicVI.eval_bvn_pdf!(elbo_vars.bvn_derivs, bvn, x)
    DeterministicVI.get_bvn_derivs!(elbo_vars.bvn_derivs, bvn, true, true)

    function bvn_function{T <: Number}(x::Vector{T}, sigma::Matrix{T})
        local_x = offset - x
        -0.5 * ((local_x' * (sigma \ local_x))[1,1] + log(det(sigma)))
    end

    x_ids = 1:2
    sig_ids = 3:5
    function wrap(x::Vector{Float64}, sigma::Matrix{Float64})
        par = zeros(Float64, length(x_ids) + length(sig_ids))
        par[x_ids] = x
        par[sig_ids] = [ sigma[1, 1], sigma[1, 2], sigma[2, 2]]
        par
    end

    function f_wrap{T <: Number}(par::Vector{T})
        x_loc = par[x_ids]
        s_vec = par[sig_ids]
        sig_loc = T[s_vec[1] s_vec[2]; s_vec[2] s_vec[3]]
        bvn_function(x_loc, sig_loc)
    end

    par = wrap(x, sigma)

    # Sanity check
    @test_approx_eq eval_bvn_log_density(elbo_vars, bvn, x) f_wrap(par)

    bvn_derivs = elbo_vars.bvn_derivs
    ad_grad = ForwardDiff.gradient(f_wrap, par)
    @test_approx_eq bvn_derivs.bvn_x_d ad_grad[x_ids]
    @test_approx_eq bvn_derivs.bvn_sig_d ad_grad[sig_ids]

    ad_hess = ForwardDiff.hessian(f_wrap, par)
    @test_approx_eq bvn_derivs.bvn_xx_h ad_hess[x_ids, x_ids]
    @test_approx_eq bvn_derivs.bvn_xsig_h ad_hess[x_ids, sig_ids]
    @test_approx_eq bvn_derivs.bvn_sigsig_h ad_hess[sig_ids, sig_ids]
end

function test_galaxy_variable_transform()
    # This is testing transform_bvn_derivs!

    # TODO: test with a real and asymmetric wcs jacobian.
    # We only need this for a psf and jacobian.
    blob, ea, three_bodies = gen_three_body_dataset()

    # Pick a single source and band for testing.
    s = 1
    b = 5

    # The pixel and world centers shouldn't matter for derivatives.
    patch = ea.patches[s, b]
    psf = patch.psf[1]

    # Pick out a single galaxy component for testing.
    gp = galaxy_prototypes[2][4]
    e_dev_dir = -1.0
    e_dev_i = 0.85

    # Test the variable transformation.
    e_angle, e_axis, e_scale = (1.1, 0.02, 4.8)

    u = Float64[5.3, 2.9]
    x = Float64[7.0, 5.0]

    # The indices in par of each variable.
    par_ids_u = [1, 2]
    par_ids_e_axis = 3
    par_ids_e_angle = 4
    par_ids_e_scale = 5
    par_ids_length = 5

    function wrap_par{T <: Number}(
            u::Vector{T}, e_angle::T, e_axis::T, e_scale::T)
        par = zeros(T, par_ids_length)
        par[par_ids_u] = u
        par[par_ids_e_angle] = e_angle
        par[par_ids_e_axis] = e_axis
        par[par_ids_e_scale] = e_scale
        par
    end

    function f_bvn_wrap{T <: Number}(par::Vector{T})
        u = par[par_ids_u]
        e_angle = par[par_ids_e_angle]
        e_axis = par[par_ids_e_axis]
        e_scale = par[par_ids_e_scale]
        u_pix = Model.linear_world_to_pix(
            patch.wcs_jacobian, patch.center, patch.pixel_center, u)

        sigma = DeterministicVI.get_bvn_cov(e_axis, e_angle, e_scale)

        function bvn_function{T <: Number}(u_pix::Vector{T}, sigma::Matrix{T})
            local_x = x - u_pix
            -0.5 * ((local_x' * (sigma \ local_x))[1,1] + log(det(sigma)))
        end

        bvn_function(u_pix, sigma)
    end

    # First just test the bvn function itself
    par = wrap_par(u, e_angle, e_axis, e_scale)
    u_pix = Model.linear_world_to_pix(
        patch.wcs_jacobian, patch.center, patch.pixel_center, u)
    sigma = DeterministicVI.get_bvn_cov(e_axis, e_angle, e_scale)
    bmc = BvnComponent{Float64}(u_pix, sigma, 1.0)
    sig_sf = DeterministicVI.GalaxySigmaDerivs(e_angle, e_axis, e_scale, sigma)
    gcc = GalaxyCacheComponent(1.0, 1.0, bmc, sig_sf)
    elbo_vars = DeterministicVI.ElboIntermediateVariables(Float64, 1, 1)
    DeterministicVI.eval_bvn_pdf!(elbo_vars.bvn_derivs, bmc, x)
    DeterministicVI.get_bvn_derivs!(elbo_vars.bvn_derivs, bmc, true, true)
    DeterministicVI.transform_bvn_derivs!(
        elbo_vars.bvn_derivs, gcc.sig_sf, patch.wcs_jacobian, true)

    f_bvn_wrap(par)

    # Check the gradient.
    ad_grad = ForwardDiff.gradient(f_bvn_wrap, par)
    bvn_derivs = elbo_vars.bvn_derivs
    @test_approx_eq ad_grad [bvn_derivs.bvn_u_d; bvn_derivs.bvn_s_d]

    ad_hess = ForwardDiff.hessian(f_bvn_wrap, par)
    @test_approx_eq ad_hess[1:2, 1:2] bvn_derivs.bvn_uu_h
    @test_approx_eq ad_hess[1:2, 3:5] bvn_derivs.bvn_us_h

    celeste_bvn_ss_h = deepcopy(bvn_derivs.bvn_ss_h)
    ad_bvn_ss_h = deepcopy(ad_hess[3:5, 3:5])
    @test_approx_eq ad_hess[3:5, 3:5] bvn_derivs.bvn_ss_h
end

function test_galaxy_cache_component()
    # TODO: eliminate some of the redundancy in these tests.

    # TODO: test with a real and asymmetric wcs jacobian.
    # We only need this for a psf and jacobian.
    blob, ea, three_bodies = gen_three_body_dataset()

    # Pick a single source and band for testing.
    s = 1
    b = 5

    # The pixel and world centers shouldn't matter for derivatives.
    patch = ea.patches[s, b]
    psf = patch.psf[1]

    # Pick out a single galaxy component for testing.
    gp = galaxy_prototypes[2][4]
    e_dev_dir = -1.0
    e_dev_i = 0.85

    # Test the variable transformation.
    e_angle, e_axis, e_scale = (1.1, 0.02, 4.8)

    u = Float64[5.3, 2.9]
    x = Float64[7.0, 5.0]

    # The indices in par of each variable.
    par_ids_u = [1, 2]
    par_ids_e_axis = 3
    par_ids_e_angle = 4
    par_ids_e_scale = 5
    par_ids_length = 5

    function f_wrap{T <: Number}(par::Vector{T})
        u = par[par_ids_u]
        e_angle = par[par_ids_e_angle]
        e_axis = par[par_ids_e_axis]
        e_scale = par[par_ids_e_scale]
        u_pix = Model.linear_world_to_pix(
            patch.wcs_jacobian, patch.center, patch.pixel_center, u)
        elbo_vars_fd = DeterministicVI.ElboIntermediateVariables(T, 1, 1)
        e_dev_i_fd = convert(T, e_dev_i)
        gcc = GalaxyCacheComponent(
                        e_dev_dir, e_dev_i_fd, gp, psf,
                        u_pix, e_axis, e_angle, e_scale, false, false)

        eval_bvn_pdf!(elbo_vars_fd.bvn_derivs, gcc.bmc, x)

        log(elbo_vars_fd.bvn_derivs.f_pre[1])
    end

    function wrap_par{T <: Number}(
            u::Vector{T}, e_angle::T, e_axis::T, e_scale::T)
        par = zeros(T, par_ids_length)
        par[par_ids_u] = u
        par[par_ids_e_angle] = e_angle
        par[par_ids_e_axis] = e_axis
        par[par_ids_e_scale] = e_scale
        par
    end

    par = wrap_par(u, e_angle, e_axis, e_scale)
    u_pix = Model.linear_world_to_pix(
        patch.wcs_jacobian, patch.center, patch.pixel_center, u)
    gcc = GalaxyCacheComponent(
                    e_dev_dir, e_dev_i, gp, psf,
                    u_pix, e_axis, e_angle, e_scale, true, true)
    elbo_vars = DeterministicVI.ElboIntermediateVariables(Float64, 1, 1)
    eval_bvn_pdf!(elbo_vars.bvn_derivs, gcc.bmc, x)
    get_bvn_derivs!(elbo_vars.bvn_derivs, gcc.bmc, true, true)
    transform_bvn_derivs!(elbo_vars.bvn_derivs, gcc.sig_sf, patch.wcs_jacobian, true)

    # Sanity check the wrapper.
    @test_approx_eq(
        -0.5 *((x - gcc.bmc.the_mean)' * gcc.bmc.precision * (x - gcc.bmc.the_mean) -
                     log(det(gcc.bmc.precision)))[1,1] - log(2pi) +
                     log(psf.alphaBar * gp.etaBar),
        f_wrap(par))

    # Check the gradient.
    ad_grad_fun = x -> ForwardDiff.gradient(f_wrap, x)
    ad_grad = ad_grad_fun(par)
    bvn_derivs = elbo_vars.bvn_derivs
    @test_approx_eq ad_grad [bvn_derivs.bvn_u_d; bvn_derivs.bvn_s_d]

    ad_hess_fun = x -> ForwardDiff.hessian(f_wrap, x)
    ad_hess = ad_hess_fun(par)

    @test_approx_eq ad_hess[1:2, 1:2] bvn_derivs.bvn_uu_h
    @test_approx_eq ad_hess[1:2, 3:5] bvn_derivs.bvn_us_h

    # I'm not sure why this requires less precision for this test.
    celeste_bvn_ss_h = deepcopy(bvn_derivs.bvn_ss_h)
    ad_bvn_ss_h = deepcopy(ad_hess[3:5, 3:5])
    @test_approx_eq ad_hess[3:5, 3:5] bvn_derivs.bvn_ss_h

end

function test_galaxy_sigma_derivs()
    # Test d sigma / d shape

    e_angle, e_axis, e_scale = (pi / 4, 0.7, 1.2)

    function wrap_par{T <: Number}(e_angle::T, e_axis::T, e_scale::T)
        par = zeros(T, length(gal_shape_ids))
        par[gal_shape_ids.e_angle] = e_angle
        par[gal_shape_ids.e_axis] = e_axis
        par[gal_shape_ids.e_scale] = e_scale
        par
    end

    for si in 1:3
        sig_i = [(1, 1), (1, 2), (2, 2)][si]
        println("Testing sigma[$(sig_i)]")
        function f_wrap{T <: Number}(par::Vector{T})
            e_angle_fd = par[gal_shape_ids.e_angle]
            e_axis_fd = par[gal_shape_ids.e_axis]
            e_scale_fd = par[gal_shape_ids.e_scale]
            this_cov = DeterministicVI.get_bvn_cov(e_axis_fd, e_angle_fd, e_scale_fd)
            this_cov[sig_i...]
        end

        par = wrap_par(e_angle, e_axis, e_scale)
        XiXi = DeterministicVI.get_bvn_cov(e_axis, e_angle, e_scale)

        gal_derivs = DeterministicVI.GalaxySigmaDerivs(e_angle, e_axis, e_scale, XiXi)

        ad_grad_fun = x -> ForwardDiff.gradient(f_wrap, x)
        ad_grad = ad_grad_fun(par)
        @test_approx_eq gal_derivs.j[si, :][:] ad_grad

        ad_hess_fun = x -> ForwardDiff.hessian(f_wrap, x)
        ad_hess = ad_hess_fun(par)
        @test_approx_eq(
            ad_hess,
            reshape(gal_derivs.t[si, :, :],
                            length(gal_shape_ids), length(gal_shape_ids)))
    end
end

function test_brightness_hessian()
    blob, ea, star_cat = gen_sample_star_dataset()
    kept_ids = [ ids.r1; ids.r2; ids.c1[:]; ids.c2[:] ]
    omitted_ids = setdiff(1:length(ids), kept_ids)
    i = 1

    for squares in [false, true], b in 1:5, i in 1:2
        squares_string = squares ? "E_G" : "E_G2"
        println("Testing brightness $(squares_string) for band $b, type $i")
        function wrap_source_brightness{NumType <: Number}(
                vp::Vector{NumType}, calculate_derivs::Bool)

            sb = DeterministicVI.SourceBrightness(vp, calculate_derivs=calculate_derivs)
            if squares
                return deepcopy(sb.E_ll_a[b, i])
            else
                return deepcopy(sb.E_l_a[b, i])
            end
        end

        function wrap_source_brightness_value{NumType <: Number}(
                bright_vp::Vector{NumType})
            vp = zeros(NumType, length(CanonicalParams))
            for b_i in 1:length(brightness_standard_alignment[i])
                vp[brightness_standard_alignment[i][b_i]] = bright_vp[b_i]
            end
            wrap_source_brightness(vp, false).v[1]
        end

        bright_vp = ea.vp[1][brightness_standard_alignment[i]]
        bright = wrap_source_brightness(ea.vp[1], true)

        @test_approx_eq bright.v[1] wrap_source_brightness_value(bright_vp)

        ad_grad = ForwardDiff.gradient(wrap_source_brightness_value, bright_vp)
        @test_approx_eq ad_grad bright.d[:, 1]

        ad_hess = ForwardDiff.hessian(wrap_source_brightness_value, bright_vp)
        @test_approx_eq ad_hess bright.h
    end
end

function test_dsiginv_dsig()
    e_angle, e_axis, e_scale = (1.1, 0.02, 4.8) # bvn_derivs.bvn_sigsig_h is large
    the_cov = DeterministicVI.get_bvn_cov(e_axis, e_angle, e_scale)
    the_mean = Float64[0., 0.]
    bvn = BvnComponent{Float64}(the_mean, the_cov, 1.0)
    sigma_vec = Float64[ the_cov[1, 1], the_cov[1, 2], the_cov[2, 2] ]

    for component_index = 1:3
        components = [(1, 1), (1, 2), (2, 2)]
        function invert_sigma{NumType <: Number}(sigma_vec::Vector{NumType})
            sigma_loc = NumType[sigma_vec[1] sigma_vec[2]; sigma_vec[2] sigma_vec[3]]
            sigma_inv = inv(sigma_loc)
            sigma_inv[components[component_index]...]
        end

        ad_grad = ForwardDiff.gradient(invert_sigma, sigma_vec)
        @test_approx_eq ad_grad bvn.dsiginv_dsig[component_index, :][:]
    end
end

########################
# Test derivatives of KL divergence functions

"""
Use Monte Carlo to check whether KL(q_dist || p_dist) matches exact_kl

Args:
    q_dist, p_dist: Distribution objects
    exact_kl: The expected exact KL
"""
function test_kl_value(q_dist, p_dist, exact_kl::Float64)
    sample_size = 2_000_000

    q_samples = rand(q_dist, sample_size)
    empirical_kl_samples =
      logpdf(q_dist, q_samples) - logpdf(p_dist, q_samples)
    empirical_kl = mean(empirical_kl_samples)
    tol = 4 * std(empirical_kl_samples) / sqrt(sample_size)
    min_diff = 1e-2 * std(empirical_kl_samples) / sqrt(sample_size)

    @test_approx_eq_eps empirical_kl exact_kl tol
end

function test_beta_kl_derivatives()
    alpha2 = 3.5
    beta2 = 4.3

    alpha1 = 4.1
    beta1 = 3.9

    par = Float64[ alpha1, beta1 ]

    beta_kl = DeterministicVI.gen_beta_kl(alpha2, beta2)
    function beta_kl_wrapper{NumType <: Number}(par::Vector{NumType})
        alpha1 = par[1]
        beta1 = par[2]
        local kl, grad, hess
        kl, grad, hess = beta_kl(alpha1, beta1, false)
        return kl
    end

    kl, grad, hess = beta_kl(alpha1, beta1, true)

    ad_grad = ForwardDiff.gradient(beta_kl_wrapper, par)
    ad_hess = ForwardDiff.hessian(beta_kl_wrapper, par)

    # Check the derivatives
    @test_approx_eq beta_kl_wrapper(par) kl
    @test_approx_eq ad_grad grad
    @test_approx_eq ad_hess hess

    # Check the value
    p1_dist = Beta(alpha1, beta1)
    p2_dist = Beta(alpha2, beta2)
    test_kl_value(p1_dist, p2_dist, kl)
end

function test_categorical_kl_derivatives()
    p1 = Float64[ 1, 2, 3, 4]
    p2 = Float64[ 5, 6, 2, 1]

    p1 = p1 / sum(p1)
    p2 = p2 / sum(p2)

    categorical_kl = DeterministicVI.gen_categorical_kl(p2)

    function categorical_kl_wrapper{NumType <: Number}(par::Vector{NumType})
        local kl, grad, hess
        kl, grad, hess = categorical_kl(par, false)
        return kl
    end

    kl, grad, hess = categorical_kl(p1, true);

    ad_grad = ForwardDiff.gradient(categorical_kl_wrapper, p1)
    ad_hess = ForwardDiff.hessian(categorical_kl_wrapper, p1)

    # Check the derivatives
    @test_approx_eq categorical_kl_wrapper(p1) kl
    @test_approx_eq ad_grad grad
    @test_approx_eq ad_hess hess

    # Check the value
    p1_dist = Categorical(p1)
    p2_dist = Categorical(p2)
    test_kl_value(p1_dist, p2_dist, kl)
end

function test_diagmvn_mvn_kl_derivatives()
    K = 4
    mean1 = rand(K)
    var1 = rand(K)
    var1 = var1 .* var1

    mean2 = rand(K)
    cov2 = rand(K, K)
    cov2 = 0.2 * cov2 * cov2' + eye(K)
    cov2 = 0.5 * (cov2 + cov2')

    diagmvn_mvn_kl = DeterministicVI.gen_diagmvn_mvn_kl(mean2, cov2)
    kl, grad_mean, grad_var, hess_mean, hess_var =
        diagmvn_mvn_kl(mean1, var1, true);

    hess = zeros(Float64, 2 * K, 2 * K)
    hess[1:K, 1:K] = hess_mean
    hess[(K + 1):(2 * K), (K + 1):(2 * K)] = hess_var

    function diagmvn_mvn_kl_wrapper{NumType <: Number}(par::Vector{NumType})
        K = Int(length(par) / 2)
        mean1 = par[1:K]
        var1 = par[(K + 1):(2 * K)]
        local kl, grad_mean, grad_var, hess_mean, hess_var
        kl, grad_mean, grad_var, hess_mean, hess_var =
            diagmvn_mvn_kl(mean1, var1, false)
        return kl
    end

    par = vcat(mean1, var1)
    ad_grad = ForwardDiff.gradient(diagmvn_mvn_kl_wrapper, par)
    ad_hess = ForwardDiff.hessian(diagmvn_mvn_kl_wrapper, par)

    # Check the derivatives
    @test_approx_eq diagmvn_mvn_kl_wrapper(par) kl
    @test_approx_eq ad_grad vcat(grad_mean, grad_var)
    @test_approx_eq ad_hess hess

    # Check the value
    p1_dist = MvNormal(mean1, diagm(var1))
    p2_dist = MvNormal(mean2, cov2)
    test_kl_value(p1_dist, p2_dist, kl)
end

function test_normal_kl_derivatives()
    mean1 = 0.5
    var1 = 2.0

    mean2 = 0.8
    var2 = 1.8

    normal_kl = DeterministicVI.gen_normal_kl(mean2, var2)

    function normal_kl_wrapper{NumType <: Number}(par::Vector{NumType})
        local kl, grad, hess
        kl, grad, hess = normal_kl(par[1], par[2], false)
        return kl
    end

    kl, grad, hess = normal_kl(mean1, var1, true)

    par = vcat(mean1, var1)
    ad_grad = ForwardDiff.gradient(normal_kl_wrapper, par)
    ad_hess = ForwardDiff.hessian(normal_kl_wrapper, par)

    # Check the derivatives
    @test_approx_eq normal_kl_wrapper(par) kl
    @test_approx_eq ad_grad grad
    @test_approx_eq ad_hess hess

    # Check the value
    p1_dist = Normal(mean1, sqrt(var1))
    p2_dist = Normal(mean2, sqrt(var2))
    test_kl_value(p1_dist, p2_dist, kl)
end

########################
# Test derivatives of sensitive float manipulations

function test_combine_sfs()
  # TODO: this test was designed for multiply_sf.  Make it more general.

  # Two sets of ids with some overlap and some disjointness.
  p = length(ids)
  S = 2

  ids1 = find((1:p) .% 2 .== 0)
  ids2 = setdiff(1:p, ids1)
  ids1 = union(ids1, 1:5)
  ids2 = union(ids2, 1:5)

  l1 = zeros(Float64, S * p);
  l2 = zeros(Float64, S * p);
  l1[ids1] = rand(length(ids1))
  l2[ids2] = rand(length(ids2))
  l1[ids1 + p] = rand(length(ids1))
  l2[ids2 + p] = rand(length(ids2))

  sigma1 = zeros(Float64, S * p, S * p);
  sigma2 = zeros(Float64, S * p, S * p);
  sigma1[ids1, ids1] = rand(length(ids1), length(ids1));
  sigma2[ids2, ids2] = rand(length(ids2), length(ids2));
  sigma1[ids1 + p, ids1 + p] = rand(length(ids1), length(ids1));
  sigma2[ids2 + p, ids2 + p] = rand(length(ids2), length(ids2));
  sigma1 = 0.5 * (sigma1 + sigma1');
  sigma2 = 0.5 * (sigma2 + sigma2');

  x = 0.1 * rand(S * p);

  function base_fun1{T <: Number}(x::Vector{T})
    (l1' * x + 0.5 * x' * sigma1 * x)[1,1]
  end

  function base_fun2{T <: Number}(x::Vector{T})
    (l2' * x + 0.5 * x' * sigma2 * x)[1,1]
  end

  function multiply_fun{T <: Number}(x::Vector{T})
    base_fun1(x) * base_fun1(x)
  end

  function combine_fun{T <: Number}(x::Vector{T})
    (base_fun1(x) ^ 2) * sqrt(base_fun2(x))
  end

  function combine_fun_derivatives{T <: Number}(x::Vector{T})
    g_d = T[2 * base_fun1(x) * sqrt(base_fun2(x)),
            0.5 * (base_fun1(x) ^ 2) / sqrt(base_fun2(x)) ]
    g_h = zeros(T, 2, 2)
    g_h[1, 1] = 2 * sqrt(base_fun2(x))
    g_h[2, 2] = -0.25 * (base_fun1(x) ^ 2) * (base_fun2(x) ^(-3/2))
    g_h[1, 2] = g_h[2, 1] = base_fun1(x) / sqrt(base_fun2(x))
    g_d, g_h
  end

  s_ind = Array(UnitRange{Int}, 2);
  s_ind[1] = 1:p
  s_ind[2] = (1:p) + p

  ret1 = zero_sensitive_float(CanonicalParams, Float64, S);
  ret1.v[1] = base_fun1(x)
  fill!(ret1.d, 0.0);
  fill!(ret1.h, 0.0);
  for s=1:S
    ret1.d[:, s] = l1[s_ind[s]] + sigma1[s_ind[s], s_ind[s]] * x[s_ind[s]];
    ret1.h[s_ind[s], s_ind[s]] = sigma1[s_ind[s], s_ind[s]];
  end

  ret2 = zero_sensitive_float(CanonicalParams, Float64, S);
  ret2.v[1] = base_fun2(x)
  fill!(ret2.d, 0.0);
  fill!(ret2.h, 0.0);
  for s=1:S
    ret2.d[:, s] = l2[s_ind[s]] + sigma2[s_ind[s], s_ind[s]] * x[s_ind[s]];
    ret2.h[s_ind[s], s_ind[s]] = sigma2[s_ind[s], s_ind[s]];
  end

  grad = ForwardDiff.gradient(base_fun1, x);
  hess = ForwardDiff.hessian(base_fun1, x);
  for s=1:S
    @test_approx_eq(ret1.d[:, s], grad[s_ind[s]])
  end
  @test_approx_eq(ret1.h, hess)

  grad = ForwardDiff.gradient(base_fun2, x);
  hess = ForwardDiff.hessian(base_fun2, x);
  for s=1:S
    @test_approx_eq(ret2.d[:, s], grad[s_ind[s]])
  end
  @test_approx_eq(ret2.h, hess)

  # Test the combinations.
  v = combine_fun(x);
  grad = ForwardDiff.gradient(combine_fun, x);
  hess = ForwardDiff.hessian(combine_fun, x);

  sf1 = deepcopy(ret1);
  sf2 = deepcopy(ret2);
  g_d, g_h = combine_fun_derivatives(x)
  combine_sfs!(sf1, sf2, sf1.v[1] ^ 2 * sqrt(sf2.v[1]), g_d, g_h);

  @test_approx_eq sf1.v[1] v
  @test_approx_eq sf1.d[:] grad
  @test_approx_eq sf1.h hess
end

function test_add_sources_sf()
  P = length(CanonicalParams)
  S = 2

  sf_all = zero_sensitive_float(CanonicalParams, Float64, S);
  sf_s = zero_sensitive_float(CanonicalParams, Float64, 1);

  function scaled_exp!{NumType <: Number}(
      sf::SensitiveFloat{CanonicalParams, NumType},
      x::Vector{NumType}, a::Vector{Float64})

    sf.v[1] = one(NumType)
    for p in 1:P
      sf.v[1] *= exp(sum(x[p] * a[p]))
    end
    if NumType == Float64
      sf.d[:, 1] = sf.v[1] * a
      sf.h[:, :] = sf.v[1] * (a * a')
    end
  end

  a1 = rand(P);
  function f1{NumType <: Number}(x::Vector{NumType})
    sf_local = zero_sensitive_float(CanonicalParams, NumType, 1);
    scaled_exp!(sf_local, x, a1)
    sf_local.v[1]
  end

  a2 = rand(P);
  function f2{NumType <: Number}(x::Vector{NumType})
    sf_local = zero_sensitive_float(CanonicalParams, NumType, 1);
    scaled_exp!(sf_local, x, a2)
    sf_local.v[1]
  end

  x1 = rand(P);
  clear!(sf_s);
  scaled_exp!(sf_s, x1, a1);
  v1 = sf_s.v[1]

  fd_grad1 = ForwardDiff.gradient(f1, x1);
  @test_approx_eq sf_s.d fd_grad1

  fd_hess1 = ForwardDiff.hessian(f1, x1);
  @test_approx_eq sf_s.h fd_hess1

  add_sources_sf!(sf_all, sf_s, 1, true)

  x2 = rand(P);
  clear!(sf_s);
  scaled_exp!(sf_s, x2, a2);
  v2 = sf_s.v

  fd_grad2 = ForwardDiff.gradient(f2, x2);
  @test_approx_eq sf_s.d fd_grad2

  fd_hess2 = ForwardDiff.hessian(f2, x2);
  @test_approx_eq sf_s.h fd_hess2

  add_sources_sf!(sf_all, sf_s, 2, true)

  @test_approx_eq (v1 + v2) sf_all.v

  @test_approx_eq (v1 + v2) sf_all.v
  @test_approx_eq fd_grad1 sf_all.d[1:P, 1]
  @test_approx_eq fd_grad2 sf_all.d[1:P, 2]
  @test_approx_eq fd_hess1 sf_all.h[1:P, 1:P]
  @test_approx_eq fd_hess2 sf_all.h[(1:P) + P, (1:P) + P]
  @test_approx_eq zeros(P, P) sf_all.h[(1:P), (1:P) + P]
  @test_approx_eq zeros(P, P) sf_all.h[(1:P) + P, (1:P)]
end

#####################################
# Transforms

function test_box_derivatives()
	blob, ea, body = gen_three_body_dataset();
	transform = Transform.get_mp_transform(ea.vp, ea.active_sources, loc_width=1.0);

	box_params = setdiff(fieldnames(ids), [:a, :k])
	vp_free = transform.from_vp(ea.vp)
	for sa = 1:length(ea.active_sources), param in box_params, ind in length(getfield(ids, param))
		# sa = 1
		# param = box_params[1]
		# ind = 1

		s = ea.active_sources[sa]
		vp_ind = getfield(ids, param)[ind]
		free_ind = [getfield(ids_free, param)[ind]]

		function wrap_transform{NumType <: Number}(vp_free_s::Vector{NumType})
			local_vp_free =
				Array{NumType, 1}[ convert(Array{NumType, 1}, vp_free[sa]) for
			                  sa = 1:length(ea.active_sources) ]
			local_vp_free[s] = vp_free_s
			vp = transform.to_vp(local_vp_free)
			vp[s][getfield(ids, param)[ind]]
		end

		ad_d  = ForwardDiff.gradient(wrap_transform, vp_free[s])[free_ind][1]
		ad_h = ForwardDiff.hessian(wrap_transform, vp_free[s])[free_ind, free_ind][1,1]

		d, h = Transform.box_derivatives(
			ea.vp[s][vp_ind][1], transform.bounds[s][param][ind])
		@test_approx_eq ad_d d
		@test_approx_eq ad_h h
	end
end

function test_box_simplex_derivatives()
	blob, ea, body = gen_three_body_dataset();
	for s = 1:ea.S
		delta = 0.01 * s # Make the parameters different for each one
		ea.vp[s][ids.a[:, 1]] = Float64[ 0.2 - delta, 0.8 + delta ]
		ea.vp[s][ids.k] = Float64[ 0.2- delta 0.2- delta; 0.8 + delta 0.8 + delta ]
	end
	transform = Transform.get_mp_transform(ea.vp, ea.active_sources, loc_width=1.0);

	simplex_params = [:a, :k]
	vp_free = transform.from_vp(ea.vp)

  for sa = 1:length(ea.active_sources), param in simplex_params
    # sa = 1
    # param = :a
    # col = 1 # For k only
    # ind = 1 # Index within the simplex

    s = ea.active_sources[sa]
    num_rows = size(getfield(ids, param))[1]
    num_cols = size(getfield(ids, param))[2]

    for col = 1:num_cols
    	vp_ind = getfield(ids, param)[:, col]

    	if length(size(getfield(ids_free, param))) == 0
    		# Hack to handle ids_free.a
    		@assert col == 1
    		free_ind = [ getfield(ids_free, param) ]
    	else
    		free_ind = getfield(ids_free, param)[:, col]
    	end

    	d, h = Transform.box_simplex_derivatives(
    		ea.vp[s][vp_ind], transform.bounds[s][param][col])

    	for row = 1:num_rows
    		# Write with a univariate output so we can take autodiff hessians.
    	  function wrap_transform{NumType <: Number}(vp_free_s::Vector{NumType})
    	  	local_vp_free =
    	  		Array{NumType, 1}[ convert(Array{NumType, 1}, vp_free[sa]) for
    	  	                     sa = 1:length(ea.active_sources) ]
    	  	local_vp_free[s] = vp_free_s
    	  	vp = transform.to_vp(local_vp_free)
    	  	vp[s][getfield(ids, param)[row, col]]
    	  end

    	  ad_d = ForwardDiff.gradient(wrap_transform, vp_free[s])[free_ind]
    	  ad_h = ForwardDiff.hessian(wrap_transform, vp_free[s])[free_ind, free_ind]
    		@test_approx_eq ad_d d[row, :][1]
    		@test_approx_eq ad_h h[row]
    	end
    end
  end
end

function test_simplex_derivatives()
	n = 4
	basic_simplex_box = Transform.SimplexBox(0, 1, n)
	z = Float64[1, 2, 4, 3]
	z /= sum(z)
	r = Transform.unsimplexify_parameter(z, basic_simplex_box)
	Transform.simplexify_parameter(r, basic_simplex_box)

	ad_d = Array(Array{Float64}, n)
	ad_h = Array(Array{Float64}, n)

	for ind = 1:n
	  function wrap_simplex{NumType <: Number}(r::Vector{NumType})
	    local z = Transform.simplexify_parameter(r, basic_simplex_box)
	    z[ind]
	  end
	  ad_d[ind] = ForwardDiff.gradient(wrap_simplex, r)
	  ad_h[ind] = ForwardDiff.hessian(wrap_simplex, r)
	end

	jacobian, hessian_vec = Transform.simplex_derivatives(z)

	@test_approx_eq jacobian' reduce(hcat, ad_d)
	for ind = 1:n
		@test_approx_eq(hessian_vec[ind], ad_h[ind])
	end
end

###################################
# Run tests

# ELBO tests:
println("Running ELBO derivative tests.")
@time test_combine_pixel_sources()
@time test_fs1m_derivatives()
@time test_fs0m_derivatives()
@time test_bvn_derivatives()
@time test_galaxy_variable_transform()
@time test_galaxy_cache_component()
@time test_galaxy_sigma_derivs()
@time test_brightness_hessian()
@time test_dsiginv_dsig()
@time test_add_log_term()
@time test_e_g_s_functions()
test_process_active_pixels()

# KL tests:
println("Running KL derivative tests.")
@time test_beta_kl_derivatives()
@time test_categorical_kl_derivatives()
@time test_diagmvn_mvn_kl_derivatives()
@time test_normal_kl_derivatives()

# SensitiveFloat tests:
println("Running SensitiveFloat derivative tests.")
@time test_combine_sfs()
@time test_add_sources_sf()

# Transform tests:
println("Running Transform derivative tests.")
@time test_box_derivatives()
@time test_box_simplex_derivatives()
@time test_simplex_derivatives()


text:
# Statistics for an array of scalars

#############################
#
#   Moments
#
#############################

# Skewness
# This is Type 1 definition according to Joanes and Gill (1998)
function skewness(v::RealArray, m::Real)
    n = length(v)
    cm2 = 0.0   # empirical 2nd centered moment (variance)
    cm3 = 0.0   # empirical 3rd centered moment
    for x_i in v
        z = x_i - m
        z2 = z * z

        cm2 += z2
        cm3 += z2 * z
    end
    cm3 /= n
    cm2 /= n
    return cm3 / sqrt(cm2 * cm2 * cm2)  # this is much faster than cm2^1.5
end

function skewness(v::RealArray, wv::WeightVec, m::Real)
    n = length(v)
    length(wv) == n || raise_dimerror()
    cm2 = 0.0   # empirical 2nd centered moment (variance)
    cm3 = 0.0   # empirical 3rd centered moment
    w = values(wv)

    @inbounds for i = 1 : n
        x_i = v[i]
        w_i = w[i]
        z = x_i - m
        z2w = z * z * w_i
        cm2 += z2w
        cm3 += z2w * z
    end
    sw = sum(wv)
    cm3 /= sw
    cm2 /= sw
    return cm3 / sqrt(cm2 * cm2 * cm2)  # this is much faster than cm2^1.5
end

skewness(v::RealArray) = skewness(v, mean(v))
skewness(v::RealArray, wv::WeightVec) = skewness(v, wv, mean(v, wv))

# (excessive) Kurtosis
# This is Type 1 definition according to Joanes and Gill (1998)
function kurtosis(v::RealArray, m::Real)
    n = length(v)
    cm2 = 0.0  # empirical 2nd centered moment (variance)
    cm4 = 0.0  # empirical 4th centered moment
    for x_i in v
        z = x_i - m
        z2 = z * z
        cm2 += z2
        cm4 += z2 * z2
    end
    cm4 /= n
    cm2 /= n
    return (cm4 / (cm2 * cm2)) - 3.0
end

function kurtosis(v::RealArray, wv::WeightVec, m::Real)
    n = length(v)
    length(wv) == n || raise_dimerror()
    cm2 = 0.0  # empirical 2nd centered moment (variance)
    cm4 = 0.0  # empirical 4th centered moment
    w = values(wv)

    @inbounds for i = 1 : n
        x_i = v[i]
        w_i = w[i]
        z = x_i - m
        z2 = z * z
        z2w = z2 * w_i
        cm2 += z2w
        cm4 += z2w * z2
    end
    sw = sum(wv)
    cm4 /= sw
    cm2 /= sw
    return (cm4 / (cm2 * cm2)) - 3.0
end

kurtosis(v::RealArray) = kurtosis(v, mean(v))
kurtosis(v::RealArray, wv::WeightVec) = kurtosis(v, wv, mean(v, wv))

#############################
#
#   Variability measurements
#
#############################

# Variation: std / mean
variation{T<:Real}(x::AbstractArray{T}, m::Real) = stdm(x, m) / m
variation{T<:Real}(x::AbstractArray{T}) = variation(x, mean(x))

# Standard error of the mean: std(a)
sem{T<:Real}(a::AbstractArray{T}) = sqrt(var(a) / length(a))

# Median absolute deviation
mad{T<:Real}(v::AbstractArray{T}, center::Real) = 1.4826 * median!(abs(v-center))

function mad!{T<:Real}(v::AbstractArray{T}, center::Real)
    for i in 1:length(v)
        v[i] = abs(v[i]-center)
    end
    1.4826 * median!(v, checknan=false)
end

mad!{T<:Real}(v::AbstractArray{T}) = mad!(v, median!(v))

mad{T<:Real}(v::AbstractArray{T}) = mad!(copy(v))
mad{T<:Real}(v::Range1{T}) = mad!([v])

#############################
#
#   min/max related
#
#############################

# middle
middle{T<:FloatingPoint}(a1::T, a2::T) = (a1 + a2) / convert(T, 2)
middle{T<:Integer}(a1::T, a2::T) = (a1 + a2) / 2

# Mid-range
midrange{T<:Real}(a::AbstractArray{T}) = middle(extrema(a)...)

# Range: xmax - xmin
range{T<:Real}(a::AbstractArray{T}) = ((m0, m1) = extrema(a); m1 - m0)

#############################
#
#   quantile and friends
#
#############################

percentile{T<:Real}(v::AbstractArray{T}, p) = quantile(v, p * 0.01)
iqr{T<:Real}(v::AbstractArray{T}) = (q = quantile(v, [.25, .75]); q[2] - q[1])

quantile{T<:Real}(v::AbstractArray{T}) = quantile(v, [.0, .25, .5, .75, 1.0])
nquantile{T<:Real}(v::AbstractArray{T}, n::Integer) = quantile(v, (0:n)/n)

#############################
#
#   mode & modes
#
#############################

# compute mode, given the range of integer values
function mode{T<:Integer}(a::AbstractArray{T}, rgn::Range1{T})
    isempty(a) && error("mode: input array cannot be empty.")

    r0 = rgn[1]
    r1 = rgn[end]
    cnts = zeros(Int, length(rgn))

    mc = 0    # maximum count
    mv = r0   # a value corresponding to maximum count

    for x in a
        if r0 <= x <= r1
            @inbounds c = (cnts[x - r0 + 1] += 1)
            if c > mc
                mc = c
                mv = x
            end
        end
    end

    return mv
end

function modes{T<:Integer}(a::AbstractArray{T}, rgn::Range1{T})
    r0 = rgn[1]
    r1 = rgn[end]
    n = length(rgn)
    cnts = zeros(Int, n)

    # find the maximum count
    mc = 0
    for x in a
        if r0 <= x <= r1
            @inbounds c = (cnts[x - r0 + 1] += 1)
            if c > mc
                mc = c
            end
        end
    end

    # find all values corresponding to maximum count
    ms = T[]
    for i = 1 : n
        @inbounds if cnts[i] == mc
            push!(ms, rgn[i])
        end
    end

    return ms
end

function mode{T}(a::AbstractArray{T})
    isempty(a) && error("mode: input array cannot be empty.")

    cnts = (T=>Int)[]

    # first element
    mc = 1
    mv = a[1]
    cnts[mv] = 1

    # find the mode along with table construction
    for i = 2 : length(a)
        @inbounds x = a[i]
        if haskey(cnts, x)
            c = (cnts[x] += 1)
            if c > mc
                mc = c
                mv = x
            end
        else
            cnts[x] = 1
            # in this case: c = 1, and thus c > mc won't happen
        end
    end

    return mv
end

function modes{T}(a::AbstractArray{T})
    isempty(a) && error("modes: input array cannot be empty.")

    cnts = (T=>Int)[]

    # first element
    mc = 1
    cnts[a[1]] = 1

    # find the mode along with table construction
    for i = 2 : length(a)
        @inbounds x = a[i]
        if haskey(cnts, x)
            c = (cnts[x] += 1)
            if c > mc
                mc = c
            end
        else
            cnts[x] = 1
            # in this case: c = 1, and thus c > mc won't happen
        end
    end

    # find values corresponding to maximum counts
    ms = T[]
    for (x, c) in cnts
        if c == mc
            push!(ms, x)
        end
    end

    return ms
end

#############################
#
#   summary
#
#############################

immutable SummaryStats{T<:FloatingPoint}
    mean::T
    min::T
    q25::T
    median::T
    q75::T
    max::T
end

function summarystats{T<:Real}(a::AbstractArray{T})
    m = mean(a)
    qs = quantile(a, [0.00, 0.25, 0.50, 0.75, 1.00])
    R = typeof(convert(FloatingPoint, zero(T)))
    SummaryStats{R}(
        convert(R, m),
        convert(R, qs[1]),
        convert(R, qs[2]),
        convert(R, qs[3]),
        convert(R, qs[4]),
        convert(R, qs[5]))
end

function Base.show(io::IO, ss::SummaryStats)
    println(io, "Summary Stats:")
    @printf(io, "Mean:         %.6f\n", ss.mean)
    @printf(io, "Minimum:      %.6f\n", ss.min)
    @printf(io, "1st Quartile: %.6f\n", ss.q25)
    @printf(io, "Median:       %.6f\n", ss.median)
    @printf(io, "3rd Quartile: %.6f\n", ss.q75)
    @printf(io, "Maximum:      %.6f\n", ss.max)
end

describe{T<:Real}(a::AbstractArray{T}) = show(summarystats(a))

ext: jl


text:
using DelayDiffEq, DiffEqProblemLibrary.DDEProblemLibrary
using Test

DDEProblemLibrary.importddeproblems()

# simple problems
@testset "simple problems" begin
  prob_ip = DDEProblemLibrary.prob_dde_constant_1delay_ip
  prob_scalar = DDEProblemLibrary.prob_dde_constant_1delay_scalar
  ts = 0:0.1:10

  # Vern6
  println("Vern6")
  @testset "Vern6" begin
    alg = MethodOfSteps(Vern6())
    sol_ip = solve(prob_ip, alg)

    @test sol_ip.errors[:l∞] < 8.7e-4
    @test sol_ip.errors[:final] < 3.9e-6
    @test sol_ip.errors[:l2] < 5.4e-4

    sol_scalar = solve(prob_scalar, alg)

    # fails due to floating point issues
    if Sys.WORD_SIZE == 32
      @test_broken sol_ip(ts, idxs=1) ≈ sol_scalar(ts)
    else
      @test sol_ip(ts, idxs=1) ≈ sol_scalar(ts)
    end
  end

  # Vern7
  println("Vern7")
  @testset "Vern7" begin
    alg = MethodOfSteps(Vern7())
    sol_ip = solve(prob_ip, alg)

    @test sol_ip.errors[:l∞] < 4.0e-4
    @test sol_ip.errors[:final] < 3.5e-7
    @test sol_ip.errors[:l2] < 1.9e-4

    sol_scalar = solve(prob_scalar, alg)

    @test sol_ip(ts, idxs=1) ≈ sol_scalar(ts)
  end

  # Vern8
  println("Vern8")
  @testset "Vern8" begin
    alg = MethodOfSteps(Vern8())
    sol_ip = solve(prob_ip, alg)

    @test sol_ip.errors[:l∞] < 2.0e-3
    @test sol_ip.errors[:final] < 1.8e-5
    @test sol_ip.errors[:l2] < 8.4e-4

    sol_scalar = solve(prob_scalar, alg)

    @test sol_ip(ts, idxs=1) ≈ sol_scalar(ts)
  end

  # Vern9
  println("Vern9")
  @testset "Vern9" begin
    alg = MethodOfSteps(Vern9())
    sol_ip = solve(prob_ip, alg)

    @test sol_ip.errors[:l∞] < 1.5e-3
    @test sol_ip.errors[:final] < 3.8e-6
    @test sol_ip.errors[:l2] < 6.5e-4

    sol_scalar = solve(prob_scalar, alg)

    @test sol_ip(ts, idxs=1) ≈ sol_scalar(ts)
  end
end

# model of Mackey and Glass
println("Mackey and Glass")
@testset "Mackey and Glass" begin
  prob = DDEProblemLibrary.prob_dde_DDETST_A1

  # Vern6
  solve(prob, MethodOfSteps(Vern6()))

  # Vern7
  solve(prob, MethodOfSteps(Vern7()))

  # Vern8
  solve(prob, MethodOfSteps(Vern8()))

  # Vern9
  solve(prob, MethodOfSteps(Vern9()))
end


text:
function angtod2hitmap(nside::Int, theta_tod::Array{T}, phi_tod::Array{T}) where {T}
    hit_map = zeros(Int64, nside2npix(nside))
    resol = Resolution(nside)
    for j in eachindex(theta_tod[1,:])
        theta_tod_jth_det = @views theta_tod[:,j]
        phi_tod_jth_det = @views phi_tod[:,j]
        @inbounds @simd for k = eachindex(theta_tod[:,1])
            ipix = @views ang2pixRing(resol, theta_tod_jth_det[k], phi_tod_jth_det[k])
            hit_map[ipix] += 1
        end
    end
    return hit_map
end

function xyztod2hitmap(nside::Int, xyz_tod::Array{T}) where {T}
    hit_map = zeros(Int64, nside2npix(nside))
    resol = Resolution(nside)
    for j in eachindex(xyz_tod[1,1,:])
        x_tod_jth_det = @views xyz_tod[1,:,j]
        y_tod_jth_det = @views xyz_tod[2,:,j]
        z_tod_jth_det = @views xyz_tod[3,:,j]
        @inbounds @simd for k = eachindex(xyz_tod[1,:,1])
            ang = @views vec2ang(x_tod_jth_det[k], y_tod_jth_det[k], z_tod_jth_det[k])
            ipix = @views ang2pixRing(resol, ang[1], ang[2])
            hit_map[ipix] += 1
        end
    end
    return hit_map
end

function pixtod2hitmap(nside::Int, pixtod::Array{T}) where {T}
    npix = nside2npix(nside)
    hit_map = zeros(Int64, npix)
    for i = eachindex(pixtod)
        hit_map[pixtod[i]] += 1
    end
    return hit_map
end

function ScanningStrategy2map(SS::ScanningStrategy, division::Int)
    resol = Resolution(SS.nside)
    npix = nside2npix(SS.nside)

    month = Int(SS.duration / division)
    ω_hwp = rpm2angfreq(SS.hwp_rpm)

    hit_map = zeros(npix)
    Cross = zeros(npix, 4, 2)
    BEGIN = 0
    p = Progress(division)
    @views @inbounds for i = 1:division
        END = i * month
        pix_tod, psi_tod, time_array = get_pointing_pixels(SS, BEGIN, END)
        @views @inbounds for j = eachindex(psi_tod[1,:])
            pix_tod_jth_det = pix_tod[:,j]
            psi_tod_jth_det = ifelse(ω_hwp == 0.0, -psi_tod[:,j], psi_tod[:,j])
            @views @inbounds @simd for k = eachindex(psi_tod[:,1])
                t = time_array[k]
                ipix = pix_tod_jth_det[k]
                psi = psi_tod_jth_det[k]
                hwp_ang = 4ω_hwp*t

                hit_map[ipix] += 1
                Cross[ipix,1,1] += sin(hwp_ang - psi)
                Cross[ipix,1,2] += cos(hwp_ang - psi)
                Cross[ipix,2,1] += sin(hwp_ang - 2psi)
                Cross[ipix,2,2] += cos(hwp_ang - 2psi)
                Cross[ipix,3,1] += sin(hwp_ang - 3psi)
                Cross[ipix,3,2] += cos(hwp_ang - 3psi)
                Cross[ipix,4,1] += sin(hwp_ang - 4psi)
                Cross[ipix,4,2] += cos(hwp_ang - 4psi)
            end
        end
        BEGIN = END
        next!(p)
    end
    link1 = @views @. (Cross[:,1,1]/hit_map)^2 + (Cross[:,1,2]/hit_map)^2
    link2 = @views @. (Cross[:,2,1]/hit_map)^2 + (Cross[:,2,2]/hit_map)^2
    link3 = @views @. (Cross[:,3,1]/hit_map)^2 + (Cross[:,3,2]/hit_map)^2
    link4 = @views @. (Cross[:,4,1]/hit_map)^2 + (Cross[:,4,2]/hit_map)^2
    out_map = @views [hit_map, link1, link2, link3, link4]
    return out_map
end

function array2map(map_array::Array)
    nside = npix2nside(length(map_array))
    m = HealpixMap{Float64, RingOrder}(nside)
    m.pixels .= map_array
    return m
end

@inline function Hᵢ(t)
    ω_HWP = @views 2 * π * (46 / 60)
    s = @views sin(4*ω_HWP*t)
    c = @views cos(4*ω_HWP*t)
    H = @views @SMatrix [
        1 0  0
        0 c  s
        0 s -c
    ]
    return H
end

@inline function HitMatrix(rho)
    s = sin(rho)
    c = cos(rho)
    W = @views (1/4).* @SMatrix [
        1 c   s
        c c^2 s*c
        s s*c s^2
    ]
    return W
end

@inline function pᵢ(pix_i, psi, t, obsmap)
    _wᵢ(psi) = @SMatrix [1.0/2.0; cos(2psi) / 2.0; sin(2psi) / 2.0]
    I = @views obsmap[1,:]
    Q = @views obsmap[2,:]
    U = @views obsmap[3,:]
    return _wᵢ(psi)' * (Hᵢ(t) * [I[pix_i]; Q[pix_i]; U[pix_i]])
end

_wᵢ(psi) = @SMatrix [1.0/2.0; cos(2psi) / 2.0; sin(2psi) / 2.0]


text:
using Documenter

using STMOZOO
using STMOZOO.Example
using STMOZOO.BeesAlgorithm

makedocs(sitename="STMO ZOO",
    format = Documenter.HTML(),
    modules=[BeesAlgorithm], # add your module
    pages=Any[
        "BeesAlgorithm" => "man/beesalgorithm.md",  # add the page to your documentation
    ])

#=
deploydocs(
            repo = "github.com/michielstock/STMOZOO.jl.git",
        )
=#

max_stars_repo_name: kirstvh/STMOZOO

max_issues_repo_name: kirstvh/STMOZOO

max_forks_repo_name: kirstvh/STMOZOO


text:
import NeurodataWithoutBorders
using FileIO
using Base.Test
using Unitful

fname = "$(tempname()).nwb"

B = Array(typeof(1.0u"mV"),1000,10)
for i in eachindex(B)
	B[i] = rand()*1.0u"mV"
end

eldata = NeurodataWithoutBorders.ElectricalSeries(B, "test", 1.0, collect(1:10), "test3",NeurodataWithoutBorders.range(0.0u"s", (1/30000.0)u"s", size(B,1)))
NeurodataWithoutBorders.save(File(format"NWB",fname), eldata)

qdata = NeurodataWithoutBorders.SpatialSeries([rand()*Unitful.cm for i in 1:1000, j in 1:2], "test", "test4", range(0.0*Unitful.s, (1.0/1000)*Unitful.s, 1000))
NeurodataWithoutBorders.save(File(format"NWB",fname), qdata)

wdata = NeurodataWithoutBorders.load(File(format"NWB", fname))

#rm(fname)
@test eldata.data == wdata[1].data
@test eldata.help == wdata[1].help
@test eldata.timestamps == wdata[1].timestamps

@test qdata.data == wdata[2].data
@test qdata.help == wdata[2].help
@test qdata.timestamps == wdata[2].timestamps

pth = NeurodataWithoutBorders.get_path(FileIO.File(format"NWB", fname))
@test "test3" in pth["ElectricalSeries"]
@test "test4" in pth["SpatialSeries"]


text:
using StatsBase
using Statistics

include("NeuralNet.jl")
foo = NeuralNet_predict(wbSGDBb,[ones(size(xvalid,1)) xvalid],nHidden)
fooR = round.(foo)
fooU = sort(unique(foor))

egg = NeuralNet_predict(wbVanb,[ones(size(xvalid,1)) xvalid],nHidden)
eggR = round.(bar2)
eggU = sort(unique(eggR))

bar1 = NeuralNet_predict(wbSGDB,[ones(size(xvalid,1)) xvalid],nHidden)
barr = round.(bar1)
barrU = sort(unique(barr))

bar2 = NeuralNet_predict(wbVan,[ones(size(xvalid,1)) xvalid],nHidden)
barr2 = round.(bar2)
barr2U = sort(unique(barr2))

x = 1
function abc()
    # x = 1
    for i in 1:5
        global x
        x = i
    end
    print(x)
end

#---------------
medyhat= NeuralNet_predict(wbSGDB,[ones(size(xvalid,1)) xvalid],nHidden)
medyhat_u = length(unique(medyhat))

medyhatv = NeuralNet_predict(bestw,[ones(size(xvalid,1)) xvalid],nHidden)
medyhatv_u = length(unique(medyhatv))


text:
using ChainRulesCore
using ZygoteRules: @adjoint

using ..ITensorAutoHOOT
using ..ITensorAutoHOOT: AbstractNetwork

struct GeneralMPSTensor <: AbstractNetwork
  mps::MPS
end

function Base.:+(A::GeneralMPSTensor, B::GeneralMPSTensor; kwargs...)
  if length(A.mps) == 1 && order(A.mps[1]) == 0
    return GeneralMPSTensor(MPS([ITensor(A[] + B[])]))
  else
    return GeneralMPSTensor(+(A.mps, B.mps; kwargs...))
  end
end

function Base.getindex(t::GeneralMPSTensor)
  @assert length(t.mps) == 1
  @assert order(t.mps[1]) == 0
  return t.mps[1][]
end

function ChainRulesCore.rrule(::typeof(getindex), x::GeneralMPSTensor)
  y = x[]
  function getindex_pullback(ȳ)
    x̄ = GeneralMPSTensor(MPS([ITensor(ȳ)]))
    return (NoTangent(), x̄)
  end
  return y, getindex_pullback
end

"""
initialization methods
"""
# TODO: general tags are not comparable
Base.isless(a::Index, b::Index) = id(a) < id(b) || (id(a) == id(b) && plev(a) < plev(b)) # && tags(a) < tags(b)

function GeneralMPSTensor(tensor::ITensor; cutoff, maxdim)
  mps_out = if size(tensor) == ()
    MPS([tensor])
  else
    MPS(tensor, sort(inds(tensor)); cutoff=cutoff, maxdim=maxdim)
  end
  return GeneralMPSTensor(mps_out)
end

@adjoint function GeneralMPSTensor(tensor::ITensor; cutoff, maxdim)
  MPSTensor_pullback(dtensor) = (ITensor(dtensor),)
  return GeneralMPSTensor(tensor; cutoff=cutoff, maxdim=maxdim), MPSTensor_pullback
end

"""
contract
"""
function general_mps_contract(mps1::MPS, mps2::MPS; cutoff, maxdim)
  ## TODO
  tensor = contract(vcat(collect(mps1), collect(mps2))...)
  return if size(tensor) == ()
    MPS([tensor])
  else
    MPS(tensor, sort(inds(tensor)); cutoff=cutoff, maxdim=maxdim)
  end
end

function ITensors.contract(t1::GeneralMPSTensor, t2::GeneralMPSTensor; cutoff, maxdim)
  mps_out = general_mps_contract(t1.mps, t2.mps; cutoff=cutoff, maxdim=maxdim)
  return GeneralMPSTensor(mps_out)
end


text:
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# ╔═╡ 0edbe34a-3714-11eb-1350-af72689d8855
md" ## _Data Structures_
In this lesson we will study how data can be collected and stored in memory. In particular we will deal with vectors, matrices, n-dimensional arrays, tuples and dictionaries.
"

# ╔═╡ 4ef2e454-3716-11eb-3014-774db1fb57d0
md" ## Vectors
A vector is a list of ordered data which share a common type (be it `Int`, `Float` or `Any`). Furthermore a vector is a one-dimensional array, and often “vector” and “array” are used a synonyms.

Contrarily to the mathematical definition of a vector, in programming a vector is simply a list of values and has no a priori geometrical meaning.

In Julia, to create a vector we use the following syntax:
"

# ╔═╡ 9833f658-3716-11eb-300e-d3b5d02b70e5
a = [1,2,3,4,5]

# ╔═╡ b2cbd7d8-3716-11eb-15cc-03c3fec2c981
b = [1.2, 3,4,5]

# ╔═╡ b2d549da-3716-11eb-0cce-b703452d7f34
c = ["Hello", "it's me", "TechyTok"]

# ╔═╡ bbe66bfa-3716-11eb-1929-ef53e271e354
md" We can access the members of an array using the indexing syntax: `array_name[element_number]`, for example, if we want to retrieve the third element of c we type `c[3]`. Contrarily to other programming languages, in **Julia vectors start at 1**, there is not much to say, it is just like that.
"

# ╔═╡ d5c03624-3716-11eb-38ed-3f395efc01f6
c[2], a[3], b[1]

# ╔═╡ e7768e7e-3716-11eb-366a-2f81166da844
md"If you want to append an element to an array, we can do it by using the `append!` function. Notice the !, this is a Julia convention to say that the function will modify the first argument given to the function, in fact `append!(a, 6)` results in `a = [1, 2, 3, 4, 5, 6]`.

Previously I have said that vector elements must share the same type, if we try to append (or add) a value to an array with a different type it will result into an error:
"

# ╔═╡ 5268b19c-3717-11eb-3854-ddd093a61057
append!(a, 3.14)

# ╔═╡ 6ced58c6-3717-11eb-1009-7909ebb6866b
md"In this case, we cannot convert a Float64 to an Int64 without losing precision, thus the error.

To check the type of an array we can type:
"

# ╔═╡ f2cb5918-3717-11eb-2b4b-77642fbeaf20
typeof(a)

# ╔═╡ ff24fcbc-3717-11eb-2ac3-6f18fde9484f
md"As we can see, `a` can only store `Int64` values or values that can be safely converted to `Int64` (such as `Int32` for example).

An alternative way to define a Vector, if we want to specify the type of the elements included, is by using a type name followed by a square bracket and the desired elements:
"

# ╔═╡ 21342288-3718-11eb-1b86-8d438637ac3a
d = Int[1,2,3,4,5]

# ╔═╡ 3fb54d2c-3718-11eb-305f-e7159394ac7b
md" ## Matrices
Matrices are two dimensional arrays. We can define a matrix with the following syntax:
"

# ╔═╡ 4c2a11c8-3718-11eb-0f59-8382e767058f
mat1 = [1 2 3; 4 5 6]

# ╔═╡ 538463b0-3718-11eb-03e7-171beded2a41
md"Rows are separated by semi colons ; and columns are separated by spaces.

When we want to access the elements of a matrix, we use the following notation `mat1[row_index, column_index]`. If we want to access the 2nd element of the first row we shall write `mat1[1, 2]`.
"

# ╔═╡ 6ed4e510-3718-11eb-3d00-9f60abc37333
md" ## N-Dimensional Arrays
Sometimes we need to create tables with more than 2 dimensions. In this case usually the tables tend to be big, so there is no explicit way to create an n-dimensional array. The suggested practice is to create an empty array first, using a new function called `zeros`, and then fill it either manually of using a loop.

For example, let’s suppose we want to create a 2x3x4 table, we would do it with `zeros(2,3,4)`. Let’s suppose we want to fill it with the product of the indexes, we can do it in the following way:

"

# ╔═╡ 9a27ee68-3718-11eb-1c7f-33eaf6e7cb3f
table = zeros(2,3,4)

# ╔═╡ ab83ddde-3718-11eb-0847-9fdeaca6d470
for k in 1:4
    for j in 1:3
        for i in 1:2
            table[i,j,k] = i*j*k
        end
    end
end

# ╔═╡ bbbc6720-3718-11eb-2c1e-978871157853
table

# ╔═╡ d3ce8154-3718-11eb-08ff-7b29208c5099
md" Please not that Julia stores values in memory differently from Python: in Julia to obtain fast loops we need to iter first over columns (which means that the first index must vary first and so on). For this reason if we plan to store, for example, 42 2x2 matrices, we need to create an array of size 2x2x42 (while in Python we would have created a 42x2x2 table).
"

# ╔═╡ df969238-3718-11eb-2196-a7a68b5a69f2
md" ## Slices
There is a convenient notation (called **slicing**) to access a subset of elements from an array. Let’s suppose we want to access the 2nd to 5th elements of an array of length 6, we can do it in the following way:

"

# ╔═╡ f82b9558-3718-11eb-0ef8-5fa509ba5b6f
begin
	e = [1,2,3,4,5,6]
	f = e[2:5]
end

# ╔═╡ 6ffa0c1a-3719-11eb-34d2-d7fc89a62e24
md"We can also use this notation to access a subset of a matrix, for example:
"

# ╔═╡ 7edb458c-3719-11eb-2bb3-4bd9173b530c
mat3 = reshape([i for i in 1:16],4,4)

# ╔═╡ 8bae5826-3719-11eb-084e-fbf852773ffd
mat4 = mat3[2:3, 2:3]

# ╔═╡ a17c0086-3719-11eb-1ab0-2381decb5422
md"On line 1 we have used a handy notation called **list comprehension**. `[i for i in 1:16]` means “create an array containing each i comprised from 1 to 16”. We then reshape it to have a size of 4x4 and store the result in `mat1`.

It is also possible to have nested list comprehensions like this one:
"

# ╔═╡ c62a2930-3719-11eb-3c28-750300385ad8
[i+j for i in 1:10 for j in 1:5]

# ╔═╡ e07fe63a-3719-11eb-05eb-81512953fb4b
md"This will create an array of length 50 containing all the sums of i+j.
"

# ╔═╡ 2655af78-371a-11eb-331c-310db7a9ccb1
md" ## Views
As in other programming languages, arrays are pointers to location in memory, thus we need to pay attention when we handle them.

If we create an array `g` and we assign a to `h` with `h=g`, the elements of `g` be modified by accessing `h`:
"

# ╔═╡ 3cdd35ea-371a-11eb-0f7f-135c3a766ecd
g=[1,2,3]

# ╔═╡ 67a3d7f2-371a-11eb-07cb-e3e233e2aab1
h=g

# ╔═╡ 67ad9e04-371a-11eb-36ad-9710ff0b0db3
g[2] = 42

# ╔═╡ 67bdc450-371a-11eb-08b2-55230f42c0e8
g

# ╔═╡ b2873584-371a-11eb-328d-37a39540905c
md"
This is particularly useful because it lets us save memory, but may have undesirable effects. If we want to make a copy of an array we need to use the function `copy`.
"

# ╔═╡ c17d2258-371a-11eb-25b5-933f413cd134
i=[1,2,3]

# ╔═╡ d902ab08-371a-11eb-31ce-8d50b3388b2b
j=copy(i)

# ╔═╡ d90d570e-371a-11eb-0e4b-075841a6a73a
j[2] = 42

# ╔═╡ dbfd0a06-371a-11eb-1bf5-e9112227d8b0
i

# ╔═╡ de4b38a0-371a-11eb-39b1-7575188e8063
j

# ╔═╡ ee1a296c-371a-11eb-0036-3b8ad05c6804
md"In some cases, when there are arrays containing other arrays, if we want to make a full copy of all the contents we need to use `deepcopy` instead."

# ╔═╡ f6fa2e36-371a-11eb-36e7-a7de898b78c5
md" ## Tuples
A tuple is a fixed size group of variables which may share a common type but don’t need to.

Unlike arrays, you cannot increase the size of a tuple once it has been created. Tuples are created using the following syntax:
"

# ╔═╡ 0e12b2be-371b-11eb-0e56-99ec2e4a224d
k =(1,2,3)

# ╔═╡ 187b73f0-371b-11eb-34e6-d7a2aeb00b58
l= 1,2,3

# ╔═╡ 1eb587ec-371b-11eb-2901-a1ebf0197aa7
md"
So tuples can be created by using regular brackets or no brackets at all! Tuples are really handy, as it is possible to “unpack” a tuple over many values:
"

# ╔═╡ 2bc6ba1e-371b-11eb-12d0-a3665fe9e1d0
tuple1 = (1, 2, 3)

# ╔═╡ 3985145c-371b-11eb-2dcc-5308e23dc5f0
m, n, o = tuple1

# ╔═╡ 39a4874c-371b-11eb-164c-692fd38680a9
("$m $n $o")

# ╔═╡ 58abede2-371b-11eb-24d3-2b16686937d5
md"
It is also possible to use tuples to emulate multiple return values from functions:
"

# ╔═╡ 6205bad0-371b-11eb-1cfb-3f1ff43c255e
function return_multiple()
    return 42, 43, 44
end

# ╔═╡ 64cd1682-371b-11eb-2097-a9267c1e63dc
p, q, r = return_multiple()

# ╔═╡ 75e813a2-371b-11eb-21d8-3d2bf2e80cac
("$p $r $q")

# ╔═╡ 827ebb9a-371b-11eb-006b-6107854f0a9c
md" ## Splatting
It is possible to “unpack” a tuple and pass its arguments to a function with the following syntax:
"

# ╔═╡ 04ea7b0a-371c-11eb-2827-7d192535d804
function splat_me(s, t, u)
    return s*t*u
end

# ╔═╡ 0e3a8240-371c-11eb-1720-7f7de97d90aa
tuple2 = (1,2,3)

# ╔═╡ 266b7540-371c-11eb-06f2-553a12a3314b
splat_me(tuple2...)

# ╔═╡ 4417da34-371c-11eb-155e-77b5eacbfc88
md" So the `...` after a tuple will unpack it! This is useful but addictive, use it only if needed as it is better for clarity (and to avoid multiple dispatch errors) to call a function with its single parameters.
"

# ╔═╡ 55d5c8a8-371c-11eb-00cf-91fa36d66a01
md" ## Named Tuples
Named tuples are like tuples but with a name identifier for a single value, for example:
"

# ╔═╡ 68bde3ec-371c-11eb-0aa7-8161a51fdbf5
namedTuple1 = (a = 1, b = "hello")

# ╔═╡ 6ba908ca-371c-11eb-1126-cfcbeb01cb0d
namedTuple1[:a]

# ╔═╡ 9b678d86-371c-11eb-0b34-2917771de661
md" or in alternative:
"

# ╔═╡ 83a9ac22-371c-11eb-09d9-c75f570c9b31
namedTuple2 = NamedTuple{(:a, :b)}((2,"hello2"))

# ╔═╡ 89ade9bc-371c-11eb-21aa-8f0831a9fff5
namedTuple2[:b]

# ╔═╡ 970399ea-371c-11eb-0cee-b54e5f3f5498
md"## Dictionaries
A dictionary is a collection of keys and values. They are unordered (which means that the order of the keys is random) and are really useful when you need to organise, for example, a dataset.

Let’s suppose we want to create an address book. A single entry should be able to store all the fundamental characteristics needed to identify a friend: the name of the contact, the phone number and the shoe size!
"

# ╔═╡ 706e79fc-371d-11eb-0894-ebb085103956
person1 = Dict("Name" => "Peppa", "Phone" => 123456789, "Shoe-size" => 40)

# ╔═╡ 7f46d19a-371d-11eb-2971-ef0eda3d5d18
person2 = Dict("Name" => "George", "Phone" => 123456789, "Shoe-size" => 36)

# ╔═╡ 9648198a-371d-11eb-3c9f-9338602af4af
md"I can even make a dictionary containing other dictionaries:
"

# ╔═╡ a25da512-371d-11eb-33ea-0f0740b919e3
addressBook = Dict("Peppa" => person1, "George" => person2)

# ╔═╡ ee8143ec-371d-11eb-3c9b-754027e67304
md"
We can add another friend to the addressBook , once it has been created, with the following syntax:
"

# ╔═╡ f94e1d2c-371d-11eb-0e5c-819b938a8d00
person3 = Dict("Name" => "Vittorio", "Phone" => 123456789, "Shoe-size" => 42)

# ╔═╡ 0560f292-371e-11eb-07de-c55444c0a3d5
addressBook["Vittorio"] = person3

# ╔═╡ 0ac0ae12-371e-11eb-3d9a-a7ef8aa4f5d9
addressBook

# ╔═╡ 19e817c2-371e-11eb-15ce-5bee1fd2f891
md" ## Conclusion
In this lesson we have learned how to operate on arrays, tuples and dictionaries. Those structures are the basics of in-memory storage for any data. Arrays are lightweight and useful solutions to pass blocks of data, so use them whenever needed! Contrarily to C++, Julia has a built in garbage collector, so you don’t have to care about freeing memory and deleting pointers, as Julia will take care of it!
"

# ╔═╡ Cell order:
# ╟─0edbe34a-3714-11eb-1350-af72689d8855
# ╟─4ef2e454-3716-11eb-3014-774db1fb57d0
# ╠═9833f658-3716-11eb-300e-d3b5d02b70e5
# ╠═b2cbd7d8-3716-11eb-15cc-03c3fec2c981
# ╠═b2d549da-3716-11eb-0cce-b703452d7f34
# ╟─bbe66bfa-3716-11eb-1929-ef53e271e354
# ╠═d5c03624-3716-11eb-38ed-3f395efc01f6
# ╟─e7768e7e-3716-11eb-366a-2f81166da844
# ╠═5268b19c-3717-11eb-3854-ddd093a61057
# ╟─6ced58c6-3717-11eb-1009-7909ebb6866b
# ╠═f2cb5918-3717-11eb-2b4b-77642fbeaf20
# ╟─ff24fcbc-3717-11eb-2ac3-6f18fde9484f
# ╠═21342288-3718-11eb-1b86-8d438637ac3a
# ╟─3fb54d2c-3718-11eb-305f-e7159394ac7b
# ╠═4c2a11c8-3718-11eb-0f59-8382e767058f
# ╟─538463b0-3718-11eb-03e7-171beded2a41
# ╟─6ed4e510-3718-11eb-3d00-9f60abc37333
# ╠═9a27ee68-3718-11eb-1c7f-33eaf6e7cb3f
# ╠═ab83ddde-3718-11eb-0847-9fdeaca6d470
# ╠═bbbc6720-3718-11eb-2c1e-978871157853
# ╟─d3ce8154-3718-11eb-08ff-7b29208c5099
# ╟─df969238-3718-11eb-2196-a7a68b5a69f2
# ╠═f82b9558-3718-11eb-0ef8-5fa509ba5b6f
# ╟─6ffa0c1a-3719-11eb-34d2-d7fc89a62e24
# ╠═7edb458c-3719-11eb-2bb3-4bd9173b530c
# ╠═8bae5826-3719-11eb-084e-fbf852773ffd
# ╟─a17c0086-3719-11eb-1ab0-2381decb5422
# ╠═c62a2930-3719-11eb-3c28-750300385ad8
# ╟─e07fe63a-3719-11eb-05eb-81512953fb4b
# ╟─2655af78-371a-11eb-331c-310db7a9ccb1
# ╠═3cdd35ea-371a-11eb-0f7f-135c3a766ecd
# ╠═67a3d7f2-371a-11eb-07cb-e3e233e2aab1
# ╠═67ad9e04-371a-11eb-36ad-9710ff0b0db3
# ╠═67bdc450-371a-11eb-08b2-55230f42c0e8
# ╟─b2873584-371a-11eb-328d-37a39540905c
# ╠═c17d2258-371a-11eb-25b5-933f413cd134
# ╠═d902ab08-371a-11eb-31ce-8d50b3388b2b
# ╠═d90d570e-371a-11eb-0e4b-075841a6a73a
# ╠═dbfd0a06-371a-11eb-1bf5-e9112227d8b0
# ╠═de4b38a0-371a-11eb-39b1-7575188e8063
# ╟─ee1a296c-371a-11eb-0036-3b8ad05c6804
# ╟─f6fa2e36-371a-11eb-36e7-a7de898b78c5
# ╠═0e12b2be-371b-11eb-0e56-99ec2e4a224d
# ╠═187b73f0-371b-11eb-34e6-d7a2aeb00b58
# ╟─1eb587ec-371b-11eb-2901-a1ebf0197aa7
# ╠═2bc6ba1e-371b-11eb-12d0-a3665fe9e1d0
# ╠═3985145c-371b-11eb-2dcc-5308e23dc5f0
# ╠═39a4874c-371b-11eb-164c-692fd38680a9
# ╟─58abede2-371b-11eb-24d3-2b16686937d5
# ╠═6205bad0-371b-11eb-1cfb-3f1ff43c255e
# ╠═64cd1682-371b-11eb-2097-a9267c1e63dc
# ╠═75e813a2-371b-11eb-21d8-3d2bf2e80cac
# ╟─827ebb9a-371b-11eb-006b-6107854f0a9c
# ╠═04ea7b0a-371c-11eb-2827-7d192535d804
# ╠═0e3a8240-371c-11eb-1720-7f7de97d90aa
# ╠═266b7540-371c-11eb-06f2-553a12a3314b
# ╟─4417da34-371c-11eb-155e-77b5eacbfc88
# ╟─55d5c8a8-371c-11eb-00cf-91fa36d66a01
# ╠═68bde3ec-371c-11eb-0aa7-8161a51fdbf5
# ╠═6ba908ca-371c-11eb-1126-cfcbeb01cb0d
# ╠═9b678d86-371c-11eb-0b34-2917771de661
# ╠═83a9ac22-371c-11eb-09d9-c75f570c9b31
# ╠═89ade9bc-371c-11eb-21aa-8f0831a9fff5
# ╟─970399ea-371c-11eb-0cee-b54e5f3f5498
# ╠═706e79fc-371d-11eb-0894-ebb085103956
# ╠═7f46d19a-371d-11eb-2971-ef0eda3d5d18
# ╟─9648198a-371d-11eb-3c9f-9338602af4af
# ╠═a25da512-371d-11eb-33ea-0f0740b919e3
# ╠═ee8143ec-371d-11eb-3c9b-754027e67304
# ╠═f94e1d2c-371d-11eb-0e5c-819b938a8d00
# ╠═0560f292-371e-11eb-07de-c55444c0a3d5
# ╠═0ac0ae12-371e-11eb-3d9a-a7ef8aa4f5d9
# ╟─19e817c2-371e-11eb-15ce-5bee1fd2f891


text:
import MPI
using Distributions

function simBer(t) return(rand(Bernoulli(t[1]))) end

function cal_optvalue(N, S, sigma, K, r, T, opt_type, n_sam, nrep)
	MPI.Init()
	comm = MPI.COMM_WORLD
	id = MPI.Comm_rank(comm)
	np = MPI.Comm_size(comm)
	# since we assume that the asset returns are log-normal with non-constant variance, u, d, and p are time-varying
	dt = T/N
	beta = 1/2*(exp(-r*dt)+exp((r+sigma^2)*dt))
	u = beta + sqrt(beta^2-1)
	d = 1/u
	p = (exp(r*dt)-d)/(u-d)
	probs = fill(p, N)
	u = fill(u, N)
	d = fill(d, N)
	M = log2(np)
	vals = zeros(nrep)
	probs_s = fill(p, convert(Int64, M))
	probs_n = fill(p, convert(Int64, (N-M)))
	probs_n = reshape(probs_n, length(probs_n), 1)
	if M%1 != 0
		error("Number of processes should be a multiple of 2")
		MPI.Finalize()
	end
	N_m = 2^(N-M)
	l_n = convert(Int64, N_m)

	if id == 0
		@printf("N: %d, np: %d\n", N, np)
		@printf("S: %f, sigma: %f, K: %f, r: %f, T: %f, %s\n", S, sigma, K, r, T, ifelse(opt_type==1, "Put", "Call"))
		tic()
	end

	v1 = 0.0
	#@printf("l_n: %d\n", l_n)
	#vt = zeros(Float64, l_n, 1)
	#p_vt = zeros(Float64, 1, l_n)
	p_vt = 0.0
	j = 0
	#for i in 0:(l_n-1)
	while j < nrep
	v1 = 0.0
	i = 0
	#srand(j)
	#smpl = sample(0:(l_n-1), n_sam, replace=true)
	while i < n_sam
		#node = smpl[i+1]
		#path = cat(2, integer_base_b(id, 2, M), integer_base_b(node, 2, N-M))
		path = cat(2, integer_base_b(id, 2, M), mapslices(simBer, probs_n, 2)')
		vt = calc_payoff(S, K, u, d, opt_type, path)
		v1 = v1 + vt
		i += 1
	end
	#v1 = exp(-r*T)*dot(vec(p_vt), vec(vt))
	#v1 = exp(-r*T)*v1*(1/n_sam)*calc_path_prob(integer_base_b(id, 2, M), probs_s)
	v1 = exp(-r*T)*v1*(1/n_sam)*(1/np)
	#@printf("id: %d, cnt: %d, v1: %f\n", id, j, v1)
	vals[j+1] = v1
	j += 1
	end
	MPI.Barrier(comm)

	#print(sprintf("P%d -- i=%d, v1: %f", id, i, v1))
	#@printf("Id: %d\n", id)
	#println(vals)
	reduced_v = MPI.Reduce(vals, MPI.SUM, 0, comm)
	if id == 0
		#@printf("Option Value: %f\n", reduced_v)
		#println(reduced_v)
		@printf("Mean of value: %f, Variance: %f\n", mean(reduced_v), var(reduced_v))
		@printf("Run time: %s\n", toc())
	end
	MPI.Finalize()
end

function calc_path_prob(path, probs)
	p1 = prod(probs[find(path -> path==1, path)])
	p2 = prod(1-probs[find(path -> path==0, path)])
	return p1*p2
end

function calc_payoff(S, K, u, d, opt_type, path)
	vt = 0.0
	if opt_type == 1
		vt = K - S*prod(u[find(path -> path==1, path)])*prod(d[find(path -> path==0, path)])
	else
		vt = S*prod(u[find(path -> path==1, path)])*prod(d[find(path -> path==0, path)]) - K
	end
	if vt < 0
		vt = 0.0
	end
	return vt
end

function integer_base_b(x, b=2, ndigits = 0)
	#xi = as.integer(x)
	if(typeof(x) != Int64)
		error("ERROR: ", x,  " is not an integer")
	end
	N = length(x)
	xMax = x
	if ndigits == 0 & x > 0
		ndigits = convert(Int64, round(floor(log(2, xMax))+1))
	end
	ndigits = convert(Int64, round(ndigits))
	#@printf("ndigits: %d\n", ndigits)
	#@printf("N: %d\n", N)
	#Base.b = array(NA, dim=c(N, ndigits))
	Base_b = zeros(Int64, 1, ndigits)

	if ndigits == 0
		return Base_b
	end
	for i in 1:ndigits
		Base_b[:,ndigits-i+1] = (x % b)
		#x = (x %/% b)
		x = div(x, b)
	end
	if N == 1
		return Base_b[1,: ]
	else
	 return Base_b
	end
end


text:
## High level

# Here is where we can add a default algorithm for computing sensitivities
# Based on problem information!
function DiffEqBase._concrete_solve_adjoint(prob,alg,sensealg::Nothing,u0,p,args...;kwargs...)
  default_sensealg = InterpolatingAdjoint()
  DiffEqBase._concrete_solve_adjoint(prob,alg,default_sensealg,u0,p,args...;kwargs...)
end

function DiffEqBase._concrete_solve_adjoint(prob,alg,
                                 sensealg::AbstractAdjointSensitivityAlgorithm,
                                 u0,p,args...;save_start=true,save_end=true,
                                 saveat = eltype(prob.tspan)[],
                                 save_idxs = nothing,
                                 kwargs...)
  _prob = remake(prob,u0=u0,p=p)

  # Force `save_start` and `save_end` in the forward pass This forces the
  # solver to do the backsolve all the way back to `u0` Since the start aliases
  # `_prob.u0`, this doesn't actually use more memory But it cleans up the
  # implementation and makes `save_start` and `save_end` arg safe.
  kwargs_fwd = NamedTuple{Base.diff_names(Base._nt_names(
  values(kwargs)), (:callback_adj,))}(values(kwargs))
  kwargs_adj = NamedTuple{Base.diff_names(Base._nt_names(values(kwargs)), (:callback_adj,:callback))}(values(kwargs))
  if haskey(kwargs, :callback_adj)
    kwargs_adj = merge(kwargs_adj, NamedTuple{(:callback,)}( [get(kwargs, :callback_adj, nothing)] ))
  end

  if ischeckpointing(sensealg)
    sol = solve(_prob,alg,args...;save_start=true,save_end=true,saveat=saveat,kwargs...)
  else
    sol = solve(_prob,alg,args...;save_start=true,save_end=true,kwargs...)
  end

  if saveat isa Number
    if _prob.tspan[2] > _prob.tspan[1]
      ts = _prob.tspan[1]:abs(saveat):_prob.tspan[2]
    else
      ts = _prob.tspan[2]:abs(saveat):_prob.tspan[1]
    end
    _out = sol(ts)
    out = save_idxs === nothing ? _out : DiffEqArray([x[save_idxs] for x in _out.u],ts)
    only_end = length(ts) == 1 && ts[1] == _prob.tspan[2]
  elseif isempty(saveat)
    no_start = !save_start
    no_end = !save_end
    sol_idxs = 1:length(sol)
    no_start && (sol_idxs = sol_idxs[2:end])
    no_end && (sol_idxs = sol_idxs[1:end-1])
    only_end = length(sol_idxs) <= 1
    _u = sol.u[sol_idxs]
    u = save_idxs === nothing ? _u : [x[save_idxs] for x in _u]
    ts = sol.t[sol_idxs]
    out = DiffEqArray(u,ts)
  else
    _saveat = saveat isa Array ? sort(saveat) : saveat # for minibatching
    ts = _saveat
    _out = sol(ts)
    out = save_idxs === nothing ? _out : DiffEqArray([x[save_idxs] for x in _out.u],ts)
    only_end = length(ts) == 1 && ts[1] == _prob.tspan[2]
  end

  function adjoint_sensitivity_backpass(Δ)
    function df(_out, u, p, t, i)
      if only_end
        _out[:] .= -vec(Δ)
      else
        _out[:] .= -adapt(typeof(u0),reshape(Δ, prod(size(Δ)[1:end-1]), size(Δ)[end])[:, i])
      end
    end

    du0, dp = adjoint_sensitivities(sol,alg,args...,df,ts; sensealg=sensealg,
                                    kwargs_adj...)

    du0 = reshape(du0,size(u0))
    dp = reshape(dp',size(p))

    (nothing,nothing,du0,dp,ntuple(_->nothing, length(args))...)
  end
  out, adjoint_sensitivity_backpass
end

# Prefer this route since it works better with callback AD
function DiffEqBase._concrete_solve_adjoint(prob,alg,sensealg::AbstractForwardSensitivityAlgorithm,
                                 u0,p,args...;
                                 kwargs...)
   _prob = ODEForwardSensitivityProblem(prob.f,u0,prob.tspan,p,sensealg)
   sol = solve(_prob,alg,args...;kwargs...)
   u,du = extract_local_sensitivities(sol, Val(true))
   function forward_sensitivity_backpass(Δ)
     adj = sum(eachindex(du)) do i
       J = du[i]
       v = @view Δ[:, i]
       J'v
     end
     (nothing,nothing,nothing,adj,ntuple(_->nothing, length(args))...)
   end
   DiffEqArray([u[:,i] for i in 1:size(u,2)],sol.t),forward_sensitivity_backpass
end

function DiffEqBase._concrete_solve_forward(prob,alg,
                                 sensealg::AbstractForwardSensitivityAlgorithm,
                                 u0,p,args...;
                                 kwargs...)
   _prob = ODEForwardSensitivityProblem(prob.f,u0,prob.tspan,p,sensealg)
   sol = solve(_prob,args...;kwargs...)
   u,du = extract_local_sensitivities(sol,Val(true))
   function _concrete_solve_pushforward(Δself, ::Nothing, ::Nothing, x3, Δp, args...)
     x3 !== nothing && error("Pushforward currently requires no u0 derivatives")
     du * Δp
   end
   DiffEqArray([u[:,i] for i in 1:size(u,2)],sol.t),_concrete_solve_pushforward
end

# Generic Fallback for ForwardDiff
function DiffEqBase._concrete_solve_adjoint(prob,alg,
                                 sensealg::ForwardDiffSensitivity,
                                 u0,p,args...;saveat=eltype(prob.tspan)[],
                                 save_idxs = nothing,
                                 kwargs...)
  save_idxs !== nothing && error("save_idxs is currently incompatible with ForwardDiffSensitivity")
  MyTag = typeof(prob.f)
  pdual = seed_duals(p,MyTag)
  u0dual = convert.(eltype(pdual),u0)

  if (convert_tspan(sensealg) === nothing && (
        (haskey(kwargs,:callback) && has_continuous_callback(kwargs.callback)) ||
        (haskey(prob.kwargs,:callback) && has_continuous_callback(prob.kwargs.callback))
        )) || (convert_tspan(sensealg) !== nothing && convert_tspan(alg))

    tspandual = convert.(eltype(pdual),prob.tspan)
  else
    tspandual = prob.tspan
  end

  _prob = remake(prob,u0=u0dual,p=pdual,tspan=tspandual)

  if saveat isa Number
    _saveat = prob.tspan[1]:saveat:prob.tspan[2]
  else
    _saveat = saveat
  end

  sol = solve(_prob,alg,args...;saveat=_saveat,kwargs...)

  u,du = extract_local_sensitivities(sol, sensealg, Val(true))
  function forward_sensitivity_backpass(Δ)
    adj = sum(eachindex(du)) do i
      J = du[i]
      v = @view Δ[:, i]
      J'v
    end
    (nothing,nothing,nothing,adj,ntuple(_->nothing, length(args))...)
  end
  DiffEqArray([u[:,i] for i in 1:size(u,2)],sol.t),forward_sensitivity_backpass
end

function DiffEqBase._concrete_solve_adjoint(prob,alg,sensealg::ZygoteAdjoint,
                                 u0,p,args...;kwargs...)
    Zygote.pullback((u0,p)->_concrete_solve(prob,alg,u0,p,args...;kwargs...),u0,p)
end

function DiffEqBase._concrete_solve_adjoint(prob,alg,sensealg::TrackerAdjoint,
                                            u0,p,args...;kwargs...)

  t = eltype(prob.tspan)[]
  u = typeof(u0)[]

  function tracker_adjoint_forwardpass(_u0,_p)
    if DiffEqBase.isinplace(prob)
      # use Array{TrackedReal} for mutation to work
      # Recurse to all Array{TrackedArray}
      _prob = remake(prob,u0=map(identity,_u0),p=_p)
    else
      # use TrackedArray for efficiency of the tape
      _f(args...) = Tracker.collect(prob.f(args...))
      if prob isa SDEProblem
        _g(args...) = Tracker.collect(prob.g(args...))
        _prob = remake(prob,f=DiffEqBase.parameterless_type(prob.f)(_f,_g),u0=_u0,p=_p)
      else
        _prob = remake(prob,f=DiffEqBase.parameterless_type(prob.f)(_f),u0=_u0,p=_p)
      end
    end
    sol = solve(_prob,alg,args...;kwargs...)
    t = sol.t
    if DiffEqBase.isinplace(prob)
      u = map.(Tracker.data,sol.u)
    else
      u = map(Tracker.data,sol.u)
    end

    if typeof(sol.u[1]) <: Array
      return adapt(typeof(u0),sol)
    else
      tmp = vec(sol.u[1])
      for i in 2:length(sol.u)
        tmp = hcat(tmp,vec(sol.u[i]))
      end
      return reshape(tmp,size(sol.u[1])...,length(sol.u))
    end
    #adapt(typeof(u0),arr)
  end

  sol,pullback = Tracker.forward(tracker_adjoint_forwardpass,u0,p)
  function tracker_adjoint_backpass(ybar)
    u0bar, pbar = pullback(ybar)
    _u0bar = u0bar isa Tracker.TrackedArray ? Tracker.data(u0bar) : Tracker.data.(u0bar)
    (nothing,nothing,_u0bar,Tracker.data(pbar),ntuple(_->nothing, length(args))...)
  end
  DiffEqArray(u,t),tracker_adjoint_backpass
end


text:
__precompile__(true)

"""
Main module for `RandomNumbers.jl` -- a random number generator package for Julia Language.

This module exports two types and four submodules:

- [`AbstractRNG`](@ref)
- [`WrappedRNG`](@ref)
- [`PCG`](@ref)
- [`MersenneTwisters`](@ref)
- [`Random123`](random123.md#Random123.Random123)
- [`Xorshifts`](@ref)
"""
module RandomNumbers

    export AbstractRNG
    export WrappedRNG
    export output_type, seed_type
    export PCG, MersenneTwisters, Random123, Xorshifts

    include("common.jl")
    include("utils.jl")

    include("wrapped_rng.jl")

    include(joinpath("PCG", "PCG.jl"))
    include(joinpath("MersenneTwisters", "MersenneTwisters.jl"))
    include(joinpath("Xorshifts", "Xorshifts.jl"))

    import Requires
    function __init__()
        # The code of Random123 has been moved to Random123.jl
        Requires.@require Random123="74087812-796a-5b5d-8853-05524746bad3" import .Random123
    end

end


text:
using StatsBase
using SparseArrays, Test

# rle & inverse_rle

z = [1, 1, 2, 2, 2, 3, 1, 2, 2, 3, 3, 3, 3]
(vals, lens) = rle(z)
@test vals == [1, 2, 3, 1, 2, 3]
@test lens == [2, 3, 1, 1, 2, 4]
@test inverse_rle(vals, lens) == z

z = [true, true, false, false, true, false, true, true, true]
vals, lens = rle(z)
@test vals == [true, false, true, false, true]
@test lens == [2, 2, 1, 1, 3]
@test inverse_rle(vals, lens) == z

# levelsmap
a = [1, 1, 2, 2, 2, 3, 1, 2, 2, 3, 3, 3, 3, 2]
b = [true, false, false, true, false, true, true, false]

@test levelsmap(a) == Dict(2=>2, 3=>3, 1=>1)
@test levelsmap(b) == Dict(false=>2, true=>1)

# indicatormat

II = [false true  false false false;
      true  false false false true;
      false false true  true  false]

x = [2, 1, 3, 3, 2]
@test indicatormat(x, 3) == II
@test Matrix(indicatormat(x, 3; sparse=true)) == II

x = ["b", "a", "c", "c", "b"]
@test indicatormat(x) == II
@test Matrix(indicatormat(x; sparse=true)) == II

io = IOBuffer()
describe(io, collect(1:10))
@test String(take!(io)) == """
                           Summary Stats:
                           Length:         10
                           Missing Count:  0
                           Mean:           5.500000
                           Minimum:        1.000000
                           1st Quartile:   3.250000
                           Median:         5.500000
                           3rd Quartile:   7.750000
                           Maximum:        10.000000
                           Type:           $Int
                           """

describe(io, Union{Float32,Missing}[1.0, 4.5, missing, missing, 33.1])
@test String(take!(io)) == """
                           Summary Stats:
                           Length:         5
                           Missing Count:  2
                           (All summary stats are missing)
                           Type:           $(Union{Float32,Missing})
                           """

describe(io, Float64[])
@test String(take!(io)) == """
                           Summary Stats:
                           Length:         0
                           Type:           Float64
                           """

describe(io, fill("s", 3))
@test String(take!(io)) == """
                           Summary Stats:
                           Length:         3
                           Type:           String
                           Number Unique:  1
                           """


text:
"""
    function frequency_ranges(
      f_min,
      f_max,
      subnet::Int64,
      directory::String,
      objective::String,
      x_axis::Array,
      y_axis::Array;
      gen_areas=Int64[],
      area_transfer=Int64[],
      gen_zones=[],
      zone_transfer=[],
      plot_vert_line::Tuple=([],""),
      plot_horiz_line::Tuple=([],""),
      xlimits::Array{Any,1}=[],
      ylimits::Array{Any,1}=[],
      output_plot_label::Tuple{String,String}=("",""),
      scopf::Bool=false,
      contingency::Int64=0,
      k_cond=[],
      k_ins=[],
      scale_load=1.0,
      scale_areas=Int64[],
      no_converter_loss=false,
      output_location_base="",
      output_results_folder=""
    )

Models and solves an OPF with frequency in specified ranges between `f_min` and `f_max`.

# Arguments
- `f_min`: lower bounds on frequency, one for each point in the frequency sweep
- `f_max`: upper bounds on frequency, one for each point in the frequency sweep. Must have the same length as `f_min`.
- `subnet::Int64`: subnetwork for which the frequency bounds are applied
- `folder::String`: the directory containing all subnetwork data, *subnetworks.csv*, and *interfaces.csv*
- `objective::String`: the objective function to use, from the following:
   - "mincost": minimize generation cost
   - "areagen": minimize generation in the areas specified in `gen_areas`
   - "zonegen": minimize generation in the zones specified in `gen_zones`
   - "minredispatch": minimize the change in generator dispatch from the initial values defined in the network data
- `x_axis::Array`: Array of Tuples identifying the x axis series for which plots should be generated over the points in the frequency sweep. A separate folder of plots is generated for each Tuple in the array. The series can be specified in the Tuple in one of three ways:
   - **results dictionary values:** A two-element Tuple, where the first element is a String matching a key in the results dictionary output from `multifrequency_opf` and the second element is an Int specifying a subnetwork. This plots the values of this key and subnetwork entry on the x axis.
   - **network data values:** A Tuple with elements corresponding to keys at each level of the network data dictionary, identifying any network variable value. This plots the values of the specified network variable on the x axis. Any key in the Tuple may be an Array, in which case a separate plot is generated for each key. For example, to generate four plots, the active and reactive power at the origin ("f") bus and destination ("t") bus for branch 1 in subnetwork 2, use the Tuple `("sn",2,"branch",1,["pt","pf","qt","qf"])`
   - **custom values:** A two-element Tuple, where the first element is a String not matching any keys in the results dictionary and the second element is an Array. This plots the values in the Array on the x axis with the label in the String.
- `y_axis::Array`: Array of Tuples identifying the y axis series for which plots should be generated over the points in the frequency sweep. A separate folder of plots is generated for each Tuple in the array. These are specified in the same way as `x_axis`.
- `gen_areas`: all areas in which generation should be minimized if `obj=="areagen"`
- `area_transfer`: two areas with power transfer between them that should be saved and plotted. Results for P, Q, S, and loss are saved for power transfer between the two areas. Must have exactly two elements.
- `gen_zones`: all zones in which generation should be minimized if `obj=="zonegen"`
- `zone_transfer`: two zones with power transfer between them that should be saved and plotted. Results for P, Q, S, and loss are saved for power transfer between the two zones. Must have exactly two elements.
- `plot_vert_line::Tuple`: x values of vertical lines to overlay on the plot. The first element is a scalar or Array specifying one or more x values to plot, and the second element is a String or Array of Strings specifying the label or labels. Default ([],"") does not add any lines to the plot.
- `plot_horiz_line::Tuple`: y values of horizontal lines to overlay on the plot. The first element is a scalar or Array specifying one or more y values to plot, and the second element is a String or Array of Strings specifying the label or labels. Default ([],"") does not add any lines to the plot.
- `xlimits::Array{Any,1}`: Array of two values specifying the min and max x axis limits to apply to the plots, overriding any other limits. Default [] does not change the plot.
- `ylimits::Array{Any,1}`: Array of two values specifying the min and max y axis limits to apply to the plots, overriding any other limits. Default [] does not change the plot.
- `output_plot_label::Tuple{String,String}`: specifies the plot to pass to the output. The first element must match the x axis label, and the second must match the y axis label.
- `scopf::Bool`: if true, model and solve the N-1 security constrained OPF for each network. Each network folder must contain a contingency specification file (_*.con_) for each subnetwork. Default false.
- `contingency::Tuple`: indices of the contingency to plot. The precontingency index is (0,). Default (0,).
- `k_cond`: conductor utilization parameter for HVDC. Only used when f==0. Default [].
- `k_ins`: insulation factor parameter for HVDC. Only used when f==0. Default [].
- `scale_load`: factor for scaling the load in the frequency sweep. Default 1.0.
- `scale_areas`: array of integer area indices for which the load scaling factor `scale_load` should be applied. Applies to all areas if this array is empty. Default Int64[].
- `no_converter_loss`: override all converter loss parameters specified in the data and replace them with the the lossless converter model.
- `output_location_base`: location in which to save the results and plots. If not specified, a folder called `results` will be created in the folder one level above the data folder.
- `output_results_folder`: specific folder in which to save the results, one level below `output_location_base`, if specified.
"""
function frequency_ranges(
      f_min,
      f_max,
      subnet::Int64,
      folder::String,
      objective::String,
      x_axis::Array,
      y_axis::Array;
      gen_areas=Int64[],
      area_transfer=Int64[],
      gen_zones=[],
      zone_transfer=[],
      plot_vert_line::Tuple=([],""),
      plot_horiz_line::Tuple=([],""),
      xlimits::Array{Any,1}=[],
      ylimits::Array{Any,1}=[],
      output_plot_label::Tuple{String,String}=("",""),
      scopf::Bool=false,
      contingency::Int64=0,
      k_cond=[],
      k_ins=[],
      scale_load=1.0,
      scale_areas=Int64[],
      no_converter_loss=false,
      output_location_base="",
      output_results_folder=""
   )
   mn_data = read_sn_data(folder, no_converter_loss=no_converter_loss)

   folder = abspath(folder)
   if length(output_location_base) > 0
      output_folder = output_location_base
   else
      folder_split = splitpath(folder)
      toplevels = folder_split[1:end-3]
      output_folder = joinpath(toplevels...,"results/$(folder_split[end-1])/$(folder_split[end])")
   end
   output_folder = joinpath(output_folder, output_results_folder)
   if !isdir(output_folder)
      mkpath(output_folder)
   end
   println("output location: $output_folder")

   (results_dict, output_plot) = frequency_ranges(
      f_min,
      f_max,
      subnet,
      mn_data,
      output_folder,
      objective,
      x_axis,
      y_axis;
      gen_areas=gen_areas,
      area_transfer=area_transfer,
      gen_zones=gen_zones,
      zone_transfer=zone_transfer,
      plot_vert_line=plot_vert_line,
      plot_horiz_line=plot_horiz_line,
      xlimits=xlimits,
      ylimits=ylimits,
      output_plot_label=output_plot_label,
      scopf=scopf,
      contingency=contingency,
      k_cond=k_cond,
      k_ins=k_ins,
      scale_load=scale_load,
      scale_areas=scale_areas,
      no_converter_loss=no_converter_loss,
      output_results_folder=""
   )
   return (results_dict, output_plot)
end

function frequency_ranges(
      f_min,
      f_max,
      subnet::Int64,
      mn_data::Dict,
      output_folder::String,
      objective::String,
      x_axis::Array,
      y_axis::Array;
      gen_areas=Int64[],
      area_transfer=Int64[],
      gen_zones=[],
      zone_transfer=[],
      plot_vert_line::Tuple=([],""),
      plot_horiz_line::Tuple=([],""),
      xlimits::Array{Any,1}=[],
      ylimits::Array{Any,1}=[],
      output_plot_label::Tuple{String,String}=("",""),
      scopf::Bool=false,
      contingency::Int64=0,
      k_cond=[],
      k_ins=[],
      scale_load=1.0,
      scale_areas=Int64[],
      no_converter_loss=false,
      output_results_folder=""
   )

   suffix=""
   if scopf
      params = ([["sn",contingency,"sn",subnet,"f_min"],["sn",contingency,"sn",subnet,"f_max"]],[f_min, f_max])
   elseif scale_load==1
      params = ([["sn",subnet,"f_min"],["sn",subnet,"f_max"]],[f_min, f_max])
   else
      println("Scaling load by $scale_load.")
      suffix = "_scaled_x$(string(scale_load))"
      scale_all_areas = false
      if length(scale_areas) == 0
         scale_all_areas = true
      end
      # get load pd and qd from areas of interest and put them in dict_filt
      dict_filt = Dict()
      for (subnet_idx,mn_subnet) in mn_data["sn"]
         for (load_idx,load) in mn_subnet["load"]
            if (mn_subnet["bus"]["$(load["load_bus"])"]["area"] in scale_areas) || (scale_all_areas)
               pd_keys = ["sn",subnet_idx,"load",load_idx,"pd"]
               pd = load["pd"]
               set_nested!(dict_filt, pd_keys, pd)
               qd_keys = ["sn",subnet_idx,"load",load_idx,"qd"]
               qd = load["qd"]
               set_nested!(dict_filt, qd_keys, qd)
            end
         end
      end

      # get the nested keys values
      (load_keys,load_values) = traverse_nested(dict_filt)

      # apply scaling
      # println("load_keys: $load_keys")
      # println("load_values: $load_values")
      # values_scaled = collect(eachrow(load_values .* scale_load'))
      values_scaled = Array{Any}(undef, length(load_values))
      for (idx,val) in enumerate(load_values)
         values_scaled[idx] = val * ones(length(f_min)) * scale_load
      end

      params = (append!([["sn",subnet,"f_min"],["sn",subnet,"f_max"]],load_keys),append!(Any[f_min, f_max], values_scaled))
   end

   if (length(k_cond)) > 0 && (length(k_ins) > 0)
      if scopf
         dc_params = ([["sn",contingency,"sn",subnet,"k_ins"],["sn",contingency,"sn",subnet,"k_cond"],["sn",contingency,"sn",subnet,"f_max"]],[k_ins, k_cond, zeros(length(k_ins))])
      else
         dc_params = ([["sn",subnet,"k_ins"],["sn",subnet,"k_cond"],["sn",subnet,"f_max"]],[k_ins, k_cond, zeros(length(k_ins))])
      end
   else
      dc_params = ()
   end

   output_folder = joinpath(output_folder, output_results_folder)
   if !isdir(output_folder)
      mkpath(output_folder)
   end
   println("output location: $output_folder")

   # println("params:")
   # println(params)

   (results_dict, output_plot) = run_multiple_params(
         mn_data,
         output_folder,
         objective,
         x_axis,
         y_axis,
         params,
         gen_areas,
         area_transfer,
         gen_zones,
         zone_transfer,
         plot_vert_line,
         plot_horiz_line,
         xlimits,
         ylimits,
         output_plot_label,
         scopf,
         contingency,
         dc_params;
         output_results_folder=""
      )
   # for config in hvdc_config
   #
   # end
   # outprint = output_plot != nothing ? output_plot : "nothing"
   # println("\n***IN frequency_ranges, output_plot = $outprint")
   return (results_dict, output_plot)
end


text:
include("DeforModelRedModule.jl")
include("MatDeforModule.jl")
include("MatDeforLinearElasticModule.jl")
include("MatDeforElastIsoModule.jl")
include("MatDeforElastOrthoModule.jl")
include("FEMMDeforLinearBaseModule.jl")
include("FEMMDeforLinearModule.jl")
include("FEMMDeforWinklerModule.jl")
include("FEMMDeforSurfaceDampingModule.jl")
include("FEMMDeforLinearIMModule.jl")
include("FEMMDeforLinearMSModule.jl")
include("FEMMDeforLinearNICEModule.jl")
include("FEMMDeforLinearESNICEModule.jl")
include("AlgoDeforLinearModule.jl")