text:
#///////////////////////////////////////
#// File Name: scene_loader.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2021/01/15
#// Description: Scene loaders to load ado positions
#///////////////////////////////////////

using PyCall

abstract type SceneLoader end

# Trajectron Scene Loader
struct TrajectronSceneParameter <: Parameter
    conf_file_name::String     # json config file name.
    test_data_name::String     # Specify "eth_test.pkl", "hotel_test.pkl", or "univ_test.pkl"
    test_scene_id::Int64       # Test Scene ID (Trajectron) (starting from 0)
    start_time_idx::Int64      # Simulation start time index (Trajectron) (starting from 0, but you should start from 2)
    incl_robot_node::Bool      # Whether or not the robot node is included in the scene.
end

mutable struct TrajectronSceneLoader <: SceneLoader
    param::TrajectronSceneParameter
    model_dir::String              # model directory to be loaded in Trajectron Predictor
    eval_env::PyObject             # evaluation Environment object (Trajectron)
    eval_scene::PyObject           # evaluation Scene object (Trajectron) (i.e. entire scene sequence)
    online_scene::PyObject         # online Scene object (Trajectron) (i.e. current scene)
    state_def::PyObject                # state definitions of agents in the scene
    dto::Float64
    curr_time_idx::Int64           # Simulation current time index (Trajectron)
    max_time_idx::Int64            # Simulation max time index (Trajectron)
end

function TrajectronSceneLoader(param::TrajectronSceneParameter;
                               verbose=true,
                               ado_id_removed=nothing)
   # note removed ado agent still exists in eval_scene and
   # will be returned by fetch_ado_positions!. You'll have
   # to manually delete the ado agent from the outputs_dict.

    @assert in(param.test_data_name, ["eth", "hotel", "univ"] .* "_test.pkl");
    file_dir = @__DIR__;
    model_dir = normpath(joinpath(file_dir, "..",
                                  "Trajectron-plus-plus", "experiments", "pedestrians",
                                  "models"));
    if param.test_data_name == "eth_test.pkl"
        if param.incl_robot_node
            model_dir = joinpath(model_dir, "models_21_Apr_2020_20_15_41_eth_ar3_rob");
        else
            model_dir = joinpath(model_dir, "eth_attention_radius_3");
        end
    elseif param.test_data_name == "hotel_test.pkl"
        if param.incl_robot_node
            model_dir = joinpath(model_dir, "models_21_Apr_2020_20_15_59_hotel_ar3_rob")
        else
            model_dir = joinpath(model_dir, "hotel_attention_radius_3");
        end
    elseif param.test_data_name == "univ_test.pkl"
        if param.incl_robot_node
            model_dir = joinpath(model_dir, "models_21_Apr_2020_20_17_14_univ_ar3_rob");
        else
            model_dir = joinpath(model_dir, "univ_attention_radius_3");
        end
    end
    hyperparams = pycall(py"load_hyperparams", PyObject,
                         model_dir, param.conf_file_name);
    maximum_history_length = pycall(hyperparams.get, PyObject, "maximum_history_length")
    state_def = pycall(hyperparams.get, PyObject, "state")

    data_dir = normpath(joinpath(@__DIR__, "..",
                                "Trajectron-plus-plus", "experiments", "processed"));

    eval_env, eval_scene, online_scene, log =
        pycall(py"load_online_scene", NTuple{4, PyObject},
               data_dir, param.test_data_name, maximum_history_length, state_def,
               param.test_scene_id, param.incl_robot_node, param.start_time_idx,
               ado_id_removed)

    dto = eval_scene.dt;

    curr_time_idx = param.start_time_idx;
    max_time_idx = pybuiltin("int")(eval_scene.timesteps) - 1;
    @assert curr_time_idx < max_time_idx "curr_time_idx must be smaller than $(max_time_idx)!"

    if verbose
        for string in log
            println(string)
        end
        println("Looking at the $(param.test_data_name) sequence, data_id $(param.test_scene_id), start_idx $(curr_time_idx)");
    end
    return TrajectronSceneLoader(param, model_dir, eval_env, eval_scene, online_scene,
                                 state_def, dto, curr_time_idx, max_time_idx)
end

function fetch_ado_positions!(loader::TrajectronSceneLoader; return_full_state=false);
    input_dict = loader.eval_scene.get_clipped_input_dict(loader.curr_time_idx,
                                                          loader.state_def)
    input_dict_returned = Dict{PyObject, Vector{Float64}}();
    for key in keys(input_dict)
        if key != loader.online_scene.robot
            if return_full_state
                # full ado states
                input_dict_returned[key] = input_dict[key][:];
            else
                # positions only
                input_dict_returned[key] = [input_dict[key][1],
                                            input_dict[key][2]]
            end
        end
    end
    loader.curr_time_idx += 1; # increment curr_time_idx by 1.
    return input_dict_returned
end

function reduce_to_positions(input_dict::Dict{T, Vector{Float64}}) where T <: Union{PyObject, String}
    input_dict_returned = Dict{T, Vector{Float64}}();
    for key in keys(input_dict)
        # positions only
        input_dict_returned[key] = [input_dict[key][1],
                                    input_dict[key][2]]
    end
    return input_dict_returned
end

function convert_nodes_to_str(input_dict::Dict)
    if typeof(input_dict) == Dict{String, Vector{Float64}}
        return input_dict
    end
    input_dict_returned = Dict{String, Vector{Float64}}();
    for key in keys(input_dict)
        input_dict_returned[pybuiltin("str")(key)] = input_dict[key][:];
    end
    return input_dict_returned
end

function get_trajectory_for_ado(scene_loader::TrajectronSceneLoader,
                                init_time::Time,
                                ado_id::String,
                                sim_horizon::Float64)
    start_time_idx = scene_loader.param.start_time_idx;
    end_time_idx = start_time_idx + Int64(round(sim_horizon/scene_loader.dto, digits=5));

    ado_trajectory = Trajectory2D();

    time = init_time;
    for ii = start_time_idx:1:end_time_idx
        input_dict_py = scene_loader.eval_scene.get_clipped_input_dict(ii, scene_loader.state_def)
        input_dict = convert_nodes_to_str(input_dict_py)
        pos_dict = input_dict;
        for key in keys(pos_dict)
            pos_dict[key] = pos_dict[key][1:2];
        end
        if in(ado_id, keys(input_dict))
            ado_trajectory[time] = pos_dict[ado_id];
        end
        time += Duration(scene_loader.dto);
    end
    if isempty(ado_trajectory)
        @error(ArgumentError("$(ado_id) does not exist in the scenes!"))
    end
    if !all(diff(to_sec.(collect(keys(ado_trajectory)))) .≈ scene_loader.dto)
        @warn("Ado trajectory is discontinuous!")
    end
    return ado_trajectory
end

# Synthetic Scene Loader
struct SyntheticSceneParameter <: Parameter
    rng::AbstractRNG
end

mutable struct SyntheticSceneLoader <: SceneLoader
    param::SyntheticSceneParameter
    ado_pos_dict::Dict{String, Vector{Float64}}
end

function fetch_ado_positions!(loader::SyntheticSceneLoader,
                              prediction_dict::Dict{String, Array{Float64, 3}})
    @assert keys(loader.ado_pos_dict) == keys(prediction_dict)
    for key in keys(prediction_dict)
        # sample next position independently per ado agent from prediction_dict
        num_samples = size(prediction_dict[key], 1);
        sample_idx = rand(loader.param.rng, 1:num_samples)
        loader.ado_pos_dict[key] = prediction_dict[key][sample_idx, 1, :];
    end
    return deepcopy(loader.ado_pos_dict)
end


text:
struct GaussProcess{A<:Array, D, F0, F1, F2, F3}
    γ::A
    β::D
    μ::F0
    ∇ₓμ::F1
    Σ::F2
    ∇ₓΣ::F3
end

function GaussProcess(γ::Array, β::Real, Σ, ∇ₓΣ)
    μ(x) = zeros(size(x, 2))
    ∇ₓμ(x) = zeros(size(x))
    return GaussProcess(γ, β, μ, ∇ₓμ, Σ, ∇ₓΣ)
end

function negloglik(GP::GaussProcess, X̄::Matrix, ȳ::Vector)
    K = GP.Σ(relu.(GP.γ), X̄, X̄)
    K += K'
    K += (GP.β^2)*I
    μ = GP.μ(X̄)
    likelihood = MvNormal(μ, K)
    nll = -1*logpdf(likelihood, ȳ)
    return nll
end

function predict(GP::GaussProcess, X̄::Matrix, ȳ::Vector, x::Matrix)
    K = GP.Σ(GP.γ, X̄, X̄) + (GP.β^2)*I
    k = GP.Σ(GP.γ, x, X̄)
    reduced_ȳ = ȳ - GP.μ(X̄)

    Cho_K = cholesky(K)
    μ = GP.μ(x) + k*(Cho_K\reduced_ȳ)
    Σp = GP.Σ(GP.γ, x, x) - k*(Cho_K\(k'))

    σ = sqrt.(diag(Σp))
    return μ, σ
end

function ∇ₓpredict(GP::GaussProcess, X̄::Matrix, ȳ::Vector, x::Matrix)
    F, N₁, N₂ = size(X̄, 1), size(X̄, 2), size(x, 2)
    K = GP.Σ(relu.(GP.γ), X̄, X̄)
    K += K'
    K += (GP.β^2)*I
    ∇ₓk = GP.∇ₓΣ(relu.(GP.γ), x, X̄)
    ∇ₓμ = convert.(eltype(GP.γ), GP.∇ₓμ(x))
    reduced_ȳ = ȳ - GP.μ(X̄)

    Cho_K = cholesky(K)
    α = Cho_K\reduced_ȳ
    for j in 1:N₁
        for i in 1:N₂
            for l in 1:F
                ∇ₓμ[l,i] += ∇ₓk[l,i,j]*α[j]
            end
        end
    end
    return ∇ₓμ
end

ext: jl

max_stars_repo_name: HamletWantToCode/DiffGPFlux

max_issues_repo_name: HamletWantToCode/DiffGPFlux


text:
using Turing
using Distributions
using AdvancedMH
using AdvancedHMC
using LinearAlgebra: I, dot
using Random
using Test

using AnnealedIS

function geomspace(start, stop, length)
    logstart = log10(start)
    logstop = log10(stop)
    points = 10 .^ range(logstart, logstop; length=length)
    points[1] = start
    points[end] = stop
    return points
end

@testset "AnnealedIS.jl" begin
    rng = MersenneTwister(42)
    Random.seed!(42)

    @testset "AnnealedISSampler" begin
        D = 5
        N = 2

        prior_sampling(rng) = rand(rng, MvNormal(D, 1.0))
        prior_density(params) = logpdf(MvNormal(D, 1.0), params)
        joint_density(params) = logpdf(MvNormal(D, 0.5), params)

        ais = AnnealedISSampler(
            prior_sampling, prior_density, joint_density, N)

        @test ais.betas == [0.0, 0.5, 1.0]
        x = ais.prior_sampling(rng)
        y = rand(rng, ais.transition_kernels[1])
        @test typeof(x) == typeof(y)

        samples = [ais.prior_sampling(rng) for _ in 1:N]

        @test AnnealedIS.logdensity(ais, 1, samples[1]) == prior_density(samples[1])
        @test AnnealedIS.logdensity(ais, 3, samples[1]) == joint_density(samples[1])
        @test AnnealedIS.logdensity(ais, 2, samples[1]) == (0.5*prior_density(samples[1]) + 0.5*joint_density(samples[1]))

        for i in 1:N
            dratio = AnnealedIS.logdensity(ais, i+1, samples[1]) - AnnealedIS.logdensity(ais, i, samples[1])
            @test AnnealedIS.logdensity_ratio(ais, i, samples[1]) == dratio
        end
    end

    @testset "Convergence test" begin
        D = 1
        N = 100
        num_samples = 1000

        y_obs = [3.0]

        prior_sampling(rng) = rand(rng, MvNormal(D, 1.0))
        prior_density(params) = logpdf(MvNormal(D, 1.0), params)
        joint_density(params) = logpdf(MvNormal(params, I), y_obs) + prior_density(params)

        ais = AnnealedISSampler(
            prior_sampling, prior_density, joint_density, N)

        samples, diagnostics = ais_sample(rng, ais, num_samples)

        posterior_mean = AnnealedIS.estimate_expectation(samples, x -> x)
        @test isapprox(posterior_mean, [1.5]; atol=3e-2)

        ess = effective_sample_size(samples)
        # TODO: What is a reasonable ESS to expect?
        @test ess > 900
        @test ess == diagnostics[:ess]
    end

    @testset "Logdensity ratio" begin
        D = 1
        N = 100
        num_samples = 1000

        y_obs = [3.0]

        prior_sampling(rng) = rand(rng, MvNormal(D, 1.0))
        prior_density(params) = logpdf(MvNormal(D, 1.0), params)
        joint_density(params) = logpdf(MvNormal(params, I), y_obs) + prior_density(params)

        betas = (geomspace(1, 1001, N+1) .- 1) ./ 1000
        anis_alg =
        ais = AnnealedISSampler(
            AnnealedIS.AnISModel(prior_sampling, prior_density, joint_density),
            AnIS(betas, (x, i) -> Normal(0, 1), SimpleRejection())
        )

        sample = [1.0]

        for i in 1:(length(ais.betas)-1)
            @test isapprox(
                AnnealedIS.logdensity_ratio(ais, i, sample),
                (AnnealedIS.logdensity(ais, i+1, sample) - AnnealedIS.logdensity(ais, i, sample));
                atol=1e-12 # TODO: What is the numerical error we would expect?
            )
        end
    end

    @testset "Sample from Turing" begin
        @model function test_model(y)
            x ~ Normal(0, 1)
            y ~ Normal(y, 1)
        end

        y_obs = 2
        tm = test_model(y_obs)

        named_tuple = sample_from_prior(rng, tm)
        named_tuple2 = sample_from_prior(rng, tm)
        @test typeof(named_tuple[:x]) == Float64
        @test named_tuple != named_tuple2

        # Test model with multiple latents and multivariate latent.
        D = 5
        @model function test_model2(y)
            a ~ Normal(0, 1)
            x ~ MvNormal(D, 1.0)
            y ~ MvNormal(x, I)
        end

        y_obs = ones(D)
        tm2 = test_model2(y_obs)

        named_tuple = sample_from_prior(rng, tm2)
        @test typeof(named_tuple[:a]) == Float64
        @test typeof(named_tuple[:x]) == typeof(rand(MvNormal(D, 1.0)))
    end

    @testset "Prior density from Turing" begin
        D = 5

        @model function test_model(y)
            a ~ MvNormal(D, 1)
            x ~ Normal(0, 1)
            y ~ Normal(x, 1)
        end

        y_obs = 2
        tm = test_model(y_obs)
        logprior_density = make_log_prior_density(tm)

        xval = 1
        aval = ones(D)
        nt = (x = xval, a = aval)
        prior_hand(nt) = logpdf(Normal(0, 1), nt[:x]) + logpdf(MvNormal(D, 1), nt[:a])
        @test logprior_density(nt) == prior_hand(nt)

        # Test that we can evaluate samples from prior
        nt = sample_from_prior(rng, tm)
        @test logprior_density(nt) == prior_hand(nt)
    end

    @testset "Joint density from Turing" begin
        @model function test_model(y)
            x ~ Normal(0, 1)
            y ~ Normal(x, 1)
        end

        y_obs = 2
        tm = test_model(y_obs)
        logjoint_density = make_log_joint_density(tm)

        xval = 1
        nt = (x = xval,)
        @test logjoint_density(nt) == logpdf(Normal(0, 1), xval) + logpdf(Normal(xval, 1), y_obs)
    end

    @testset "Turing AnnealedISSampler" begin
        N = 3
        num_samples = 10

        @model function test_model(y)
            x ~ Normal(0, 1)
            z ~ MvNormal(zeros(2))
            y ~ Normal(x, 1)
        end

        y_obs = 2
        tm = test_model(y_obs)

        ais = AnnealedISSampler(tm, N)
        samples, diagnostics = ais_sample(rng, ais, num_samples)

        chains = sample(rng, tm, AnIS(N), num_samples)

        @test typeof(chains) <: MCMCChains.Chains
        @test typeof(chains.logevidence) <: Real
        @test haskey(chains.info, :ess)
        @test haskey(get(chains, :log_weight), :log_weight)
        @test haskey(get(chains, :x), :x)

        # Test intermediate samples
        chains = sample(
            rng, tm, AnIS(20), num_samples; store_intermediate_samples=true
        )
        @test haskey(chains.info, :intermediate_samples)

        inter_samples = hcat(chains.info[:intermediate_samples]...)
        inter_samples = permutedims(inter_samples, [2, 1])
        @test size(inter_samples) == (num_samples, 11)

        @test haskey(chains.info, :intermediate_log_weights)
        inter_weights = hcat(chains.info[:intermediate_log_weights]...)
        inter_weights = permutedims(inter_weights, [2, 1])
        @test size(inter_weights) == (num_samples, 11)
    end

    @testset "Transition kernels" begin
        prior_sample = (a = 0.0, b = [1; 2])
        expected_transition_kernel = (a = Normal(0, 1), b = MvNormal(2, 1.0))
        transition_kernel = AnnealedIS.get_normal_transition_kernel(prior_sample, 1)
        @test expected_transition_kernel == transition_kernel
    end

    @testset "Intermediate Samples" begin
        N = 20
        num_samples = 10

        @model function test_model(y)
            x ~ Normal(0, 1)
            y ~ Normal(x, 1)
        end

        y_obs = 2
        tm = test_model(y_obs)

        ais = AnnealedISSampler(tm, N)
        samples, diagnostics = ais_sample(
            rng, ais, num_samples; store_intermediate_samples=true)

        inter_samples = hcat(diagnostics[:intermediate_samples]...)
        inter_samples = permutedims(inter_samples, [2, 1])
        @test size(inter_samples) == (num_samples, 11)
        @test map(x->x.params, samples) == inter_samples[:,end]

        @test haskey(diagnostics, :intermediate_log_weights)
    end

    @testset "Resampling Tests" begin
        D = 1
        N = 10
        num_samples = 100

        y_obs = [3.0]

        # Roughly half the samples should get rejected.
        prior_sampling(rng) = rand(rng, MvNormal(D, 1.0))
        prior_density(params) = logpdf(MvNormal(D, 1.0), params)
        joint_density(params) = logpdf(MvNormal(params, I), y_obs) + prior_density(params) + log(max(params[1],0))

        ais_simple_rejection = AnnealedISSampler(
            prior_sampling,
            prior_density,
            joint_density,
            N
        )
        samples_no_resampling, no_resampling_diag = ais_sample(
            rng,
            ais_simple_rejection,
            num_samples
        )

        ais_resampling = AnnealedISSampler(
            prior_sampling,
            prior_density,
            joint_density,
            N,
            RejectionResample()
        )
        samples_resampling, resampling_diag = ais_sample(rng, ais_resampling, num_samples)

        non_zero_weights(samples) = sum(samples) do weighted_sample
            Int(weighted_sample.log_weight != -Inf)
        end

        @test non_zero_weights(samples_resampling) == num_samples
        @test sum(resampling_diag[:num_rejected]) > 0

        num_accepted_samples = non_zero_weights(samples_no_resampling)
        @test num_accepted_samples < num_samples
        @test (num_samples - sum(no_resampling_diag[:num_rejected])) == num_accepted_samples
        @test maximum(no_resampling_diag[:num_rejected]) == 1
    end

    @testset "Out of distribution catch" begin
        @model function model()
            x ~ LogNormal(1, 1)
        end

        prior_density = make_log_prior_density(model())
        joint_density = make_log_joint_density(model())

        @test prior_density((x = -1,)) == -Inf
        @test joint_density((x = -1,)) == -Inf
    end

    @testset "AnISHMC" begin
        betas = [0.0, 0.5, 1.0]
        proposal = AdvancedHMC.StaticTrajectory(AdvancedHMC.Leapfrog(0.05), 10)
        anis = AnISHMC(
            betas,
            proposal,
            10,
            SimpleRejection()
        )
        @test anis.betas[1] == 0.0
    end

    @testset "AnnealedISSamplerHMC" begin
        N = 100
        num_samples = 10

        @model function test_model(y)
            x ~ Normal(0, 1)
            y ~ Normal(x, 1)
        end

        yval = 3
        tm = test_model(yval)

        #betas = collect(range(0, 1, length=100))
        betas = (geomspace(1, 1001, 101) .- 1) ./ 1000
        proposal = AdvancedHMC.StaticTrajectory(AdvancedHMC.Leapfrog(0.05), 10)
        anis = AnISHMC(
            betas,
            proposal,
            5,
            SimpleRejection()
        )
        spl = Turing.DynamicPPL.Sampler(anis, tm)
        anis_sampler = AnnealedISSamplerHMC(tm, anis, spl)

        # Check that hamiltonians are correct. The first hamiltonian should
        # have the target density specified by the first non-zero beta.
        prior_sample = anis_sampler.prior_sampling()
        for i in 2:(length(betas)-1)
            h = anis_sampler.hamiltonians[i-1]
            d = AnnealedIS.logdensity(anis_sampler, i, prior_sample)
            @test h.ℓπ(prior_sample) == d

            dr = AnnealedIS.logdensity_ratio(anis_sampler, i, prior_sample)
            hp1 = anis_sampler.hamiltonians[i]
            hdr = hp1.ℓπ(prior_sample) - h.ℓπ(prior_sample)
            @test isapprox(dr, hdr, atol=1e-12)
        end

        chain = sample(
            Random.GLOBAL_RNG,
            test_model(yval),
            anis,
            num_samples;
            store_intermediate_samples=true
        )

        acceptance_rates = chain.info[:intermediate_acceptance_rate]
        @test haskey(chain.info, :intermediate_acceptance_rate)
        @test haskey(chain.info, :acceptance_rate)

        #samples = Array(chain[:x])[:,1]
        #log_weights = Array(chain[:log_weight])[:,1]

        #posterior_mean = dot(samples, exp.(log_weights)) / (num_samples * exp(chain.logevidence))
        #@test isapprox(posterior_mean, 1.5; atol=3e-2)
    end

    # @testset "Constrained Variables" begin
    #     @model function test_model(y)
    #         x ~ Beta(1, 1)
    #         y ~ Normal(x, 1)
    #     end
    #     yval = 3

    #     betas = [0.0, 0.5, 1.0]
    #     proposal = AdvancedHMC.StaticTrajectory(AdvancedHMC.Leapfrog(0.05), 10)
    #     anis = AnISHMC(
    #         betas,
    #         proposal,
    #         10,
    #         SimpleRejection()
    #     )

    #     #chain = sample(Random.GLOBAL_RNG, test_model(yval), anis, 10)
    #     #xs = Array(chain[:x])
    #     #@test all(0 .< xs) && all(xs .<= 1)

    #     true_logjoint(x) = logpdf(Beta(1, 1), x) + logpdf(Normal(x, 1), yval)
    #     tm = test_model(yval)
    #     spl = Turing.DynamicPPL.Sampler(anis, tm)
    #     Turing.DynamicPPL.link!(spl.state.vi, spl)
    #     logjoint = AnnealedIS.gen_logjoint(spl.state.vi, tm, spl)

    #     @test logjoint([2.0]) == true_logjoint(2)
    #     @show logjoint([2.0])
    #     @test logjoint([0.5]) == true_logjoint(0.5)
    # end

    @testset "Density closures" begin
        @model function test_model(y)
            x ~ Normal(0, 1)
            y ~ Normal(x, 1)
        end

        betas = [0.0, 0.5, 1.0]
        proposal = AdvancedHMC.StaticTrajectory(AdvancedHMC.Leapfrog(0.05), 10)
        anis = AnISHMC{Turing.Core.ReverseDiffAD{false}}(
            betas,
            proposal,
            10,
            SimpleRejection()
        )

        yval = 3
        tm = test_model(yval)
        spl = Turing.DynamicPPL.Sampler(anis, tm)

        vi = spl.state.vi
        logjoint = AnnealedIS.gen_logjoint(vi, tm, spl)
        logprior = AnnealedIS.gen_logprior(vi, tm, spl)

        true_logjoint(x) = logpdf(Normal(0, 1), x[1]) + logpdf(Normal(x[1], 1), yval)
        @test logjoint([2.0]) == true_logjoint([2])

        true_logprior(x) = logpdf(Normal(0, 1), x[1])
        @test logprior([2.0]) == true_logprior([2])

        # Tests to check whether functions are ThreadSafe.
        num_vals = 20
        test_vals = rand(num_vals)
        test_vals = [[x] for x in test_vals]

        vi = Turing.VarInfo(tm)
        logjoint = AnnealedIS.gen_logjoint(vi, tm, spl)
        logjoint_grad = AnnealedIS.gen_logjoint_grad(vi, tm, spl)

        logprior = AnnealedIS.gen_logprior(vi, tm, spl)
        logprior_grad = AnnealedIS.gen_logprior_grad(vi, tm, spl)

        for fn in [logjoint, logjoint_grad, logprior, logprior_grad]
            vals_seq = map(fn, test_vals)
            vals_parallel = similar(vals_seq)
            Threads.@threads for i in 1:num_vals
                vals_parallel[i] = fn(test_vals[i])
            end
            @test vals_seq == vals_parallel

            vals_broadcast = fn.(test_vals)
            @test vals_seq == vals_broadcast
        end
    end

    @testset "AD backends" begin
        @model function test_model(y)
            x ~ Normal(0, 1)
            y ~ Normal(x, 1)
        end

        yval = 3
        tm = test_model(yval)

        num_vals = 20
        test_vals = rand(num_vals)
        test_vals = [[x] for x in test_vals]

        # Forward mode AD
        betas = [0.0, 0.5, 1.0]
        proposal = AdvancedHMC.StaticTrajectory(AdvancedHMC.Leapfrog(0.05), 10)
        anis_fw = AnISHMC{Turing.ForwardDiffAD{1}}(
            betas,
            proposal,
            10,
            SimpleRejection()
        )
        spl = Turing.DynamicPPL.Sampler(anis_fw, tm)
        vi = spl.state.vi

        logjoint_grad_fw = AnnealedIS.gen_logjoint_grad(vi, tm, spl)
        logprior_grad_fw = AnnealedIS.gen_logprior_grad(vi, tm, spl)

        joint_grads_fw = logjoint_grad_fw.(test_vals)
        prior_grads_fw = logprior_grad_fw.(test_vals)

        # ReverseDiff AD
        anis_rv = AnISHMC{Turing.Core.ReverseDiffAD{false}}(
            betas,
            proposal,
            10,
            SimpleRejection()
        )
        spl = Turing.DynamicPPL.Sampler(anis_rv, tm)
        vi = spl.state.vi

        logjoint_grad_rv = AnnealedIS.gen_logjoint_grad(vi, tm, spl)
        logprior_grad_rv = AnnealedIS.gen_logprior_grad(vi, tm, spl)

        joint_grads_rv = logjoint_grad_rv.(test_vals)
        prior_grads_rv = logprior_grad_rv.(test_vals)

        @test joint_grads_fw == joint_grads_rv
        @test prior_grads_fw == prior_grads_rv

        # Some sanity check that the joint and prior are different
        @test prior_grads_rv != joint_grads_rv
        @test prior_grads_fw != joint_grads_fw
    end
end


text:
module MathOptSetDistances

import MathOptInterface
const MOI = MathOptInterface
using BlockDiagonals: BlockDiagonal
import FillArrays
using LinearAlgebra
import ChainRulesCore
const CRC = ChainRulesCore
import FiniteDifferences

export distance_to_set, projection_on_set, projection_gradient_on_set

include("utils.jl")
include("distances.jl")
include("distance_sets.jl")
include("projections.jl")
include("chainrules.jl")

end # module


text:
function LoadTestAgaricusData(Oracles)

	(features, labels) = readData("data/TestAgaricus/TestAgaricus", (8124,126), 0)

	(classLabels, numClasses) = createClassLabels(labels)

	(trf, trl, trl2, numTP, tef, tel, tel2) = trainTestRandomSeparate(features, MinusPlusOneVector(labels, Set([1.0])), classLabels)

	L2regs=[
	false,
	true
	]

	for L2reg in L2regs
		(gradientOracle, numVars, numTP, csDataType, LossFunctionString, myOutputter, L2reg) = createBLOracles(trf, trl, numTP, tef, tel, L2reg; outputLevel = 0)

		push!(Oracles, (gradientOracle, numVars, numTP, csDataType, LossFunctionString, myOutputter, L2reg,  "TestAgaricus", (t)-> Problem(L2reg, "TestAgaricus", LossFunctionString,numTP, t,(j)->getSampleFunctionAt(j,gradientOracle))))

		(gradientOracle2, numVars2, numTP2, csDataType2, LossFunctionString2, myOutputter2, L2reg2) = createMLOracles(trf, trl2, numTP, numClasses, tef, tel2, L2reg; outputLevel = 0)

		push!(Oracles,(gradientOracle2, numVars2, numTP2, csDataType2, LossFunctionString2, myOutputter2, L2reg2, "TestAgaricus", (t)-> Problem(L2reg2, "TestAgaricus", LossFunctionString2,numTP2, t,(j)->getSampleFunctionAt(j,gradientOracle2))))
	end
end

LoadTestAgaricusData(Oracles)


text:
#
# This file is part of the Simulate.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#
# Benchmark adopted from the dicegame.ipynb example
#

using DiscreteEvents, Distributions, Random, BenchmarkTools
import Dates.now

mutable struct Worker
    nr::Int64              # worker number
    input::Channel
    output::Channel
    dist::Distribution     # distribution of processing time
    retard::Float64        # worker retard factor, 1: no retardation, >1: retardation
    done::Int64            # number of finished items

    Worker(nr, input, output, dist, perform) = new(nr, input, output, dist, 1/perform, 0)
end

function work(clk::Clock, w::Worker)
    job = take!(w.input)
    delay!(clk, rand(w.dist) * w.retard)
    put!(w.output, job)
    w.done += 1
end

function setup( n::Int64, mw::Int64,
                    vp::Distribution, vw::Distribution;
                    d=1000, seed=1234, jobs=1000)
    global clk = Clock()
    Random.seed!(seed)                  # seed random number generator
    global C = [Channel{Int64}(mw) for i in 1:n+1] # create n+1 channels with given buffer sizes
    C[1] = Channel{Int64}(Inf)                 # unlimited sizes for channels 1 and n+1
    C[n+1] = Channel{Int64}(Inf)
    j = reverse(Array(1:(n-1)*2))
    for i in n:-1:2                     # seed channels 2:(n-1) each with 2 inventory items
        C[i].sz_max > 0 ? put!(C[i], j[(i-1)*2]) : nothing
        C[i].sz_max > 1 ? put!(C[i], j[(i-1)*2-1]) : nothing
    end
    for i in ((n-1)*2+1):jobs           # put other jobs into channel 1
        put!(C[1], i)
    end

    wp = rand(vw, n)                    # calculate worker performance
    global W = [Worker(i, C[i], C[i+1], vp, wp[i]) for i in 1:n]
    for i in 1:n
        process!(clk, Prc(i, work, W[i]))
    end
end

# first run to compile
setup(5, 5, Gamma(10,1/10), Normal(1,0), jobs=1200)
run!(clk, 1000)

# then take measurements
println(now())
@time onthread(2) do; setup(5, 5, Gamma(10,1/10), Normal(1,0), jobs=1200); end
println(@time onthread(2) do; run!(clk, 1000); end)
println("events=", clk.evcount, " throughput:", length(C[end].data))

t = run(@benchmarkable onthread(2) do; run!(clk, 1000); end setup=onthread(2) do; setup(5, 5, Gamma(10,1/10), Normal(1,0), jobs=1200); end  evals=1 seconds=15.0 samples=50)


text:
module Pyramids

using Images, Interpolations, Colors

export ImagePyramid, PyramidType, ComplexSteerablePyramid, LaplacianPyramid, GaussianPyramid
export subband, toimage, update_subband, update_subband!, test

"Abstract supertype for the variety of pyramids that can be constructed."
abstract PyramidType
abstract SimplePyramid <: PyramidType

"""Type that indicates a complex steerable pyramid. [1]

[1] http://www.cns.nyu.edu/~eero/steerpyr/"""
type ComplexSteerablePyramid <: PyramidType
end

"""Type that indicates a Laplacian pyramid. [1]

[1] persci.mit.edu/pub_pdfs/pyramid83.pdf"""
type LaplacianPyramid <: SimplePyramid
end

"""Type that indicates a Gaussian pyramid. [1]

[1] http://persci.mit.edu/pub_pdfs/RCA84.pdf"""
type GaussianPyramid <: SimplePyramid
end

"""Type that represents a concrete pyramidal representation of a given input image. Each type of pyramid has its own parameters. The basic construction method is, for example

```pyramid = ImagePyramid(im, ComplexSteerablePyramid(), scale=0.5^0.25)```

See the code for more information on optional arguments."""
type ImagePyramid
    scale::Real
    num_levels::Int
    num_orientations::Int
    t::PyramidType
    pyramid_bands::Dict

    function ImagePyramid(pyr::ImagePyramid)
        return deepcopy(pyr)
    end

    function ImagePyramid(im::Array, t::ComplexSteerablePyramid; scale=0.5, min_size=15, num_orientations=8, max_levels=23, twidth=1)
        this = new()

        this.t = ComplexSteerablePyramid()
        this.scale = scale
        this.num_orientations = num_orientations

        im_dims = size(im)
        h = im_dims[1]
        w = im_dims[2]

        this.num_levels = min(ceil(log2(minimum([h w]))/log2(1/scale) - (log2(min_size)/log2(1/scale))),max_levels);

        pyramid_bands, mtx, harmonics = build_complex_steerable_pyramid(im, this.num_levels, this.num_levels, order=num_orientations-1, twidth=twidth, scale=scale)

        this.pyramid_bands = pyramid_bands

        return this
    end

    function ImagePyramid(im::Image, t::ComplexSteerablePyramid; scale=0.5, min_size=15, num_orientations=8, max_levels=23, twidth=1)
        im_arr = image_to_array(im)

        return ImagePyramid(im_arr, t, scale=scale, min_size=min_size, num_orientations=num_orientations, max_levels=23, twidth=1)
    end

    function ImagePyramid(im::Array, t::GaussianPyramid; min_size=15, max_levels=23, filter=[0.0625; 0.25; 0.375; 0.25; 0.0625])
        this = new()

        pyramid_bands, num_levels = generate_gaussian_pyramid(im, min_size=min_size, max_levels=max_levels, filter=filter)

        this.scale = 0.5
        this.pyramid_bands = pyramid_bands
        this.num_orientations = 1
        this.num_levels = num_levels
        this.t = t

        return this
    end

    function ImagePyramid(im::Array, t::LaplacianPyramid; min_size=15, max_levels=23, filter=[0.0625; 0.25; 0.375; 0.25; 0.0625])
        this = new()

        pyramid_bands, num_levels = generate_laplacian_pyramid(im, min_size=min_size, max_levels=max_levels, filter=filter)

        this.scale = 0.5
        this.pyramid_bands = pyramid_bands
        this.num_orientations = 1
        this.num_levels = num_levels
        this.t = t

        return this
    end

    function ImagePyramid(im::Image, t::SimplePyramid; min_size=15, max_levels=23, filter=[0.0625; 0.25; 0.375; 0.25; 0.0625])
        im_arr = image_to_array(im)

        return ImagePyramid(im_arr, t, min_size=min_size, max_levels=max_levels, filter=filter)
    end

    function ImagePyramid(pyramid_bands::Dict, scale, t, num_levels, num_orientations)
        this = new()

        this.scale = scale
        this.pyramid_bands = pyramid_bands
        this.num_orientations = num_orientations
        this.num_levels = num_levels
        this.t = t

        return this
    end

end

##############################
# Functions

"""
    subband(pyramid::ImagePyramid, level; orientation=())

Returns the sub-band of the `pyramid` at `level` and at `orientation.` If orientation is not provided, `subband` assumes that the pyramid is not oriented.

Level 0 is always the high frequency residual."""
function subband(pyramid::ImagePyramid, level; orientation = ())
    if isempty(orientation)
        return pyramid.pyramid_bands[level]
    else
        return pyramid.pyramid_bands[level][orientation]
    end
end

function update_subband!(pyramid::ImagePyramid, level, new_subband; orientation = ())
    if isempty(orientation)
        pyramid.pyramid_bands[level] = copy(new_subband)
    else
        pyramid.pyramid_bands[level][orientation] = copy(new_subband)
    end

    return pyramid
end

function update_subband(pyramid::ImagePyramid, level, new_subband; orientation = ())
    newpyramid = ImagePyramid(pyramid)
    return update_subband!(newpyramid, level, new_subband, orientation=orientation)
end

"Converts a pyramid to a 2-D array (not of type `Image`!)"
function toimage(pyramid::ImagePyramid)
    if typeof(pyramid.t) <: ComplexSteerablePyramid
        im = reconstruct_complex_steerable_pyramid(pyramid, scale=pyramid.scale)
    elseif typeof(pyramid.t) <: LaplacianPyramid
        im = reconstruct_laplacian_pyramid(pyramid)
    elseif typeof(pyramid.t) <: GaussianPyramid
        im = subband(pyramid, 0)
    else
        error("Unsupported pyramid type $(typeof(pyramid.T))")
    end

    return im
end

##############################
# Private helper functions.

function image_to_array(im::Image)
    # Please someone tell me a better way of checking for image type.
    if typeof(im.data[1]) <: Gray
        im_arr = separate(convert(Array{Float64,2}, im))
    else
        warn("ImagePyramid only supports 2D images currently. Returning luminance channel.")
        im_lab = convert(Image{Lab}, im)
        im_arr = convert(Array{Float64,3}, separate(im))[:,:,1]
    end

    return im_arr
end

##############################
# Private functions for building Gaussian and Laplacian pyramids.

function convolve_reflect(im, filter)
    filter_len = length(filter)
    filter_offset::Int = (filter_len-1)/2

    padded_im = zeros((size(im,1) + filter_len*2), (size(im,1) + filter_len*2))

    padded_im[(filter_len+1):(end-filter_len), (filter_len+1):(end-filter_len)] = im

    padded_im[(filter_len+1):(end-filter_len), 1:filter_len] = flipdim(padded_im[(filter_len+1):(end-filter_len), (filter_len+2):(2*filter_len+1)], 2)
    padded_im[(filter_len+1):(end-filter_len), (end-filter_len+1):end] = flipdim(padded_im[(filter_len+1):(end-filter_len), (end-2*filter_len):(end-filter_len-1)], 2)
    padded_im[1:filter_len, (filter_len+1):(end-filter_len),] = flipdim(padded_im[(filter_len+2):(1+2*filter_len), (filter_len+1):(end-filter_len)], 1)
    padded_im[(end-filter_len+1):end, (filter_len+1):(end-filter_len)] = flipdim(padded_im[(end-2*filter_len):(end-filter_len-1), (filter_len+1):(end-filter_len)], 1)

    new_im = conv2(filter, filter, padded_im)

    new_im = new_im[(1+filter_len+filter_offset):(end-filter_len-filter_offset), (1+filter_len+filter_offset):(end-filter_len-filter_offset)]
    return new_im
end

function reduce(im; filter=[0.0625; 0.25; 0.375; 0.25; 0.0625])
    new_im = convolve_reflect(im, filter)
    return new_im[1:2:end, 1:2:end]
end

function expand(im; filter=[0.0625; 0.25; 0.375; 0.25; 0.0625])
    new_im = zeros((size(im, 1)*2, size(im, 2)*2))
    new_im[1:2:end, 1:2:end] = im

    new_im = convolve_reflect(new_im, filter)
    return new_im
end

function generate_gaussian_pyramid(im; min_size=15, max_levels=23, filter=[0.0625; 0.25; 0.375; 0.25; 0.0625])
    im = convert(Array{Float64}, copy(im))

    pyramid_bands = Dict{Integer, Array}()

    im_dims = collect(size(im))
    num_levels = min(max_levels, ceil(Int, log2(minimum(im_dims)) - log2(min_size)))

    for i = 1:num_levels
        pyramid_bands[i-1] = im
        im = reduce(im, filter=filter)
    end

    return (pyramid_bands, num_levels)
end

function generate_laplacian_pyramid(im; min_size=15, max_levels=23, filter=[0.0625; 0.25; 0.375; 0.25; 0.0625])
    im = convert(Array{Float64}, copy(im))

    pyramid_bands = Dict{Integer, Array}()

    im_dims = collect(size(im))

    num_levels = min(max_levels, ceil(Int, log2(minimum(im_dims)) - log2(min_size)))

    filter_offset::Int = (length(filter)-1)/2

    for i = 1:num_levels
        reduced_im = reduce(im, filter=filter)
        next_im = expand(reduced_im)
        diff_im = im - next_im

        pyramid_bands[i-1] = diff_im

        im = reduced_im
    end

    pyramid_bands[num_levels] = im

    return (pyramid_bands, num_levels)
end

function reconstruct_laplacian_pyramid(pyramid::ImagePyramid)
    output_im::Array{Float64} = subband(pyramid, pyramid.num_levels)

    for i = (pyramid.num_levels-1):-1:0
        output_im = expand(output_im)
        output_im += subband(pyramid, i)
    end

    return output_im
end

##############################
# Private functions for building complex steerable pyramids.

function construct_steering_matrix(harmonics, angles; even = true)
    numh = 2*length(harmonics) - any(harmonics .== 0)

    imtx = zeros(length(angles), numh)
    col = 1
    for h in harmonics
        args = h * angles

        if h == 0
            imtx[:,col] = ones(angles)
            col += 1
        elseif !even
            imtx[:,col] = sin(args)
            imtx[:,col+1] = -cos(args)
            col += 2
        else
            imtx[:,col] = cos(args)
            imtx[:,col+1] = sin(args)
            col += 2
        end
    end

    r = rank(imtx)

    if (r != numh) && (r != length(angles))
        warning("Matrix is not full rank")
    end

    return pinv(imtx)
end

function raisedcosine(width=1, position=0, values=[0,1])
    sz = 256 # arbitrary
    X = pi * (-sz-1:1) / (2 * sz)
    Y = values[1] + (values[2]-values[1]) * cos(X).^2

    Y[1] = Y[2]
    Y[sz+3] = Y[sz+2]
    X = position + (2*width/pi) * (X + pi/4)

    return (X, Y)
end

function build_complex_steerable_pyramid(im, height, nScales; order=3, twidth=1, scale=0.5)
    pyramid_bands = Dict{Integer, Union{Array, Dict{Integer, Array}}}()

    num_orientations = order + 1

    # generate steering matrix and harmonics info
    if mod(num_orientations, 2) == 0
        harmonics = ((0:((num_orientations/2)-1))*2 + 1)'
    else
        harmonics = ((0:((num_orientations-1)/2))*2)'
    end

    steeringmatrix = construct_steering_matrix(harmonics, pi*(0:(num_orientations-1))/num_orientations, even=true)

    imdft = fftshift(fft(im))
    im_dims = collect(size(im))
    ctr = ceil(Int, (im_dims+0.5)/2)

    angle = broadcast(atan2, (((1:im_dims[1]) - ctr[1]) ./ (im_dims[1]/2)), (((1:im_dims[2]) - ctr[2]) ./ (im_dims[2]/2))') #SLOOW
    log_rad = broadcast((x,y) -> log2(sqrt(x.^2 + y.^2)), (((1:im_dims[1]) - ctr[1]) ./ (im_dims[1]/2)), (((1:im_dims[2]) - ctr[2]) ./ (im_dims[2]/2))') #SLOOW
    log_rad[ctr[1], ctr[2]] = log_rad[ctr[1], ctr[2]-1]

    Xrcos, Yrcos = raisedcosine(twidth, (-twidth/2), [0 1])
    Yrcos = sqrt(Yrcos)
    YIrcos = sqrt(1 - Yrcos.^2)

    # generate high frequency residual
    Yrcosinterpolant = interpolate((Xrcos,), Yrcos, Gridded(Linear()))
    hi0mask = reshape(Yrcosinterpolant[log_rad], size(log_rad))
    hi0dft =  imdft .* hi0mask;

    pyramid_bands[0] = ifft(ifftshift(hi0dft));

    YIrcosinterpolant = interpolate((Xrcos,), YIrcos, Gridded(Linear()))
    lo0mask = reshape(YIrcosinterpolant[log_rad], size(log_rad))

    lo0dft = imdft .* lo0mask

    for ht = height:-1:1
        Xrcos = Xrcos - log2(1/scale)
        order = num_orientations - 1

        cnst = (2^(2*order))*(factorial(order)^2)/(num_orientations*factorial(2*order))

        Yrcosinterpolant = interpolate((Xrcos,), Yrcos, Gridded(Linear()))
        himask = reshape(Yrcosinterpolant[log_rad], size(log_rad))

        # loop through each orientation band
        pyramid_level = Dict{Integer, Array}()

        for b in 1:num_orientations
            banddft = zeros(lo0dft)

            for iter in eachindex(lo0dft)
                ang = angle[iter]-pi*(b-1)/num_orientations
                a = (abs(mod(pi+ang, 2*pi) - pi) .< pi/2) .* (2*sqrt(cnst) * (cos(ang).^order))
                banddft[iter] = ((complex(0,-1)).^(num_orientations-1)) * lo0dft[iter] * himask[iter] * a
            end

            pyramid_level[b] = ifft(ifftshift(banddft))
        end

        pyramid_bands[height-ht+1] = pyramid_level

        dims = collect(size(lo0dft))
        ctr = ceil(Int, (dims+0.5)/2)

        lodims = round(Int, im_dims[1:2]*scale^(nScales-ht+1))

        loctr = ceil(Int, (lodims+0.5)/2)
        lostart = ctr - loctr+1
        loend = lostart + lodims - 1

        log_rad = log_rad[lostart[1]:loend[1], lostart[2]:loend[2]]
        angle = angle[lostart[1]:loend[1], lostart[2]:loend[2]]
        lodft = lo0dft[lostart[1]:loend[1], lostart[2]:loend[2]]
        YIrcos = abs(sqrt(1 - Yrcos.^2))

        YIrcosinterpolant = interpolate((Xrcos,), YIrcos, Gridded(Linear()))
        lomask = reshape(YIrcosinterpolant[log_rad], size(log_rad))

        lo0dft = lomask .* lodft
    end

    pyramid_bands[height+1] = real(ifft(ifftshift(lo0dft)))

    return (pyramid_bands, steeringmatrix, harmonics)
end

function make_angle_grid(sz, phase=0, origin=-1)
    if length(sz) == 1
        sz = [sz, sz]
    end

    if origin == -1
        origin = (sz + 1)/2
    end

    xramp = ones(round(Int, sz[1]), 1) * collect((1:sz[2]) - origin[2])'
    yramp = collect((1:sz[1]) - origin[1]) * ones(1, round(Int, sz[2]))

    res = atan2(yramp, xramp)

    res = mod(res+(pi-phase), 2*pi) - pi

    return res
end

function reconstruct_steerable_pyramid(pyr::ImagePyramid; levs="all", bands="all", twidth=1, scale=0.5)
    dims = collect(size(subband(pyr, 0)))
    im_dft = zeros(Complex{Float64}, size(subband(pyr, 0)))

    ctr = ceil(Int, (dims + 0.5)/2)

    angle = broadcast(atan2, (((1:dims[1]) - ctr[1]) ./ (dims[1]/2)), (((1:dims[2]) - ctr[2]) ./ (dims[2]/2))') #SLOOW
    log_rad = broadcast((x,y) -> log2(sqrt(x.^2 + y.^2)), (((1:dims[1]) - ctr[1]) ./ (dims[1]/2)), (((1:dims[2]) - ctr[2]) ./ (dims[2]/2))') #SLOOW
    log_rad[ctr[1], ctr[2]] = log_rad[ctr[1], ctr[2]-1]
    log_rad0 = log_rad
    angle0 = angle

    Xrcos, Yrcos = raisedcosine(twidth, (-twidth/2), [0 1])
    Yrcos = sqrt(Yrcos)
    YIrcos = sqrt(1 - Yrcos.^2)

    # Start with low frequency residual
    low_dft = fftshift(fft(subband(pyr, pyr.num_levels+1)))

    lodims = collect(size(low_dft))
    loctr = ceil(Int, (lodims+0.5)/2)
    lostart = ctr - loctr+1
    loend = lostart + lodims - 1

    log_rad = log_rad0[lostart[1]:loend[1], lostart[2]:loend[2]]
    angle = angle0[lostart[1]:loend[1], lostart[2]:loend[2]]

    Xrcos = Xrcos - log2(1/scale)*pyr.num_levels
    YIrcosinterpolant = interpolate((Xrcos,), YIrcos, Gridded(Linear()))
    lomask = reshape(YIrcosinterpolant[log_rad], size(log_rad))

    im_dft[lostart[1]:loend[1], lostart[2]:loend[2]] += low_dft .* lomask

    # Accumulate mid-bamds
    for level = pyr.num_levels:-1:1
        lodims = collect(size(subband(pyr, level, orientation=1)))
        loctr = ceil(Int, (lodims+0.5)/2)
        lostart = ctr - loctr+1
        loend = lostart + lodims - 1

        log_rad = copy(log_rad0[lostart[1]:loend[1], lostart[2]:loend[2]])
        angle = copy(angle0[lostart[1]:loend[1], lostart[2]:loend[2]])

        Yrcosinterpolant = interpolate((Xrcos,), Yrcos, Gridded(Linear()))
        himask = reshape(Yrcosinterpolant[log_rad], size(log_rad))

        Xrcos = Xrcos + log2(1/scale)
        order = pyr.num_orientations - 1

        cnst = ((complex(0,1)).^(pyr.num_orientations-1)) * sqrt((2^(2*order))*(factorial(order)^2)/(pyr.num_orientations*factorial(2*order)))

        for orientation = 1:pyr.num_orientations
            band_dft = fftshift(fft(subband(pyr, level, orientation=orientation)))
            ang = angle-pi*(orientation-1)/pyr.num_orientations
            angle_mask = 2 *(abs(mod(pi+ang, 2*pi) - pi) .< pi/2) .* (cnst * (cos(ang).^order))

            im_dft[lostart[1]:loend[1], lostart[2]:loend[2]] += band_dft .* himask .* angle_mask
        end

        YIrcosinterpolant = interpolate((Xrcos,), YIrcos, Gridded(Linear()))
        lomask = reshape(YIrcosinterpolant[log_rad], size(log_rad))

        # Everything must be scaled by the low-frequency mask
        im_dft[lostart[1]:loend[1], lostart[2]:loend[2]] .*= lomask
    end

    # Add high frequency residual
    Yrcosinterpolant = interpolate((Xrcos,), Yrcos, Gridded(Linear()))
    hi0mask = reshape(Yrcosinterpolant[log_rad0], size(log_rad))
    im_dft += fftshift(fft(subband(pyr, 0))) .* hi0mask;

    return real(ifft(ifftshift(im_dft)))
end

# make the complex steerable pyramid real
function convert_complex_steerable_pyramid_to_real(pyramid::ImagePyramid; levs="all", bands="all", twidth=1, scale=0.5)
    pyramid = ImagePyramid(pyramid)
    num_levels = pyramid.num_levels
    num_orientations = pyramid.num_orientations

    for nsc in 1:num_levels
        dims = collect(size(subband(pyramid, nsc, orientation=1)))
        ctr = ceil(Int, (dims+0.5)/2)
        ang = make_angle_grid(dims, 0, ctr)
        ang[ctr[1], ctr[2]] = -pi/2

        for nor = 1:num_orientations
            ch = subband(pyramid, nsc, orientation=nor)

            ang0 = pi*(nor-1)/num_orientations
            xang = mod(ang-ang0+pi, 2*pi) - pi

            # this creates an angular mask
            amask = 2*(abs(xang) .< pi/2) + (abs(xang) .== pi/2)
            amask[ctr[1], ctr[2]] = 1
            amask[:,1] = 1
            amask[1,:] = 1

            # and masks the fft by it
            amask = fftshift(amask)
            ch = ifft(amask.*fft(ch))
            ch = 0.5 * real(ch)

            # then creates a new pyramid
            update_subband!(pyramid, nsc, ch, orientation=nor)
        end
    end

    # and returns it as a real
    return pyramid
end

function reconstruct_complex_steerable_pyramid(pyramid::ImagePyramid; levs="all", bands="all", twidth=1, scale=0.5)
    real_pyramid = convert_complex_steerable_pyramid_to_real(pyramid, levs=levs, bands=bands, twidth=twidth, scale=scale)
    return reconstruct_steerable_pyramid(real_pyramid, levs=levs, bands=bands, twidth=twidth, scale=scale)
end

end


text:
export jacobi, jacobi_weight

"""
    jacobi_weight(α, β, x)

Compute the weighting function (1-x)^α * (1-x)^β .
"""
function jacobi_weight(α, β, x)
	return (1-x)^α * (1-x )^β
end

"""
    abc(n, α, β, x)

Compute the three coefficients in the recurrence relation for the Jacobi polynomial.

n is the order of the polynomial.

α,β are the two weighting parameters.

x is the argument or point at which the polynomial will be evaluated.
"""
function abc(n, α, β, x)
	# the parens emphasize the terms, not PEMDAS
	a = (2n) * (n + α + β) * (2n + α + β - 2)
	b1 = (2n + α + β -1)
	b2 = (2n + α + β)
	b2 = b2 * (b2 - 2)
	b = b1 * (b2 * x + α^2 - β^2)
	c = 2 * (n + α - 1) * (n + β - 1) * (2n + α + β)
	return a, b, c
end

"""
    jacobi(n, α, β, x)

Compute the Jacobi polynomial of order n and weights α,β at point x.

This function uses a recurrence relation and is numerical stable to very high
order.  The computation time is linear w.r.t. the order n.  jacobi.(n, a, b, x)
is also linear w.r.t. the size of argument x.
"""
function jacobi(n, α, β, x)
	if n == 0
		return 1.0
	elseif n == 1
		return (α + 1) + (α + β + 2) * ((x-1)/2)
	end
	# this uses a loop to avoid recursion
	# we init P of n-2, n-1, and n=2.
	# the latter is the first recursive term.
	# then loop from 3 up to n, updating
	# Pnm2, Pnm1, a, b, c, Pn
	Pnm2 = jacobi(0, α, β, x)
	Pnm1 = jacobi(1, α, β, x)
	a, b, c = abc(2, α, β, x)
	Pn = ((b * Pnm1) - (c * Pnm2)) / a
	if n == 2
		return Pn
	end
	for i = 3:n
		Pnm2 = Pnm1
		Pnm1 = Pn
		a, b, c = abc(n, α, β, x)
		Pn = ((b * Pnm1) - (c * Pnm2)) / a
	end
	return Pn
end


text:
@testset "Julia Reserved Words Constant Tests" begin

  @test isdefined(Julz, :julia_reserved_words) == true

end


text:
# RELOG: Reverse Logistics Optimization
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using CRC
using CSV
using DataFrames
using Shapefile
using Statistics
using ZipFile
using ProgressBars
using OrderedCollections

import Downloads: download
import Base: parse

crc32 = crc(CRC_32)

struct GeoPoint
    lat::Float64
    lon::Float64
end

struct GeoRegion
    centroid::GeoPoint
    population::Int
    GeoRegion(; centroid, population) = new(centroid, population)
end

DB_CACHE = Dict{String,Dict{String,GeoRegion}}()

function centroid(geom::Shapefile.Polygon)::GeoPoint
    x_max, x_min, y_max, y_min = -Inf, Inf, -Inf, Inf
    for p in geom.points
        x_max = max(x_max, p.x)
        x_min = min(x_min, p.x)
        y_max = max(y_max, p.y)
        y_min = min(y_min, p.y)
    end
    x_center = (x_max + x_min) / 2.0
    y_center = (y_max + y_min) / 2.0
    return GeoPoint(round(y_center, digits = 5), round(x_center, digits = 5))
end

function _download_file(url, output, expected_crc32)::Nothing
    if isfile(output)
        return
    end
    mkpath(dirname(output))
    @info "Downloading: $url"
    fname = download(url)
    actual_crc32 = open(crc32, fname)
    expected_crc32 == actual_crc32 || error("CRC32 mismatch")
    cp(fname, output)
    return
end

function _download_zip(url, outputdir, expected_output_file, expected_crc32)::Nothing
    if isfile(expected_output_file)
        return
    end
    mkpath(outputdir)
    @info "Downloading: $url"
    zip_filename = download(url)
    actual_crc32 = open(crc32, zip_filename)
    expected_crc32 == actual_crc32 || error("CRC32 mismatch")
    open(zip_filename) do zip_file
        zr = ZipFile.Reader(zip_file)
        for file in zr.files
            open(joinpath(outputdir, file.name), "w") do output_file
                write(output_file, read(file))
            end
        end
    end
    return
end

function _geodb_load_gov_census(;
    db_name,
    extract_cols,
    shp_crc32,
    shp_filename,
    shp_url,
    population_url,
    population_crc32,
    population_col,
    population_preprocess,
    population_join,
)::Dict{String,GeoRegion}
    basedir = joinpath(dirname(@__FILE__), "..", "..", "data", db_name)
    csv_filename = "$basedir/locations.csv"
    if !isfile(csv_filename)
        # Download required files
        _download_zip(shp_url, basedir, joinpath(basedir, shp_filename), shp_crc32)
        _download_file(population_url, "$basedir/population.csv", population_crc32)

        # Read shapefile
        @info "Processing: $shp_filename"
        table = Shapefile.Table(joinpath(basedir, shp_filename))
        geoms = Shapefile.shapes(table)

        # Build empty dataframe
        df = DataFrame()
        cols = extract_cols(table, 1)
        for k in keys(cols)
            df[!, k] = []
        end
        df[!, "latitude"] = Float64[]
        df[!, "longitude"] = Float64[]

        # Add regions to dataframe
        for (i, geom) in tqdm(enumerate(geoms))
            c = centroid(geom)
            cols = extract_cols(table, i)
            push!(df, [values(cols)..., c.lat, c.lon])
        end
        sort!(df)

        # Join with population data
        population = DataFrame(CSV.File("$basedir/population.csv"))
        population_preprocess(population)
        population = population[:, [population_join, population_col]]
        rename!(population, population_col => "population")
        df = leftjoin(df, population, on = population_join)

        # Write output
        CSV.write(csv_filename, df)
    end
    if db_name ∉ keys(DB_CACHE)
        csv = CSV.File(csv_filename)
        DB_CACHE[db_name] = Dict(
            string(row.id) => GeoRegion(
                centroid = GeoPoint(row.latitude, row.longitude),
                population = (row.population === missing ? 0 : row.population),
            ) for row in csv
        )
    end
    return DB_CACHE[db_name]
end

# 2018 US counties
# -----------------------------------------------------------------------------
function _extract_cols_2018_us_county(
    table::Shapefile.Table,
    i::Int,
)::OrderedDict{String,Any}
    return OrderedDict(
        "id" => table.STATEFP[i] * table.COUNTYFP[i],
        "statefp" => table.STATEFP[i],
        "countyfp" => table.COUNTYFP[i],
        "name" => table.NAME[i],
    )
end

function _population_preprocess_2018_us_county(df)
    df[!, "id"] = [@sprintf("%02d%03d", row.STATE, row.COUNTY) for row in eachrow(df)]
end

function _geodb_load_2018_us_county()::Dict{String,GeoRegion}
    return _geodb_load_gov_census(
        db_name = "2018-us-county",
        extract_cols = _extract_cols_2018_us_county,
        shp_crc32 = 0x83eaec6d,
        shp_filename = "cb_2018_us_county_500k.shp",
        shp_url = "https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_county_500k.zip",
        population_url = "https://www2.census.gov/programs-surveys/popest/datasets/2010-2019/counties/totals/co-est2019-alldata.csv",
        population_crc32 = 0xf85b0405,
        population_col = "POPESTIMATE2019",
        population_join = "id",
        population_preprocess = _population_preprocess_2018_us_county,
    )
end

# US States
# -----------------------------------------------------------------------------
function _extract_cols_us_state(table::Shapefile.Table, i::Int)::OrderedDict{String,Any}
    return OrderedDict(
        "id" => table.STUSPS[i],
        "statefp" => parse(Int, table.STATEFP[i]),
        "name" => table.NAME[i],
    )
end

function _population_preprocess_us_state(df)
    rename!(df, "STATE" => "statefp")
end

function _geodb_load_us_state()::Dict{String,GeoRegion}
    return _geodb_load_gov_census(
        db_name = "us-state",
        extract_cols = _extract_cols_us_state,
        shp_crc32 = 0x9469e5ca,
        shp_filename = "cb_2018_us_state_500k.shp",
        shp_url = "https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_state_500k.zip",
        population_url = "http://www2.census.gov/programs-surveys/popest/datasets/2010-2019/national/totals/nst-est2019-alldata.csv",
        population_crc32 = 0x191cc64c,
        population_col = "POPESTIMATE2019",
        population_join = "statefp",
        population_preprocess = _population_preprocess_us_state,
    )
end

function geodb_load(db_name::AbstractString)::Dict{String,GeoRegion}
    db_name == "2018-us-county" && return _geodb_load_2018_us_county()
    db_name == "us-state" && return _geodb_load_us_state()
    error("Unknown database: $db_name")
end

function geodb_query(name)::GeoRegion
    db_name, id = split(name, ":")
    return geodb_load(db_name)[id]
end

max_stars_repo_name: ANL-CEEESA/RELOG

max_issues_repo_name: ANL-CEEESA/RELOG

max_forks_repo_name: ANL-CEEESA/RELOG


text:
# Single line comments start with a hash (pound) symbol.
#= Multiline comments can be written
   by putting '#=' before the text  and '=#'
   after the text. They can also be nested.
=#

####################################################
## 1. Primitive Datatypes and Operators
####################################################

# Everything in Julia is an expression.

# There are several basic types of numbers.
3 # => 3 (Int64)
3.2 # => 3.2 (Float64)
2 + 1im # => 2 + 1im (Complex{Int64})
2//3 # => 2//3 (Rational{Int64})

# All of the normal infix operators are available.
1 + 1 # => 2
8 - 1 # => 7
10 * 2 # => 20
35 / 5 # => 7.0
5 / 2 # => 2.5 # dividing an Int by an Int always results in a Float
div(5, 2) # => 2 # for a truncated result, use div
5 \ 35 # => 7.0
2 ^ 2 # => 4 # power, not bitwise xor
12 % 10 # => 2

# Enforce precedence with parentheses
(1 + 3) * 2 # => 8

# Bitwise Operators
~2 # => -3   # bitwise not
3 & 5 # => 1 # bitwise and
2 | 4 # => 6 # bitwise or
2 $ 4 # => 6 # bitwise xor
2 >>> 1 # => 1 # logical shift right
2 >> 1  # => 1 # arithmetic shift right
2 << 1  # => 4 # logical/arithmetic shift left

# You can use the bits function to see the binary representation of a number.
bits(12345)
# => "0000000000000000000000000000000000000000000000000011000000111001"
bits(12345.0)
# => "0100000011001000000111001000000000000000000000000000000000000000"

# Boolean values are primitives
true
false

# Boolean operators
!true # => false
!false # => true
1 == 1 # => true
2 == 1 # => false
1 != 1 # => false
2 != 1 # => true
1 < 10 # => true
1 > 10 # => false
2 <= 2 # => true
2 >= 2 # => true
# Comparisons can be chained
1 < 2 < 3 # => true
2 < 3 < 2 # => false

# Strings are created with "
"This is a string."

# Julia has several types of strings, including ASCIIString and UTF8String.
# More on this in the Types section.

# Character literals are written with '
'a'

# Some strings can be indexed like an array of characters
"This is a string"[1] # => 'T' # Julia indexes from 1
# However, this is will not work well for UTF8 strings,
# so iterating over strings is recommended (map, for loops, etc).

# $ can be used for string interpolation:
"2 + 2 = $(2 + 2)" # => "2 + 2 = 4"
# You can put any Julia expression inside the parentheses.

# Another way to format strings is the printf macro.
@printf "%d is less than %f" 4.5 5.3 # 4.5 is less than 5.300000

# Printing is easy
println("I'm Julia. Nice to meet you!")

# String can be compared lexicographically
"good" > "bye" # => true
"good" == "good" # => true
"1 + 2 = 3" == "1 + 2 = $(1+2)" # => true

####################################################
## 2. Variables and Collections
####################################################

# You don't declare variables before assigning to them.
some_var = 5 # => 5
some_var # => 5

# Accessing a previously unassigned variable is an error
try
    some_other_var # => ERROR: some_other_var not defined
catch e
    println(e)
end

# Variable names start with a letter or underscore.
# After that, you can use letters, digits, underscores, and exclamation points.
SomeOtherVar123! = 6 # => 6

# You can also use certain unicode characters
☃ = 8 # => 8
# These are especially handy for mathematical notation
2 * π # => 6.283185307179586

# A note on naming conventions in Julia:
#
# * Word separation can be indicated by underscores ('_'), but use of
#   underscores is discouraged unless the name would be hard to read
#   otherwise.
#
# * Names of Types begin with a capital letter and word separation is shown
#   with CamelCase instead of underscores.
#
# * Names of functions and macros are in lower case, without underscores.
#
# * Functions that modify their inputs have names that end in !. These
#   functions are sometimes called mutating functions or in-place functions.

# Arrays store a sequence of values indexed by integers 1 through n:
a = Int64[] # => 0-element Int64 Array

# 1-dimensional array literals can be written with comma-separated values.
b = [4, 5, 6] # => 3-element Int64 Array: [4, 5, 6]
b = [4; 5; 6] # => 3-element Int64 Array: [4, 5, 6]
b[1] # => 4
b[end] # => 6

# 2-dimensional arrays use space-separated values and semicolon-separated rows.
matrix = [1 2; 3 4] # => 2x2 Int64 Array: [1 2; 3 4]

# Arrays of a particular Type
b = Int8[4, 5, 6] # => 3-element Int8 Array: [4, 5, 6]

# Add stuff to the end of a list with push! and append!
push!(a,1)     # => [1]
push!(a,2)     # => [1,2]
push!(a,4)     # => [1,2,4]
push!(a,3)     # => [1,2,4,3]
append!(a,b) # => [1,2,4,3,4,5,6]

# Remove from the end with pop
pop!(b)        # => 6 and b is now [4,5]

# Let's put it back
push!(b,6)   # b is now [4,5,6] again.

a[1] # => 1 # remember that Julia indexes from 1, not 0!

# end is a shorthand for the last index. It can be used in any
# indexing expression
a[end] # => 6

# we also have shift and unshift
shift!(a) # => 1 and a is now [2,4,3,4,5,6]
unshift!(a,7) # => [7,2,4,3,4,5,6]

# Function names that end in exclamations points indicate that they modify
# their argument.
arr = [5,4,6] # => 3-element Int64 Array: [5,4,6]
sort(arr) # => [4,5,6]; arr is still [5,4,6]
sort!(arr) # => [4,5,6]; arr is now [4,5,6]

# Looking out of bounds is a BoundsError
try
    a[0] # => ERROR: BoundsError() in getindex at array.jl:270
    a[end+1] # => ERROR: BoundsError() in getindex at array.jl:270
catch e
    println(e)
end

# Errors list the line and file they came from, even if it's in the standard
# library. If you built Julia from source, you can look in the folder base
# inside the julia folder to find these files.

# You can initialize arrays from ranges
a = [1:5;] # => 5-element Int64 Array: [1,2,3,4,5]

# You can look at ranges with slice syntax.
a[1:3] # => [1, 2, 3]
a[2:end] # => [2, 3, 4, 5]

# Remove elements from an array by index with splice!
arr = [3,4,5]
splice!(arr,2) # => 4 ; arr is now [3,5]

# Concatenate lists with append!
b = [1,2,3]
append!(a,b) # Now a is [1, 2, 3, 4, 5, 1, 2, 3]

# Check for existence in a list with in
in(1, a) # => true

# Examine the length with length
length(a) # => 8

# Tuples are immutable.
tup = (1, 2, 3) # => (1,2,3) # an (Int64,Int64,Int64) tuple.
tup[1] # => 1
try:
    tup[1] = 3 # => ERROR: no method setindex!((Int64,Int64,Int64),Int64,Int64)
catch e
    println(e)
end

# Many list functions also work on tuples
length(tup) # => 3
tup[1:2] # => (1,2)
in(2, tup) # => true

# You can unpack tuples into variables
a, b, c = (1, 2, 3) # => (1,2,3)  # a is now 1, b is now 2 and c is now 3

# Tuples are created even if you leave out the parentheses
d, e, f = 4, 5, 6 # => (4,5,6)

# A 1-element tuple is distinct from the value it contains
(1,) == 1 # => false
(1) == 1 # => true

# Look how easy it is to swap two values
e, d = d, e  # => (5,4) # d is now 5 and e is now 4

# Dictionaries store mappings
empty_dict = Dict() # => Dict{Any,Any}()

# You can create a dictionary using a literal
filled_dict = Dict("one"=> 1, "two"=> 2, "three"=> 3)
# => Dict{ASCIIString,Int64}

# Look up values with []
filled_dict["one"] # => 1

# Get all keys
keys(filled_dict)
# => KeyIterator{Dict{ASCIIString,Int64}}(["three"=>3,"one"=>1,"two"=>2])
# Note - dictionary keys are not sorted or in the order you inserted them.

# Get all values
values(filled_dict)
# => ValueIterator{Dict{ASCIIString,Int64}}(["three"=>3,"one"=>1,"two"=>2])
# Note - Same as above regarding key ordering.

# Check for existence of keys in a dictionary with in, haskey
in(("one" => 1), filled_dict) # => true
in(("two" => 3), filled_dict) # => false
haskey(filled_dict, "one") # => true
haskey(filled_dict, 1) # => false

# Trying to look up a non-existent key will raise an error
try
    filled_dict["four"] # => ERROR: key not found: four in getindex at dict.jl:489
catch e
    println(e)
end

# Use the get method to avoid that error by providing a default value
# get(dictionary,key,default_value)
get(filled_dict,"one",4) # => 1
get(filled_dict,"four",4) # => 4

# Use Sets to represent collections of unordered, unique values
empty_set = Set() # => Set{Any}()
# Initialize a set with values
filled_set = Set([1,2,2,3,4]) # => Set{Int64}(1,2,3,4)

# Add more values to a set
push!(filled_set,5) # => Set{Int64}(5,4,2,3,1)

# Check if the values are in the set
in(2, filled_set) # => true
in(10, filled_set) # => false

# There are functions for set intersection, union, and difference.
other_set = Set([3, 4, 5, 6]) # => Set{Int64}(6,4,5,3)
intersect(filled_set, other_set) # => Set{Int64}(3,4,5)
union(filled_set, other_set) # => Set{Int64}(1,2,3,4,5,6)
setdiff(Set([1,2,3,4]),Set([2,3,5])) # => Set{Int64}(1,4)

####################################################
## 3. Control Flow
####################################################

# Let's make a variable
some_var = 5

# Here is an if statement. Indentation is not meaningful in Julia.
if some_var > 10
    println("some_var is totally bigger than 10.")
elseif some_var < 10    # This elseif clause is optional.
    println("some_var is smaller than 10.")
else                    # The else clause is optional too.
    println("some_var is indeed 10.")
end
# => prints "some var is smaller than 10"

# For loops iterate over iterables.
# Iterable types include Range, Array, Set, Dict, and AbstractString.
for animal=["dog", "cat", "mouse"]
    println("$animal is a mammal")
    # You can use $ to interpolate variables or expression into strings
end
# prints:
#    dog is a mammal
#    cat is a mammal
#    mouse is a mammal

# You can use 'in' instead of '='.
for animal in ["dog", "cat", "mouse"]
    println("$animal is a mammal")
end
# prints:
#    dog is a mammal
#    cat is a mammal
#    mouse is a mammal

for a in Dict("dog"=>"mammal","cat"=>"mammal","mouse"=>"mammal")
    println("$(a[1]) is a $(a[2])")
end
# prints:
#    dog is a mammal
#    cat is a mammal
#    mouse is a mammal

for (k,v) in Dict("dog"=>"mammal","cat"=>"mammal","mouse"=>"mammal")
    println("$k is a $v")
end
# prints:
#    dog is a mammal
#    cat is a mammal
#    mouse is a mammal

# While loops loop while a condition is true
x = 0
while x < 4
    println(x)
    x += 1  # Shorthand for x = x + 1
end
# prints:
#   0
#   1
#   2
#   3

# Handle exceptions with a try/catch block
try
   error("help")
catch e
   println("caught it $e")
end
# => caught it ErrorException("help")

####################################################
## 4. Functions
####################################################

# The keyword 'function' creates new functions
#function name(arglist)
#  body...
#end
function add(x, y)
    println("x is $x and y is $y")

    # Functions return the value of their last statement
    x + y
end

add(5, 6) # => 11 after printing out "x is 5 and y is 6"

# Compact assignment of functions
f_add(x, y) = x + y # => "f (generic function with 1 method)"
f_add(3, 4) # => 7

# Function can also return multiple values as tuple
f(x, y) = x + y, x - y
f(3, 4) # => (7, -1)

# You can define functions that take a variable number of
# positional arguments
function varargs(args...)
    return args
    # use the keyword return to return anywhere in the function
end
# => varargs (generic function with 1 method)

varargs(1,2,3) # => (1,2,3)

# The ... is called a splat.
# We just used it in a function definition.
# It can also be used in a function call,
# where it will splat an Array or Tuple's contents into the argument list.
add([5,6]...) # this is equivalent to add(5,6)

x = (5,6)     # => (5,6)
add(x...)     # this is equivalent to add(5,6)

# You can define functions with optional positional arguments
function defaults(a,b,x=5,y=6)
    return "$a $b and $x $y"
end

defaults('h','g') # => "h g and 5 6"
defaults('h','g','j') # => "h g and j 6"
defaults('h','g','j','k') # => "h g and j k"
try
    defaults('h') # => ERROR: no method defaults(Char,)
    defaults() # => ERROR: no methods defaults()
catch e
    println(e)
end

# You can define functions that take keyword arguments
function keyword_args(;k1=4,name2="hello") # note the ;
    return Dict("k1"=>k1,"name2"=>name2)
end

keyword_args(name2="ness") # => ["name2"=>"ness","k1"=>4]
keyword_args(k1="mine") # => ["k1"=>"mine","name2"=>"hello"]
keyword_args() # => ["name2"=>"hello","k1"=>4]

# You can combine all kinds of arguments in the same function
function all_the_args(normal_arg, optional_positional_arg=2; keyword_arg="foo")
    println("normal arg: $normal_arg")
    println("optional arg: $optional_positional_arg")
    println("keyword arg: $keyword_arg")
end

all_the_args(1, 3, keyword_arg=4)
# prints:
#   normal arg: 1
#   optional arg: 3
#   keyword arg: 4

# Julia has first class functions
function create_adder(x)
    adder = function (y)
        return x + y
    end
    return adder
end

# This is "stabby lambda syntax" for creating anonymous functions
(x -> x > 2)(3) # => true

# This function is identical to create_adder implementation above.
function create_adder(x)
    y -> x + y
end

# You can also name the internal function, if you want
function create_adder(x)
    function adder(y)
        x + y
    end
    adder
end

add_10 = create_adder(10)
add_10(3) # => 13

# There are built-in higher order functions
map(add_10, [1,2,3]) # => [11, 12, 13]
filter(x -> x > 5, [3, 4, 5, 6, 7]) # => [6, 7]

# We can use list comprehensions for nicer maps
[add_10(i) for i=[1, 2, 3]] # => [11, 12, 13]
[add_10(i) for i in [1, 2, 3]] # => [11, 12, 13]

####################################################
## 5. Types
####################################################

# Julia has a type system.
# Every value has a type; variables do not have types themselves.
# You can use the `typeof` function to get the type of a value.
typeof(5) # => Int64

# Types are first-class values
typeof(Int64) # => DataType
typeof(DataType) # => DataType
# DataType is the type that represents types, including itself.

# Types are used for documentation, optimizations, and dispatch.
# They are not statically checked.

# Users can define types
# They are like records or structs in other languages.
# New types are defined using the `type` keyword.

# type Name
#   field::OptionalType
#   ...
# end
type Tiger
  taillength::Float64
  coatcolor # not including a type annotation is the same as `::Any`
end

# The default constructor's arguments are the properties
# of the type, in the order they are listed in the definition
tigger = Tiger(3.5,"orange") # => Tiger(3.5,"orange")

# The type doubles as the constructor function for values of that type
sherekhan = typeof(tigger)(5.6,"fire") # => Tiger(5.6,"fire")

# These struct-style types are called concrete types
# They can be instantiated, but cannot have subtypes.
# The other kind of types is abstract types.

# abstract Name
abstract Cat # just a name and point in the type hierarchy

# Abstract types cannot be instantiated, but can have subtypes.
# For example, Number is an abstract type
subtypes(Number) # => 2-element Array{Any,1}:
                 #     Complex{T<:Real}
                 #     Real
subtypes(Cat) # => 0-element Array{Any,1}

# AbstractString, as the name implies, is also an abstract type
subtypes(AbstractString)    # 8-element Array{Any,1}:
                            #  Base.SubstitutionString{T<:AbstractString}
                            #  DirectIndexString
                            #  RepString
                            #  RevString{T<:AbstractString}
                            #  RopeString
                            #  SubString{T<:AbstractString}
                            #  UTF16String
                            #  UTF8String

# Every type has a super type; use the `super` function to get it.
typeof(5) # => Int64
super(Int64) # => Signed
super(Signed) # => Integer
super(Integer) # => Real
super(Real) # => Number
super(Number) # => Any
super(super(Signed)) # => Real
super(Any) # => Any
# All of these type, except for Int64, are abstract.
typeof("fire") # => ASCIIString
super(ASCIIString) # => DirectIndexString
super(DirectIndexString) # => AbstractString
# Likewise here with ASCIIString

# <: is the subtyping operator
type Lion <: Cat # Lion is a subtype of Cat
  mane_color
  roar::AbstractString
end

# You can define more constructors for your type
# Just define a function of the same name as the type
# and call an existing constructor to get a value of the correct type
Lion(roar::AbstractString) = Lion("green",roar)
# This is an outer constructor because it's outside the type definition

type Panther <: Cat # Panther is also a subtype of Cat
  eye_color
  Panther() = new("green")
  # Panthers will only have this constructor, and no default constructor.
end
# Using inner constructors, like Panther does, gives you control
# over how values of the type can be created.
# When possible, you should use outer constructors rather than inner ones.

####################################################
## 6. Multiple-Dispatch
####################################################

# In Julia, all named functions are generic functions
# This means that they are built up from many small methods
# Each constructor for Lion is a method of the generic function Lion.

# For a non-constructor example, let's make a function meow:

# Definitions for Lion, Panther, Tiger
function meow(animal::Lion)
  animal.roar # access type properties using dot notation
end

function meow(animal::Panther)
  "grrr"
end

function meow(animal::Tiger)
  "rawwwr"
end

# Testing the meow function
meow(tigger) # => "rawwr"
meow(Lion("brown","ROAAR")) # => "ROAAR"
meow(Panther()) # => "grrr"

# Review the local type hierarchy
issubtype(Tiger,Cat) # => false
issubtype(Lion,Cat) # => true
issubtype(Panther,Cat) # => true

# Defining a function that takes Cats
function pet_cat(cat::Cat)
  println("The cat says $(meow(cat))")
end

pet_cat(Lion("42")) # => prints "The cat says 42"
try
    pet_cat(tigger) # => ERROR: no method pet_cat(Tiger,)
catch e
    println(e)
end

# In OO languages, single dispatch is common;
# this means that the method is picked based on the type of the first argument.
# In Julia, all of the argument types contribute to selecting the best method.

# Let's define a function with more arguments, so we can see the difference
function fight(t::Tiger,c::Cat)
  println("The $(t.coatcolor) tiger wins!")
end
# => fight (generic function with 1 method)

fight(tigger,Panther()) # => prints The orange tiger wins!
fight(tigger,Lion("ROAR")) # => prints The orange tiger wins!

# Let's change the behavior when the Cat is specifically a Lion
fight(t::Tiger,l::Lion) = println("The $(l.mane_color)-maned lion wins!")
# => fight (generic function with 2 methods)

fight(tigger,Panther()) # => prints The orange tiger wins!
fight(tigger,Lion("ROAR")) # => prints The green-maned lion wins!

# We don't need a Tiger in order to fight
fight(l::Lion,c::Cat) = println("The victorious cat says $(meow(c))")
# => fight (generic function with 3 methods)

fight(Lion("balooga!"),Panther()) # => prints The victorious cat says grrr
try
  fight(Panther(),Lion("RAWR")) # => ERROR: no method fight(Panther,Lion)
catch
end

# Also let the cat go first
fight(c::Cat,l::Lion) = println("The cat beats the Lion")
# => Warning: New definition
#    fight(Cat,Lion) at none:1
# is ambiguous with
#    fight(Lion,Cat) at none:2.
# Make sure
#    fight(Lion,Lion)
# is defined first.
#fight (generic function with 4 methods)

# This warning is because it's unclear which fight will be called in:
fight(Lion("RAR"),Lion("brown","rarrr")) # => prints The victorious cat says rarrr
# The result may be different in other versions of Julia

fight(l::Lion,l2::Lion) = println("The lions come to a tie")
fight(Lion("RAR"),Lion("brown","rarrr")) # => prints The lions come to a tie

# Under the hood
# You can take a look at the llvm  and the assembly code generated.

square_area(l) = l * l      # square_area (generic function with 1 method)

square_area(5) #25

# What happens when we feed square_area an integer?
code_native(square_area, (Int32,))
	#	    .section    __TEXT,__text,regular,pure_instructions
	#	Filename: none
	#	Source line: 1              # Prologue
	#	    push    RBP
	#	    mov RBP, RSP
	#	Source line: 1
	#	    movsxd  RAX, EDI        # Fetch l from memory?
	#	    imul    RAX, RAX        # Square l and store the result in RAX
	#	    pop RBP                 # Restore old base pointer
	#	    ret                     # Result will still be in RAX

code_native(square_area, (Float32,))
	#	    .section    __TEXT,__text,regular,pure_instructions
	#	Filename: none
	#	Source line: 1
	#	    push    RBP
	#	    mov RBP, RSP
	#	Source line: 1
	#	    vmulss  XMM0, XMM0, XMM0  # Scalar single precision multiply (AVX)
	#	    pop RBP
	#	    ret

code_native(square_area, (Float64,))
	#	    .section    __TEXT,__text,regular,pure_instructions
	#	Filename: none
	#	Source line: 1
	#	    push    RBP
	#	    mov RBP, RSP
	#	Source line: 1
	#	    vmulsd  XMM0, XMM0, XMM0 # Scalar double precision multiply (AVX)
	#	    pop RBP
	#	    ret
	#
# Note that julia will use floating point instructions if any of the
# arguments are floats.
# Let's calculate the area of a circle
circle_area(r) = pi * r * r     # circle_area (generic function with 1 method)
circle_area(5)                  # 78.53981633974483

code_native(circle_area, (Int32,))
	#	    .section    __TEXT,__text,regular,pure_instructions
	#	Filename: none
	#	Source line: 1
	#	    push    RBP
	#	    mov RBP, RSP
	#	Source line: 1
	#	    vcvtsi2sd   XMM0, XMM0, EDI          # Load integer (r) from memory
	#	    movabs  RAX, 4593140240              # Load pi
	#	    vmulsd  XMM1, XMM0, QWORD PTR [RAX]  # pi * r
	#	    vmulsd  XMM0, XMM0, XMM1             # (pi * r) * r
	#	    pop RBP
	#	    ret
	#

code_native(circle_area, (Float64,))
	#	    .section    __TEXT,__text,regular,pure_instructions
	#	Filename: none
	#	Source line: 1
	#	    push    RBP
	#	    mov RBP, RSP
	#	    movabs  RAX, 4593140496
	#	Source line: 1
	#	    vmulsd  XMM1, XMM0, QWORD PTR [RAX]
	#	    vmulsd  XMM0, XMM1, XMM0
	#	    pop RBP
	#	    ret
	#

max_stars_repo_name: MateoCerquetella/dark-plus-syntax

max_issues_repo_name: MateoCerquetella/dark-plus-syntax

max_forks_repo_name: MateoCerquetella/dark-plus-syntax


text:
using Base.Test
using JLab
import PyPlot.plt

# #############################################################################
#Testing "getfit"
println("########## Testing \"getfit\" ##########")

#Sample data
x = linspace(0,10,100)
y = 1 + 3*x
y0 = y + 0.1*randn(100)
y0err = fill(0.1,length(y))
#y1 = y + (0.01*x+0.1).*randn(100)
y1 = [0.960833039028134,1.407268029451956,1.7041362562849804,1.7439408410837736,2.408365436088803,2.3558604670692187,2.6911766100562473,2.985346937471822,3.4325552973746216,3.6518536870040856,3.8198571085348556,4.249638029398669,4.5191875756320465,4.869643189063813,5.200098855165331,5.6028720774344825,5.8669985679112155,6.13446639319173,6.6589983913598445,6.583599635749586,7.0643398456691475,7.22881263766044,7.766208155795485,7.9279248867527095,8.031876149276476,8.622312209495371,8.91270431877012,9.253120846258275,9.46051424213136,9.864708345572803,10.196800366860534,10.262262734504358,10.703336679844554,11.159367909410706,11.448384037553025,11.419407340631926,11.882087020095623,12.123486147998683,12.589137110943783,12.915894183324713,13.195867464426314,13.381066468422823,13.498128075316496,13.946806258004903,14.140720490641081,14.65645590142797,14.92980075431186,15.422243674966104,15.606744487089607,15.664750441796148,16.035908118733797,16.4953299199874,16.958066212026043,16.994285363697312,17.354783816298067,17.700581921115532,18.0306766566351,18.382659996933242,18.231343132691336,19.063032842217314,19.255377262197005,19.81437282587234,19.956484257580982,19.97657224631401,20.360730031763094,20.65882754494313,21.106841540004787,21.252496222358445,21.766760201336467,22.19490679290051,22.6047978108247,22.529362909631008,22.915135996360842,23.17387598157397,23.43698104498548,23.598722726534746,23.73428978953588,24.241320075552377,24.65076910198924,24.78716650508645,25.05422998129448,25.58384488593317,25.93224573773362,25.843980681757863,26.45408533700855,26.888994042344866,27.25415441864504,27.83797148019114,27.550899599039617,28.03717936818957,28.18093241086242,28.548423792990217,28.591866793393848,29.208822396144086,29.49443477158012,29.793649140093873,29.867114283987977,30.526713526135133,30.7660738194201,30.89530630541283]
y1err = 0.01*x + 0.1

#Model
f(x,p) = p[1] + p[2]*x

#Get fit
fit = getfit(f,x,y,ones(Float64,length(y)),[0.,0.])
fit0 = getfit(f,x,y0,y0err,[0.,0.])
fit1 = getfit(f,x,y1,y1err,[0.,0.])

#Fit fields
println("Testing JLabFitresult fields ...")
@test fit0.xdata == x
@test fit0.ydata == y0
@test fit0.yerrors == y0err
@test fit0.model == f
@test fit0.dof == 98

#Fit accuracy
println("Testing accuracy ...")
@test_approx_eq_eps(fit.params[1],1,1e-10)
@test_approx_eq_eps(fit.params[2],3,1e-10)

p0 = sgolaycoeff(x,1)*y0
@test_approx_eq_eps(fit0.params[1],p0[1],1e-7)
@test_approx_eq_eps(fit0.params[2],p0[2],1e-10)

@test_approx_eq_eps(fit1.params[1],1,fit1.perrors[1])
@test_approx_eq_eps(fit1.params[2],3,fit1.perrors[2])

println("Function \"getfit\" passes tests.")


text:
using Clang.cindex

DAILOC="/home/bana/GSP/research/samc/samcnet/deps/libdai/include"

function wrap_dai_file(
	fname::ASCIIString, namespace::ASCIIString, class::ASCIIString, ofname::ASCIIString,
	funcs::Vector{Symbol}, exclude::Vector{Symbol}=Symbol[]
	)
	topcu = cindex.parse_header("$DAILOC/dai/$fname"; cplusplus = true, args = ["-I$DAILOC"])

	namespace = cindex.search(topcu, namespace)
	cls = mapreduce(x->cindex.search(x, class), vcat, namespace)

	of = open("$(ofname).cc", "w")
	ofjl = open("$(ofname).jl", "w")
	println(of, "//c++ -shared -fPIC $(ofname).cc -o $(ofname)")
	println(of, "#include <dai/$fname>")
	println(of, "extern \"C\" {")
	#println(map(cindex.spelling,children(topcu)))
	for c in cls
		length(children(c)) == 0 && continue
		println(name(c))
		#println(length(children(c)))
		#println(cindex.spelling(children(c)[1]))

		MethodCount = Dict{Symbol, Int}()
		for decl in children(c)
			const declname = symbol(strip(spelling(decl)))
			wrapped = Any[]
			if (declname in funcs) && ( (isa(decl, cindex.CXXMethod)) || isa(decl, cindex.Constructor))

				#check for public
				Clang.cindex.getCXXAccessSpecifier(decl)==1 || continue

			    id = (MethodCount[declname] = get(MethodCount, declname, 0) + 1)
			    if (symbol("$(declname)$(id)") in exclude)
			    	println("excluding $(name(c))::$(declname)::$(id)")
			    	continue
			    end
			    Clang.wt.wrap(of, decl, id)
			    Clang.wt.wrapjl(ofjl, decl, id)
					#Clang.wt.wrapjl(ofjl, ":libdai",decl, id)
			end
		end
	end
	println(of, "} //extern C")
	close(of)
	close(ofjl)
end

mkpath("gen")
wrap_dai_file("var.h", "dai", "Var", "gen/var", [:states, :label])


text:
struct ScaledArray{T<:AbstractFloat, N} <: AbstractArray{T, N}
    A::AbstractArray{T,N}
    s::T
end

Base.size(SA::ScaledArray{T,N}) where {T,N} = size(SA.A)
Base.getindex(SA::ScaledArray,i...) = getindex(SA.A,i...)
Base.setindex!(SA::ScaledArray,x::Real,i...) = setindex!(SA.A,x,i...)

ScaledArray(A::AbstractArray) = ScaledArray(A,one(eltype(A)))
ScaledArray(A::AbstractArray{T},s::Real) where T = ScaledArray(A,T(s))

unscale(SA::ScaledArray{T,N}) where {T,N} = ScaledArray(SA.A*SA.s,one(T))
rescale(SA::ScaledArray{T,N},s::Real) where {T,N} = ScaledArray(SA.A*(SA.s/T(s)),T(s))

Base.:(*)(X::ScaledArray{T,N},Y::ScaledArray{T,N}) where {T,N} = ScaledArray(X.A*Y.A,X.s*Y.s)


text:
module Antennae
using Swizzles.Properties

export Antenna

using Base.Broadcast: broadcasted

struct Antenna{F}
    f::F
end

(a::Antenna)(args...) = broadcasted(a.f, args...)

@inline Properties.initial(a::Antenna{F}, T) where {F} = Ref(initial(a.f, eltype(T)))

end


text:
#__precompile__()
module AnalyticDomains

# Curve parametrization
struct CurveDesc
    tau
    dtau
    d2tau
    normal
    tangent
end

export starfish, CurveDesc

using SymPy

# Use any curve parametrized in t in [0, 2 pi), described as string
# Ex. desc = paramcurve("(1 + 0.3*cos(5*t))*exp(I*t)")
# !Note that imaginary unit is I
function paramcurve(str::String)
    desc = Sym(str)
    t = symbols("t")
    # Differentiate curve
    diff1 = diff(desc, t)
    diff2 = diff(diff1, t)
    tang = diff1 / abs(diff1)
    normal = 1im*tang
    return CurveDesc(
        lambdify(desc),
        lambdify(diff1),
        lambdify(diff2),
        lambdify(normal),
        lambdify(tang)
    )
end

function starfish(;n_arms = 5, amplitude = 0.3, radius = 1.0, center=(0,0), rotation=0.0, interior=false)
    if interior
        s = -1.0
    else
        s = 1.0
    end
    R = exp(1im*rotation)
    desc(t) = R*radius*((1 + amplitude*cos(n_arms*t)).*exp(1im*s*t)) + (center[1] + 1im*center[2])
    diff1(t) = R*radius*(exp(t*1im*s)*(amplitude*cos(n_arms*t) + 1)*1im*s - amplitude*n_arms*exp(t*1im*s)*sin(n_arms*t))
    diff2(t) = R*radius*(- exp(t*1im*s)*(amplitude*cos(n_arms*t) + 1) - amplitude*n_arms*exp(t*1im*s)*sin(n_arms*t)*2*1im*s - amplitude*n_arms^2*exp(t*1im*s)*cos(n_arms*t))
    tang(t) = diff1(t) / abs(diff1(t))
    normal(t) = 1im*tang(t)
    return CurveDesc(
        desc,
        diff1,
        diff2,
        normal,
        tang
    )
end

function rectangle(;width = 1, height = 1, order = 10, interior=false)
    @assert iseven(order)
    a = width
    b = height
    p = order;
    desc(t) = (a*cos(t) + b*sin(t)*1im)/(cos(t)^p + sin(t)^p)^(1/p)
    _diff1(t) = -(a*cos(t)*sin(t)^p - b*cos(t)^p*sin(t)*1im)/(cos(t)*sin(t)*(cos(t)^p + sin(t)^p)^((p + 1)/p))
    _diff2(t) = (-(a*cos(t)^(p + 1)*sin(t)^p + b*cos(t)^(2*p + 2)*sin(t)*2im - 2*a*cos(t)^3*sin(t)^(2*p) - 2*a*cos(t)^(p + 3)*sin(t)^p - b*cos(t)^p*sin(t)^(p + 1)*1im - b*cos(t)^(2*p)*sin(t)*2im + b*cos(t)^(p + 2)*sin(t)^(p + 1)*2im + a*p*cos(t)^(p + 1)*sin(t)^p + b*p*cos(t)^p*sin(t)^(p + 1)*1im)/(cos(t)^2*sin(t)^2*(cos(t)^p + sin(t)^p)^((2*p + 1)/p)))

    # Fudge because derivative expression is NaN at t=0
    diff1(t) = _diff1(t + eps(t))
    diff2(t) = _diff2(t + eps(t))

    tang(t) = diff1(t) / abs(diff1(t))
    normal(t) = 1im*tang(t)

    return CurveDesc(
        desc,
        diff1,
        diff2,
        normal,
        tang
    )
end

end


text:
# file = "clusterIds_temporal-(perc_50%,1980-2015).tif"
# x = GeoArrays.read(file);

# ds_raster = GDAL.gdalopen(raster_file, GDAL.GA_ReadOnly)
# band = GDAL.gdalgetrasterband(ds, band)
# GDAL.gdalsetrasternodatavalue(band.ptr, -1)
file = "OUTPUT/clusterIds_temporal-(perc_20%,1980-2015).tif"
@time x = GeoArrays.read(file, 1:1)
# @time gdal_polygonize(file, 1, "OUTPUT/shp_UrbanSpatioTemporalCluster/urban.shp"; fieldname = "urbanId")


text:
using OrdinaryDiffEq
using Flux
using Flux: logitcrossentropy
using MLDatasets: MNIST
using MLDataUtils
using DiffEqFlux
using CuArrays
using NNlib

function loadmnist(batchsize=bs)
	# Use MLDataUtils LabelEnc for natural onehot conversion
    onehot(labels_raw) = convertlabel(LabelEnc.OneOfK, labels_raw, LabelEnc.NativeLabels(collect(0:9)))
	# Load MNIST
	imgs, labels_raw = MNIST.traindata();
	# Process images into (H,W,C,BS) batches
	x_train = reshape(imgs,size(imgs,1),size(imgs,2),1,size(imgs,3))|>gpu
	x_train = batchview(x_train,batchsize);
	# Onehot and batch the labels
	y_train = onehot(labels_raw)|>gpu
	y_train = batchview(y_train,batchsize)
	return x_train, y_train
end

# Main
const bs = 128
x_train, y_train = loadmnist(bs)

down = Chain(
             x->reshape(x,(28*28,:)),
             Dense(784,20,tanh)
            )|>gpu
nfe=0
nn = Chain(
           #= x->(global nfe+=1;x), =#
           Dense(20,10,tanh),
           Dense(10,10,tanh),
           Dense(10,20,tanh)
          ) |>gpu
fc = Chain(
           Dense(20,10)
          )|>gpu

CuArrays.allowscalar(false)
nn_ode(x) = neural_ode(nn,x,gpu((0.f0,1.f0)), Tsit5(),save_everystep=false,reltol=1e-3,abstol=1e-3, save_start=false)
m = Chain(down,nn_ode,fc)
m_no_ode = Chain(down,nn,fc)

x_d = down(x_train[1])
nn_ode(x_d)

## Stop here

# Showing this works
x_m = m(x_train[1])

classify(x) = argmax.(eachcol(x))

function accuracy(model,data; n_batches=100)
    total_correct = 0
    total = 0
    for (x,y) in collect(data)[1:n_batches]
        target_class = classify(cpu(y))
        predicted_class = classify(cpu(model(x)))
        total_correct += sum(target_class .== predicted_class)
        total += length(target_class)
    end
    return total_correct/total
end

accuracy(m, zip(x_train,y_train))

function loss(x,y)
    y_hat = m(x)
    return logitcrossentropy(y_hat,y)
end

loss(x_train[1],y_train[1])

opt=ADAM()
iter = 0
cb() = begin
    global iter += 1
    @show iter
    @show nfe
    @show loss(x_train[1],y_train[1])
    @show cpu(down)[2].W[1]
    if iter%10 == 0
        @show accuracy(m, zip(x_train,y_train))
    end
    global nfe=0
end

Flux.train!(loss,params(down,nn,fc),zip(x_train,y_train),opt, cb=cb)

# Saving doesn't work yet
using BSON: @save, @load

m_cpu = cpu(m)
@save "saved_models/mnist_node.bson" m_cpu

@load "saved_models/mnist_node.bson" m_cpu

max_stars_repo_name: jessebett/model-zoo

max_issues_repo_name: jessebett/model-zoo

max_forks_repo_name: jessebett/model-zoo


text:
function ht21(A,debug=false;itmax=20)
    n,m = size(A)
    T = eltype(A)
    Tr = real(T)
    if n != m
        error("norm 1 estimation can only be applied to square matrices")
    end
    x = ones(T,n)/n
    y, z, ξ = similar(x), similar(x), similar(x)
    γ = 0
    itcount = itmax
    for k in 1:itmax
        y = mul!(y, A, x) #y = A*x
        ξ .= sign.(y) #ξ = sign(y)
        z = mul!(z, A', ξ) #z = A'*ξ
        z_dot_x = dot(z, x)
        if debug
            # z_dot_x should be real to within floating point errors
            @show z_dot_x
            @assert abs(imag(z_dot_x)) < 10*eps(real(z_dot_x))
        end
        if norm(z,Inf) <= real(z_dot_x) && k > 1
            γ = norm(y, 1)
            itcount = k
            break
        end
        x = fill!(x, 0) #x = zeros(Float64,n);
        x[lazyfindmax(abs, z)[2]] = 1 #x[indmax(abs(z))] = 1.
    end

    # b = Float64[(-1)^(i+1) * (1+(i-1)/(n-1)) for i in 1:n]
    # out = max(γ,norm(A*b,1)/norm(b,1))
    for i in 1:n
        sgn = ifelse(iseven(i+1), 1, -1)
        x[i] = sgn * (1 + Tr(i-1) / Tr(n-1))
    end
    y = mul!(y, A, x)
    est = max(γ, norm(y,1) / norm(x,1))

    return est, itcount
end

function ht22(A,t::Int=2,debug=false;itmax=10)
    n,m = size(A)
    T = eltype(A)
    Tr = real(T)
    if n!=m
        error("norm 1 estimation can only be applied to square matrices")
    end
    # start with a random matrix of normalized columns
    X = randn(T,n,t)
    Y, Z = similar(X), similar(X)
    for col=1:t
        # scale!(slice(X,1:n,col),1/norm(slice(X,1:n,col),1))
        @views X[:,col] .*= 1/norm(X[:,col],1)
    end
    # Id = Array(I,n,n) #Id = eye(n)

    h = zeros(Tr,n)
    itcount = itmax
    for k in 1:itmax
        Y = mul!(Y,A,X) #Y = A*X
        @views g = Tr[norm(Y[:,j],1) for j in 1:t]
        _, ind_best = findmax(g)
        g = sort!(g; rev=true)
        if debug
            @show g[1]
        end
        S = sign.(Y)
        Z = mul!(Z,A',S) #Z = A'*S
        @views h = Tr[maximum(abs, Z[i,:]) for i in 1:n] #Float64[maximum(abs(vec(Z[i,:]))) for i in 1:n]
        ind = [1:n;]
        h_ind = sortslices(hcat(h,typeof(h)(ind)); dims=1, by=first, rev=true) #sortrows(hcat(h,ind),by=x->x[1],rev=true)
        h = copyto!(h, h_ind[:,1])
        ind = round.(Int, h_ind[:,2])

        @views Z_dot_X = dot(Z[:,ind_best], X[:,ind_best])
        if debug
            # Z_dot_X should be real to within floating point errors
            @show Z_dot_X
            @assert abs(imag(Z_dot_X)) < 10*eps(real(Z_dot_X))
        end
        if maximum(h) <= real(Z_dot_X)
            itcount = k
            break
        end
        if debug
            @show h[1]
        end
        for j in 1:t
            @views fill!(X[:,j], 0) # X[:,j] = Id[:,ind[j]]
            X[ind[j],j] = 1
        end
    end
    est = h[1]

    return est, itcount
end

function ht23(A,t::Int=2;itmax=2)
    n,m = size(A)
    T = eltype(A)
    Tr = real(T)
    if n!=m
        error("norm 1 estimation can only be applied to square matrices")
    end
    # start with a random matrix of normalized columns
    X = randn(T,n,t)
    Y, Z = similar(X), similar(X)
    for col=1:t
        @views X[:,col] .*= 1/norm(X[:,col],1)
    end
    # Id = eye(n)
    est_old = zero(Tr) #0 # initial old estimate
    est = zero(Tr) #0 # initial estimate
    # ind = zeros(Int,n)
    S = zeros(T,n,t) #S = zeros(Int8,n,t)
    S_old = zeros(T,n,t) #S_old = zeros(n,t)
    ind_hist = [] # integer vector recording indices of used unit vectors
    itcount = itmax
    for k=1:itmax+1
        Y = mul!(Y,A,X) #Y = A*X # TODO: make Y preallocated, and in place multiplication
        # est, est_indx = findmax([norm(slice(Y,1:n,col),1) for col in 1:t])
        est, est_indx = lazyfindmax(y->norm(y,1), view(Y,:,col) for col in 1:t)
        if est > est_old || k == 2
            ind_best = est_indx
            # w = Y[:,ind_best]
        end
        if est <= est_old && k >= 2
            est = est_old
            itcount = k
            break
        end
        est_old = est
        S_old = copyto!(S_old, S)
        if k > itmax
            break
        end
        S .= sign.(Y) #copy!(S,sign(Y))
        # If every column of S is parallel to a column of S_old, break
        if parallel_cols(S, S_old)
            itcount = k
            break
        end
        if t>1
            randomize_parallel!(S,S_old)
        end
        Z = mul!(Z,A',S) #Z = transpose(A)*S
        @views h = Tr[norm(Z[i,:], Inf) for i in 1:n]
        if k>=2 && maximum(h) == h[ind_best]
            itcount = k
            break
        end
        # Sort h in decreasing order, and reorder ind correspondingly
        ind = [1:n;]
        h_ind = sortslices(hcat(h,typeof(h)(ind)); dims=1, by=first, rev=true) #sortrows(hcat(h,ind),by=x->x[1],rev=true)
        h = h_ind[:,1]
        #ind = integer(h_ind[:,2])
        ind = round.(Int, h_ind[:,2]) #ind = @compat round(Integer,h_ind[:,2])
        if t>1
            if ind[1:t] ⊆ ind_hist
                itcount = k
                break
            end
            # replace ind[1:t] with the first t indices in
            # ind[1:n] that are not in ind_hist
            ind[1:t] = filter(x->!(x ⊆ ind_hist),ind)[1:t]
        end
        for j in 1:t
            @views fill!(X[:,j], 0) # X[:,j] = Id[:,ind[j]]
            X[ind[j],j] = 1
        end
        ind_hist = [ind_hist; ind[1:t]]
    end

    return est, itcount
end

function parallel_cols(S,So)
    n = size(S,1)
    # the assumption is that S,So are matrices with elements ±1,
    # so being parallel would translate to abs value of inner product ≈ n
    return any(map(x->isapprox(x,n), abs.(So'*S)); dims = 1) |> all
end

function randomize_parallel!(S,So)
    n = size(S,1)
    # Ensure no column of S is parallel to a column of S
    ips = map(x->isapprox(x,n), abs.(S'*S))
    for i = 1:size(ips,1) #1:n
        ips[i,i] = false
        if any(ips[:,i])
            @views S[:,i] .= rand([-1,1], n)
        end
    end
    # Ensure no column of S is parallel to a column of So
    ips = map(x->isapprox(x,n), abs.(So'*S))
    for i = 1:size(ips,1) #1:n
        ips[i,i] = false
        if any(ips[:,i])
            @views S[:,i] .= rand([-1,1], n)
        end
    end
end

# using BenchmarkTools
# using LinearMaps
#
# function print_ests(A,t,debug=false)
#     # println("ht21: $(ht21(A,debug)[1])")
#     # println("ht22: $(ht22(A,t,debug)[1])")
#     # println("ht23: $(ht23(A,t)[1])")
#     # println("norm: $(opnorm(Array(A),1))")
#     println("ht21: $(@btime ht21($A,$debug)[1] evals = 3)")
#     println("ht22: $(@btime ht22($A,$t,$debug)[1] evals = 3)")
#     println("ht23: $(@btime ht23($A,$t)[1] evals = 3)")
#     println("norm: $(@btime opnorm(Array($A),1) evals = 3)")
#     return nothing
# end
#
# function test_normest1(t=2; d=100, p=6/d)
#     println("Full...")
#     A = randn(Float64,d,d)
#     ests = print_ests(A,t)
#
#     println("Complex Full...")
#     A = randn(ComplexF64,d,d)
#     ests = print_ests(A,t)
#
#     println("Sparse...")
#     A = sprandn(Float64,d,d,p)
#     ests = print_ests(A,t)
#
#     println("Complex Sparse...")
#     A = sprandn(ComplexF64,d,d,p)
#     ests = print_ests(A,t)
#
#     println("Complex LinearMaps...")
#     A = randn(ComplexF64,d,d) |> LinearMap
#     ests = print_ests(A,t)
#
#     return nothing
# end


text:
module FactoriesTests

using Test
using UnstructuredGrids.Core
using UnstructuredGrids.Kernels
using UnstructuredGrids.Factories
using UnstructuredGrids.RefCellGallery
using UnstructuredGrids.VTK

grid = UGrid(domain=(0,1,-1,0),partition=(2,2))

@test isa(grid,UGrid)

p = [0.0 1.0 2.0 0.0 1.0 2.0 0.0 1.0 2.0;
     -1.0 -1.0 -1.0 -0.5 -0.5 -0.5 0.0 0.0 0.0]

clist = [1, 2, 4, 5, 2, 3, 5, 6, 4, 5, 7, 8, 5, 6, 8, 9]
cptrs = [1, 5, 9, 13, 17]

@test coordinates(grid) ≈ p
@test list(connections(grid)) == clist
@test ptrs(connections(grid)) == cptrs
@test celltypes(grid) == ones(Int,length(cptrs)-1)

grid = UGrid(domain=(0,1,-1,0,2,4),partition=(1,2,3))

p = [0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0;
     -1.0 -1.0 -0.5 -0.5 0.0 0.0 -1.0 -1.0 -0.5 -0.5 0.0 0.0 -1.0 -1.0 -0.5 -0.5 0.0 0.0 -1.0 -1.0 -0.5 -0.5 0.0 0.0;
     2.0 2.0 2.0 2.0 2.0 2.0 2.5 2.5 2.5 2.5 2.5 2.5 3.0 3.0 3.0 3.0 3.0 3.0 3.5 3.5 3.5 3.5 3.5 3.5]

clist = [
  1, 2, 3, 4, 7, 8, 9, 10, 3, 4, 5, 6, 9, 10, 11,
  12, 7, 8, 9, 10, 13, 14, 15, 16, 9, 10, 11, 12,
  15, 16, 17, 18, 13, 14, 15, 16, 19, 20, 21, 22,
  15, 16, 17, 18, 21, 22, 23, 24]

cptrs = [1, 9, 17, 25, 33, 41, 49]

@test coordinates(grid) ≈ p
@test list(connections(grid)) == clist
@test ptrs(connections(grid)) == cptrs
@test celltypes(grid) == ones(Int,length(cptrs)-1)

grid = UGrid(domain=(0,1,0,1),partition=(2,2))

ltcell_to_lpoints = [[1,2,3],[4,3,2]]

tgrid = UGrid(grid,ltcell_to_lpoints,TRIANGLE)

@test coordinates(grid) === coordinates(tgrid)

clist = [1, 2, 4, 5, 4, 2, 2, 3, 5, 6, 5, 3, 4, 5, 7, 8, 7, 5, 5, 6, 8, 9, 8, 6]
cptrs = [1, 4, 7, 10, 13, 16, 19, 22, 25]

@test list(connections(tgrid)) == clist
@test ptrs(connections(tgrid)) == cptrs
@test celltypes(tgrid) == ones(Int,length(cptrs)-1)

ltcell_to_lpoints = [
 [7,3,2,1], [7,5,2,1], [7,4,3,2], [7,4,8,2], [7,6,5,2], [7,6,8,2]]

grid = UGrid(domain=(0,1,0,1,0,1),partition=(2,2,3))

tgrid = UGrid(grid,ltcell_to_lpoints,TETRAHEDRON)

D = 1

grid = UGrid(domain=(0,1,0,1),partition=(2,2))

ltcell_to_lnodes = [[1,2,3],[4,3,2]]
tgrid = UGrid(grid,ltcell_to_lnodes,TRIANGLE)

cell_to_faces = generate_cell_to_faces(D,grid)
tcell_to_tfaces = generate_cell_to_faces(D,tgrid)

refcell = grid.refcells[1]
reftcell = tgrid.refcells[1]

ltface_to_ltnodes = reftcell.faces[D+1]
lface_to_lnodes = refcell.faces[D+1]

ntfaces = maximum(tcell_to_tfaces.list)

tface_to_face = generate_tface_to_face(
  cell_to_faces.list,
  cell_to_faces.ptrs,
  tcell_to_tfaces.list,
  tcell_to_tfaces.ptrs,
  ltcell_to_lnodes,
  ltface_to_ltnodes,
  lface_to_lnodes,
  ntfaces)

@test tface_to_face == [1, 0, 3, 2, 4, 5, 0, 6, 7, 0, 9, 8, 10, 0, 11, 12]

end # module FactoriesTests