text:
module SpecialGraphs

import LightGraphs
const LG = LightGraphs

using LightGraphs: nv, ne, outneighbors,
                   inneighbors, vertices, edges, Edge,
                   has_vertex, has_edge

export WheelGraph, PathGraph, CompleteGraph

include("WheelGraph.jl")
include("PathGraph.jl")
include("CompleteGraph.jl")

end # module


text:
cd(@__DIR__)

using Pkg

using Test, TestSetExtensions, SafeTestsets
using Validator

@testset ExtendedTestSet "validator tests" begin
  @includetests ARGS #[(endswith(t, ".jl") && t[1:end-3]) for t in ARGS]
end

ext: jl


text:
#=
  Aircrat position.
  From BLOG example/aircraft-position.blog

  Distributions of variable numAircraft (num:0)
  3.00000 =>    5857  (0.585700)
  4.00000 =>    3178  (0.317800)
  5.00000 =>     833  (0.083300)
  6.00000 =>     114  (0.011400)
  7.00000 =>      18  (0.001800)

  Also see ~/webppl/aircraft_position.wppl
=#

using Turing, StatsPlots, DataFrames
# using ReverseDiff, Zygote, Tracker
# Turing.setadbackend(:reversediff)
# Turing.setadbackend(:zygote)
# Turing.setadbackend(:tracker)

include("jl_utils.jl")

@model function aircraft_position()
    epsilon = 0.05

    # Poisson might get 0 which is not good for this model
    pp ~ Poisson(5)
    numAircraft ~ Dirac(pp == 0 ? 1 : pp)

    blips = tzeros(numAircraft)
    positions = tzeros(numAircraft)
    obsPos = tzeros(numAircraft)
    for b in 1:numAircraft
        blips[b] ~ flip(0.9)
        positions[b] ~ Normal(0,10)
        obsPos[b] ~ Normal(positions[b],1)
    end

    sumBlips ~ Dirac(sum(blips))

    function inRange(x,y,epsilon)
        Dirac(x > y - epsilon && x < y + epsilon)
    end

    true ~ Dirac(sumBlips == 3)

    inRange(obsPos[1],5.0,epsilon)
    numAircraftPositionLt5 ~ Dirac(sum([positions[a] > 5 ? 1 : 0 for a in 1:numAircraft]))
end

model = aircraft_position()

num_chns = 4
# chns = sample(model, Prior(), 1000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, IS(), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, SMC(), MCMCThreads(), 10_000, num_chns)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns, :numAircraft)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
#!/usr/bin/env julia

using MPI
using PencilArrays
using Random
using Test

MPI.Init()

comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
Nproc = MPI.Comm_size(comm)

let dev_null = @static Sys.iswindows() ? "nul" : "/dev/null"
    rank == 0 || redirect_stdout(open(dev_null, "w"))
end

function test_pencil(pen)
    A = PencilArray{Float64}(undef, pen)
    G = global_view(A)
    randn!(A)
    perm = Tuple(permutation(A))

    @testset "Broadcast $(nameof(typeof(x)))" for x in (A, G)
        test_broadcast(x)
        let y = similar(x)
            broadcast!(+, y, x, x, 3)  # precompile before measuring allocations
            alloc = @allocated broadcast!(+, y, x, x, 3)
            if VERSION ≥ v"1.5"  # there are small allocations in Julia 1.3
                @test alloc == 0
            end
            @test y ≈ 2x .+ 3
        end
    end

    @testset "Combinations" begin
        # Combine with regular Array
        P = parent(A) :: Array
        P′ = perm === nothing ? P : PermutedDimsArray(P, perm)
        @test typeof(P′ .+ A) == typeof(A)

        # Combine PencilArray and GlobalPencilArray
        @test_throws ArgumentError A .+ G

        # Combine Array and GlobalPencilArray
        @test_throws ArgumentError P .+ G
    end
end

function test_broadcast(A)
    @test typeof(2A) == typeof(A)
    @test typeof(A .+ A) == typeof(A)
    @test typeof(A .+ A .+ 3) == typeof(A)
    @test parent(2A) == 2parent(A)
    nothing
end

topo = MPITopology(comm, (Nproc, ))

pencils = (
    "Non-permuted" => Pencil(topo, (11, 12), (2, )),
    "Permuted" => Pencil(topo, (11, 12), (2, ), permute=Permutation(2, 1)),
)

@testset "$s" for (s, pen) in pencils
    test_pencil(pen)
end


text:
module SpikeNetOpt

using PyCall
using OrderedCollections
using LinearAlgebra
using SpikeSynchrony
using SpikingNeuralNetworks
SNN = SpikingNeuralNetworks
SNN.@load_units
using Evolutionary
import DataStructures
using JLD
using UnicodePlots
using SpikeSynchrony
using Statistics
using JLD
using Distributed
using Plots
using UnicodePlots
using Distributions
using LightGraphs
using ClearStacktrace
using SharedArrays
using WaspNet
using SpikingNN
using Reexport

#using Plots

#using Metaheuristics

# Export Function names and struct names found in file utils
export get_vm
function get_vm end
function checkmodel end
include("utils.jl")
include("current_search.jl")
#include("sciunit.jl")

end


text:
## --- Load the StatGeochem package which has the resampling functions we'll want

    using StatGeochem
    using Plots; gr();
    using Statistics, DelimitedFiles, SpecialFunctions

## --- Download and unzip Keller and Schoene (2012) dataset

    if ~isfile("ign.h5") # Unless it already exists
        download("https://storage.googleapis.com/statgeochem/ign.h5.gz","./ign.h5.gz")
        download("https://storage.googleapis.com/statgeochem/err2srel.csv","./err2srel.csv")
        run(`gunzip -f ign.h5.gz`) # Unzip file
    end

    # Read HDF5 file
    using HDF5
    ign = h5read("ign.h5","vars")

## --- Compute proximity coefficients (inverse weights)

    # # Compute inverse weights
    # k = invweight(ign["Latitude"] .|> Float32, ign["Longitude"] .|> Float32, ign["Age"] .|> Float32)

    # Since this is somewhat computatually intensive, let's load a precomputed version instead
    k = ign["k"]

    # Probability of keeping a given data point when sampling
    p = 1.0 ./ ((k .* median(5.0 ./ k)) .+ 1.0) # Keep rougly one-fith of the data in each resampling
    p[vec(ign["Elevation"].<-100)] .= 0 # Consider only continental crust

    # Set absolute uncertainties for each element where possible, using errors defined inerr2srel.csv
    err2srel = importdataset("err2srel.csv", ',')
    for e in ign["elements"]
        # If there's an err2srel for this variable, create a "_sigma" if possible
        if haskey(err2srel, e) && !haskey(ign, e*"_sigma")
            ign[e*"_sigma"] = ign[e] .* (err2srel[e] / 2);
        end
    end

    # Special cases: age uncertainty
    ign["Age_sigma"] = (ign["Age_Max"]-ign["Age_Min"])/2;
    t = (ign["Age_sigma"] .< 50) .| isnan.(ign["Age_sigma"]) # Find points with < 50 Ma absolute uncertainty
    ign["Age_sigma"][t] .= 50 # Set 50 Ma minimum age uncertainty (1-sigma)

    # Special cases: location uncertainty
    ign["Latitude_sigma"] = ign["Loc_Prec"]
    ign["Longitude_sigma"] = ign["Loc_Prec"]

## --- Single element differentiation example

    xelem = "SiO2"
    xmin = 45
    xmax = 75
    nbins = 8
    elem = "K2O"

    h = plot(xlabel=xelem, ylabel="$(elem)",xlims=(xmin,xmax),framestyle=:box,grid=:off,fg_color_legend=:white) # Format plot

    rt = [0,1,2,3,4] # Time range (Ga)
    colors = reverse(resize_colormap(viridis[1:end-20],length(rt)-1))
    for i=1:length(rt)-1
        t = rt[i]*1000 .< ign["Age"] .< rt[i+1]*1000

        # Resample, returning binned means and uncertainties
        # (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)
        (c,m,el,eu) = bin_bsr_means(ign[xelem][t],ign[elem][t],xmin,xmax,nbins, p=p[t],
                        x_sigma=ign[xelem*"_sigma"][t], y_sigma=ign[elem*"_sigma"][t])

        # Plot results
        plot!(h, c,m,yerror=(el,eu),color=colors[i],mscolor=colors[i],seriestype=:scatter,label="$(rt[i])-$(rt[i+1]) Ga")
        plot!(h, c,m,style=:dot,color=colors[i],mscolor=colors[i],label="")
    end
    display(h)
## --- Ratio differentiation example

    xelem = "SiO2"
    xmin = 45
    xmax = 75
    nbins = 8
    num = "Rb" # Numerator
    denom = "Sr" # Denominator

    h = plot(xlabel=xelem, ylabel="$(num) / $(denom)",xlims=(xmin,xmax),framestyle=:box,grid=:off,legend=:topleft,fg_color_legend=:white) # Format plot

    rt = [0,1,2,3,4]
    colors = reverse(resize_colormap(viridis[1:end-20],length(rt)-1))
    for i=1:length(rt)-1
        t = rt[i]*1000 .< ign["Age"] .< rt[i+1]*1000

        # Resample, returning binned means and uncertainties
        # (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)
        (c,m,el,eu) = bin_bsr_ratios(nanbinmedian!, ign[xelem][t],ign[num][t],ign[denom][t],xmin,xmax,nbins, p=p[t],
                        x_sigma=ign[xelem*"_sigma"][t], num_sigma=ign[num*"_sigma"][t], denom_sigma=ign[denom*"_sigma"][t])

        # Plot results
        plot!(h, c,m,yerror=(el,eu),color=colors[i],mscolor=colors[i],seriestype=:scatter,label="$(rt[i])-$(rt[i+1]) Ga")
        plot!(h, c,m,style=:dot,color=colors[i],mscolor=colors[i],label="")
    end
    display(h)

## --- High-level functions for calculating combined averages over a set or silica ranges

    function constprop(binbsrfunction::Function, dataset::Dict, elem, xmin, xmax, nbins, p; xelem="Age", norm_by="SiO2", norm_bins=[43,55,65,78], nresamplings=1000)
        c = zeros(nbins)
        m = zeros(nbins)
        el = zeros(nbins)
        eu = zeros(nbins)
        for i=1:length(norm_bins)-1
            # Find the samples we're looking for
            t = (norm_bins[i] .< dataset[norm_by] .< norm_bins[i+1]) .& (dataset[elem] .> 0)

            # See what proportion of the modern crust falls in this norm_bin
            prop = sum((norm_bins[i] .< dataset[norm_by] .< norm_bins[i+1]) .& (p .> 0)) / sum((norm_bins[1] .< dataset[norm_by] .< norm_bins[end]) .& (p .> 0))

            # Resample, returning binned means and uncertainties
            # (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)
            (c[:],m1,el1,eu1) = binbsrfunction(dataset[xelem][t], dataset[elem][t], xmin, xmax, nbins, x_sigma=dataset["$(xelem)_sigma"][t], nresamplings=nresamplings, p=p[t])

            m .+= prop.*m1
            el .+= prop.*el1
            eu .+= prop.*eu1
        end
        el ./= sqrt(length(norm_bins)-1) # Standard error
        eu ./= sqrt(length(norm_bins)-1) # Standard error

        return c, m, el, eu
    end

    function constprop(binbsrfunction::Function, dataset::Dict, num, denom, xmin, xmax, nbins, p; xelem="Age", norm_by="SiO2", norm_bins=[43,55,65,78], nresamplings=1000)
        c = zeros(nbins)
        m = zeros(nbins)
        el = zeros(nbins)
        eu = zeros(nbins)
        for i=1:length(norm_bins)-1
            # Find the samples we're looking for
            t = (norm_bins[i] .< dataset[norm_by] .< norm_bins[i+1]) .& (dataset[num].>0) .& (dataset[denom].>0)

            # See what proportion of the modern crust falls in this norm_bin
            prop = sum((norm_bins[i] .< dataset[norm_by] .< norm_bins[i+1]) .& (p .> 0)) / sum((norm_bins[1] .< dataset[norm_by] .< norm_bins[end]) .& (p .> 0))

            # Resample, returning binned means and uncertainties
            # (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)
            (c[:],m1,el1,eu1) = binbsrfunction(dataset[xelem][t],dataset[num][t],dataset[denom][t],xmin,xmax,nbins,x_sigma=dataset["$(xelem)_sigma"][t],num_sigma=dataset["$(num)_sigma"][t],denom_sigma=dataset["$(denom)_sigma"][t],nresamplings=nresamplings,p=p[t])

            m .+= prop.*m1
            el .+= prop.*el1
            eu .+= prop.*eu1
        end
        el ./= sqrt(length(norm_bins)-1) # Standard error
        eu ./= sqrt(length(norm_bins)-1) # Standard error

        return c, m, el, eu
    end

## --- Single element constant-silica reference model

    tmin = 0
    tmax = 3900
    nbins = 39

    plotmin = 0
    plotmax = 4000

    elem = "MgO"

    (c, m, el, eu) = constprop(bin_bsr_means, ign, elem, tmin, tmax, nbins, p)

    # Plot results
    h = plot(c,m,yerror=(el,eu),seriestype=:scatter,color=:darkblue,mscolor=:darkblue,label="")
    plot!(h, c,m,style=:dot,color=:darkblue,mscolor=:darkblue,label="")
    plot!(h, xlabel="Age (Ma)", ylabel="$(elem)",xlims=(plotmin,plotmax),framestyle=:box,grid=:off,xflip=true) # Format plot
    savefig(h,"Constant Silica $(elem).pdf")
    display(h)

## --- Ratio constant-silica reference model

    tmin = 0
    tmax = 4000
    nbins = 8

    plotmin = 0
    plotmax = 4000

    num = "Rb"
    denom = "Sr"

    (c, m, el, eu) = constprop(bin_bsr_ratio_medians, ign, num, denom, tmin, tmax, nbins, p)

    # Plot results
    h = plot(c,m,yerror=(el,eu),seriestype=:scatter,color=:darkblue,mscolor=:darkblue,label="")
    plot!(h, c,m,style=:dot,color=:darkblue,mscolor=:darkblue,label="")
    plot!(h, xlabel="Age (Ma)", ylabel="$(num) / $(denom)",xlims=(plotmin,plotmax),framestyle=:box,grid=:off,xflip=true) # Format plot
    savefig(h,"Constant Silica $(num)_$(denom).pdf")
    display(h)

## --- End of File


text:
# Test Core data accessors
@testset "Core Data Accessers" begin
    # Setup data
    m = InfiniteModel()
    ind_idx = IndependentParameterIndex(1)
    fin_idx = FiniteParameterIndex(1)
    set = IntervalSet(0, 1)
    supps_dict = SortedDict{Float64, Set{DataType}}(0. => Set{DataType}([UserDefined]))
    method = InfiniteOpt.DefaultDerivativeMethod
    ind_param = IndependentParameter(set, supps_dict, 5, method)
    fin_param = FiniteParameter(42)
    ind_object = ScalarParameterData(ind_param, 1, 1, "ind")
    fin_object = ScalarParameterData(fin_param, -1, -1, "fin")
    ind_pref = IndependentParameterRef(m, ind_idx)
    fin_pref = FiniteParameterRef(m, fin_idx)
    ind_gvref = GeneralVariableRef(m, 1, IndependentParameterIndex)
    fin_gvref = GeneralVariableRef(m, 1, FiniteParameterIndex)
    bad_ind_pref = IndependentParameterRef(m, IndependentParameterIndex(-1))
    bad_fin_pref = FiniteParameterRef(m, FiniteParameterIndex(-1))
    # test dispatch_variable_ref
    @testset "dispatch_variable_ref" begin
        @test dispatch_variable_ref(m, ind_idx) == ind_pref
        @test dispatch_variable_ref(ind_gvref) == ind_pref
        @test dispatch_variable_ref(m, fin_idx) == fin_pref
        @test dispatch_variable_ref(fin_gvref) == fin_pref
    end
    # test _add_data_object
    @testset "_add_data_object" begin
        @test InfiniteOpt._add_data_object(m, ind_object) == ind_idx
        @test InfiniteOpt._add_data_object(m, fin_object) == fin_idx
        @test InfiniteOpt._param_object_indices(m)[end] == ind_idx
    end
    # test _data_dictionary
    @testset "_data_dictionary" begin
        @test InfiniteOpt._data_dictionary(ind_pref) === m.independent_params
        @test InfiniteOpt._data_dictionary(ind_gvref) === m.independent_params
        @test InfiniteOpt._data_dictionary(fin_pref) === m.finite_params
        @test InfiniteOpt._data_dictionary(fin_gvref) === m.finite_params
        @test InfiniteOpt._data_dictionary(m, IndependentParameter) === m.independent_params
        @test InfiniteOpt._data_dictionary(m, FiniteParameter) === m.finite_params
    end
    # test _data_object
    @testset "_data_object" begin
        @test InfiniteOpt._data_object(ind_pref) == ind_object
        @test InfiniteOpt._data_object(ind_gvref) == ind_object
        @test InfiniteOpt._data_object(fin_pref) == fin_object
        @test InfiniteOpt._data_object(fin_gvref) == fin_object
        @test_throws ErrorException InfiniteOpt._data_object(bad_ind_pref)
        @test_throws ErrorException InfiniteOpt._data_object(bad_fin_pref)
    end
    # test _core_variable_object
    @testset "_core_variable_object" begin
        @test InfiniteOpt._core_variable_object(ind_pref) == ind_param
        @test InfiniteOpt._core_variable_object(ind_gvref) == ind_param
        @test InfiniteOpt._core_variable_object(fin_pref) == fin_param
        @test InfiniteOpt._core_variable_object(fin_gvref) == fin_param
    end
    # test _parameter_number
    @testset "_parameter_number" begin
        @test InfiniteOpt._parameter_number(ind_pref) == 1
        @test InfiniteOpt._parameter_number(ind_gvref) == 1
    end
    # test _parameter_numbers
    @testset "_parameter_numbers" begin
        @test InfiniteOpt._parameter_numbers(ind_pref) == [1]
        @test InfiniteOpt._parameter_numbers(ind_gvref) == [1]
    end
    # test _object_number
    @testset "_object_number" begin
        @test InfiniteOpt._object_number(ind_pref) == 1
        @test InfiniteOpt._object_number(ind_gvref) == 1
    end
    # test _object_numbers
    @testset "_object_numbers" begin
        @test InfiniteOpt._object_numbers(ind_pref) == [1]
        @test InfiniteOpt._object_numbers(ind_gvref) == [1]
    end
    # test _adaptive_data_update
    @testset "_adaptive_data_update" begin
        # test change of same type
        p = IndependentParameter(IntervalSet(0, 2), supps_dict, 12, method)
        data = InfiniteOpt._data_object(ind_pref)
        @test InfiniteOpt._adaptive_data_update(ind_pref, p, data) isa Nothing
        @test InfiniteOpt._core_variable_object(ind_pref) == p
        # test change of different types
        p = IndependentParameter(IntervalSet(0, 2), supps_dict, 12, TestMethod())
        data = InfiniteOpt._data_object(ind_pref)
        @test InfiniteOpt._adaptive_data_update(ind_pref, p, data) isa Nothing
        @test InfiniteOpt._core_variable_object(ind_pref) == p
    end
    # test _set_core_variable_object
    @testset "_set_core_variable_object" begin
        @test InfiniteOpt._set_core_variable_object(ind_pref, ind_param) isa Nothing
        @test InfiniteOpt._set_core_variable_object(fin_pref, fin_param) isa Nothing
    end
    # test _delete_data_object
    @testset "_delete_data_object" begin
        @test is_valid(m, ind_pref)
        @test is_valid(m, fin_pref)
        @test InfiniteOpt._delete_data_object(ind_pref) isa Nothing
        @test InfiniteOpt._delete_data_object(fin_pref) isa Nothing
        @test !is_valid(m, ind_pref)
        @test !is_valid(m, fin_pref)
    end
end

# Test macro methods
@testset "Macro Helpers" begin
    @testset "Symbol Methods" begin
        @test InfiniteOpt._is_set_keyword(:(lower_bound = 0))
    end
    # test ParameterInfoExpr datatype
    @testset "_ParameterInfoExpr" begin
        @test InfiniteOpt._ParameterInfoExpr isa DataType
        @test InfiniteOpt._ParameterInfoExpr(ones(Bool, 8)...).has_lb
        @test !InfiniteOpt._ParameterInfoExpr().has_lb
    end
    # test InfiniteOpt._set_lower_bound_or_error
    @testset "InfiniteOpt._set_lower_bound_or_error" begin
        info = InfiniteOpt._ParameterInfoExpr()
        # test normal operation
        @test isa(InfiniteOpt._set_lower_bound_or_error(error, info, 0), Nothing)
        @test info.has_lb && info.lower_bound == 0
        # test double/lack of input errors
        @test_throws ErrorException InfiniteOpt._set_lower_bound_or_error(error,
                                                                   info, 0)
        info.has_lb = false; info.has_dist = true
        @test_throws ErrorException InfiniteOpt._set_lower_bound_or_error(error,
                                                                   info, 0)
        info.has_dist = false; info.has_set = true
        @test_throws ErrorException InfiniteOpt._set_lower_bound_or_error(error,
                                                                   info, 0)
    end
    # test InfiniteOpt._set_upper_bound_or_error
    @testset "InfiniteOpt._set_upper_bound_or_error" begin
        info = InfiniteOpt._ParameterInfoExpr()
        # test normal operation
        @test isa(InfiniteOpt._set_upper_bound_or_error(error, info, 0), Nothing)
        @test info.has_ub && info.upper_bound == 0
        # test double/lack of input errors
        @test_throws ErrorException InfiniteOpt._set_upper_bound_or_error(error,
                                                                   info, 0)
        info.has_ub = false; info.has_dist = true
        @test_throws ErrorException InfiniteOpt._set_upper_bound_or_error(error,
                                                                   info, 0)
        info.has_dist = false; info.has_set = true
        @test_throws ErrorException InfiniteOpt._set_upper_bound_or_error(error,
                                                                   info, 0)
    end
    # _dist_or_error
    @testset "_dist_or_error" begin
        info = InfiniteOpt._ParameterInfoExpr()
        # test normal operation
        @test isa(InfiniteOpt._dist_or_error(error, info, 0), Nothing)
        @test info.has_dist && info.distribution == 0
        # test double/lack of input errors
        @test_throws ErrorException InfiniteOpt._dist_or_error(error, info, 0)
        info.has_dist = false; info.has_lb = true
        @test_throws ErrorException InfiniteOpt._dist_or_error(error, info, 0)
        info.has_lb = false; info.has_set = true
        @test_throws ErrorException InfiniteOpt._dist_or_error(error, info, 0)
    end
    # _set_or_error
    @testset "_set_or_error" begin
        info = InfiniteOpt._ParameterInfoExpr()
        # test normal operation
        @test isa(InfiniteOpt._set_or_error(error, info, 0), Nothing)
        @test info.has_set && info.set == 0
        # test double/lack of input errors
        @test_throws ErrorException InfiniteOpt._set_or_error(error, info, 0)
        info.has_set = false; info.has_lb = true
        @test_throws ErrorException InfiniteOpt._set_or_error(error, info, 0)
        info.has_lb = false; info.has_dist = true
        @test_throws ErrorException InfiniteOpt._set_or_error(error, info, 0)
    end
    # _constructor_set
    @testset "InfiniteOpt._constructor_set" begin
        info = InfiniteOpt._ParameterInfoExpr()
        @test_throws ErrorException InfiniteOpt._constructor_set(error, info)
        info.has_lb = true; info.lower_bound = 0
        @test_throws ErrorException InfiniteOpt._constructor_set(error, info)
        info.has_ub = true; info.upper_bound = 1
        check = :(isa($(info.lower_bound), Real))
        expected = :($(check) ? IntervalSet($(info.lower_bound), $(info.upper_bound)) : error("Bounds must be a number."))
        @test InfiniteOpt._constructor_set(error, info) == expected
        info = InfiniteOpt._ParameterInfoExpr(distribution = Normal())
        check = :(isa($(info.distribution), Distributions.UnivariateDistribution))
        expected = :($(check) ? UniDistributionSet($(info.distribution)) : error("Distribution must be a Distributions.UnivariateDistribution."))
        @test InfiniteOpt._constructor_set(error, info) == expected
        info = InfiniteOpt._ParameterInfoExpr(set = IntervalSet(0, 1))
        check1 = :(isa($(info.set), InfiniteScalarSet))
        check2 = :(isa($(info.set), Distributions.UnivariateDistribution))
        expected = :($(check1) ? $(info.set) : ($(check2) ? UniDistributionSet($(info.set)) : error("Set must be a subtype of InfiniteScalarSet.")))
        @test InfiniteOpt._constructor_set(error, info) == expected
    end

    # _parse_one_operator_parameter
    @testset "_parse_one_operator_parameter" begin
        info = InfiniteOpt._ParameterInfoExpr()
        @test isa(InfiniteOpt._parse_one_operator_parameter(error, info,
                                                          Val(:<=), 0), Nothing)
        @test info.has_ub && info.upper_bound == 0
        @test isa(InfiniteOpt._parse_one_operator_parameter(error, info,
                                                          Val(:>=), 0), Nothing)
        @test info.has_lb && info.lower_bound == 0
        info = InfiniteOpt._ParameterInfoExpr()
        @test isa(InfiniteOpt._parse_one_operator_parameter(error, info,
                                                          Val(:in), esc(0)), Nothing)
        @test info.has_set && info.set == esc(0)
        info = InfiniteOpt._ParameterInfoExpr()
        @test isa(InfiniteOpt._parse_one_operator_parameter(error, info,
                                                          Val(:in), esc(:([0, 1]))), Nothing)
        @test info.has_lb && info.has_ub
        @test_throws ErrorException InfiniteOpt._parse_one_operator_parameter(error, info,
                                                                              Val(:d), 0)
    end
    # _parse_ternary_parameter
    @testset "_parse_ternary_parameter" begin
        info = InfiniteOpt._ParameterInfoExpr()
        @test isa(InfiniteOpt._parse_ternary_parameter(error, info, Val(:<=), 0,
                                                       Val(:<=), 0), Nothing)
        @test info.has_ub && info.upper_bound == 0
        @test info.has_lb && info.lower_bound == 0
        info = InfiniteOpt._ParameterInfoExpr()
        @test isa(InfiniteOpt._parse_ternary_parameter(error, info, Val(:>=), 0,
                                                       Val(:>=), 0), Nothing)
        @test info.has_ub && info.upper_bound == 0
        @test info.has_lb && info.lower_bound == 0
        @test_throws ErrorException InfiniteOpt._parse_ternary_parameter(error,
                                                 info, Val(:<=), 0, Val(:>=), 0)
    end
    # _parse_parameter
    @testset "_parse_parameter" begin
        info = InfiniteOpt._ParameterInfoExpr()
        @test InfiniteOpt._parse_parameter(error, info, :<=, :x, 0) == :x
        @test info.has_ub && info.upper_bound == 0
        @test InfiniteOpt._parse_parameter(error, info, :<=, 0, :x) == :x
        @test info.has_lb && info.lower_bound == 0
        info = InfiniteOpt._ParameterInfoExpr()
        @test InfiniteOpt._parse_parameter(error, info, 0, :<=, :x, :<=, 0) == :x
        @test info.has_ub && info.upper_bound == 0
        @test info.has_lb && info.lower_bound == 0
    end
end

# Test parameter definition methods
@testset "Definition" begin
    # _check_supports_in_bounds
    @testset "_check_supports_in_bounds" begin
        set = IntervalSet(0, 1)
        @test isa(InfiniteOpt._check_supports_in_bounds(error, 0, set), Nothing)
        @test_throws ErrorException InfiniteOpt._check_supports_in_bounds(error,
                                                                        -1, set)
        @test_throws ErrorException InfiniteOpt._check_supports_in_bounds(error,
                                                                         2, set)
        set = UniDistributionSet(Uniform())
        @test isa(InfiniteOpt._check_supports_in_bounds(error, 0, set), Nothing)
        @test_throws ErrorException InfiniteOpt._check_supports_in_bounds(error,
                                                                        -1, set)
        @test_throws ErrorException InfiniteOpt._check_supports_in_bounds(error,
                                                                         2, set)
    end
    # build_independent_parameter
    @testset "build_parameter (IndependentParameter)" begin
        set = IntervalSet(0, 1)
        supps = 0.
        supps_dict = SortedDict{Float64, Set{DataType}}(0. => Set([UserDefined]))
        method = TestMethod()
        @test build_parameter(error, set, supports = supps).set == set
        @test build_parameter(error, set, supports = supps).supports == supps_dict
        @test_throws ErrorException build_parameter(error, set, bob = 42)
        warn = "Ignoring num_supports since supports is not empty."
        @test_logs (:warn, warn) build_parameter(error, set,
                                            supports = [0, 1], num_supports = 2)
        repeated_supps = [1, 1]
        expected = IndependentParameter(set, SortedDict{Float64, Set{DataType}}(1. => Set{DataType}()), 5, method)
        warn = "Support points are not unique, eliminating redundant points."
        @test_logs (:warn, warn) build_parameter(error, set, supports = repeated_supps,
                                                 derivative_method = method) == expected
        set = UniDistributionSet(Normal())
        @test length(build_parameter(error, set, num_supports = 5).supports) == 5
        @test build_parameter(error, set, derivative_method = method).derivative_method == method
    end
    # build_finite_parameter
    @testset "build_parameter (FiniteParameter)" begin
        @test_throws ErrorException build_parameter(error, 1, bob = 42)
        expected = FiniteParameter(1)
        @test build_parameter(error, 1) == expected
    end

    # add_parameter
    @testset "add_parameter" begin
        m = InfiniteModel()
        method = InfiniteOpt.DefaultDerivativeMethod
        param = IndependentParameter(IntervalSet(0, 1),
                                    SortedDict{Float64, Set{DataType}}(), 5, method)
        expected = GeneralVariableRef(m, 1, IndependentParameterIndex, -1)
        @test add_parameter(m, param) == expected
        @test InfiniteOpt._core_variable_object(expected) == param
        @test InfiniteOpt._param_object_indices(m)[InfiniteOpt._object_number(expected)] == index(expected)
        param = FiniteParameter(1.5)
        expected = GeneralVariableRef(m, 1, FiniteParameterIndex, -1)
        @test add_parameter(m, param) == expected
        @test InfiniteOpt._core_variable_object(expected) == param
    end
end

# Test Reference Queries
@testset "Basic Reference Queries" begin
    m = InfiniteModel()
    p = build_parameter(error, IntervalSet(0, 1), derivative_method = TestMethod())
    pref = add_parameter(m, p)
    dpref = dispatch_variable_ref(pref)
    # JuMP.index
    @testset "JuMP.index" begin
        @test JuMP.index(pref) == IndependentParameterIndex(1)
        @test JuMP.index(dpref) == IndependentParameterIndex(1)
    end
    # JuMP.owner_model
    @testset "JuMP.owner_model" begin
        @test owner_model(pref) === m
        @test owner_model(dpref) === m
    end
    # has_derivative_supports
    @testset "has_derivative_supports" begin
        @test !has_derivative_supports(pref)
        InfiniteOpt._data_object(pref).has_derivative_supports = true
        @test has_derivative_supports(dpref)
    end
    # _set_has_derivative_supports
    @testset "_set_has_derivative_supports" begin
        @test InfiniteOpt._set_has_derivative_supports(pref, true) isa Nothing
        @test has_derivative_supports(pref)
        @test InfiniteOpt._set_has_derivative_supports(dpref, false) isa Nothing
        @test !has_derivative_supports(pref)
    end
    # derivative_method
    @testset "derivative_method" begin
        @test derivative_method(pref) isa TestMethod
        @test derivative_method(dpref) isa TestMethod
    end
    # has_internal_supports
    @testset "has_internal_supports" begin
        @test !has_internal_supports(pref)
        InfiniteOpt._data_object(pref).has_internal_supports = true
        @test has_internal_supports(dpref)
    end
    # _set_has_internal_supports
    @testset "_set_has_internal_supports" begin
        @test InfiniteOpt._set_has_internal_supports(pref, true) isa Nothing
        @test has_internal_supports(pref)
        @test InfiniteOpt._set_has_internal_supports(dpref, false) isa Nothing
        @test !has_internal_supports(pref)
    end
    # has_derivative_constraints
    @testset "has_derivative_constraints" begin
        @test !has_derivative_constraints(pref)
        InfiniteOpt._data_object(pref).has_deriv_constrs = true
        @test has_derivative_constraints(dpref)
    end
    # _set_has_derivative_constraints
    @testset "_set_has_derivative_constraints" begin
        @test InfiniteOpt._set_has_derivative_constraints(pref, true) isa Nothing
        @test has_derivative_constraints(pref)
        @test InfiniteOpt._set_has_derivative_constraints(dpref, false) isa Nothing
        @test !has_derivative_constraints(pref)
    end
end

# Test name methods
@testset "Name" begin
    m = InfiniteModel()
    param = build_parameter(error, IntervalSet(0, 1))
    pref = add_parameter(m, param, "test")
    dpref = dispatch_variable_ref(pref)
    bad = TestVariableRef(m, TestIndex(-1))
    bad_pref = FiniteParameterRef(m, FiniteParameterIndex(-1))
    # JuMP.name
    @testset "JuMP.name" begin
        @test_throws ArgumentError name(bad)
        @test name(pref) == "test"
        @test name(dpref) == "test"
        @test name(bad_pref) == ""
    end
    # JuMP.set_name
    @testset "JuMP.set_name" begin
        @test_throws ArgumentError set_name(bad, "test")
        @test isa(set_name(pref, "new"), Nothing)
        @test name(pref) == "new"
        @test isa(set_name(dpref, "test"), Nothing)
        @test name(dpref) == "test"
    end
    # _make_parameter_ref
    @testset "_make_parameter_ref" begin
        @test InfiniteOpt._make_parameter_ref(m, IndependentParameterIndex(1)) == pref
    end
    # _param_name_dict
    @testset "_param_name_dict" begin
        @test isa(InfiniteOpt._param_name_dict(m), Nothing)
    end
    # _update_param_name_dict
    @testset "_update_param_name_dict" begin
        m.name_to_param = Dict{String, AbstractInfOptIndex}()
        @test InfiniteOpt._update_param_name_dict(m, m.independent_params) isa Nothing
        @test InfiniteOpt._update_param_name_dict(m, m.dependent_params) isa Nothing
        @test m.name_to_param["test"] == IndependentParameterIndex(1)
        m.name_to_param = nothing
    end
    # parameter_by_name
    @testset "parameter_by_name" begin
        @test parameter_by_name(m, "test") == pref
        @test isa(parameter_by_name(m, "test2"), Nothing)
        pref = add_parameter(m, param, "test")
        m.name_to_param = nothing
        @test_throws ErrorException parameter_by_name(m, "test")
    end
end

# Test the parameter macro
@testset "Macro" begin
    m = InfiniteModel()
    # single parameter
    @testset "Single" begin
        pref = GeneralVariableRef(m, 1, IndependentParameterIndex, -1)
        @test @independent_parameter(m, 0 <= a <= 1) == pref
        @test InfiniteOpt._core_variable_object(pref).set == IntervalSet(0, 1)
        @test name(pref) == "a"
        pref = GeneralVariableRef(m, 2, IndependentParameterIndex, -1)
        @test @independent_parameter(m, b in Normal(), supports = [1; 2]) == pref
        @test InfiniteOpt._core_variable_object(pref).set == UniDistributionSet(Normal())
        @test InfiniteOpt._core_variable_object(pref).supports == SortedDict(i => Set{DataType}([UserDefined]) for i in [1,2])
        pref = GeneralVariableRef(m, 3, IndependentParameterIndex, -1)
        @test @independent_parameter(m, c in IntervalSet(0, 1)) == pref
        @test InfiniteOpt._core_variable_object(pref).set == IntervalSet(0, 1)
        pref = GeneralVariableRef(m, 4, IndependentParameterIndex, -1)
        @test @independent_parameter(m, set = IntervalSet(0, 1),
                                  base_name = "d") == pref
        @test name(pref) == "d"
        pref = GeneralVariableRef(m, 5, IndependentParameterIndex, -1)
        @test @independent_parameter(m, set = IntervalSet(0, 1)) == pref
        @test name(pref) == ""
        pref = GeneralVariableRef(m, 6, IndependentParameterIndex, -1)
        @test @independent_parameter(m, z in [0, 1], derivative_method = TestMethod()) == pref
        @test InfiniteOpt._core_variable_object(pref).set == IntervalSet(0, 1)
        @test name(pref) == "z"
        @test derivative_method(pref) isa TestMethod
        @test InfiniteOpt._param_object_indices(m)[InfiniteOpt._object_number(pref)] == index(pref)
    end
    # multiple parameters
    @testset "Array" begin
        prefs = [GeneralVariableRef(m, 7, IndependentParameterIndex, -1),
                 GeneralVariableRef(m, 8, IndependentParameterIndex, -1)]
        @test @independent_parameter(m, 0 <= e[1:2] <= 1) == prefs
        @test InfiniteOpt._core_variable_object(prefs[1]).set == IntervalSet(0, 1)
        @test InfiniteOpt._core_variable_object(prefs[2]).set == IntervalSet(0, 1)
        prefs = [GeneralVariableRef(m, 9, IndependentParameterIndex, -1),
                 GeneralVariableRef(m, 10, IndependentParameterIndex, -1)]
        @test @independent_parameter(m, [1:2], set = IntervalSet(0, 1)) == prefs
        @test InfiniteOpt._core_variable_object(prefs[1]).set == IntervalSet(0, 1)
        @test InfiniteOpt._core_variable_object(prefs[2]).set == IntervalSet(0, 1)
        prefs = [GeneralVariableRef(m, 11, IndependentParameterIndex, -1),
                 GeneralVariableRef(m, 12, IndependentParameterIndex, -1)]
        sets = [IntervalSet(0, 1), IntervalSet(-1, 2)]
        @test @independent_parameter(m, f[i = 1:2], set = sets[i]) == prefs
        @test InfiniteOpt._core_variable_object(prefs[1]).set == IntervalSet(0, 1)
        @test InfiniteOpt._core_variable_object(prefs[2]).set == IntervalSet(-1, 2)
        prefs = [GeneralVariableRef(m, 13, IndependentParameterIndex, -1),
                 GeneralVariableRef(m, 14, IndependentParameterIndex, -1)]
        @test @independent_parameter(m, [i = 1:2], set = sets[i]) == prefs
        @test InfiniteOpt._core_variable_object(prefs[1]).set == IntervalSet(0, 1)
        @test InfiniteOpt._core_variable_object(prefs[2]).set == IntervalSet(-1, 2)
        prefs = [GeneralVariableRef(m, 15, IndependentParameterIndex, -1),
                 GeneralVariableRef(m, 16, IndependentParameterIndex, -1)]
        @test @independent_parameter(m, [0, -1][i] <= g[i = 1:2] <= [1, 2][i]) == prefs
        @test InfiniteOpt._core_variable_object(prefs[1]).set == IntervalSet(0, 1)
        @test InfiniteOpt._core_variable_object(prefs[2]).set == IntervalSet(-1, 2)
        prefs = [GeneralVariableRef(m, 17, IndependentParameterIndex, -1),
                 GeneralVariableRef(m, 18, IndependentParameterIndex, -1)]
        prefs = convert(JuMP.Containers.SparseAxisArray, prefs)
        @test @independent_parameter(m, 0 <= i[1:2] <= 1,
                                  container = SparseAxisArray) == prefs
        @test InfiniteOpt._core_variable_object(prefs[1]).set == IntervalSet(0, 1)
        @test InfiniteOpt._core_variable_object(prefs[2]).set == IntervalSet(0, 1)
        @test InfiniteOpt._param_object_indices(m)[InfiniteOpt._object_number(prefs[2])] == index(prefs[2])
    end
    # test for errors
    @testset "Errors" begin
        @test_macro_throws ErrorException @independent_parameter(m, 0 <= [1:2] <= 1)
        @test_macro_throws ErrorException @independent_parameter(m, 0 <= "bob" <= 1)
        @test_macro_throws ErrorException @independent_parameter(m, 0 <= a <= 1)
        @test_macro_throws ErrorException @independent_parameter(m, 0 <= j)
        @test_macro_throws ErrorException @independent_parameter(m, j)
        @test_macro_throws ErrorException @independent_parameter(m, j, foo = 42)
        @test_macro_throws ErrorException @independent_parameter(m, j in Multinomial(3, [1/3, 1/3]))
        @test_macro_throws ErrorException @independent_parameter(m, 0 <= k <= 1, Int)
    end
end

# Test if used
@testset "Used" begin
    m = InfiniteModel()
    @independent_parameter(m, pref1 in [0, 1])
    @finite_parameter(m, pref2, 1)
    dpref1 = dispatch_variable_ref(pref1)
    dpref2 = dispatch_variable_ref(pref2)
    bad_pref = IndependentParameterRef(m, IndependentParameterIndex(-1))
    # _infinite_variable_dependencies
    @testset "_infinite_variable_dependencies" begin
        @test InfiniteOpt._infinite_variable_dependencies(pref1) == InfiniteVariableIndex[]
        @test InfiniteOpt._infinite_variable_dependencies(dpref1) == InfiniteVariableIndex[]
        @test InfiniteOpt._infinite_variable_dependencies(pref2) == InfiniteVariableIndex[]
        @test InfiniteOpt._infinite_variable_dependencies(dpref2) == InfiniteVariableIndex[]
        @test_throws ErrorException InfiniteOpt._infinite_variable_dependencies(bad_pref)
    end
    # _parameter_function_dependencies
    @testset "_parameter_function_dependencies" begin
        @test InfiniteOpt._parameter_function_dependencies(pref1) == ParameterFunctionIndex[]
        @test InfiniteOpt._parameter_function_dependencies(dpref1) == ParameterFunctionIndex[]
        @test InfiniteOpt._parameter_function_dependencies(pref2) == ParameterFunctionIndex[]
        @test InfiniteOpt._parameter_function_dependencies(dpref2) == ParameterFunctionIndex[]
        @test_throws ErrorException InfiniteOpt._parameter_function_dependencies(bad_pref)
    end
    # _derivative_dependencies
    @testset "_derivative_dependencies" begin
        @test InfiniteOpt._derivative_dependencies(pref1) == DerivativeIndex[]
        @test InfiniteOpt._derivative_dependencies(dpref1) == DerivativeIndex[]
        @test InfiniteOpt._derivative_dependencies(pref2) == DerivativeIndex[]
        @test InfiniteOpt._derivative_dependencies(dpref2) == DerivativeIndex[]
        @test_throws ErrorException InfiniteOpt._derivative_dependencies(bad_pref)
    end
    # _measure_dependencies
    @testset "_measure_dependencies" begin
        @test InfiniteOpt._measure_dependencies(pref1) == MeasureIndex[]
        @test InfiniteOpt._measure_dependencies(dpref1) == MeasureIndex[]
        @test InfiniteOpt._measure_dependencies(pref2) == MeasureIndex[]
        @test InfiniteOpt._measure_dependencies(dpref2) == MeasureIndex[]
        @test_throws ErrorException InfiniteOpt._measure_dependencies(bad_pref)
    end
    # _constraint_dependencies
    @testset "_constraint_dependencies" begin
        @test InfiniteOpt._constraint_dependencies(pref1) == ConstraintIndex[]
        @test InfiniteOpt._constraint_dependencies(dpref1) == ConstraintIndex[]
        @test InfiniteOpt._constraint_dependencies(pref2) == ConstraintIndex[]
        @test InfiniteOpt._constraint_dependencies(dpref2) == ConstraintIndex[]
        @test_throws ErrorException InfiniteOpt._constraint_dependencies(bad_pref)
    end
    # used_by_constraint
    @testset "used_by_constraint" begin
        @test !used_by_constraint(pref1)
        @test !used_by_constraint(pref2)
        @test !used_by_constraint(dpref1)
        @test !used_by_constraint(dpref2)
        push!(InfiniteOpt._constraint_dependencies(dpref1), ConstraintIndex(1))
        push!(InfiniteOpt._constraint_dependencies(dpref2), ConstraintIndex(1))
        @test used_by_constraint(pref1)
        @test used_by_constraint(pref2)
        @test used_by_constraint(dpref1)
        @test used_by_constraint(dpref2)
        popfirst!(InfiniteOpt._constraint_dependencies(dpref1))
        popfirst!(InfiniteOpt._constraint_dependencies(dpref2))
    end
    # used_by_measure
    @testset "used_by_measure" begin
        @test !used_by_measure(pref1)
        @test !used_by_measure(pref2)
        @test !used_by_measure(dpref1)
        @test !used_by_measure(dpref2)
        push!(InfiniteOpt._measure_dependencies(dpref1), MeasureIndex(1))
        push!(InfiniteOpt._measure_dependencies(dpref2), MeasureIndex(1))
        @test used_by_measure(pref1)
        @test used_by_measure(pref2)
        @test used_by_measure(dpref1)
        @test used_by_measure(dpref2)
        popfirst!(InfiniteOpt._measure_dependencies(dpref1))
        popfirst!(InfiniteOpt._measure_dependencies(dpref2))
    end
    # used_by_infinite_variable
    @testset "used_by_infinite_variable" begin
        @test !used_by_infinite_variable(pref1)
        @test !used_by_infinite_variable(pref2)
        @test !used_by_infinite_variable(dpref1)
        @test !used_by_infinite_variable(dpref2)
        push!(InfiniteOpt._infinite_variable_dependencies(dpref1), InfiniteVariableIndex(1))
        @test used_by_infinite_variable(pref1)
        @test used_by_infinite_variable(dpref1)
        popfirst!(InfiniteOpt._infinite_variable_dependencies(dpref1))
    end
    # used_by_parameter_function
    @testset "used_by_parameter_function" begin
        @test !used_by_parameter_function(pref1)
        @test !used_by_parameter_function(pref2)
        @test !used_by_parameter_function(dpref1)
        @test !used_by_parameter_function(dpref2)
        push!(InfiniteOpt._parameter_function_dependencies(dpref1), ParameterFunctionIndex(1))
        @test used_by_parameter_function(pref1)
        @test used_by_parameter_function(dpref1)
        popfirst!(InfiniteOpt._parameter_function_dependencies(dpref1))
    end
    # used_by_derivative
    @testset "used_by_derivative" begin
        @test !used_by_derivative(pref1)
        @test !used_by_derivative(pref2)
        @test !used_by_derivative(dpref1)
        @test !used_by_derivative(dpref2)
        push!(InfiniteOpt._derivative_dependencies(dpref1), DerivativeIndex(1))
        @test used_by_derivative(pref1)
        @test used_by_derivative(dpref1)
        popfirst!(InfiniteOpt._derivative_dependencies(dpref1))
    end
    # used_by_objective
    @testset "used_by_objective" begin
        @test !used_by_objective(pref1)
        @test !used_by_objective(dpref1)
        @test !used_by_objective(pref2)
        @test !used_by_objective(dpref2)
        InfiniteOpt._data_object(pref2).in_objective = true
        @test used_by_objective(pref2)
        @test used_by_objective(dpref2)
        InfiniteOpt._data_object(pref2).in_objective = false
    end
    # is_used
    @testset "is_used" begin
        @test !is_used(pref1)
        @test !is_used(pref2)
        @test !is_used(dpref1)
        @test !is_used(dpref2)
        push!(InfiniteOpt._constraint_dependencies(dpref1), ConstraintIndex(1))
        push!(InfiniteOpt._constraint_dependencies(dpref2), ConstraintIndex(1))
        @test is_used(pref1)
        @test is_used(pref2)
        @test is_used(dpref1)
        @test is_used(dpref2)
        popfirst!(InfiniteOpt._constraint_dependencies(dpref1))
        popfirst!(InfiniteOpt._constraint_dependencies(dpref2))
        push!(InfiniteOpt._measure_dependencies(dpref1), MeasureIndex(1))
        push!(InfiniteOpt._measure_dependencies(dpref2), MeasureIndex(1))
        @test is_used(pref1)
        @test is_used(pref2)
        @test is_used(dpref1)
        @test is_used(dpref2)
        popfirst!(InfiniteOpt._measure_dependencies(dpref1))
        popfirst!(InfiniteOpt._measure_dependencies(dpref2))
        push!(InfiniteOpt._infinite_variable_dependencies(dpref1), InfiniteVariableIndex(1))
        @test is_used(pref1)
        @test is_used(dpref1)
        popfirst!(InfiniteOpt._infinite_variable_dependencies(dpref1))
        push!(InfiniteOpt._parameter_function_dependencies(dpref1), ParameterFunctionIndex(1))
        @test is_used(pref1)
        @test is_used(dpref1)
        popfirst!(InfiniteOpt._parameter_function_dependencies(dpref1))
    end
end

# Test derivative methods
@testset "Derivative Methods" begin
    m = InfiniteModel()
    @independent_parameter(m, pref in [0, 1])
    dpref = dispatch_variable_ref(pref)
    func = (x) -> NaN
    num = 0.
    info = VariableInfo{Float64, Float64, Float64, Function}(true, num, true,
                                        num, true, num, false, func, false, false)
    d = Derivative(info, true, pref, pref) # this is wrong but that is ok
    object = VariableData(d)
    idx = InfiniteOpt._add_data_object(m, object)
    push!(InfiniteOpt._derivative_dependencies(pref), idx)
    dref = DerivativeRef(m, idx)
    gdref = GeneralVariableRef(m, idx.value, DerivativeIndex)
    cref = @constraint(m, gdref == 0)
    # test _reset_derivative_evaluations
    @testset "_reset_derivative_evaluations" begin
        # test empty
        @test InfiniteOpt._reset_derivative_evaluations(dpref) isa Nothing
        # test warning
        InfiniteOpt._set_has_derivative_constraints(pref, true)
        @test push!(InfiniteOpt._derivative_constraint_dependencies(dref), index(cref)) isa Vector
        warn = "Support/method changes will invalidate existing derivative evaluation " *
               "constraints that have been added to the InfiniteModel. Thus, " *
               "these are being deleted."
        @test_logs (:warn, warn) InfiniteOpt._reset_derivative_evaluations(dpref) isa Nothing
        @test !has_derivative_constraints(pref)
        @test !is_valid(m, cref)
        # test has derivative supports to deal with
        supps = SortedDict{Float64, Set{DataType}}(42 => Set([InternalLabel]))
        param = IndependentParameter(IntervalSet(0, 1), supps, 12, TestGenMethod())
        InfiniteOpt._set_core_variable_object(dpref, param)
        InfiniteOpt._set_has_derivative_supports(pref, true)
        @test InfiniteOpt._reset_derivative_evaluations(dpref) isa Nothing
        @test !has_derivative_supports(dpref)
        @test length(InfiniteOpt._core_variable_object(dpref).supports) == 0
    end
    # test set_derivative_method
    @testset "set_derivative_method" begin
        push!(InfiniteOpt._constraint_dependencies(dpref), ConstraintIndex(1))
        @test set_derivative_method(pref, TestMethod()) isa Nothing
        @test derivative_method(dpref) isa TestMethod
    end
end

# Test parameter set methods
@testset "Infinite Set" begin
    m = InfiniteModel()
    @independent_parameter(m, pref_gen in [0, 1])
    pref_disp = dispatch_variable_ref(pref_gen)
    bad = Bad()
    bad_pref = IndependentParameterRef(m, IndependentParameterIndex(-1))
    # _parameter_set
    @testset "_parameter_set" begin
        @test InfiniteOpt._parameter_set(pref_disp) == IntervalSet(0, 1)
        @test_throws ErrorException InfiniteOpt._parameter_set(bad_pref)
    end
    # _update_parameter_set
    @testset "_update_parameter_set " begin
        push!(InfiniteOpt._constraint_dependencies(pref_disp), ConstraintIndex(1))
        @test isa(InfiniteOpt._update_parameter_set(pref_disp,
                                                    IntervalSet(1, 2)), Nothing)
        @test InfiniteOpt._parameter_set(pref_disp) == IntervalSet(1, 2)
    end
    # infinite_set
    @testset "infinite_set" begin
        @test_throws ArgumentError infinite_set(bad)
        @test infinite_set(pref_disp) == IntervalSet(1, 2)
        @test infinite_set(pref_gen) == IntervalSet(1, 2)
        @test_throws ErrorException infinite_set(bad_pref)
    end
    # set_infinite_set
    @testset "set_infinite_set" begin
        @test_throws ArgumentError set_infinite_set(bad, IntervalSet(0, 1))
        @test isa(set_infinite_set(pref_disp, IntervalSet(2, 3)), Nothing)
        @test infinite_set(pref_disp) == IntervalSet(2, 3)
        @test isa(set_infinite_set(pref_gen, IntervalSet(1, 3)), Nothing)
        @test infinite_set(pref_gen) == IntervalSet(1, 3)
        @test set_infinite_set(pref_gen, UniDistributionSet(Normal())) isa Nothing
        @test infinite_set(pref_disp) isa UniDistributionSet
        push!(InfiniteOpt._data_object(pref_gen).measure_indices, MeasureIndex(1))
        @test_throws ErrorException set_infinite_set(pref_gen, IntervalSet(1, 3))
    end
end

# Test parameter support methods
@testset "Supports" begin
    m = InfiniteModel()
    @independent_parameter(m, pref in [0, 1], sig_digits = 5)
    @independent_parameter(m, pref2 in [0, 1])
    pref_disp = dispatch_variable_ref(pref)
    bad = Bad()
    push!(InfiniteOpt._data_object(pref).constraint_indices, ConstraintIndex(1))
    # _parameter_supports
    @testset "_parameter_supports" begin
        @test InfiniteOpt._parameter_supports(pref_disp) == SortedDict{Float64, Set{DataType}}()
    end
    @testset "_parameter_support_values" begin
        @test InfiniteOpt._parameter_support_values(pref_disp) == Float64[]
    end
    # _update_parameter_supports
    @testset "_update_parameter_supports " begin
        dict = SortedDict{Float64, Set{DataType}}(1. => Set{DataType}([MCSample]))
        @test isa(InfiniteOpt._update_parameter_supports(pref_disp, dict), Nothing)
        @test InfiniteOpt._parameter_support_values(pref_disp) == [1.]
    end
    # significant_digits
    @testset "significant_digits" begin
        @test significant_digits(pref) == 5
    end
    # num_supports
    @testset "num_supports" begin
        @test_throws ArgumentError num_supports(bad)
        @test num_supports(pref_disp) == 1
        @test num_supports(pref_disp, label = UserDefined) == 0
        @test num_supports(pref) == 1
        @test num_supports(pref, label = UserDefined) == 0
        @test num_supports(pref, label = SampleLabel) == 1
        @test num_supports(pref, label = All) == 1
    end
    # has_supports
    @testset "has_supports" begin
        @test_throws ArgumentError has_supports(bad)
        @test has_supports(pref_disp)
        @test has_supports(pref)
        InfiniteOpt._update_parameter_supports(pref_disp, SortedDict{Float64, Set{DataType}}())
        @test !has_supports(pref_disp)
        @test !has_supports(pref)
    end
    # supports
    @testset "supports" begin
        @test supports(pref_disp) == []
        @test supports(pref) == []
        dict = SortedDict{Float64, Set{DataType}}(1. => Set{DataType}([MCSample]))
        InfiniteOpt._update_parameter_supports(pref_disp, dict)
        @test supports(pref_disp) == [1.]
        @test supports(pref) == [1.]
        @test supports(pref, label = MCSample) == [1.]
        @test supports(pref, label = All) == [1.]
    end
    # supports (vector)
    @testset "supports (vector)" begin
        # test simple case
        @test supports([pref_disp], label = SampleLabel) == ones(1, 1)
        # test typical combinatorial case
        supps = [[-1, 0, 1], [-1, 1]]
        @infinite_parameter(m, x[i = 1:2] in [-1, 1], supports = supps[i], independent = true)
        @test supports(x) == Float64[-1 0 1 -1 0 1; -1 -1 -1 1 1 1]
        # test non-combinatorial case
        @test_throws ErrorException supports(x, use_combinatorics = false)
        @test set_supports(x[2], supps[1], force = true) isa Nothing
        @test supports(x, use_combinatorics = false) == Float64[-1 0 1; -1 0 1]
    end
    # set_supports
    @testset "set_supports" begin
        @test_throws ArgumentError set_supports(bad, [0, 1])
        @test isa(set_supports(pref_disp, [0, 1], force = true), Nothing)
        @test isa(set_supports(pref, [0, 1], force = true), Nothing)
        @test supports(pref) == [0., 1.]
        @test_throws ErrorException set_supports(pref, [2, 3])
        warn = "Support points are not unique, eliminating redundant points."
        @test_logs (:warn, warn) set_supports(pref, [1, 1], force = true)
        @test_throws ErrorException set_supports(pref, [0.5])
        @test !has_internal_supports(pref)
    end
    # add_supports
    @testset "add_supports" begin
        @test_throws ArgumentError add_supports(bad, 0.5)
        @test isa(add_supports(pref_disp, 0.25), Nothing)
        @test isa(add_supports(pref, 0.5), Nothing)
        @test supports(pref) == [0.25, 0.5, 1.]
        @test isa(add_supports(pref, [0, 0.25, 1], check = false), Nothing)
        @test supports(pref) == [0, 0.25, 0.5, 1.]
        @test add_supports(pref, 0.2, label = InternalLabel) isa Nothing
        @test supports(pref) == [0, 0.25, 0.5, 1.]
        @test supports(pref, label = All) == [0, 0.2, 0.25, 0.5, 1.]
        @test has_internal_supports(pref)
    end
    # delete_supports
    @testset "delete_supports" begin
        # test bad parameter
        @test_throws ArgumentError delete_supports(bad)
        # test label deletion
        set_derivative_method(pref, TestGenMethod())
        InfiniteOpt._set_has_derivative_supports(pref, true)
        add_supports(pref, 0.1, label = UniformGrid)
        @test delete_supports(pref_disp, label = UniformGrid) isa Nothing
        @test !has_derivative_supports(pref)
        @test !has_internal_supports(pref)
        @test supports(pref, label = All) == [0, 0.25, 0.5, 1.]
        # test total deletion
        @test isa(delete_supports(pref), Nothing)
        @test !has_derivative_supports(pref)
        @test !has_internal_supports(pref)
        @test supports(pref) == []
        # test array input
        @test isa(delete_supports([pref]), Nothing)
        # prepare to test derivative constraints
        func = (x) -> NaN
        num = 0.
        info = VariableInfo{Float64, Float64, Float64, Function}(true, num, true,
                                            num, true, num, false, func, false, false)
        deriv = Derivative(info, true, pref, pref)
        object = VariableData(deriv)
        idx = InfiniteOpt._add_data_object(m, object)
        push!(InfiniteOpt._derivative_dependencies(pref), idx)
        dref = DerivativeRef(m, idx)
        gdref = GeneralVariableRef(m, 1, DerivativeIndex)
        cref = @constraint(m, gdref == 0)
        InfiniteOpt._set_has_derivative_constraints(pref, true)
        @test push!(InfiniteOpt._derivative_constraint_dependencies(dref), index(cref)) isa Vector
        # test derivative constraint warning
        warn = "Deleting supports invalidated derivative evaluations. Thus, these " *
               "are being deleted as well."
        @test_logs (:warn, warn) delete_supports(pref) isa Nothing
        @test !has_derivative_constraints(pref)
        @test !is_valid(m, cref)
        # test measure error
        push!(InfiniteOpt._data_object(pref).measure_indices, MeasureIndex(1))
        @test_throws ErrorException delete_supports(pref)
    end
end

# Test lower bound functions
@testset "Lower Bound" begin
    m = InfiniteModel()
    @independent_parameter(m, pref1 in [0, 1])
    @independent_parameter(m, pref2 in Normal())
    @independent_parameter(m, pref3 in BadScalarSet())
    bad = TestVariableRef(m, TestIndex(-1))
    # JuMP.has_lower_bound
    @testset "JuMP.has_lower_bound" begin
        @test_throws ArgumentError has_lower_bound(bad)
        @test has_lower_bound(dispatch_variable_ref(pref1))
        @test has_lower_bound(pref1)
        @test has_lower_bound(dispatch_variable_ref(pref2))
        @test has_lower_bound(pref2)
    end
    # JuMP.lower_bound
    @testset "JuMP.lower_bound" begin
        @test_throws ArgumentError lower_bound(bad)
        @test lower_bound(dispatch_variable_ref(pref1)) == 0
        @test lower_bound(dispatch_variable_ref(pref2)) == -Inf
        @test lower_bound(pref1) == 0
        @test lower_bound(pref2) == -Inf
        @test_throws ErrorException lower_bound(pref3)
    end
    # JuMP.set_lower_bound
    @testset "JuMP.set_lower_bound" begin
        @test_throws ArgumentError set_lower_bound(bad, 2)
        @test_throws ErrorException set_lower_bound(dispatch_variable_ref(pref1), 2)
        @test_throws ErrorException set_lower_bound(dispatch_variable_ref(pref3), 2)
        @test_throws ErrorException set_lower_bound(pref1, 2)
        @test_throws ErrorException set_lower_bound(pref3, 2)
        @test isa(set_lower_bound(dispatch_variable_ref(pref1), -1), Nothing)
        @test lower_bound(pref1) == -1
        @test isa(set_lower_bound(pref1, -2), Nothing)
        @test lower_bound(pref1) == -2
    end
end

# Test upper bound functions
@testset "Upper Bound" begin
    m = InfiniteModel()
    @independent_parameter(m, pref1 in [0, 1])
    @independent_parameter(m, pref2 in Normal())
    @independent_parameter(m, pref3 in BadScalarSet())
    bad = TestVariableRef(m, TestIndex(-1))
    # JuMP.has_upper_bound
    @testset "JuMP.has_upper_bound" begin
        @test_throws ArgumentError has_upper_bound(bad)
        @test has_upper_bound(dispatch_variable_ref(pref1))
        @test has_upper_bound(dispatch_variable_ref(pref2))
        @test has_upper_bound(pref1)
        @test has_upper_bound(pref2)
    end
    # JuMP.lower_bound
    @testset "JuMP.upper_bound" begin
        @test_throws ArgumentError upper_bound(bad)
        @test upper_bound(dispatch_variable_ref(pref1)) == 1
        @test upper_bound(dispatch_variable_ref(pref2)) == Inf
        @test upper_bound(pref1) == 1
        @test upper_bound(pref2) == Inf
        @test_throws ErrorException upper_bound(pref3)
    end
    # JuMP.set_lower_bound
    @testset "JuMP.set_upper_bound" begin
        @test_throws ArgumentError set_upper_bound(bad, -1)
        @test_throws ErrorException set_upper_bound(dispatch_variable_ref(pref1), -1)
        @test_throws ErrorException set_upper_bound(dispatch_variable_ref(pref3), -1)
        @test_throws ErrorException set_upper_bound(pref1, -1)
        @test_throws ErrorException set_upper_bound(pref3, -1)
        @test isa(set_upper_bound(dispatch_variable_ref(pref1), 2), Nothing)
        @test upper_bound(pref1) == 2
        @test isa(set_upper_bound(pref1, 3), Nothing)
        @test upper_bound(pref1) == 3
    end
end

# Test methods for finite parameters
@testset "Finite Parameters" begin
    # initialize the model
    m = InfiniteModel()
    bad = TestVariableRef(m, TestIndex(-1))
    # test @finite_parameter
    @testset "@finite_parameter" begin
        m2 = Model()
        # test errors
        @test_macro_throws ErrorException @finite_parameter(m)
        @test_macro_throws ErrorException @finite_parameter(m, a, 2, 3)
        @test_macro_throws ErrorException @finite_parameter(m, (2, 3, 4), 2)
        @test_macro_throws ErrorException @finite_parameter(m2, 2)
        @test_macro_throws ErrorException @finite_parameter(m, "bob")
        @test_macro_throws ErrorException @finite_parameter(m2, test, 2)
        @test_macro_throws ErrorException @finite_parameter(m, test, 2, bob = 2)
        # test anonymous definition
        pref = GeneralVariableRef(m, 1, FiniteParameterIndex)
        @test @finite_parameter(m, 42) == pref
        @test InfiniteOpt._core_variable_object(pref).value == 42
        # test vector anonymous definition
        prefs = [GeneralVariableRef(m, 2, FiniteParameterIndex),
                 GeneralVariableRef(m, 3, FiniteParameterIndex)]
        @test @finite_parameter(m, [1:2], 42, base_name = "a") == prefs
        @test InfiniteOpt._core_variable_object(prefs[1]).value == 42
        @test name.(prefs) == ["a[1]", "a[2]"]
        # test named definition
        pref = GeneralVariableRef(m, 4, FiniteParameterIndex)
        @test @finite_parameter(m, b, 42) == pref
        @test InfiniteOpt._core_variable_object(pref).value == 42
        @test name(pref) == "b"
        # test named vector definition
        prefs = [GeneralVariableRef(m, 5, FiniteParameterIndex),
                 GeneralVariableRef(m, 6, FiniteParameterIndex)]
        prefs = convert(JuMPC.SparseAxisArray, prefs)
        @test @finite_parameter(m, c[i = 1:2], [3, 7][i],
                                container = SparseAxisArray) == prefs
        @test InfiniteOpt._core_variable_object(prefs[2]).value == 7
        @test name(prefs[2]) == "c[2]"
    end
    # initialize the model
    m = InfiniteModel()
    # test parameter_value
    @testset "parameter_value" begin
        pref = GeneralVariableRef(m, 1, FiniteParameterIndex)
        dpref = dispatch_variable_ref(pref)
        @test_throws ArgumentError parameter_value(bad)
        @test @finite_parameter(m, g, 1) == pref
        @test parameter_value(dpref) == 1
        @test parameter_value(pref) == 1
    end
    # test JuMP.set_value
    @testset "JuMP.set_value" begin
        pref = GeneralVariableRef(m, 1, FiniteParameterIndex)
        dpref = dispatch_variable_ref(pref)
        push!(InfiniteOpt._constraint_dependencies(dpref), ConstraintIndex(1))
        @test_throws ArgumentError set_value(bad, 42)
        @test isa(set_value(dpref, 42), Nothing)
        @test parameter_value(pref) == 42
        @test isa(set_value(pref, 41), Nothing)
        @test parameter_value(pref) == 41
    end
end

# Test support flll-in and geneartion functions
@testset "Support Fill-in and Generation" begin
    @testset "generate_and_add_supports! (AbstractInfiniteSet)" begin
        m = InfiniteModel()
        gvref1 = @independent_parameter(m, 0 <= a <= 1)
        pref1 = dispatch_variable_ref(gvref1)
        set1 = infinite_set(pref1)
        dist = Normal(0., 1.)
        gvref2 = @independent_parameter(m, c in dist)
        pref2 = dispatch_variable_ref(gvref2)
        set2 = infinite_set(pref2)
        @test generate_and_add_supports!(pref1, set1, num_supports = 4) isa Nothing
        @test generate_and_add_supports!(pref2, set2, num_supports = 4) isa Nothing
        @test length(supports(pref1)) == 4
        @test length(supports(pref2)) == 4
    end
    # fill_in_supports! (ParameterRef)
    @testset "fill_in_supports! (ParameterRef)" begin
        m = InfiniteModel()
        pref1 = @independent_parameter(m, 0 <= a <= 1)
        pref2 = @independent_parameter(m, 0 <= b[1:2] <= 1)
        dist = Normal(0., 1.)
        pref3 = @independent_parameter(m, c in dist, supports = [-0.5, 0.5])
        @test fill_in_supports!(pref1, num_supports = 11) isa Nothing
        @test fill_in_supports!.(pref2, num_supports = 11) isa Array{Nothing}
        @test fill_in_supports!(pref3, num_supports = 11) isa Nothing
        @test length(supports(pref1)) == 11
        @test length(supports(pref2[1])) == 11
        @test length(supports(pref2[2])) == 11
        @test length(supports(pref3)) == 11
        @test -0.5 in supports(pref3)
        @test 0.5 in supports(pref3)
        @test fill_in_supports!(pref1, num_supports = 20) isa Nothing
        @test length(supports(pref1)) == 20
    end
end


text:
#This script tests the time it takes to compute a projection vs model size in
#3D. Tests serial PARSDMM, parallel PARSDMM, multilevel serial PARSDMM, and
#multilevel parallel PARSDMM. Requires 5 Julia workers (julia -p 5). timings
#in the paper use export JULIA_NUM_THREADS=4

using Distributed
@everywhere using SetIntersectionProjection
@everywhere using LinearAlgebra

using HDF5

@everywhere mutable struct compgrid
  d :: Tuple
  n :: Tuple
end

width=[50 100 200 340]

test = Vector{Any}(undef,length(width))
N    = Vector{Any}(undef,length(width))

#PARSDMM options:
options              = PARSDMM_options()
options.FL           = Float32
options.evol_rel_tol = 10*eps(options.FL)
set_zero_subnormals(true)

#select working precision
# if options.FL==Float64
#   TF = Float64
#   TI = Int64
# elseif options.FL==Float32
#   TF = Float32
#   TI = Int32
# end
TF = Float32

#define constraints
constraint = Vector{SetIntersectionProjection.set_definitions}() #initialize

#bound constraints
m_min     = 1500.0
m_max     = 6000.0
set_type  = "bounds"
TD_OP     = "identity"
app_mode  = ("tensor","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#monotonicity
m_min     = 0.0
m_max     = 1e6
set_type  = "bounds"
TD_OP     = "D_z"
app_mode  = ("tensor","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#lateral smoothness via a bound on the gradient
m_min     = -1.0
m_max     = 1.0
set_type  = "bounds"
TD_OP     = "D_x"
app_mode  = ("tensor","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

m_min     = -1.0
m_max     = 1.0
set_type  = "bounds"
TD_OP     = "D_y"
app_mode  = ("tensor","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

# Load velocity model
if ~isfile("overthrust_3D_true_model.h5")
  run(`wget ftp://slim.gatech.edu/data/SoftwareRelease/WaveformInversion.jl/3DFWI/overthrust_3D_true_model.h5`)
end
n, d, o, m_full = read(h5open("overthrust_3D_true_model.h5","r"), "n", "d", "o", "m")

m_full .= 1000.0 ./ sqrt.(m_full);

#initialize arrays to save timings
log_T_serial   = Vector{Any}(undef,length(width))
T_tot_serial   = Vector{Any}(undef,length(width))
log_T_parallel = Vector{Any}(undef,length(width))
T_tot_parallel = Vector{Any}(undef,length(width))

log_T_serial_multilevel   = Vector{Any}(undef,length(width))
T_tot_serial_multilevel   = Vector{Any}(undef,length(width))
log_T_parallel_multilevel = Vector{Any}(undef,length(width))
T_tot_parallel_multilevel = Vector{Any}(undef,length(width))

options.rho_ini = [1.0;1000.0;1000.0;1000.0;1.0]

for i=length(width):-1:1
  print(i)

  m         = m_full[1:width[i],1:width[i],:];
  comp_grid = compgrid( (TF(d[1]), TF(d[2]), TF(d[3])),( size(m,1), size(m,2), size(m,3) ) )
  m         = convert(Vector{TF},vec(m))

  N[i]=prod(size(m));

  #parallel
  @everywhere GC.gc()
  println("")
  println("parallel")
  options.parallel=true
  @everywhere BLAS.set_num_threads(2)
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  (TD_OP,AtA,l,y)        = PARSDMM_precompute_distribute(TD_OP,set_Prop,comp_grid,options)
  (x,log_PARSDMM)        = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);
  val, t, bytes, gctime, memallocs = @timed (x,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);
  println(t)
  log_T_parallel[i] = log_PARSDMM
  T_tot_parallel[i] = t

  #serial
  @everywhere GC.gc()
  println("")
  println("serial")
  options.parallel=false
  #options.rho_ini = [10.0]
  BLAS.set_num_threads(4)
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  (TD_OP,AtA,l,y)        = PARSDMM_precompute_distribute(TD_OP,set_Prop,comp_grid,options)
  (x,log_PARSDMM)        = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);
  val, t, bytes, gctime, memallocs = @timed (x,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);
  println(t)
  log_T_serial[i] = log_PARSDMM
  T_tot_serial[i] = t

  #serial multilevel
  @everywhere GC.gc()
  println("")
  println("serial multilevel")
  BLAS.set_num_threads(4)
  #options.rho_ini = [1000.0]
  options.parallel  = false
  n_levels          = 3
  coarsening_factor = 2
  (TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,constraint_level)=setup_multi_level_PARSDMM(m,n_levels,coarsening_factor,comp_grid,constraint,options)
  (x,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);
  val, t, bytes, gctime, memallocs = @timed (x,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);
  println(t)
  log_T_serial_multilevel[i] = log_PARSDMM
  T_tot_serial_multilevel[i] = t

  #parallel multilevel
  @everywhere GC.gc()
  println("")
  println("parallel multilevel")
  @everywhere BLAS.set_num_threads(2)
  options.parallel  = true
  n_levels          = 3
  coarsening_factor = 2
  (TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,constraint_level)=setup_multi_level_PARSDMM(m,n_levels,coarsening_factor,comp_grid,constraint,options)
  (x,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);
  val, t, bytes, gctime, memallocs = @timed (x,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);
  println(t)
  log_T_parallel_multilevel[i] = log_PARSDMM
  T_tot_parallel_multilevel[i] = t

end

ENV["MPLBACKEND"]="qt5agg"
using PyPlot
#plot results
fig, ax = subplots()
ax[:loglog](N, T_tot_serial, marker="o", markersize=10, label="serial",linewidth=5)
ax[:loglog](N, T_tot_parallel, marker="x", markersize=10, label="parallel",linewidth=5)
ax[:loglog](N, T_tot_serial_multilevel, marker="v", markersize=10, label="serial multilevel",linewidth=5)
ax[:loglog](N, T_tot_parallel_multilevel, marker="^", markersize=10, label="parallel multilevel",linewidth=5)
ax[:legend]()
title(string("time 3D vs grid size, JuliaThreads=",Threads.nthreads(),", BLAS threads=",ccall((:openblas_get_num_threads64_, Base.libblas_name), Cint, ())), fontsize=12)
xlabel("N gridpoints", fontsize=15)
ylabel("time [seconds]", fontsize=15)
savefig("projection_intersection_timings3D_1.eps",bbox_inches="tight",dpi=1200)

#######################################################################################
#######################################################################################

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
module PredictorCorrector
    using ..HomConBase
    using ..Homotopy
    import ..HomConBase: HomConAlgorithm, AbstractHomotopy, solve

    export AbstractPredictorCorrectorAlgorithm
    abstract type AbstractPredictorCorrectorAlgorithm{Projective} <: HomConAlgorithm end

    const APCA{P} = AbstractPredictorCorrectorAlgorithm{P}

    """
        predict(alg, H, J_H, ∂H∂t, x, t, Δt)

    Make a prediction step for the algorithm `alg` and the homotopy `H` with the jacobian
    `J_H`, the derivative w.r.t t `∂H∂t` at `x` to the time `t` with a stepwidth `Δt`.
    """
    function predict end

    """
        correct!(u, alg, H, J_H, x, t, tol, max_iterations)::Bool

    Make a correction step for the algorithm `alg` and the homotopy `H` with the jacobian
    `J_H` at `x` to the time `t`. Stores the result in `u` and returns a boolean indicating
    whether ``|u-x|<tol`` withing `max_iterations` iterations.
    """
    function correct! end

    export correct!, predict

    include("spherical.jl")
    include("affine.jl")

    include("prepare.jl")
    include("trackpath.jl")
    include("result_types.jl")
    include("cauchyendgame.jl")
    include("solve.jl")

    export Spherical, Affine, solve, trackpath, cauchyendgame
    export Result, CauchyEndgameResult, PathResult, ConvergentCluster
end


text:
@recipe function f(ga::AbstractDimensionalArray{T,2,<:Tuple{<:AbstractDimension,<:Time}}) where T
    ylabel --> label(ga)
    xlabel --> label(dims(ga)[1])
    legendtitle --> label(dims(ga)[1])
    title --> label(refdims(ga))
    parent(ga)
end

@recipe function f(ga::AbstractDimensionalArray{T,2,<:Tuple{<:Time,<:AbstractDimension}}) where T
    permutedims(ga)
end

@recipe function f(ga::AbstractDimensionalArray{T,1,<:Tuple{<:AbstractDimension}}) where T
    ylabel --> label(ga)
    xlabel --> label(dims(ga)[1])
    # legend --> false
    title --> label(refdims(ga))
    val(dims(ga)[1]), parent(ga)
end


text:
function readYahooAdjClose(dates::StepRange,
                          ticker::Array{String, 1},
                          freq=:d)

    ## parallelized version for multiple stocks

    nStocks = length(ticker)
    allStocks = @distributed(joinSortedIdx_outer) for ii=1:nStocks
        getAdjClose(dates, ticker[ii], freq)
    end

    return allStocks
end

function getAdjClose(dates::StepRange,
                            ticker::String,
                            freq=:d)
    stock = readYahooFinance(dates, ticker, freq)

    ## refine stock
    if isequal(stock, NA)
        adjClose = Timenum(DataFrame(Adj_Close = NA), [dates[end]])
    else
        adjClose = stock[:Adj_Close]
    end
    names!(adjClose.vals, [symbol(ticker)])

    return adjClose
end

function readYahooFinance(dates::StepRange,
                          ticker::String,
                          freq=:d)

    #############################
    ## get download parameters ##
    #############################

    ## get download urls
    urls = getUrls(dates, [ticker], freq)

    ###################
    ## download data ##
    ###################

    try
        filepath = download(urls[1])
        td = readTimedata(filepath)
        vals = flipud(td)
    catch
        vals = NA
    end
end

function getUrls(dates::StepRange,
                 ticker::Array{String, 1},
                 freq)
    nStocks = length(ticker)

    urls = Array(String, nStocks)
    (bd, bm, by, ed, em, ey) = getDates(dates)

    for ii=1:nStocks
        urls[ii] = string("http://real-chart.finance.yahoo.com/table.csv?s=",
                          ticker[ii],
                          "&a=",bm,"&b=",bd,"&c=",by,
                          "&d=",em,"&e=",ed,"&f=",ey,
                          "&g=",freq,"&ignore=.csv")
    end
    urls
end

function getDates(dates::StepRange)
    ## get start and end date
    startDate = dates[1]
    endDate = dates[end]

    ## decompose dates
    (bd, bm, by) = (day(startDate), month(startDate), year(startDate))
    (ed, em, ey) = (day(endDate), month(endDate), year(endDate))

    (lpad(string(bd-1), 2, "0"),
     lpad(string(bm-1), 2, "0"),
     string(by),
     lpad(string(ed-1), 2, "0"),
     lpad(string(em-1), 2, "0"),
     string(ey))
end

################
## test cases ##
################

## one function that returns complete data for single ticker
## this function can be applied to each element of ticker array
##
## another function that returns adj_closing prices for array of
## ticker symbols in one large


text:
using Distributions

function test_gaussian_recon_loss_layer(backend::Backend, T, eps)
    println("-- Testing GaussianReconLossLayer on $(typeof(backend)){$T}...")

    ############################################################
    # Prepare Data for Testing
    ############################################################
    dims = (3, 4, 5)
    tensor_dim = length(dims)
    p, n = prod(dims[1:end-1]), dims[end]
    println("        > $dims")

    ############################################################
    # Setup
    ############################################################
    layer = GaussianReconLossLayer(bottoms=[:mu, :sigma, :x])

    mu = rand(p, n)  # sigmoid output is in [0, 1]
    sigma = max(1e-2, 0.2rand(p, n)) # PReLU output is > 0
    x = max(1e-2, rand(dims))  # normalized pixel values are in [0, 1]

    inputs = Blob[
        make_blob(backend, convert(Array{T}, mu)),
        make_blob(backend, convert(Array{T}, sigma)),
        make_blob(backend, convert(Array{T}, x))]

    diffs = Blob[
        make_blob(backend, rand(T, size(inputs[1]))),  # gets overwritten
        make_blob(backend, rand(T, size(inputs[2]))),  # gets overwritten
        NullBlob()]  # the data is fixed

    state = setup(backend, layer, inputs, diffs)

    ############################################################
    # Forward Propagation
    ############################################################
    forward(backend, state, inputs)

    loss = 0.
    for i in 1:n
        Xi_given_zi = MvNormal(mu[:, i], sigma[:, i])  # a random var
        xi = slicedim(x, tensor_dim, i)[:]  # an instance
        loss -= logpdf(Xi_given_zi, xi)
    end
    loss /= n
    @test loss >= 0

    info("expected $loss; got $(state.loss)")
    @test -eps < loss - state.loss < eps

    ############################################################
    # Backward Propagation
    ############################################################
    backward(backend, state, inputs, diffs)

    for i in 1:2
        input_i = Array(T, size(inputs[i]))
        copy!(input_i, inputs[i])
        diff_i = Array(T, size(diffs[i]))
        copy!(diff_i, diffs[i])

        for k in 1:p, j in 1:n
            x0 = input_i[k, j]
            function f(x::FloatingPoint)
                input_i[k, j] = x
                copy!(inputs[i], input_i)
                forward(backend, state, inputs)
                input_i[k, j] = x0
                copy!(inputs[i], input_i)
                state.loss
            end
            test_deriv(f, x0, diff_i[k, j])
        end
    end

    ############################################################
    # Shutdown
    ############################################################
    shutdown(backend, state)
end

function test_gaussian_recon_loss_layer(backend::Backend)
    test_gaussian_recon_loss_layer(backend, Float64, 1e-8)
end

if test_cpu
    test_gaussian_recon_loss_layer(backend_cpu)
end
if test_gpu
    test_gaussian_recon_loss_layer(backend_gpu)
end


text:
struct d1_coloring <: AbstractColoring
    colpack_coloring::String
end

struct d2_coloring <: AbstractColoring
    colpack_coloring::String
end

struct acyclic_coloring <: AbstractColoring
    colpack_coloring::String
end

struct star_coloring <: AbstractColoring
    colpack_coloring::String
end

d1_coloring() = d1_coloring("DISTANCE_ONE")
d2_coloring()= d2_coloring("DISTANCE_TWO")
acyclic_coloring() = acyclic_coloring("ACYCLIC")
star_coloring() = star_coloring("STAR")

COLORINGS = [
    d1_coloring(),
    d2_coloring(),
    acyclic_coloring(),
    star_coloring()
]

export COLORINGS
export d1_coloring, d2_coloring, acyclic_coloring, star_coloring


text:
using KhepriUnreal
using Test

@testset "KhepriUnreal.jl" begin
    # Write your tests here.
end


text:
function check_simplex(sx)
    if !(eltype(sx) <: Integer)
        throw(ArgumentError("Invalid simplex '$sx'!" *
                            "Labels should <: Integer!"))
    end
    if isa(sx, AbstractSet)
        sort!(collect(sx))
    else
        if !allunique(sx)
            throw(ArgumentError("Invalid simplex '$sx'!" *
                                "Vertices in simplex should be unique."))
        end
        if !issorted(sx)
            throw(ArgumentError("Invalid simplex '$sx'!" *
                                "The labels in simplex should be sorted."))
        end
        sx
    end
end

function check_simplex(sxtree, sx)
    check_simplex(sx)
    if any(sx .> num_children(sxtree))
        lbl = filter(l -> l > num_children(sxtree), sx)
        throw(ArgumentError("Invalid simplex '$sx'!" *
                            "Label $lbl too large for this tree."))
    end
end


text:
@safetestset "Create new config file" begin

#===============================================================================

  using SearchLight
  using YAML

  workdir = pwd()
  testdir = mktempdir()
  cd(testdir)

  @testset "Create config without explicit arguments" begin
    SearchLight.Generator.newconfig()
    @test isdir(joinpath(testdir, SearchLight.DB_PATH)) == true
    @test isfile(joinpath(testdir, SearchLight.DB_PATH, SearchLight.SEARCHLIGHT_DB_CONFIG_FILE_NAME)) == true
  end;

  @testset "Create config with custom file name" begin
    filename = "db.yml"
    SearchLight.Generator.newconfig(filename = filename)
    @test isfile(joinpath(testdir, SearchLight.DB_PATH, filename)) == true
  end;

  @testset "Create config with custom path and filename" begin
    filename = "db.yml"
    filepath = joinpath("foo", "bar", "baz")
    SearchLight.Generator.newconfig(filepath, filename = filename)
    @test isfile(joinpath(testdir, filepath, filename)) == true
  end;

  @testset "Create config in current dir" begin
    filename = "db.yml"
    filepath = "."
    SearchLight.Generator.newconfig(filepath, filename = filename)
    @test isfile(joinpath(testdir, filepath, filename)) == true
  end;

  cd(workdir)
end;

@safetestset "Adapter config properties" begin
  workdir = pwd()
  testdir = mktempdir()
  cd(testdir)

  @safetestset "Default adapter is MySQL" begin
    using SearchLight
    using YAML

    SearchLight.Generator.newconfig(".")
    db_conn_data = YAML.load(open(SearchLight.SEARCHLIGHT_DB_CONFIG_FILE_NAME))

    @test isa(db_conn_data, Dict) == true
    @test haskey(db_conn_data, SearchLight.config.app_env) == true
    @test haskey(db_conn_data[SearchLight.config.app_env], "adapter") == true
    @test db_conn_data[SearchLight.config.app_env]["adapter"] == "MySQL"
    @test haskey(db_conn_data[SearchLight.config.app_env], "port") == true
    @test db_conn_data[SearchLight.config.app_env]["port"] == 3306
  end;

  cd(workdir)

=#

end


text:
## Step2 is iteration to find the stationary distribution

#module Step2

using LinearAlgebra, Random, Distributions, Statistics, Plots
using BenchmarkTools, Interpolations
include("Step1Grid.jl")

"""
    Relies on step1 (either with simple grid search or endogenous) and finds
    stationary distribution.
"""

# Define function for inverse of policy function a^{-1}
function pol_inv(A, pol_func, aprime, eprime)

    """
    Function that finds the inverse of the optimal policy function.
    Inputs:
        - A: set of assets
        - pol_func: optimal policy function
        - aprime: assets in the next period
        - eprime: exogenous state variable in the next period
    Output: inverse of the policy function
    """

    # if a' < minimum value of policy function, then a must be a_lowerbar
    if minimum(pol_func[:, eprime]) > A[aprime]

        return 1

    # when e'=2, distinguish pre vs post cross 45° line (endogenous a_UpperBar)
    elseif eprime == 2

        return min(findlast(pol_func[:, eprime] .<= A[aprime]),
            findfirst((pol_func[:, eprime] - A) .< 0.01))

    # any other case, find element of A s.t. policy function is closest to a'
    else

        return findlast(pol_func[:, eprime] .<= A[aprime])

    end

end

# Define function for Step2
function step2(A, E, Π, pol_func; tol = 0.01, maxT = 600)

    """
    Function that finds the stationary distribution Ψ.
    Relies on pol_inv in the code.
    Inputs:
        - A: set of assets
        - E: set of exogenous shocks
        - Π: transition matrix
        - pol_func: optimal policy function
    Output:
        - F_update: stationary distribution Ψ
    """

    # retrieve dimensions
    a_size = length(A)
    e_size = size(E)[1]

    # Check that everything works
    @assert size(pol_func)[1] == a_size
    @assert size(pol_func)[2] == e_size

    # Find endogenous position where a' = a
    aBar = findlast((pol_func[:,2] - A) .< tol)
    cut = max(aBar, 2)           # if aBar = 1 then range will not work

    aux1 = range(zero(eltype(A)), one(eltype(A))/e_size, length = cut)
    aux2 = ones(eltype(A), a_size - cut)/e_size
    aux = [aux1; aux2]
    F0 = repeat(aux, 1, e_size)

    F_update = copy(F0)

    # Iteration of distribution
    normdiff = Inf
    t = 1

    while normdiff > tol && t <= maxT

        for a = 1:a_size
            for e = 1:e_size

                F_update[a,e] = Π[1, e] * F0[pol_inv(A, pol_func, a, 1), 1] +
                          Π[2, e] * F0[pol_inv(A, pol_func, a, 2), 2]

            end
        end

        t += 1
        normdiff = maximum(abs.(F_update - F0))

        F0 = copy(F_update)

    end

    return F_update

end

#export pol_inv, step2

#end


text:
# HalidePerovskites
# Uses the newly forked off FeynmanKadanoffOsakaHellwarth.jl package
# Codes by Jarvist Moore Frost, 2017
# These codes were developed with Julia 0.5.0, and requires the Optim and Plots packages.

# This file, when run under Julia, should regenerate polaron data associated with Arxiv paper:
# https://arxiv.org/abs/1708.04158
# Slow cooling of hot polarons in halide perovskite solar cells
# Jarvist Moore Frost, Lucy D. Whalley and Aron Walsh
# (Submitted on 14 Aug 2017 [v1])

push!(LOAD_PATH,"../src/") # load module from local directory

using PolaronMobility
using PlotPolaron # Plots dependency

##### load in library routines... #####
# Plot figures with Plots, which defaults to Pyplot backend
using Plots
pyplot() # PyPlot (matplotlib) backend, to Plots
#gr() # GR backend, to Plots
default(grid=false) # No silly dotted grid lines
default(size=(400,300)) # A good small size for two-column EPS output
#default(size=(800,600)) # Nice size for small-ish PNGs for slides

# Physical constants
const hbar = const ħ = 1.05457162825e-34;          # kg m2 / s
const eV = const q = const ElectronVolt = 1.602176487e-19;                         # kg m2 / s2
const me=MassElectron = 9.10938188e-31;                          # kg
const Boltzmann = const kB =  1.3806504e-23;                  # kg m2 / K s2
const ε_0 = 8.854E-12 #Units: C2N−1m−2, permittivity of free space
# Units
Å=1E-10

# --> 0.12 for electrons and 0.15 for holes, in MAPI. See 2014 PRB.
# MAPI  4.5, 24.1, 2.25THz - 75 cm^-1 ; α=
MAPIe=polaronmobility(10:20:1000, 4.5, 24.1, 2.25E12, 0.12)
savepolaron("MAPI-electron",MAPIe)
plotpolaron("MAPI-electron", MAPIe)
MAPIh=polaronmobility(10:20:1000, 4.5, 24.1, 2.25E12, 0.15)
plotpolaron("MAPI-hole", MAPIh)
savepolaron("MAPI-hole",MAPIh)

## Polaron radius vs. Temperature
# More complex figure for the Thermal Pathways paper;
# https://github.com/WMD-group/2017-03-ThermalPathways/
p=MAPIe

plot(p.T,p.rfsi./Å, markersize=2,marker=:rect,
    label="Polaron radius",xlab="Temperature (K)",ylab="Polaron Radius (Angstrom)",ylims=(0,Inf))
plot!(p.T,p.rfsmallalpha./Å,label="T=0 Schultz small alpha polaron radius")
savefig("ThermalPathways-radius.png")
savefig("ThermalPathways-radius.pdf")

plot(p.T,p.v./p.w,label="v/w",markersize=2,marker=:circle,
    xlab="Temperature (K)", ylab="\hbar\omega")
savefig("ThermalPathways-vwratio.png")
savefig("ThermalPathways-vwratio.pdf")

    plot!(p.T,p.v,label="v",markersize=2,marker=:uptriangle)
plot!(p.T,p.w,label="w",markersize=2,marker=:downtriangle)
savefig("ThermalPathways-vw.png")
savefig("ThermalPathways-vw.pdf")

println("That's me!")


text:
function generate_point(rng = Random.GLOBAL_RNG)
    thres = rand(rng)
    if thres >= 0.5
        rand(rng, LogNormal())
    else
        rand(rng, Uniform(2.0, 3.0))
    end
end

mt = Random.MersenneTwister(42)
xs = [generate_point(mt) for _ in 1:5000]
bandwidths = [0.05, 0.1, 0.5]
densities = [KernelDensity.kde(xs, bandwidth = bw) for bw in bandwidths]

p = Plots.plot()
for (b,d) in zip(bandwidths, densities)
    Plots.plot!(p, d.x, d.density, labels = "bw = $b")
end
Plots.xlims!(p, 0.0, 8.0)
Plots.title!("KDE with Gaussian Kernel")
Plots.savefig(joinpath(base_img, "KDE.pdf"))


text:
using RecurrentNN
# reload("RecurrentNN.jl")

# # global settings
# const generator = "rnn" # can be 'rnn' or 'lstm'
srand(12345)
const generator = "lstm" # can be 'rnn' or 'lstm' or 'gru' or 'gflstm' or 'gfgru'
const hiddensizes = [100,100]
const lettersize = 7 # size of letter embeddings

# optimization
const regc = 0.000001 # L2 regularization strength
const learning_rate = 0.001 # Initial learning rate for lstm
const clipval = 5.0 # clip gradients at this value

function initVocab(inpath::String)

    f = open(inpath,"r")
    sents = [string(l[1:end-1]) for l in readlines(f)] # split(str,"\r\n") # array of sentences
    str = ""
    for i in 1:length(sents)
        s = sents[i]
        if contains(s,"\r")
            idx = findfirst(s,'\r')
            s = s[1:idx-1]
            sents[i] = s
        end
        str = "$str $(s[1:end-1])"
    end
    vocab = sort(setdiff(unique(str),['\r','\n'])) # unique characters in data
    inputsize = length(vocab) + 1 # 1 additional token (zero) in used for beginning and end tokens
    outputsize = length(vocab) + 1
    epochsize = length(sents) # nmber of sentence in sample

    # build char <-> index lookups
    letterToIndex =  Dict{Char,Int}()
    indexToLetter =  Dict{Int,Char}()
    [letterToIndex[vocab[i]] = i for i in 1:length(vocab)]
    [indexToLetter[i] = vocab[i] for i in 1:length(vocab)]

    return sents, vocab, letterToIndex, indexToLetter, inputsize, outputsize, epochsize
end

function initModel(inputsize::Int, lettersize::Int, hiddensizes::Array{Int,1},outputsize::Int)
    wil = randNNMat(inputsize,lettersize,0.08)
    nn = if generator == "rnn"
        RNN(lettersize,hiddensizes,outputsize)
    elseif generator == "gru"
        GRU(lettersize,hiddensizes,outputsize)
    elseif generator == "lstm"
        LSTM(lettersize,hiddensizes,outputsize)
    elseif generator == "gflstm"
        GFLSTM(lettersize,hiddensizes,outputsize)
    else
        GFGRU(lettersize,hiddensizes,outputsize)
    end
    # println((typeof(wil),typeof(nn)))
    return wil, nn
end

#########################################
#          initialize the model         #
#########################################

solver = Solver() # RMSProp optimizer

# init the text source
sents, vocab, letterToIndex, indexToLetter, inputsize, outputsize, epochsize =
    initVocab(joinpath(dirname(@__FILE__),"samples.txt"))

# init the rnn/lstm
wil, model = initModel(inputsize, lettersize, hiddensizes, outputsize)

tickiter = 0
pplmedian = Inf # perplexity will slowly move towards 0
pplcurve = Array(FloatingPoint,0) # track perplexity

#########################################
#             run the model             #
#########################################

function predictsentence(model::Model, wil::NNMatrix, samplei::Bool=false, temp::FloatingPoint=1.0)

    g = Graph(false) # backprop not needed
    s = ""
    max_chars_gen = 100
    prevhd   = Array(NNMatrix,0) # final hidden layer of the recurrent model after each forward step
    prevcell = Array(NNMatrix,0) # final cell output of the LSTM model after each forward step
    prevout  = NNMatrix(outputsize,1) # output of the recurrent model after each forward step
    prev = (prevhd, prevcell, prevout)
    cnt = 0
    while cnt < 100

        # RNN tick
        ix = length(s) == 0 ? 0 : letterToIndex[s[end]]
        x = rowpluck!(g, wil, ix+1) # get letter's embedding (vector)

        # returns a 2-tuple (RNN) or 3-tuple (LSTM). Last part is always outputNNMatrix
        prev = forwardprop!(g, model, x, prev)

        # sample predicted letter
        logprobs =  prev[end] # interpret output (last position in tuple) as logprobs
        if temp != 1.0 && samplei
            # scale log probabilities by temperature and renormalize
            # if temp is high, logprobs will go towards zero
            # and the softmax outputs will be more diffuse. if temp is
            # very low, the softmax outputs will be more peaky
            logprobs.w ./= temp
        end

        probs = softmax(logprobs)
        rndIX = rand()
        ix = samplei? findfirst(p-> p >= rndIX, cumsum(probs.w)) : indmax(probs.w);

        # break on out to the other side!
        if ix-1 == 0 break end # END token predicted, break out
        if length(s) > max_chars_gen break end # something is wrong

        letter = indexToLetter[ix-1]
        s = "$(s)$(letter)"
        cnt = cnt + 1
    end

    return s
end

function costfunc(model:: Model, wil::NNMatrix, sent::String)

    # takes a model and a sentence and
    # calculates the loss. Also returns the Graph
    # object which can be used to do backprop
    n = length(sent)
    g =  Graph()
    log2ppl = 0.0
    cost = 0.0
    prevhd   = Array(NNMatrix,0) # final hidden layer of the recurrent model after each forward step
    prevcell = Array(NNMatrix,0) # final cell output of the LSTM model after each forward step
    prevout  = NNMatrix(outputsize,1) # output of the recurrent model after each forward step
    prev = (prevhd, prevcell, prevout)
    for i= 0:length(sent)

        # start and end tokens are zeros
        ix_source = i == 0 ? 0 : letterToIndex[sent[i]] # first step: start with START token
        ix_target = i == n ? 0 : letterToIndex[sent[i+1]] # last step: end with END token

        # get the letter embbeding of the char
        x = rowpluck!(g, wil, ix_source+1)

        # returns a 2-tuple (RNN) or 3-tuples (LSTM). Last part of tuple is always the outputNNMatrix
        prev = forwardprop!(g, model, x, prev)

        # set gradients into logprobabilities
        logprobs =  prev[end] # interpret output (last position in tuple) as logprobs
        probs = softmax(logprobs) # compute the softmax probabilities

        log2ppl += -log2(probs.w[ix_target+1]) # accumulate base 2 log prob and do smoothing
        cost += -log(probs.w[ix_target+1])

        # write gradients into log probabilities
        logprobs.dw = probs.w
        logprobs.dw[ix_target+1] -= 1
    end
    ppl = (log2ppl/(n))^2
    return g, ppl, cost
end

function tick(model::Model, wil::NNMatrix, sents::Array, solver::Solver, tickiter::Int, pplcurve::Array{FloatingPoint,1}, pplmedian)

    # sample sentence fromd data
#     sent = sents[rand(1:21)] # use this if just kicking tires (faster)
    sent = sents[rand(1:length(sents))] # switch to this for a proper model (1-3hrs)

    t1 = time_ns() # log start timestamp

    # evaluate cost function on a sentence
    g, ppl, cost = costfunc(model,wil, sent)

    # use built up graph of backprop functions to compute backprop
    # i.e. set .dw fields in matrices
    backprop!(g)

    # perform param update ( learning_rate, regc, clipval are global constants)
    solverstats = step(solver, model, learning_rate, regc, clipval)

    tm = (time_ns()-t1)/1e9

    push!(pplcurve, ppl) #keep track of perplexity

    tickiter += 1;

    if tickiter % 50 == 0 # output sample sentences every Xth iteration
        pplmedian = round(median(pplcurve),4)
        println("Perplexity = $(pplmedian) @ $tickiter")
        push!(pplcurve, pplmedian)
        # draw samples to see how we're doing
        for i =1:2
            pred = predictsentence(model, wil, true, 0.7)
            println("   $(i). $pred ")
        end
        # greedy argmax (i.e if we were to select the most likely letter at each point)
        println("   Argmax: $(predictsentence(model, wil, false))")
    end

    return model, wil, solver, tickiter, pplcurve, pplmedian
end

maxIter = 100 # make this about 100_000 to run full model
trgppl = 1.1 # stop if this perplexity score is reached
@time while tickiter < maxIter && pplmedian > trgppl
    model, wil, solver, tickiter, pplcurve, pplmedian = tick(model, wil, sents, solver, tickiter, pplcurve, pplmedian)
end

# Profile.print(format=:flat)
#using ProfileView
#ProfileView.view()
#readline()

max_forks_repo_name: haroldsoh/RecurrentNN