text:
using Gen
using PyPlot
using Printf: @sprintf
import Random
import Distributions
using Statistics: median, mean
using JLD

include("scenes.jl")
include("path_planner.jl")
include("piecewise_normal.jl")

#####################################
# lightweight implementation of hmm #
#####################################

@gen function lightweight_hmm(step::Int, path::Path, distances_from_start::Vector{Float64},
                              times::Vector{Float64}, speed::Float64,
                              noise::Float64, dist_slack::Float64)
    @assert step >= 1

    # walk path
    locations = Vector{Point}(undef, step)
    dist = @trace(normal(speed * times[1], dist_slack), (:dist, 1))
    locations[1] = walk_path(path, distances_from_start, dist)
    for t=2:step
        dist = @trace(normal(dist + speed * (times[t] - times[t-1]), dist_slack), (:dist, t))
        locations[t] = walk_path(path, distances_from_start, dist)
    end

    # generate noisy observations
    for t=1:step
        point = locations[t]
        @trace(normal(point.x, noise), (:x, t))
        @trace(normal(point.y, noise), (:y, t))
    end

    return locations
end

################################################################
# lightweight and static implementations of hmm, with unfold #
################################################################

struct KernelState
    dist::Float64
    loc::Point
end

struct KernelParams
    times::Vector{Float64}
    path::Path
    distances_from_start::Vector{Float64}
    speed::Float64
    dist_slack::Float64
    noise::Float64
end

function dist_mean(prev_state::KernelState, params::KernelParams, t::Int)
    if t > 1
        prev_state.dist + (params.times[t] - params.times[t-1]) * params.speed
    else
        params.speed * params.times[1]
    end
end

@gen function lightweight_hmm_kernel(t::Int, prev_state::Any, params::KernelParams)
    # NOTE: if t = 1 then prev_state will have all NaNs
    dist = @trace(normal(dist_mean(prev_state, params, t), params.dist_slack), :dist)
    loc = walk_path(params.path, params.distances_from_start, dist)
    @trace(normal(loc.x, params.noise), :x)
    @trace(normal(loc.y, params.noise), :y)
    return KernelState(dist, loc)
end

lightweight_hmm_with_markov = Unfold(lightweight_hmm_kernel)

@gen (static) function static_hmm_kernel(t::Int, prev_state::KernelState, params::KernelParams)
    # NOTE: if t = 1 then prev_state will have all NaNs
    dist::Float64 = @trace(normal(dist_mean(prev_state, params, t), params.dist_slack), :dist)
    loc::Point = walk_path(params.path, params.distances_from_start, dist)
    @trace(normal(loc.x, params.noise), :x)
    @trace(normal(loc.y, params.noise), :y)
    ret::KernelState = KernelState(dist, loc)
    return ret
end

# construct a hidden markov model generator
# args: (num_steps::Int, init_state::KernelState, params::KernelParams)
static_hmm = Unfold(static_hmm_kernel)

####################
# custom proposals #
####################

function compute_custom_proposal_params(dt::Float64, prev_dist::Float64, noise::Float64, obs::Point,
                                        posterior_var_d::Float64, posterior_covars::Vector{Matrix{Float64}},
                                        path::Path, distances_from_start::Vector{Float64},
                                        speed::Float64, dist_slack::Float64)

    N = length(path.points)

    # Initialize parameters for truncated normal
    unnormalized_log_segment_probs = Vector{Float64}(undef, N+1)
    mus  = Vector{Float64}(undef, N+1)
    stds = Vector{Float64}(undef, N+1)

    # First segment
    mus[1] = prev_dist + dt * speed
    stds[1] = dist_slack
    unnormalized_log_segment_probs[1] = Distributions.logcdf(Distributions.Normal(mus[1], stds[1]), 0) + Distributions.logpdf(Distributions.Normal(path.start.x, noise), obs.x) + Distributions.logpdf(Distributions.Normal(path.start.y, noise), obs.y)

    # Last segment
    mus[N+1] = prev_dist + dt * speed
    stds[N+1] = dist_slack
    unnormalized_log_segment_probs[N+1] = Distributions.logccdf(Distributions.Normal(mus[N+1], stds[N+1]), distances_from_start[end]) +  Distributions.logpdf(Distributions.Normal(path.goal.x, noise), obs.x) + Distributions.logpdf(Distributions.Normal(path.goal.y, noise), obs.y)

    # Middle segments
    for i=2:N
        dx = path.points[i].x - path.points[i-1].x
        dy = path.points[i].y - path.points[i-1].y
        dd = distances_from_start[i] - distances_from_start[i-1]

        prior_mu_d = mus[1] - distances_from_start[i-1]
        x_obs = obs.x - path.points[i-1].x
        y_obs = obs.y - path.points[i-1].y

        posterior_mu_d = posterior_var_d * ((dx * x_obs + dy * y_obs) / (dd * noise^2) + prior_mu_d / dist_slack^2)

        mu_xy = prior_mu_d/dd .* [dx, dy]

        mus[i] = posterior_mu_d + distances_from_start[i-1]
        stds[i] = sqrt(posterior_var_d)
        unnormalized_log_segment_probs[i] = Distributions.logpdf(Distributions.MvNormal(mu_xy, posterior_covars[i]), [x_obs, y_obs]) + log(Distributions.cdf(Distributions.Normal(posterior_mu_d, stds[i]), dd) - Distributions.cdf(Distributions.Normal(posterior_mu_d, stds[i]), 0))
    end

    # You can think of the piecewise truncated normal as doing (1) categorical
    # draw, and (2) within that, a truncated normal draw.
    log_total_weight = logsumexp(unnormalized_log_segment_probs)
    log_normalized_weights = unnormalized_log_segment_probs .- log_total_weight
    probabilities = exp.(log_normalized_weights)

    return (probabilities, mus, stds, distances_from_start)
end

@gen function lightweight_fancy_step_proposal(step::Int, prev_dist::Float64, noise :: Float64, obs :: Point,
                                        posterior_var_d :: Float64, posterior_covars :: Vector{Matrix{Float64}},
                                        path :: Path, distances_from_start :: Vector{Float64},
                                        speed::Float64, times::Vector{Float64}, dist_slack::Float64)

    dt = step == 1 ? times[step] : times[step] - times[step-1]
    (probabilities, mus, stds, distances_from_start) = compute_custom_proposal_params(dt, prev_dist, noise, obs,
                                                            posterior_var_d, posterior_covars, path, distances_from_start,
                                                            speed, dist_slack)

    @trace(piecewise_normal(probabilities, mus, stds, distances_from_start), (:dist, step))
end

@gen function markov_fancy_proposal_inner(dt::Float64, prev_dist::Float64, noise :: Float64, obs :: Point,
                                        posterior_var_d :: Float64, posterior_covars :: Vector{Matrix{Float64}},
                                        path :: Path, distances_from_start :: Vector{Float64},
                                        speed::Float64, dist_slack::Float64)

    (probabilities, mus, stds, distances_from_start) = compute_custom_proposal_params(dt, prev_dist, noise, obs,
                                                            posterior_var_d, posterior_covars, path, distances_from_start,
                                                            speed, dist_slack)

    @trace(piecewise_normal(probabilities, mus, stds, distances_from_start), :dist)
end

lightweight_markov_fancy_proposal = at_dynamic(markov_fancy_proposal_inner, Int)

@gen (static) function static_markov_fancy_proposal_inner(dt::Float64, prev_dist::Float64, noise :: Float64, obs :: Point,
                                        posterior_var_d :: Float64, posterior_covars :: Vector{Matrix{Float64}},
                                        path :: Path, distances_from_start :: Vector{Float64},
                                        speed::Float64, dist_slack::Float64)

    dist_params = compute_custom_proposal_params(dt, prev_dist, noise, obs,
                                                 posterior_var_d, posterior_covars, path, distances_from_start,
                                                 speed, dist_slack)

    @trace(piecewise_normal(dist_params[1], dist_params[2], dist_params[3], dist_params[4]), :dist)
end

static_fancy_proposal = at_dynamic(static_markov_fancy_proposal_inner, Int)

#############################
# define some fixed context #
#############################

function make_scene()
    scene = Scene(0, 1, 0, 1)
    add!(scene, Tree(Point(0.30, 0.20), size=0.1))
    add!(scene, Tree(Point(0.83, 0.80), size=0.1))
    add!(scene, Tree(Point(0.80, 0.40), size=0.1))
    horiz = 1
    vert = 2
    wall_height = 0.30
    wall_thickness = 0.02
    walls = [
        Wall(Point(0.20, 0.40), horiz, 0.40, wall_thickness, wall_height)
        Wall(Point(0.60, 0.40), vert, 0.40, wall_thickness, wall_height)
        Wall(Point(0.60 - 0.15, 0.80), horiz, 0.15 + wall_thickness, wall_thickness, wall_height)
        Wall(Point(0.20, 0.80), horiz, 0.15, wall_thickness, wall_height)
        Wall(Point(0.20, 0.40), vert, 0.40, wall_thickness, wall_height)]
    for wall in walls
        add!(scene, wall)
    end
    return scene
end

const scene = make_scene()
const times = collect(range(0, stop=1, length=20))
const start_x = 0.1
const start_y = 0.1
const stop_x = 0.5
const stop_y = 0.5
const speed = 0.5
const noise = 0.02
const dist_slack = 0.2
const start = Point(start_x, start_y)
const stop = Point(stop_x, stop_y)

####################
# trace renderings #
####################

# the only difference between rendering the lightweight and static versions
# how the locations are extracted from the trace, and the addresses of the
# measurements.

function render_lightweight_hmm_trace(scene::Scene, start::Point, stop::Point,
                maybe_path::Nullable{Path},
                times::Vector{Float64}, speed::Float64,
                noise::Float64, dist_slack::Float64, trace, ax;
                show_measurements=true,
                show_start=true, show_stop=true,
                show_path=true, show_noise=true,
                start_alpha=1., stop_alpha=1., path_alpha=1.)

    # set current axis
    sca(ax)

    # render obstacles
    render(scene, ax)

    # plot start and stop
    if show_start
        scatter([start.x], [start.y], color="blue", s=100, alpha=start_alpha)
    end
    if show_stop
        scatter([stop.x], [stop.y], color="red", s=100, alpha=stop_alpha)
    end

    # plot path lines
    if !isnull(maybe_path) && show_path
        path = get(maybe_path)
        for i=1:length(path.points)-1
            prev = path.points[i]
            next = path.points[i+1]
            plot([prev.x, next.x], [prev.y, next.y], color="black", alpha=0.5, linewidth=5)
        end
    end

    # plot locations with measurement noise around them
    locations = get_call_record(trace).retval
    scatter([pt.x for pt in locations], [pt.y for pt in locations],
        color="orange", alpha=path_alpha, s=25)
    if show_noise
        for pt in locations
            circle = patches[:Circle]((pt.x, pt.y), noise, facecolor="purple", alpha=0.2)
            ax[:add_patch](circle)
        end
    end

    # plot measured locations
    assignment = get_choices(trace)
    if show_measurements
        measured_xs = [assignment[(:x, i)] for i=1:length(locations)]
        measured_ys = [assignment[(:y, i)] for i=1:length(locations)]
        scatter(measured_xs, measured_ys, marker="x", color="black", alpha=1., s=25)
    end
end

function render_static_hmm_trace(scene::Scene, start::Point, stop::Point,
                maybe_path::Nullable{Path},
                times::Vector{Float64}, speed::Float64,
                noise::Float64, dist_slack::Float64, trace, ax;
                show_measurements=true,
                show_start=true, show_stop=true,
                show_path=true, show_noise=true,
                start_alpha=1., stop_alpha=1., path_alpha=1.)

    # set current axis
    sca(ax)

    # render obstacles
    render(scene, ax)

    # plot start and stop
    if show_start
        scatter([start.x], [start.y], color="blue", s=100, alpha=start_alpha)
    end
    if show_stop
        scatter([stop.x], [stop.y], color="red", s=100, alpha=stop_alpha)
    end

    # plot path lines
    if !isnull(maybe_path) && show_path
        path = get(maybe_path)
        for i=1:length(path.points)-1
            prev = path.points[i]
            next = path.points[i+1]
            plot([prev.x, next.x], [prev.y, next.y], color="black", alpha=0.5, linewidth=5)
        end
    end

    # plot locations with measurement noise around them
    states = get_call_record(trace).retval
    locations = Point[state.loc for state in states]
    scatter([pt.x for pt in locations], [pt.y for pt in locations],
        color="orange", alpha=path_alpha, s=25)
    if show_noise
        for pt in locations
            circle = patches[:Circle]((pt.x, pt.y), noise, facecolor="purple", alpha=0.2)
            ax[:add_patch](circle)
        end
    end

    # plot measured locations
    assignment = get_choices(trace)
    if show_measurements
        measured_xs = [assignment[i => :x] for i=1:length(locations)]
        measured_ys = [assignment[i => :y] for i=1:length(locations)]
        scatter(measured_xs, measured_ys, marker="x", color="black", alpha=1., s=25)
    end
end

######################
# show prior samples #
######################

function show_prior_samples()
    Random.seed!(0)

    # generate a path
    maybe_path = plan_path(start, stop, scene, PlannerParams(300, 3.0, 2000, 1.))
    @assert !isnull(maybe_path)
    path = get(maybe_path)
    distances_from_start = compute_distances_from_start(path)

    # show samples from the lightweight model
    model_args = (length(times), path, distances_from_start, times, speed, noise, dist_slack)
    figure(figsize=(32, 32))
    for i=1:15
        subplot(4, 4, i)
        ax = gca()
        trace = simulate(lightweight_hmm, model_args)
        render_lightweight_hmm_trace(scene, start, stop, maybe_path, times, speed, noise, dist_slack, trace, ax)
    end
    savefig("lightweight_hmm_prior_samples.png")

    # show samples from the static model
    kernel_params = KernelParams(times, path, distances_from_start, speed, dist_slack, noise)
    model_args_rest = (KernelState(NaN, Point(NaN, NaN)), kernel_params)
    figure(figsize=(32, 32))
    for i=1:15
        subplot(4, 4, i)
        ax = gca()
        trace = simulate(static_hmm, (length(times), model_args_rest...))
        render_static_hmm_trace(scene, start, stop, maybe_path, times, speed, noise, dist_slack, trace, ax)
    end
    savefig("static_hmm_prior_samples.png")

end

##################################
# particle filtering experiments #
##################################

struct Params
    times::Vector{Float64}
    speed::Float64
    dist_slack::Float64
    noise::Float64
    path::Path
end

struct PrecomputedPathData
    posterior_var_d::Float64
    posterior_covars::Vector{Matrix{Float64}}
    distances_from_start::Vector{Float64}
end

function PrecomputedPathData(params::Params)
    times = params.times
    speed = params.speed
    dist_slack = params.dist_slack
    noise = params.noise
    path = params.path

    distances_from_start = compute_distances_from_start(path)

    # posterior_var_d is the variance of the posterior on d' given x and y.
    # posterior_covars is a vector of 2x2 covariance matrices, representing the
    # covariance of x and y when d' has been marginalized out.
    posterior_var_d = dist_slack^2 * noise^2 / (dist_slack^2 + noise^2)
    posterior_covars = Vector{Matrix{Float64}}(undef, length(distances_from_start))
    for i = 2:length(distances_from_start)
        dd = distances_from_start[i] - distances_from_start[i-1]
        dx = path.points[i].x - path.points[i-1].x
        dy = path.points[i].y - path.points[i-1].y
        posterior_covars[i] = [noise 0; 0 noise] .+ (dist_slack^2/dd^2 .* [dx^2 dx*dy; dy*dx dy^2])
    end

    PrecomputedPathData(posterior_var_d, posterior_covars, distances_from_start)
end

### (1) static, with markov, default proposal ###

function particle_filtering_static_hmm_default_proposal(params::Params,
        measured_xs::Vector{Float64}, measured_ys::Vector{Float64},
        num_particles_list::Vector{Int}, num_reps::Int)

    println("static hmm, default proposal")

    precomputed = PrecomputedPathData(params)
    times = params.times
    speed = params.speed
    dist_slack = params.dist_slack
    noise = params.noise
    path = params.path

    function obs(step::Int)
        assignment = choicemap()
        assignment[step => :x] = measured_xs[step]
        assignment[step => :y] = measured_ys[step]
        return (assignment, UnfoldCustomArgDiff(true, false, false))
    end

    max_steps = length(times)
    kernel_params = KernelParams(times, path,
                                 precomputed.distances_from_start,
                                 speed, dist_slack, noise)
    model_args_rest = (KernelState(NaN, Point(NaN, NaN)), kernel_params)
    results = Dict{Int, Tuple{Vector{Float64},Vector{Float64}}}()
    for num_particles in num_particles_list
        ess_threshold = num_particles / 2
        elapsed = Vector{Float64}(undef, num_reps)
        lmls = Vector{Float64}(undef, num_reps)
        for rep=1:num_reps
            start = time_ns()
            (traces, _, lml) = particle_filter(static_hmm, model_args_rest, max_steps,
                                               num_particles, ess_threshold, obs)
            elapsed[rep] = Int(time_ns() - start) / 1e9
            println("num_particles: $num_particles, lml estimate: $lml, elapsed: $(elapsed[rep])")
            lmls[rep] = lml
        end
        results[num_particles] = (lmls, elapsed)
    end

    return results
end

### (2) static, with markov, custom proposal ###

function particle_filtering_static_hmm_custom_proposal(params::Params,
        measured_xs::Vector{Float64}, measured_ys::Vector{Float64},
        num_particles_list::Vector{Int}, num_reps::Int)

    println("static hmm, custom proposal")

    precomputed = PrecomputedPathData(params)
    times = params.times
    speed = params.speed
    dist_slack = params.dist_slack
    noise = params.noise
    path = params.path

    function init()
        observations = choicemap()
        observations[1 => :x] = measured_xs[1]
        observations[1 => :y] = measured_ys[1]
        proposal_args = (1, (times[1], Float64(0.0), noise,
                         Point(measured_xs[1], measured_ys[1]),
                         precomputed.posterior_var_d, precomputed.posterior_covars, path,
                         precomputed.distances_from_start, speed, dist_slack))
        return (observations, proposal_args)
    end

    function step(step::Int, trace)
        @assert step > 1
        observations = choicemap()
        observations[step => :x] = measured_xs[step]
        observations[step => :y] = measured_ys[step]
        prev_dist = get_choices(trace)[step-1 => :dist]
        dt = step==1 ? times[step] : times[step] - times[step-1]
        proposal_args = (step, (dt, prev_dist, noise,
                         Point(measured_xs[step], measured_ys[step]),
                         precomputed.posterior_var_d, precomputed.posterior_covars,  path,
                         precomputed.distances_from_start, speed, dist_slack))
        return (observations, proposal_args, UnfoldCustomArgDiff(true, false, false))
    end

    max_steps = length(times)
    kernel_params = KernelParams(times, params.path,
                                 precomputed.distances_from_start,
                                 speed, dist_slack, noise)
    model_args_rest = (KernelState(NaN, Point(NaN, NaN)), kernel_params)
    results = Dict{Int, Tuple{Vector{Float64},Vector{Float64}}}()
    for num_particles in num_particles_list
        ess_threshold = num_particles / 2
        elapsed = Vector{Float64}(undef, num_reps)
        lmls = Vector{Float64}(undef, num_reps)
        for rep=1:num_reps
            start = time_ns()
            (traces, _, lml) = particle_filter(static_hmm, model_args_rest, max_steps,
                                               num_particles, ess_threshold,
                                               init, step,
                                               static_fancy_proposal,
                                               static_fancy_proposal)
            elapsed[rep] = Int(time_ns() - start) / 1e9
            println("num_particles: $num_particles, lml estimate: $lml, elapsed: $(elapsed[rep])")
            lmls[rep] = lml
        end
        results[num_particles] = (lmls, elapsed)
    end

    return results
end

### (3) lightweight, default proposal ###

function particle_filtering_lightweight_hmm_default_proposal(params::Params,
        measured_xs::Vector{Float64}, measured_ys::Vector{Float64},
        num_particles_list::Vector{Int}, num_reps::Int)

    println("lightweight hmm, default proposal")

    precomputed = PrecomputedPathData(params)
    times = params.times
    speed = params.speed
    dist_slack = params.dist_slack
    noise = params.noise
    path = params.path

    function obs(step::Int)
        observations = choicemap()
        observations[(:x, step)] = measured_xs[step]
        observations[(:y, step)] = measured_ys[step]
        return (observations, unknownargdiff)
    end

    max_steps = length(times)
    model_args_rest = (path, precomputed.distances_from_start, times, speed, noise, dist_slack)
    results = Dict{Int, Tuple{Vector{Float64},Vector{Float64}}}()
    for num_particles in num_particles_list
        ess_threshold = num_particles / 2
        elapsed = Vector{Float64}(undef, num_reps)
        lmls = Vector{Float64}(undef, num_reps)
        for rep=1:num_reps
            start = time_ns()
            (traces, _, lml) = particle_filter(lightweight_hmm, model_args_rest, max_steps,
                                               num_particles, ess_threshold, obs)
            elapsed[rep] = Int(time_ns() - start) / 1e9
            println("num_particles: $num_particles, lml estimate: $lml, elapsed: $(elapsed[rep])")
            lmls[rep] = lml
        end
        results[num_particles] = (lmls, elapsed)
    end

    return results
end

### (4) lightweight, custom proposal ###

function particle_filtering_lightweight_hmm_custom_proposal(params::Params,
        measured_xs::Vector{Float64}, measured_ys::Vector{Float64},
        num_particles_list::Vector{Int}, num_reps::Int)

    println("lightweight hmm, custom proposal")

    precomputed = PrecomputedPathData(params)
    times = params.times
    speed = params.speed
    dist_slack = params.dist_slack
    noise = params.noise
    path = params.path

    function init()
        observations = choicemap()
        observations[(:x, 1)] = measured_xs[1]
        observations[(:y, 1)] = measured_ys[1]
        proposal_args = (1, Float64(0.0), noise,
                         Point(measured_xs[1], measured_ys[1]),
                         precomputed.posterior_var_d, precomputed.posterior_covars,  path,
                         precomputed.distances_from_start, speed, times, dist_slack)
        return (observations, proposal_args)
    end

    function step(step::Int, trace)
        @assert step > 1
        observations = choicemap()
        observations[(:x, step)] = measured_xs[step]
        observations[(:y, step)] = measured_ys[step]
        prev_dist = get_choices(trace)[(:dist, step-1)]
        proposal_args = (step, prev_dist, noise,
                         Point(measured_xs[step], measured_ys[step]),
                         precomputed.posterior_var_d, precomputed.posterior_covars,  path,
                         precomputed.distances_from_start, speed, times, dist_slack)
        return (observations, proposal_args, unknownargdiff)
    end

    max_steps = length(times)
    kernel_params = KernelParams(times, params.path,
                                 precomputed.distances_from_start,
                                 speed, dist_slack, noise)
    model_args_rest = (path, precomputed.distances_from_start, times, speed, noise, dist_slack)
    results = Dict{Int, Tuple{Vector{Float64},Vector{Float64}}}()
    for num_particles in num_particles_list
        ess_threshold = num_particles / 2
        elapsed = Vector{Float64}(undef, num_reps)
        lmls = Vector{Float64}(undef, num_reps)
        for rep=1:num_reps
            start = time_ns()
            (traces, _, lml) = particle_filter(lightweight_hmm, model_args_rest, max_steps,
                                               num_particles, ess_threshold,
                                               init, step,
                                               lightweight_fancy_step_proposal,
                                               lightweight_fancy_step_proposal)
            elapsed[rep] = Int(time_ns() - start) / 1e9
            println("num_particles: $num_particles, lml estimate: $lml, elapsed: $(elapsed[rep])")
            lmls[rep] = lml
        end
        results[num_particles] = (lmls, elapsed)
    end

    return results
end

### (5) lightweight, with markov, default proposal ###

function particle_filtering_lightweight_markov_hmm_default_proposal(params::Params,
        measured_xs::Vector{Float64}, measured_ys::Vector{Float64},
        num_particles_list::Vector{Int}, num_reps::Int)

    println("lightweight hmm with markov, default proposal")

    precomputed = PrecomputedPathData(params)
    times = params.times
    speed = params.speed
    dist_slack = params.dist_slack
    noise = params.noise
    path = params.path

    function obs(step::Int)
        assignment = choicemap()
        assignment[step => :x] = measured_xs[step]
        assignment[step => :y] = measured_ys[step]
        return (assignment, UnfoldCustomArgDiff(true, false, false))
    end

    max_steps = length(times)
    kernel_params = KernelParams(times, path,
                                 precomputed.distances_from_start,
                                 speed, dist_slack, noise)
    model_args_rest = (KernelState(NaN, Point(NaN, NaN)), kernel_params)
    results = Dict{Int, Tuple{Vector{Float64},Vector{Float64}}}()
    for num_particles in num_particles_list
        ess_threshold = num_particles / 2
        elapsed = Vector{Float64}(undef, num_reps)
        lmls = Vector{Float64}(undef, num_reps)
        for rep=1:num_reps
            start = time_ns()
            (traces, _, lml) = particle_filter(lightweight_hmm_with_markov, model_args_rest, max_steps,
                                               num_particles, ess_threshold, obs)
            elapsed[rep] = Int(time_ns() - start) / 1e9
            println("num_particles: $num_particles, lml estimate: $lml, elapsed: $(elapsed[rep])")
            lmls[rep] = lml
        end
        results[num_particles] = (lmls, elapsed)
    end

    return results
end

### (6) lightweight, with markov, custom proposal ###

function particle_filtering_lightweight_markov_hmm_custom_proposal(params::Params,
        measured_xs::Vector{Float64}, measured_ys::Vector{Float64},
        num_particles_list::Vector{Int}, num_reps::Int)

    println("lightweight hmm with markov, custom proposal")

    precomputed = PrecomputedPathData(params)
    times = params.times
    speed = params.speed
    dist_slack = params.dist_slack
    noise = params.noise
    path = params.path

    function init()
        observations = choicemap()
        observations[1 => :x] = measured_xs[1]
        observations[1 => :y] = measured_ys[1]
        proposal_args = (1, (times[1], Float64(0.0), noise,
                         Point(measured_xs[1], measured_ys[1]),
                         precomputed.posterior_var_d, precomputed.posterior_covars, path,
                         precomputed.distances_from_start, speed, dist_slack))
        return (observations, proposal_args)
    end

    function step(step::Int, trace)
        @assert step > 1
        observations = choicemap()
        observations[step => :x] = measured_xs[step]
        observations[step => :y] = measured_ys[step]
        prev_dist = get_choices(trace)[step-1 => :dist]
        dt = step==1 ? times[step] : times[step] - times[step-1]
        proposal_args = (step, (dt, prev_dist, noise,
                         Point(measured_xs[step], measured_ys[step]),
                         precomputed.posterior_var_d, precomputed.posterior_covars,  path,
                         precomputed.distances_from_start, speed, dist_slack))
        return (observations, proposal_args, UnfoldCustomArgDiff(true, false, false))
    end

    max_steps = length(times)
    kernel_params = KernelParams(times, params.path,
                                 precomputed.distances_from_start,
                                 speed, dist_slack, noise)
    model_args_rest = (KernelState(NaN, Point(NaN, NaN)), kernel_params)
    results = Dict{Int, Tuple{Vector{Float64},Vector{Float64}}}()
    for num_particles in num_particles_list
        ess_threshold = num_particles / 2
        elapsed = Vector{Float64}(undef, num_reps)
        lmls = Vector{Float64}(undef, num_reps)
        for rep=1:num_reps
            start = time_ns()
            (traces, _, lml) = particle_filter(lightweight_hmm_with_markov, model_args_rest, max_steps,
                                               num_particles, ess_threshold,
                                               init, step,
                                               lightweight_markov_fancy_proposal,
                                               lightweight_markov_fancy_proposal)
            elapsed[rep] = Int(time_ns() - start) / 1e9
            println("num_particles: $num_particles, lml estimate: $lml, elapsed: $(elapsed[rep])")
            lmls[rep] = lml
        end
        results[num_particles] = (lmls, elapsed)
    end

    return results
end

function plot_results(results::Dict, num_particles_list::Vector{Int}, label::String, color::String)
    median_elapsed = [median(results[num_particles][2]) for num_particles in num_particles_list]
    mean_lmls = [mean(results[num_particles][1]) for num_particles in num_particles_list]
    plot(median_elapsed, mean_lmls, label=label, color=color)
end

function experiment()

    Random.seed!(0)

    # generate a path
    maybe_path = plan_path(start, stop, scene, PlannerParams(300, 3.0, 2000, 1.))
    @assert !isnull(maybe_path)
    path = get(maybe_path)

    println("path:")
    println(path)

    # precomputation
    params = Params(times, speed, dist_slack, noise, path)
    precomputed = PrecomputedPathData(params)

    # generate ground truth locations and observations
    model_args = (length(times), path, precomputed.distances_from_start, times, speed, noise, dist_slack)
    trace = simulate(lightweight_hmm, model_args)
    assignment = get_choices(trace)
    measured_xs = [assignment[(:x, i)] for i=1:length(times)]
    measured_ys = [assignment[(:y, i)] for i=1:length(times)]
    actual_dists = [assignment[(:dist, i)] for i=1:length(times)]

    println("measured_xs:")
    println(measured_xs)

    println("measured_ys:")
    println(measured_ys)

    # parameters for particle filtering
    num_particles_list = [1, 10, 30]#[1, 2, 3, 5, 7, 10, 15, 20, 25, 30, 35, 40, 45, 50, 60, 70, 80, 90, 100, 200, 300]
    #num_particles_list = [1, 2, 3, 5, 7, 10, 15, 20, 25, 30, 35, 40, 45, 50, 60, 70, 80, 90, 100, 200, 300]
    #num_reps = 50
    num_reps = 20

    # experiments with static model
    results_static_default_proposal = particle_filtering_static_hmm_default_proposal(params,
        measured_xs, measured_ys, num_particles_list, num_reps)
    results_static_custom_proposal = particle_filtering_static_hmm_custom_proposal(params,
        measured_xs, measured_ys, num_particles_list, num_reps)

    # experiments with lightweight model (no markov)
    results_lightweight_default_proposal = particle_filtering_lightweight_hmm_default_proposal(params,
        measured_xs, measured_ys, num_particles_list, num_reps)
    results_lightweight_custom_proposal = particle_filtering_lightweight_hmm_custom_proposal(params,
        measured_xs, measured_ys, num_particles_list, num_reps)

    # experiments with markov
    results_lightweight_markov_default_proposal = particle_filtering_lightweight_markov_hmm_default_proposal(params,
        measured_xs, measured_ys, num_particles_list, num_reps)
    results_lightweight_markov_custom_proposal = particle_filtering_lightweight_markov_hmm_custom_proposal(params,
        measured_xs, measured_ys, num_particles_list, num_reps)

    save("results.jld",
        "results_static_default_proposal", results_static_default_proposal,
        "results_static_custom_proposal", results_static_custom_proposal,
        "results_lightweight_default_proposal", results_lightweight_default_proposal,
        "results_lightweight_custom_proposal", results_lightweight_custom_proposal,
        "results_lightweight_markov_default_proposal", results_lightweight_markov_default_proposal,
        "results_lightweight_markov_custom_proposal", results_lightweight_markov_custom_proposal)
end

# NOTE: you have to separate static gen function definitions from calling API
# methods on them with a top-level call to this function
Gen.load_generated_functions()

# show_prior_samples()

experiment()


text: print("\nPlease install ZeroMQ using your system's package manager\n")

ext: jl


text:
const maxMantissa = (UInt64(1) << 53) - 1
todouble(sign, exp, mant) = Core.bitcast(Float64, (UInt64(sign) << 63) | (UInt64(exp) << 52) | (UInt64(mant)))

@testset "Parsers.writeshortest" begin

@testset "Float64" begin

@testset "Basic" begin
    @test Parsers.writeshortest(0.0) == "0.0"
    @test Parsers.writeshortest(-0.0) == "-0.0"
    @test Parsers.writeshortest(1.0) == "1.0"
    @test Parsers.writeshortest(-1.0) == "-1.0"
    @test Parsers.writeshortest(NaN) == "NaN"
    @test Parsers.writeshortest(Inf) == "Inf"
    @test Parsers.writeshortest(-Inf) == "-Inf"
end

@testset "SwitchToSubnormal" begin
    @test "2.2250738585072014e-308" == Parsers.writeshortest(2.2250738585072014e-308)
end

@testset "MinAndMax" begin
    @test "1.7976931348623157e308" == Parsers.writeshortest(Core.bitcast(Float64, 0x7fefffffffffffff))
    @test "5.0e-324" == Parsers.writeshortest(Core.bitcast(Float64, Int64(1)))
end

@testset "LotsOfTrailingZeros" begin
    @test "2.9802322387695312e-8" == Parsers.writeshortest(2.98023223876953125e-8)
end

@testset "Regression" begin
    @test "-2.109808898695963e16" == Parsers.writeshortest(-2.109808898695963e16)
    @test "4.940656e-318" == Parsers.writeshortest(4.940656e-318)
    @test "1.18575755e-316" == Parsers.writeshortest(1.18575755e-316)
    @test "2.989102097996e-312" == Parsers.writeshortest(2.989102097996e-312)
    @test "9.0608011534336e15" == Parsers.writeshortest(9.0608011534336e15)
    @test "4.708356024711512e18" == Parsers.writeshortest(4.708356024711512e18)
    @test "9.409340012568248e18" == Parsers.writeshortest(9.409340012568248e18)
    @test "1.2345678" == Parsers.writeshortest(1.2345678)
end

@testset "LooksLikePow5" begin
    # These numbers have a mantissa that is a multiple of the largest power of 5 that fits,
    # and an exponent that causes the computation for q to result in 22, which is a corner
    # case for Ryu.
    @test "5.764607523034235e39" == Parsers.writeshortest(Core.bitcast(Float64, 0x4830F0CF064DD592))
    @test "1.152921504606847e40" == Parsers.writeshortest(Core.bitcast(Float64, 0x4840F0CF064DD592))
    @test "2.305843009213694e40" == Parsers.writeshortest(Core.bitcast(Float64, 0x4850F0CF064DD592))
end

@testset "OutputLength" begin
    @test "1.0" == Parsers.writeshortest(1.0) # already tested in Basic
    @test "1.2" == Parsers.writeshortest(1.2)
    @test "1.23" == Parsers.writeshortest(1.23)
    @test "1.234" == Parsers.writeshortest(1.234)
    @test "1.2345" == Parsers.writeshortest(1.2345)
    @test "1.23456" == Parsers.writeshortest(1.23456)
    @test "1.234567" == Parsers.writeshortest(1.234567)
    @test "1.2345678" == Parsers.writeshortest(1.2345678) # already tested in Regressi
    @test "1.23456789" == Parsers.writeshortest(1.23456789)
    @test "1.234567895" == Parsers.writeshortest(1.234567895) # 1.234567890 would be trimm
    @test "1.2345678901" == Parsers.writeshortest(1.2345678901)
    @test "1.23456789012" == Parsers.writeshortest(1.23456789012)
    @test "1.234567890123" == Parsers.writeshortest(1.234567890123)
    @test "1.2345678901234" == Parsers.writeshortest(1.2345678901234)
    @test "1.23456789012345" == Parsers.writeshortest(1.23456789012345)
    @test "1.234567890123456" == Parsers.writeshortest(1.234567890123456)
    @test "1.2345678901234567" == Parsers.writeshortest(1.2345678901234567)

  # Test 32-bit chunking
    @test "4.294967294" == Parsers.writeshortest(4.294967294) # 2^32 -
    @test "4.294967295" == Parsers.writeshortest(4.294967295) # 2^32 -
    @test "4.294967296" == Parsers.writeshortest(4.294967296) # 2^
    @test "4.294967297" == Parsers.writeshortest(4.294967297) # 2^32 +
    @test "4.294967298" == Parsers.writeshortest(4.294967298) # 2^32 +
end

# Test min, max shift values in shiftright128
@testset "MinMaxShift" begin
    # 32-bit opt-size=0:  49 <= dist <= 50
    # 32-bit opt-size=1:  30 <= dist <= 50
    # 64-bit opt-size=0:  50 <= dist <= 50
    # 64-bit opt-size=1:  30 <= dist <= 50
    @test "1.7800590868057611e-307" == Parsers.writeshortest(todouble(false, 4, 0))
    # 32-bit opt-size=0:  49 <= dist <= 49
    # 32-bit opt-size=1:  28 <= dist <= 49
    # 64-bit opt-size=0:  50 <= dist <= 50
    # 64-bit opt-size=1:  28 <= dist <= 50
    @test "2.8480945388892175e-306" == Parsers.writeshortest(todouble(false, 6, maxMantissa))
    # 32-bit opt-size=0:  52 <= dist <= 53
    # 32-bit opt-size=1:   2 <= dist <= 53
    # 64-bit opt-size=0:  53 <= dist <= 53
    # 64-bit opt-size=1:   2 <= dist <= 53
    @test "2.446494580089078e-296" == Parsers.writeshortest(todouble(false, 41, 0))
    # 32-bit opt-size=0:  52 <= dist <= 52
    # 32-bit opt-size=1:   2 <= dist <= 52
    # 64-bit opt-size=0:  53 <= dist <= 53
    # 64-bit opt-size=1:   2 <= dist <= 53
    @test "4.8929891601781557e-296" == Parsers.writeshortest(todouble(false, 40, maxMantissa))

    # 32-bit opt-size=0:  57 <= dist <= 58
    # 32-bit opt-size=1:  57 <= dist <= 58
    # 64-bit opt-size=0:  58 <= dist <= 58
    # 64-bit opt-size=1:  58 <= dist <= 58
    @test "1.8014398509481984e16" == Parsers.writeshortest(todouble(false, 1077, 0))
    # 32-bit opt-size=0:  57 <= dist <= 57
    # 32-bit opt-size=1:  57 <= dist <= 57
    # 64-bit opt-size=0:  58 <= dist <= 58
    # 64-bit opt-size=1:  58 <= dist <= 58
    @test "3.6028797018963964e16" == Parsers.writeshortest(todouble(false, 1076, maxMantissa))
    # 32-bit opt-size=0:  51 <= dist <= 52
    # 32-bit opt-size=1:  51 <= dist <= 59
    # 64-bit opt-size=0:  52 <= dist <= 52
    # 64-bit opt-size=1:  52 <= dist <= 59
    @test "2.900835519859558e-216" == Parsers.writeshortest(todouble(false, 307, 0))
    # 32-bit opt-size=0:  51 <= dist <= 51
    # 32-bit opt-size=1:  51 <= dist <= 59
    # 64-bit opt-size=0:  52 <= dist <= 52
    # 64-bit opt-size=1:  52 <= dist <= 59
    @test "5.801671039719115e-216" == Parsers.writeshortest(todouble(false, 306, maxMantissa))

    # https:#github.com/ulfjack/ryu/commit/19e44d16d80236f5de25800f56d82606d1be00b9#commitcomment-30146483
    # 32-bit opt-size=0:  49 <= dist <= 49
    # 32-bit opt-size=1:  44 <= dist <= 49
    # 64-bit opt-size=0:  50 <= dist <= 50
    # 64-bit opt-size=1:  44 <= dist <= 50
    @test "3.196104012172126e-27" == Parsers.writeshortest(todouble(false, 934, 0x000FA7161A4D6E0C))
end

@testset "SmallIntegers" begin
    @test "9.007199254740991e15" == Parsers.writeshortest(9007199254740991.0)
    @test "9.007199254740992e15" == Parsers.writeshortest(9007199254740992.0)

    @test "1.0" == Parsers.writeshortest(1.0e+0)
    @test "12.0" == Parsers.writeshortest(1.2e+1)
    @test "123.0" == Parsers.writeshortest(1.23e+2)
    @test "1234.0" == Parsers.writeshortest(1.234e+3)
    @test "12345.0" == Parsers.writeshortest(1.2345e+4)
    @test "123456.0" == Parsers.writeshortest(1.23456e+5)
    @test "1.234567e6" == Parsers.writeshortest(1.234567e+6)
    @test "1.2345678e7" == Parsers.writeshortest(1.2345678e+7)
    @test "1.23456789e8" == Parsers.writeshortest(1.23456789e+8)
    @test "1.23456789e9" == Parsers.writeshortest(1.23456789e+9)
    @test "1.234567895e9" == Parsers.writeshortest(1.234567895e+9)
    @test "1.2345678901e10" == Parsers.writeshortest(1.2345678901e+10)
    @test "1.23456789012e11" == Parsers.writeshortest(1.23456789012e+11)
    @test "1.234567890123e12" == Parsers.writeshortest(1.234567890123e+12)
    @test "1.2345678901234e13" == Parsers.writeshortest(1.2345678901234e+13)
    @test "1.23456789012345e14" == Parsers.writeshortest(1.23456789012345e+14)
    @test "1.234567890123456e15" == Parsers.writeshortest(1.234567890123456e+15)

  # 10^i
    @test "1.0" == Parsers.writeshortest(1.0e+0)
    @test "10.0" == Parsers.writeshortest(1.0e+1)
    @test "100.0" == Parsers.writeshortest(1.0e+2)
    @test "1000.0" == Parsers.writeshortest(1.0e+3)
    @test "10000.0" == Parsers.writeshortest(1.0e+4)
    @test "100000.0" == Parsers.writeshortest(1.0e+5)
    @test "1.0e6" == Parsers.writeshortest(1.0e+6)
    @test "1.0e7" == Parsers.writeshortest(1.0e+7)
    @test "1.0e8" == Parsers.writeshortest(1.0e+8)
    @test "1.0e9" == Parsers.writeshortest(1.0e+9)
    @test "1.0e10" == Parsers.writeshortest(1.0e+10)
    @test "1.0e11" == Parsers.writeshortest(1.0e+11)
    @test "1.0e12" == Parsers.writeshortest(1.0e+12)
    @test "1.0e13" == Parsers.writeshortest(1.0e+13)
    @test "1.0e14" == Parsers.writeshortest(1.0e+14)
    @test "1.0e15" == Parsers.writeshortest(1.0e+15)

  # 10^15 + 10^i
    @test "1.000000000000001e15" == Parsers.writeshortest(1.0e+15 + 1.0e+0)
    @test "1.00000000000001e15" == Parsers.writeshortest(1.0e+15 + 1.0e+1)
    @test "1.0000000000001e15" == Parsers.writeshortest(1.0e+15 + 1.0e+2)
    @test "1.000000000001e15" == Parsers.writeshortest(1.0e+15 + 1.0e+3)
    @test "1.00000000001e15" == Parsers.writeshortest(1.0e+15 + 1.0e+4)
    @test "1.0000000001e15" == Parsers.writeshortest(1.0e+15 + 1.0e+5)
    @test "1.000000001e15" == Parsers.writeshortest(1.0e+15 + 1.0e+6)
    @test "1.00000001e15" == Parsers.writeshortest(1.0e+15 + 1.0e+7)
    @test "1.0000001e15" == Parsers.writeshortest(1.0e+15 + 1.0e+8)
    @test "1.000001e15" == Parsers.writeshortest(1.0e+15 + 1.0e+9)
    @test "1.00001e15" == Parsers.writeshortest(1.0e+15 + 1.0e+10)
    @test "1.0001e15" == Parsers.writeshortest(1.0e+15 + 1.0e+11)
    @test "1.001e15" == Parsers.writeshortest(1.0e+15 + 1.0e+12)
    @test "1.01e15" == Parsers.writeshortest(1.0e+15 + 1.0e+13)
    @test "1.1e15" == Parsers.writeshortest(1.0e+15 + 1.0e+14)

  # Largest power of 2 <= 10^(i+1)
    @test "8.0" == Parsers.writeshortest(8.0)
    @test "64.0" == Parsers.writeshortest(64.0)
    @test "512.0" == Parsers.writeshortest(512.0)
    @test "8192.0" == Parsers.writeshortest(8192.0)
    @test "65536.0" == Parsers.writeshortest(65536.0)
    @test "524288.0" == Parsers.writeshortest(524288.0)
    @test "8.388608e6" == Parsers.writeshortest(8388608.0)
    @test "6.7108864e7" == Parsers.writeshortest(67108864.0)
    @test "5.36870912e8" == Parsers.writeshortest(536870912.0)
    @test "8.589934592e9" == Parsers.writeshortest(8589934592.0)
    @test "6.8719476736e10" == Parsers.writeshortest(68719476736.0)
    @test "5.49755813888e11" == Parsers.writeshortest(549755813888.0)
    @test "8.796093022208e12" == Parsers.writeshortest(8796093022208.0)
    @test "7.0368744177664e13" == Parsers.writeshortest(70368744177664.0)
    @test "5.62949953421312e14" == Parsers.writeshortest(562949953421312.0)
    @test "9.007199254740992e15" == Parsers.writeshortest(9007199254740992.0)

  # 1000 * (Largest power of 2 <= 10^(i+1))
    @test "8000.0" == Parsers.writeshortest(8.0e+3)
    @test "64000.0" == Parsers.writeshortest(64.0e+3)
    @test "512000.0" == Parsers.writeshortest(512.0e+3)
    @test "8.192e6" == Parsers.writeshortest(8192.0e+3)
    @test "6.5536e7" == Parsers.writeshortest(65536.0e+3)
    @test "5.24288e8" == Parsers.writeshortest(524288.0e+3)
    @test "8.388608e9" == Parsers.writeshortest(8388608.0e+3)
    @test "6.7108864e10" == Parsers.writeshortest(67108864.0e+3)
    @test "5.36870912e11" == Parsers.writeshortest(536870912.0e+3)
    @test "8.589934592e12" == Parsers.writeshortest(8589934592.0e+3)
    @test "6.8719476736e13" == Parsers.writeshortest(68719476736.0e+3)
    @test "5.49755813888e14" == Parsers.writeshortest(549755813888.0e+3)
    @test "8.796093022208e15" == Parsers.writeshortest(8796093022208.0e+3)
end

end # Float64

@testset "Float32" begin

@testset "Basic" begin
    @test "0.0" == Parsers.writeshortest(Float32(0.0))
    @test "-0.0" == Parsers.writeshortest(Float32(-0.0))
    @test "1.0" == Parsers.writeshortest(Float32(1.0))
    @test "-1.0" == Parsers.writeshortest(Float32(-1.0))
    @test "NaN" == Parsers.writeshortest(Float32(NaN))
    @test "Inf" == Parsers.writeshortest(Float32(Inf))
    @test "-Inf" == Parsers.writeshortest(Float32(-Inf))
end

@testset "SwitchToSubnormal" begin
    @test "1.1754944e-38" == Parsers.writeshortest(1.1754944f-38)
end

@testset "MinAndMax" begin
    @test "3.4028235e38" == Parsers.writeshortest(Core.bitcast(Float32, 0x7f7fffff))
    @test "1.0e-45" == Parsers.writeshortest(Core.bitcast(Float32, Int32(1)))
end

# Check that we return the exact boundary if it is the shortest
# representation, but only if the original floating point number is even.
@testset "BoundaryRoundeven" begin
    @test "3.355445e7" == Parsers.writeshortest(3.355445f7)
    @test "9.0e9" == Parsers.writeshortest(8.999999f9)
    @test "3.436672e10" == Parsers.writeshortest(3.4366717f10)
end

# If the exact value is exactly halfway between two shortest representations,
# then we round to even. It seems like this only makes a difference if the
# last two digits are ...2|5 or ...7|5, and we cut off the 5.
@testset "exactValueRoundeven" begin
    @test "305404.12" == Parsers.writeshortest(3.0540412f5)
    @test "8099.0312" == Parsers.writeshortest(8.0990312f3)
end

@testset "LotsOfTrailingZeros" begin
    # Pattern for the first test: 00111001100000000000000000000000
    @test "0.00024414062" == Parsers.writeshortest(2.4414062f-4)
    @test "0.0024414062" == Parsers.writeshortest(2.4414062f-3)
    @test "0.0043945312" == Parsers.writeshortest(4.3945312f-3)
    @test "0.0063476562" == Parsers.writeshortest(6.3476562f-3)
end

@testset "Regression" begin
    @test "4.7223665e21" == Parsers.writeshortest(4.7223665f21)
    @test "8.388608e6" == Parsers.writeshortest(8388608f0)
    @test "1.6777216e7" == Parsers.writeshortest(1.6777216f7)
    @test "3.3554436e7" == Parsers.writeshortest(3.3554436f7)
    @test "6.7131496e7" == Parsers.writeshortest(6.7131496f7)
    @test "1.9310392e-38" == Parsers.writeshortest(1.9310392f-38)
    @test "-2.47e-43" == Parsers.writeshortest(-2.47f-43)
    @test "1.993244e-38" == Parsers.writeshortest(1.993244f-38)
    @test "4103.9004" == Parsers.writeshortest(4103.9003f0)
    @test "5.3399997e9" == Parsers.writeshortest(5.3399997f9)
    @test "6.0898e-39" == Parsers.writeshortest(6.0898f-39)
    @test "0.0010310042" == Parsers.writeshortest(0.0010310042f0)
    @test "2.882326e17" == Parsers.writeshortest(2.8823261f17)
    @test "7.038531e-26" == Parsers.writeshortest(7.0385309f-26)
    @test "9.223404e17" == Parsers.writeshortest(9.2234038f17)
    @test "6.710887e7" == Parsers.writeshortest(6.7108872f7)
    @test "1.0e-44" == Parsers.writeshortest(1.0f-44)
    @test "2.816025e14" == Parsers.writeshortest(2.816025f14)
    @test "9.223372e18" == Parsers.writeshortest(9.223372f18)
    @test "1.5846086e29" == Parsers.writeshortest(1.5846085f29)
    @test "1.1811161e19" == Parsers.writeshortest(1.1811161f19)
    @test "5.368709e18" == Parsers.writeshortest(5.368709f18)
    @test "4.6143166e18" == Parsers.writeshortest(4.6143165f18)
    @test "0.007812537" == Parsers.writeshortest(0.007812537f0)
    @test "1.0e-45" == Parsers.writeshortest(1.4f-45)
    @test "1.18697725e20" == Parsers.writeshortest(1.18697724f20)
    @test "1.00014165e-36" == Parsers.writeshortest(1.00014165f-36)
    @test "200.0" == Parsers.writeshortest(200f0)
    @test "3.3554432e7" == Parsers.writeshortest(3.3554432f7)
end

@testset "LooksLikePow5" begin
    # These numbers have a mantissa that is the largest power of 5 that fits,
    # and an exponent that causes the computation for q to result in 10, which is a corner
    # case for Ryu.
    @test "6.7108864e17" == Parsers.writeshortest(Core.bitcast(Float32, 0x5D1502F9))
    @test "1.3421773e18" == Parsers.writeshortest(Core.bitcast(Float32, 0x5D9502F9))
    @test "2.6843546e18" == Parsers.writeshortest(Core.bitcast(Float32, 0x5E1502F9))
end

@testset "OutputLength" begin
    @test "1.0" == Parsers.writeshortest(Float32(1.0))
    @test "1.2" == Parsers.writeshortest(Float32(1.2))
    @test "1.23" == Parsers.writeshortest(Float32(1.23))
    @test "1.234" == Parsers.writeshortest(Float32(1.234))
    @test "1.2345" == Parsers.writeshortest(Float32(1.2345))
    @test "1.23456" == Parsers.writeshortest(Float32(1.23456))
    @test "1.234567" == Parsers.writeshortest(Float32(1.234567))
    @test "1.2345678" == Parsers.writeshortest(Float32(1.2345678))
    @test "1.23456735e-36" == Parsers.writeshortest(Float32(1.23456735e-36))
end

end # Float32

@testset "Float16" begin

@testset "Basic" begin
    @test "0.0" == Parsers.writeshortest(Float16(0.0))
    @test "-0.0" == Parsers.writeshortest(Float16(-0.0))
    @test "1.0" == Parsers.writeshortest(Float16(1.0))
    @test "-1.0" == Parsers.writeshortest(Float16(-1.0))
    @test "NaN" == Parsers.writeshortest(Float16(NaN))
    @test "Inf" == Parsers.writeshortest(Float16(Inf))
    @test "-Inf" == Parsers.writeshortest(Float16(-Inf))
end

let x=floatmin(Float16)
    while x <= floatmax(Float16)
        @test parse(Float16, Parsers.writeshortest(x)) == x
        x = nextfloat(x)
    end
end

# function testfloats(T)
#     x = floatmin(T)
#     i = 0
#     fails = 0
#     success = 0
#     while x < floatmax(T)
#         test = parse(T, Parsers.writeshortest(x)) == x
#         if !test

#             fails += 1
#         else
#             success += 1
#         end
#         x = nextfloat(x)
#         i += 1

#     end
#     return fails / (fails + success)
# end

end # Float16

end


text:
using SimpleDiffEq, Test

prob = DiscreteProblem(0.5,(0.0,1.0))
sol =solve(prob,SimpleFunctionMap())

@test sol[1] == sol[end]

integrator = init(prob,SimpleFunctionMap())
step!(integrator)

@test integrator.u == integrator.uprev
@test integrator.t == 1.0

prob2 = DiscreteProblem(rand(4,2),(0.0,1.0))
sol = solve(prob2,SimpleFunctionMap())

@test sol[1] == sol[end]

integrator = init(prob,SimpleFunctionMap())
step!(integrator)

@test integrator.u == integrator.uprev
@test integrator.t == 1.0

prob = DiscreteProblem((u,p,t)->1.01u,0.5,(0.0,1.0))
sol = solve(prob,SimpleFunctionMap())

@test sol[end] ≈ .505

integrator = init(prob,SimpleFunctionMap())
step!(integrator)

@test integrator.u ≈ .505
@test integrator.uprev ≈ 0.5
@test integrator.t == 1.0

prob = DiscreteProblem((du,u,p,t)->du .= 1.01.*u,rand(4,2),(0.0,1.0))
sol = solve(prob,SimpleFunctionMap())

@test sol[end]./sol[1] ≈ fill(1.01,4,2)

integrator = init(prob,SimpleFunctionMap())
step!(integrator)

@test integrator.u ≈ prob.u0*1.01
@test integrator.uprev ≈ prob.u0
@test integrator.t == 1.0

@test_broken step!(integrator)


text:
using JuMP
using CPLEX
include("LSP_PLNE.jl")
include("InstanceLoader.jl")
include("ResolvePlne.jl")
include("Helper.jl")
include("VRP_Heuristic.jl")
include("VRP_PLNE.jl")

#=
Initialise la résolution heuristique du PDI :
 --> créer le modèle (PLNE) LSP
 --> y ajoute les SC dans la fonction objectif
 --> y ajoute les nouvelles variables de décision z_i,t binaire
 --> créer les contraintes liant les z_i,t et q_i,t
=#
function initialisation_PDI_heuristique(INSTANCE_PATH)
        #Récuperer l'instance
    params, nodes, demands, costs = readPRP(INSTANCE_PATH)
        #Créer le PLNE LSP associé à l'instance
    lsp_model = createLSP(params, nodes, demands, costs)
        #Creer les SC
    SC=Array{Int,2}(undef,params["n"],params["l"])
        #Prendre la fonction objectif
    fonctionObj=objective_function(lsp_model)
    fonctionObjInitial=objective_function(lsp_model)

    #####Initialisation des variables binaires z_i,t#####
    @variable(lsp_model,z[1:params["n"], 1:params["l"]],Bin)

    for revendeur in 1:params["n"]
        for pasDeTemps in 1:params["l"]
    ######Couplage entre q_i,t et z_i,t#####
            #si z=1, on peut produire pour le revendeur autant que params["C]
            #qui est la capacité de production du fournisseur, sinon on ne produit pas
            @constraint(lsp_model,variable_by_name(lsp_model, "q[$revendeur,$pasDeTemps]")<=z[revendeur,pasDeTemps]*params["C"])
            #si on ne produit pas alors z est forcément nul.
            @constraint(lsp_model,z[revendeur,pasDeTemps]<=variable_by_name(lsp_model, "q[$revendeur,$pasDeTemps]"))

            #initialiser les SC_i,t
            SC[revendeur,pasDeTemps]=costs[0,revendeur]+costs[revendeur,0]

            #Changer la fonction objectif pour ajouter les coûts de visite
            fonctionObj+=SC[revendeur,pasDeTemps]*z[revendeur,pasDeTemps]
        end
    end

    #mettre la nouvelle fonction objectif
    set_objective_function(lsp_model,fonctionObj)

    return lsp_model, params, nodes, demands, costs, SC, fonctionObjInitial
end

#=
    Crée un modèle (PLNE) de VRP
=#
function getVRPfromLSP(lsp_model,params, nodes, demands, costs)
    qteALivrer = Array{Float64, 2}(undef, params["n"], params["l"]) #qteALivrer[i,t] = qté à livrer au revendeur i au pas de temps t
	notEmptyIndexAtT = [[] for i in 1:params["l"]] #Pour stocker les revendeurs qui ont été livré au pas de temps t

	for i in 1:params["n"]
		for t in 1:params["l"]
			qteALivrer[i, t] = value(variable_by_name(lsp_model, "q[$i,$t]")) #stocker les qté à livrer pour tous les revendeurs, pour tous les pas de temps
			if qteALivrer[i, t] != 0 #si on doit livrer quelque chose i au pas de temps t, on l'ajoute dans la liste
				push!(notEmptyIndexAtT[t], i)
			end

		end
	end

	tToVrp = Any[]
	for t in 1:params["l"]
		#On copie les paramètres mais n prend la valeur du nombre de noeud avec une demande supérieur à 0
		copyParams = copy(params)
		copyParams["n"] = size(notEmptyIndexAtT[t], 1)

		#On copie les informations des noeuds mais dans sans les noeud avec une demande 0
		copyNodes = Array{Dict, 1}(undef, size(notEmptyIndexAtT[t], 1) + 1) # copyNodes[i]= le ième noeud ayant une qté à livrer non nulle

		#On copie les demande à chaque temps mais sans les noeud avec une demande à 0
		copyQteALivrer = Array{Float64, 2}(undef, size(notEmptyIndexAtT[t], 1), params["l"]) #copyQteALivrer[i,t] = la qté à livrer du ième noeud ayant une qté à livrer non nulle

		#Initialise le noeud 0 à part
		copyNodes[1] = nodes[0]

		for i in eachindex(notEmptyIndexAtT[t])
			copyNodes[i+1] = nodes[notEmptyIndexAtT[t][i]]# copyNodes[i]= le ième noeud ayant une qté à livrer non nulle
			copyQteALivrer[i, :] = qteALivrer[notEmptyIndexAtT[t][i], :]#copyQteALivrer[i,t] = la qté à livrer du ième noeud ayant une qté à livrer non nulle
		end

		#On remet les index en commençant à 0
		copyNodes = OffsetVector(copyNodes, 0:(size(copyNodes, 1) - 1))

		#Un dictionnaire faisant le lien entre les anciens index et les nouveaux (Dû au décalage d'index dans les array)
		isPresent = Dict{Int, Int}()

		for (index, elem) in enumerate(notEmptyIndexAtT[t])
			isPresent[elem] = index
		end

		isPresent[0] = 0

		#On copie les coûts de transports mais sans les arêtes passant par un noeud avec une demande à 0
		copyCost = Dict{Tuple{Int, Int}, Float64}()
		for (edge, edgeCost) in costs

			if in(edge[1], keys(isPresent)) && in(edge[2], keys(isPresent)) #si les deux noeuds formant une arête sont présent dans la liste des noeuds ayant une qté à livrer non nulle
				copyCost[(isPresent[edge[1]], isPresent[edge[2]])] = edgeCost
			end

		end

		push!(tToVrp, (copyParams, copyNodes, copyQteALivrer, copyCost))
        #=
        On fait un VRP avec un nombre de noeuds n réduit (< au n initial) où:
        On ne garde que les noeuds dont on doit livrer à un certain pas de temps, on a gardé dans le dictionnaire des noeuds l'index d'origine pour mieux se retrouver
        On ne garde que les qté à livrer non nulle pour chaque pas de temps
        on ne garde que les arcs et les coûts de ces arcs que s'ils ne sont reliés à deux noeuds auxquels on doit livrer (sous graphe induit par les noeuds qu'on doit livrer)
        =#

	end

    #On renvoi les modèles de VRP à résoudre
    vrp_models_et_parametres_en_fonction_du_pasDeTemps=Dict(t => (createVRP_MTZ(tToVrp[t][1],tToVrp[t][2],tToVrp[t][3],tToVrp[t][4],t),tToVrp[t]) for t in 1:params["l"])
	return vrp_models_et_parametres_en_fonction_du_pasDeTemps

end

function getSCfromVRPCircuits(vrp_circuits, SC, pasDeTemps, params, costs)
    revendeur_visites=Dict{Int,Tuple}() #nodeID => (predecesseur, successeur)=(i-,i+)
    for circuit in vrp_circuits #on trouve tous les noeuds visités ainsi que ses successeur et predecesseur dans son circuit
        for (index,revendeur) in enumerate(circuit)
            if index==1 #le revendeur est 0 (c'est le fournisseur), rappel : on a toujour circuit[1]=0
                revendeur_visites[revendeur]=(circuit[length(circuit)],circuit[index+1])
            elseif index==length(circuit)
                revendeur_visites[revendeur]=(circuit[index-1],circuit[1])
            else
                revendeur_visites[revendeur]=(circuit[index-1],circuit[index+1])
            end
        end
    end
    for (revendeur,(predecesseur, successeur)) in revendeur_visites #on met à jour les SC pour les visités
        SC[revendeur,pasDeTemps]=costs[predecesseur,revendeur]+costs[revendeur,successeur]-costs[predecesseur,successeur]
    end
    for nodeToInsert in 1:params["n"] #on met à jour les SC pour les non visités
        if !(nodeToInsert in keys(revendeur_visites))#si le revendeur n'est pas visité
            #coût si on inserait nodeToInsert après revendeur
            coutInsertion=[cost[revendeur,nodeToInsert]+cost[nodeToInsert,successeur]-cost[revendeur,successeur] for (revendeur,(predecesseur, successeur)) in revendeur_visites]
            #SC[nodeToInsert,pasDeTemps]=min([coutInsertion]) #marche pas ... chelou
            min=coutInsertion[1]
            for cout in coutInsertion
                min=min<cout ? min : cout
            end
            SC[nodeToInsert,pasDeTemps]=min
        end
    end
    return SC
end

function update_LSP(lsp_model,SC,fonctionObjInitial,params)
    fonctionObj=fonctionObjInitial
    for revendeur in 1:params["n"]
        for pasDeTemps in 1:params["l"]
            fonctionObj+=SC[revendeur,pasDeTemps]*variable_by_name(lsp_model, "z[$revendeur,$pasDeTemps]")
        end
    end

    #mettre la nouvelle fonction objectif
    set_objective_function(lsp_model,fonctionObj)

end

function PDI_heuristique(lsp_model, params, nodes, demands, costs, SC , fonctionObjInitial, nbMaxIte=10, resoudreVRPwithHeuristic=true)
    circuits=Dict()
    for ite in 1:nbMaxIte
        #résolution LSP
        println("============================================================================SOLVING============================================================================")
        println("============================================================================SOLVING============================================================================")
        println("============================================================================SOLVING============================================================================")
        println("============================================================================SOLVING============================================================================\n\n\n")
        #lsp_model=resolvePlne(lsp_model,2,"")
        optimize!(lsp_model)

        #récuperer les VRP pour chaque pas de temps
        vrp_models_et_parametres_en_fonction_du_pasDeTemps=getVRPfromLSP(lsp_model,params, nodes, demands, costs)

        #résolution VRP à chaque pas de temps
        allCircuits=[]
        for pasDeTemps in 1:params["l"]

            #résoudre VRP
            if(resoudreVRPwithHeuristic) #si on fait avec l'heuristique clark wright
                vrp_model,(copyParams, copyNodes, copyQteALivrer, copyCost)=vrp_models_et_parametres_en_fonction_du_pasDeTemps[pasDeTemps]
                #vrp_circuits_temp=clark_wright(copyParams, copyNodes, copyQteALivrer, copyCost, pasDeTemps)
                vrp_circuits_temp=mixMetaheuristic(clark_wright(copyParams, copyNodes, copyQteALivrer, copyCost, pasDeTemps),copyParams,copyCost)
            else #si on fait une résolution exacte
                #récuperer le VRP pour le pas de temps correspondant ainsi les params, nodes, demands, costs de ce VRP
                vrp_model,(copyParams, copyNodes, copyQteALivrer, copyCost)=vrp_models_et_parametres_en_fonction_du_pasDeTemps[pasDeTemps]
                #résolution exacte de VRP
                vrp_model=resolvePlne(vrp_model,0,"")
                #prendre les circuits à partir du modèle de VRP, ils sont décrit avec les index locaux
                vrp_circuits_temp=vrpToCircuit(vrp_model, copyParams)
            end

            #retrouver les indexOriginaux
            originalIndexes=Dict(i-1 => e["initial_index"] for (i,e) in enumerate(copyNodes))
            #revenir aux index initiaux
            vrp_circuits=[]
            for circuit in vrp_circuits_temp
                circuit_with_Original_index=[originalIndexes[k] for k in circuit]
                push!(vrp_circuits,circuit_with_Original_index)
            end
            #prendre les SC_i,pasDeTemps depuis les circuits calculé dans le VRP
            if(empty!(vrp_circuits)==false)
                SC=getSCfromVRPCircuits(vrp_circuits, SC, pasDeTemps, params, costs)
            end
            circuits[pasDeTemps]=vrp_circuits
        end

        #mettre à jour les SC dans le modèle du LSP
        if(ite<nbMaxIte)
            update_LSP(lsp_model,SC,fonctionObjInitial,params)
        end

    end
    return lsp_model,circuits
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -16.424759279035918 <= q <= 20.23452089871623)

                     add_NL_constraint(m, :(gelu(0.8420310557117703 + 0.34198141401583015*gelu(-0.5389217244927917 + 0.5798025520034438*gelu(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])) + 0.4428632765655749*gelu(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])) + 0.9033224385695577*gelu(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))) + -0.1887519896617813*gelu(-0.23921485861490321 + 0.5030319708107833*gelu(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])) + -0.2081481501163327*gelu(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])) + -0.2939216839670169*gelu(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))) + -0.2882490310085508*gelu(-0.6545847926105979 + 0.9774628615021363*gelu(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])) + -0.21180887804098658*gelu(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])) + -0.7013494059895882*gelu(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))) + gelu(-0.35188345522271236 + -0.11415509244934263*gelu(-0.5389217244927917 + 0.5798025520034438*gelu(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])) + 0.4428632765655749*gelu(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])) + 0.9033224385695577*gelu(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))) + -0.8088728193024375*gelu(-0.23921485861490321 + 0.5030319708107833*gelu(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])) + -0.2081481501163327*gelu(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])) + -0.2939216839670169*gelu(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))) + 0.8271426964720727*gelu(-0.6545847926105979 + 0.9774628615021363*gelu(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])) + -0.21180887804098658*gelu(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])) + -0.7013494059895882*gelu(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))) + gelu(0.4798469380896484 + 0.8134725364327324*gelu(-0.5389217244927917 + 0.5798025520034438*gelu(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])) + 0.4428632765655749*gelu(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])) + 0.9033224385695577*gelu(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))) + -0.12585836968848696*gelu(-0.23921485861490321 + 0.5030319708107833*gelu(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])) + -0.2081481501163327*gelu(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])) + -0.2939216839670169*gelu(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))) + 0.613011911779838*gelu(-0.6545847926105979 + 0.9774628615021363*gelu(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])) + -0.21180887804098658*gelu(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])) + -0.7013494059895882*gelu(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
struct DielectricDistribution{T, N, S} <: AbstractArray{T, N}
    data::Array{T, N}
    grid::Grid{T, N, S}
end

@inline size(ϵ::DielectricDistribution{T, N, S}) where {T, N, S} = size(ϵ.data)
@inline length(ϵ::DielectricDistribution{T, N, S}) where {T, N, S} = length(ϵ.data)
@inline getindex(ϵ::DielectricDistribution{T, N, S}, I::Vararg{Int, N}) where {T, N, S} = getindex(ϵ.data, I...)
@inline getindex(ϵ::DielectricDistribution{T, N, S}, i::Int) where {T, N, S} = getindex(ϵ.data, i)
@inline getindex(ϵ::DielectricDistribution{T, N, S}, s::Symbol) where {T, N, S} = getindex(ϵ.grid, s)

function DielectricDistribution(fss::PotentialSimulationSetup{T, N, S})::DielectricDistribution{T, N, S} where {T, N, S}
    return DielectricDistribution{T, N, S}( fss.ϵ, fss.grid )
end

function NamedTuple(ϵ::DielectricDistribution{T}) where {T <: SSDFloat}
    return (
        grid = NamedTuple(ϵ.grid),
        values = ϵ.data * internal_voltage_unit,
    )
end
Base.convert(T::Type{NamedTuple}, x::DielectricDistribution) = T(x)

function DielectricDistribution(nt::NamedTuple)
    grid = Grid(nt.grid)
    T = typeof(ustrip(nt.values[1]))
    S = get_coordinate_type(grid)
    N = get_number_of_dimensions(grid)
    DielectricDistribution{T, N, S}( ustrip.(uconvert.(internal_voltage_unit, nt.values)), grid)
end
Base.convert(T::Type{DielectricDistribution}, x::NamedTuple) = T(x)

@recipe function f( ϵ::DielectricDistribution{T, 3, :cylindrical};
                    r = missing,
                    φ = missing,
                    z = missing ) where {T}
    g::Grid{T, 3, :cylindrical} = ϵ.grid

    # seriescolor --> :viridis
    st --> :heatmap
    aspect_ratio --> 1
    foreground_color_border --> nothing
    tick_direction --> :out

    cross_section::Symbol, idx::Int = if ismissing(φ) && ismissing(r) && ismissing(z)
        :φ, 1
    elseif !ismissing(φ) && ismissing(r) && ismissing(z)
        φ_rad::T = T(deg2rad(φ))
        while !(g[:φ].interval.left <= φ_rad <= g[:φ].interval.right)
            if φ_rad > g[:φ].interval.right
                φ_rad -= g[:φ].interval.right - g[:φ].interval.left
            elseif φ_rad < g[:φ].interval.left
                φ_rad += g[:φ].interval.right - g[:φ].interval.left
            end
        end
        :φ, searchsortednearest(g[:φ], φ_rad)
    elseif ismissing(φ) && !ismissing(r) && ismissing(z)
        :r, searchsortednearest(g[:r], T(r))
    elseif ismissing(φ) && ismissing(r) && !ismissing(z)
        :z, searchsortednearest(g[:z], T(z))
    else
        error(ArgumentError, ": Only one of the keywords `r, φ, z` is allowed.")
    end
    value::T = if cross_section == :φ
        g[:φ][idx]
    elseif cross_section == :r
        g[:r][idx]
    elseif cross_section == :z
        g[:z][idx]
    end

    @series begin
        if cross_section == :φ
            title --> "Dielectric Distribution @$(cross_section) = $(round(rad2deg(value), sigdigits = 2))"
            xlabel --> "r / m"
            ylabel --> "z / m"
            size --> ( 400, 350 / (g[:r][end] - g[:r][1]) * (g[:z][end] - g[:z][1]) )
            g[:r], g[:z], ϵ.data[:, idx,:]'
        elseif cross_section == :r
            title --> "Dielectric Distribution @$(cross_section) = $(round(value, sigdigits = 2))"
            g[:φ], g[:z], ϵ.data[idx,:,:]'
        elseif cross_section == :z
            title --> "Dielectric Distribution @$(cross_section) = $(round(value, sigdigits = 2))"
            proj --> :polar
            g[:φ], g[:r], ϵ.data[:,:,idx]
        end
    end
end

@recipe function f( ϵ::DielectricDistribution{T, 3, :cartesian};
                    x = missing,
                    y = missing,
                    z = missing ) where {T}
    g::Grid{T, 3, :cartesian} = ϵ.grid

    seriescolor --> :viridis
    st --> :heatmap
    aspect_ratio --> 1
    foreground_color_border --> nothing
    tick_direction --> :out

    cross_section::Symbol, idx::Int = if ismissing(x) && ismissing(y) && ismissing(z)
        :x, 1
    elseif !ismissing(x) && ismissing(y) && ismissing(z)
        :x, searchsortednearest(g[:x], T(x))
    elseif ismissing(x) && !ismissing(y) && ismissing(z)
        :y, searchsortednearest(g[:y], T(y))
    elseif ismissing(x) && ismissing(y) && !ismissing(z)
        :z, searchsortednearest(g[:z], T(z))
    else
        error(ArgumentError, ": Only one of the keywords `r, y, z` is allowed.")
    end
    value::T = if cross_section == :x
        g[:x][idx]
    elseif cross_section == :y
        g[:y][idx]
    elseif cross_section == :z
        g[:z][idx]
    end
    @series begin
            title --> "Dielectric Distribution @$(cross_section) = $(round(value, sigdigits = 2))"
        if cross_section == :x
            xlabel --> "y / m"
            ylabel --> "z / m"
            xlims --> (g[:y][1], g[:y][end])
            ylims --> (g[:z][1], g[:z][end])
            g[:y], g[:z], ϵ.data[idx, :, :]'
        elseif cross_section == :y
            xlabel --> "x / m"
            ylabel --> "z / m"
            xlims --> (g[:x][1], g[:x][end])
            ylims --> (g[:z][1], g[:z][end])
            g[:x], g[:z], ϵ.data[:, idx, :]'
        elseif cross_section == :z
            xlabel --> "x / m"
            ylabel --> "y / m"
            xlims --> (g[:x][1], g[:x][end])
            ylims --> (g[:y][1], g[:y][end])
            g[:x], g[:y], ϵ.data[:, :, idx]
        end
    end
end


text:
module SkeletonTriangulationTests

using Gridap

model = CartesianDiscreteModel(partition=(3,4,2))

tags = [27,]
grid = SkeletonGrid(model,tags)

trian = Triangulation(grid)
descr1 = grid.descriptor1
descr2 = grid.descriptor2

strian = SkeletonTriangulation(trian,descr1,descr2)

test_triangulation(strian)

strian = SkeletonTriangulation(model)
#writevtk(strian,"strian")

end # module


text:
using JSON
using uJSON
using Base.Test

function print_object_summary(ob, max_depth::Int, depth=0)
	if isa(ob, Union(Dict, Array))
		println(string(typeof(ob)), "(len ", length(ob), "): ")
		if depth*4 + 1 < max_depth*4
			d2 = depth + 1
			if isa(ob, Dict)
				for (key, ob2) in ob
					print(" "^4d2, key, ": ")
					print_object_summary(ob2, max_depth, d2)
				end
			elseif isa(ob,Array)
				for ob2 in ob
					print(" "^4d2)
					print_object_summary(ob2, max_depth, d2)
				end
			end
		end
	else
		println("(", string(typeof(ob)), ") ", string(ob))
	end
end

function generate_jagged(l1, l2)
	data = Any[]
	arrays_only = true
	randrand(maxlen) = rand(int(rand()*maxlen))
	for r in rand(int(l1))
	    if r < 0.25
	    	if arrays_only
	        	push!(data, string(randrand(l2)))
	        else
	        	push!(data, {string(i)=>int(i*10) for i = randrand(l2)})
	    	end
	    elseif r < 0.5
	        push!(data, randrand(l2))
	    elseif r < 0.75
	        push!(data, int(10*randrand(l2)))
	    else
	    	if arrays_only
	        	push!(data, [nothing for i = randrand(l2)])
	    	else
	        	push!(data, {string(int(i))=>nothing for i = 10*randrand(l2)})
	        end
	    end
	end
	data
end
data = generate_jagged(1e5, 500)
json_str = json(data)
fn = "testing.json"
f=open(fn, "w")
write(f, json_str)
close(f)

print("parsing with uJSON:      ")
@time uJSON.parse(json_str)
print("parsing with uJSON:      ")
@time uJSON.parse(json_str)
print("parsing with JSON:       ")
@time JSON.parse(json_str)
print("parsing with JSON:       ")
@time JSON.parse(json_str)
# @test uresult == result

fn = Pkg.dir("uJSON", "test", "sample_u4c.json")

print("JSON running deep test:  ")
@time result = JSON.parse(open(fn, "r"))
print("uJSON running deep test: ")
@time uresult = uJSON.parse(open(fn, "r"))

# @profile result = uJSON.parse(json_str)
# # Profile.print()
# using ProfileView
# ProfileView.view()
# readline(STDIN)


text:
export IndiaNSECalendar

abstract IndiaCalendar <: WesternCalendar

type IndiaNSECalendar <: IndiaCalendar; end

string(::IndiaNSECalendar) = "National Stock Exchange of India (Equities Trading)"

function isWorkingDay(c::IndiaNSECalendar, dt::CalendarTime)
	d=day(dt)
	dd=dayofyear(dt)
	m=month(dt)
	w=dayofweek(dt)
	y=year(dt)
	em=easterMonday(c,y)

	if (isWeekend(c,w)
            # Republic Day
            || (d == 26 && m == January)
            # Good Friday
            || (dd == em-3)
            # Ambedkar Jayanti
            || (d == 14 && m == April)
            # Independence Day
            || (d == 15 && m == August)
            # Gandhi Jayanti
            || (d == 2 && m == October)
            # Christmas
            || (d == 25 && m == December))

            return false
    end

         if (y == 2005)
            # Moharram, Holi, Maharashtra Day, and Ramzan Id fall
            # on Saturday or Sunday in 2005
            if (# Bakri Id
                (d == 21 && m == January)
                # Ganesh Chaturthi
                || (d == 7 && m == September)
                # Dasara
                || (d == 12 && m == October)
                # Laxmi Puja
                || (d == 1 && m == November)
                # Bhaubeej
                || (d == 3 && m == November)
                # Guru Nanak Jayanti
                || (d == 15 && m == November))

                return false
            end
        end

        if (y == 2006)
            if (# Bakri Id
                (d == 11 && m == January)
                # Moharram
                || (d == 9 && m == February)
                # Holi
                || (d == 15 && m == March)
                # Ram Navami
                || (d == 6 && m == April)
                # Mahavir Jayanti
                || (d == 11 && m == April)
                # Maharashtra Day
                || (d == 1 && m == May)
                # Bhaubeej
                || (d == 24 && m == October)
                # Ramzan Id
                || (d == 25 && m == October))

                return false
            end
        end

        if (y == 2007)
            if (# Bakri Id
                (d == 1 && m == January)
                # Moharram
                || (d == 30 && m == January)
                # Mahashivratri
                || (d == 16 && m == February)
                # Ram Navami
                || (d == 27 && m == March)
                # Maharashtra Day
                || (d == 1 && m == May)
                # Buddha Pournima
                || (d == 2 && m == May)
                # Laxmi Puja
                || (d == 9 && m == November)
                # Bakri Id (again)
                || (d == 21 && m == December))

                return false
            end
        end

        if (y == 2008)
            if (# Mahashivratri
                (d == 6 && m == March)
                # Id-E-Milad
                || (d == 20 && m == March)
                # Mahavir Jayanti
                || (d == 18 && m == April)
                # Maharashtra Day
                || (d == 1 && m == May)
                # Buddha Pournima
                || (d == 19 && m == May)
                # Ganesh Chaturthi
                || (d == 3 && m == September)
                # Ramzan Id
                || (d == 2 && m == October)
                # Dasara
                || (d == 9 && m == October)
                # Laxmi Puja
                || (d == 28 && m == October)
                # Bhau bhij
                || (d == 30 && m == October)
                # Gurunanak Jayanti
                || (d == 13 && m == November)
                # Bakri Id
                || (d == 9 && m == December))

                return false
            end
        end

        if (y == 2009)
            if (# Moharram
                (d == 8 && m == January)
                # Mahashivratri
                || (d == 23 && m == February)
                # Id-E-Milad
                || (d == 10 && m == March)
                # Holi
                || (d == 11 && m == March)
                # Ram Navmi
                || (d == 3 && m == April)
                # Mahavir Jayanti
                || (d == 7 && m == April)
                # Maharashtra Day
                || (d == 1 && m == May)
                # Ramzan Id
                || (d == 21 && m == September)
                # Dasara
                || (d == 28 && m == September)
                # Bhau Bhij
                || (d == 19 && m == October)
                # Gurunanak Jayanti
                || (d == 2 && m == November)
                # Moharram (again)
                || (d == 28 && m == December))

                return false
            end
        end

        if (y == 2010)
            if (# New Year's Day
                (d == 1 && m == January)
                # Mahashivratri
                || (d == 12 && m == February)
                # Holi
                || (d == 1 && m == March)
                # Ram Navmi
                || (d == 24 && m == March)
                # Ramzan Id
                || (d == 10 && m == September)
                # Laxmi Puja
                || (d == 5 && m == November)
                # Bakri Id
                || (d == 17 && m == November)
                # Moharram
                || (d == 17 && m == December))

                return false
            end
        end

        if (y == 2011)
            if (# Mahashivratri
                (d == 2 && m == March)
                # Ram Navmi
                || (d == 12 && m == April)
                # Ramzan Id
                || (d == 31 && m == August)
                # Ganesh Chaturthi
                || (d == 1 && m == September)
                # Dasara
                || (d == 6 && m == October)
                # Laxmi Puja
                || (d == 26 && m == October)
                # Diwali - Balipratipada
                || (d == 27 && m == October)
                # Bakri Id
                || (d == 7 && m == November)
                # Gurunanak Jayanti
                || (d == 10 && m == November)
                # Moharram
                || (d == 6 && m == December))

                return false
             end
        end

        if (y == 2012)
            if (# Mahashivratri
                (d == 20 && m == February)
                # Holi
                || (d == 8 && m == March)
                # Mahavir Jayanti
                || (d == 5 && m == April)
                # May Day
                || (d == 1 && m == May)
                # Ramzan Id
                || (d == 20 && m == August)
                # Ganesh Chaturthi
                || (d == 19 && m == September)
                # Dasara
                || (d == 24 && m == October)
                # Diwali - Balipratipada
                || (d == 14 && m == November)
                # Gurunanak Jayanti
                || (d == 28 && m == November))

                return false
             end
        end

        return true

end


text:
"""Lazy Taylor series in 1 variable.
`T` is the type of the coefficients.
`coeffs` are the stored (memoized) coefficients.
`f` is the function that calculates coefficients.

Note: `f` must take *two* variables;
the first is used as an explicit reference to the current object when necessary, e.g. for `exp`.

Each `Taylor1` object keeps a list (`Set`) of other `Taylor1` objects
that depend on it, in `children`, and those on which it depends, in `parents`.
"""

using StaticArrays

const ContainerType = Union{Vector{T}, Dict{Int,T}, MVector{N,T}} where {T,N}

# similar to eltype:
coeff_type(::Type{Vector{T}}) where {T} = T
coeff_type(::Type{Dict{Int,T}}) where {T} = T
coeff_type(::Type{MVector{N,T}}) where {N,T} = T

struct Taylor1{T,F,C<:ContainerType}
    f::F
    coeffs::C  # C is the type of the container
    parents:: Set{Taylor1{T,G,C} where {G}}  # objects this one depends on
    children::Set{Taylor1{T,G,C} where {G}}  # other Taylor1 objects that use this one
end

init(::Type{Vector{T}}) where {T} = T[]
init(::Type{Dict{Int,T}}) where {T} = Dict{Int,T}()
init(::Type{MVector{N,T}}) where {N,T} = zero(MVector{N,T})

Taylor1(f::F, coeffs::C) where {F,C} = Taylor1{coeff_type(C),F,C}(f, coeffs, Set(), Set())

Taylor1(f::F, coeffs::C, parents) where {F,C} = Taylor1{coeff_type(C),F,C}(f, coeffs, parents, Set())

Base.haskey(v::AbstractVector, i::Int) = i in keys(v)   # type piracy

# Taylor1(T, f::F, C) where {F} = Taylor1{T,F,C}(f, init(C))

# Taylor1(f::F, C) where {F} = Taylor1(coeff_type(C), f, init(C))
# function Taylor1(f::F) where F
#     t = Taylor1{Float64,F}(f, Float64[])
#     dummy = t[0]  # compile getindex by calculating first coefficient
#     return t
# end

# function Taylor1(T, f::F) where {F}
#     t = Taylor1(f, T[])
#     #dummy = t[0]  # compile getindex by calculating first coefficient
#     return t
# end

# Taylor1(T, f::Function) = Taylor1(T, f)

# Memoized; use NaN to indicate value not yet calculated
# function getindex(t::Taylor1{T,F,C}, i::Int) where {T, F, C}
#     j = i + 1
#     coeffs = t.coeffs
#
#     @inbounds if j <= length(coeffs)
#         if isnan(coeffs[j])
#             coeffs[j] = (t.f)(t, i)  # pass in the object as the first argument to the function for those functions that are recursive
#         end
#
#         return coeffs[j]
#
#     else  # too short
#         current_length = length(coeffs)
#         resize!(coeffs, j)
#         @inbounds coeffs[current_length+1:end] .= NaN
#         @inbounds coeffs[end] = (t.f)(t, i)
#
#         @inbounds return coeffs[end]
#     end
# end

function Base.setindex!(t::Taylor1{T,F,<:AbstractVector{T}}, val, i::Int) where {T,F}
    j = i + 1
    coeffs = t.coeffs

    current_length = length(coeffs)

    if j > current_length
        resize!(coeffs, j)  # fills with junk
    end
    # coeffs[(current_length+1):(end-1)] .= NaN
    @inbounds coeffs[j] = val
    return val
end

function Base.setindex!(t::Taylor1{T,F,Dict{Int,T}}, val, i::Int) where {T,F}
    @inbounds t.coeffs[i+1] = val
    @inbounds return val
end

function Base.getindex(t::Taylor1{T,F,C}, i::Int) where {T, F, C}
    j = i + 1  # offset so that numbering is 0-based
    coeffs = t.coeffs

    if !haskey(coeffs, j)
        @inbounds t[i] = (t.f)(t, i)
    end

    return @inbounds coeffs[j]
end

# tt is the independent variable; non-memoized (nothing stored in memory)
# tt = Taylor1( i::Int -> (i == 1) * 1.0, false )

# tt = Taylor1(Interval{Float64}, (t, i) -> (i == 1) ? (1..1) : (0..0), true)
# use as constant(3); also non-memoized
constant(c::Real) = Taylor1( i::Int -> (i == 0) ? c : 0.0, false )
constant(T, c::Real) = Taylor1(T, i::Int -> (i == 0) ? c : 0.0, false )

# these are memoized, but should look at performance without memoizing perhaps

# use promotion!

import Base: +, -, *

"Add a child to a parent and its parents recursively"
function add_child!(parent, child)
    push!(parent.children, child)

    for p in parent.parents
        add_child!(p, child)
    end
end

function make_taylor(f, C, parents)
    obj = Taylor1( f, init(C), Set(parents) )

    for p in parents
        add_child!(p, obj)
    end

    return obj
end

+(f::Taylor1{T,F1,C}, g::Taylor1{T,F2,C}) where {T,F1,F2,C} = make_taylor( (t, i) -> f[i] + g[i], C, [f, g] )
-(f::Taylor1{T,F1,C}, g::Taylor1{T,F2,C}) where {T,F1,F2,C} = make_taylor( (t, i) -> f[i] - g[i], C, [f, g] )

-(f::Taylor1{T,F,C}) where {T,F,C} = make_taylor( (t, i) -> -f[i], C, [f] )

+(a::Real, f::Taylor1{T,F,C}) where {T,F,C} = make_taylor( (t, i) -> (i == 0) ? a+f[0] : f[i], C, [f] )
-(a::Real, f::Taylor1) = a + (-f)

+(f::Taylor1, a::Real) = a + f
-(f::Taylor1, a::Real) = (-a) + f

# formulas from Warwick Tucker, *Validated Numerics*

*(f::Taylor1{T,F1,C}, g::Taylor1{T,F2,C}) where {T,F1,F2,C} = make_taylor( (t, i) -> sum(f[k] * g[i-k] for k in 0:i), C, [f, g] )

*(a::Real, f::Taylor1{T,F,C}) where {T,F,C} = make_taylor( (t, i) -> a*f[i], C, [f] )
*(f::Taylor1, a::Real) = a * f

# self is a reference to the object exp(g), that is used recursively
# function exp(g::Taylor1{T}) where {T}
#     function f(self, k)
#         k == 0 && return exp(g[0])
#         # dummy = g[k]  # preallocate g
#         @inbounds return sum(i * g[i] * self[k-i] for i in 1:k) / k
#     end
#
#     return Taylor1(T, f)
#
# end

"""
Evaluate using Horner rule
"""
function (f::Taylor1)(x)
    total = f.coeffs[end]

    for i in length(f.coeffs)-1 : -1 : 1
        @inbounds total = x * total + f.coeffs[i]
    end

    return total
end

Base.literal_pow(::typeof(^), t::Taylor1, n::Integer) = t^n

import Base: ^
^(t::Taylor1, n::Integer) = power_by_squaring(t, n)  # uses modified version in powers.jl

"""
Reset a Taylor series by emptying its own coefficients and the
coefficients of all other Taylor objects that depend on it.
"""
function reset!(t::Taylor1)
    empty!(t.coeffs)

    for c in t.children
        empty!(c.coeffs)
    end
end

max_stars_repo_name: dpsanders/LazyTaylorSeries

max_issues_repo_name: dpsanders/LazyTaylorSeries

max_forks_repo_name: dpsanders/LazyTaylorSeries


text:
"""Reference for parameters update when present"""
struct UpdateRef{T}
    access_ref::Symbol
end

const DRDA = Dict{UpdateRef, JuMP.Containers.DenseAxisArray}


text:
using ADCME
using PyPlot
encoding_dim = 32
x = placeholder(Float64, shape=[nothing, 784])
variable_scope("encoder", reuse=AUTO_REUSE) do
    net = x
    net = dense(net, 200, activation="relu")
    net = dense(net, 200, activation="relu")
    global encoded = dense(net, encoding_dim, activation="relu")
end
variable_scope("decoder", reuse=AUTO_REUSE) do
    net = dense(encoded, 200, activation="relu")
    net = dense(net, 200, activation="relu")
    global decoded = dense(net, 784, activation="sigmoid")
end

x_train, _, x_test, = mnist(Float64)
x_train /= 255
x_test /= 255
x_train = reshape(x_train, size(x_train,1), length(x_train[1,:,:]))
x_test = reshape(x_test, size(x_test,1), length(x_test[1,:,:]))
loss = sum(tf.keras.backend.binary_crossentropy(x, decoded))
variable_scope("nn2", reuse=AUTO_REUSE) do
    opt = AdamOptimizer()
    global train_op = minimize(opt, loss)
end
sess = Session()
init(sess)
batchsize=256
for i = 1:10000
    II = (1+(i-1)*batchsize):i*batchsize
    II = mod.(II, size(x_train,1)).+1
    data = x_train[II, :]
    _, los = run(sess, [train_op, loss], feed_dict=Dict(x=>data))
    if mod(i,100)==0
        println("#iter=$i, loss=$los")
    end
end

p = rand(1:size(x_train,1),6)
V = run(sess, decoded, feed_dict=Dict(x=>x_train[p,:]))
# V = Array{Int64}(V .> 0.5)
figure()
for i = 1:5
    subplot(2,5,i)
    imshow(reshape(V[i,:],28,28))
    subplot(2,5,i+5)
    imshow(reshape(x_train[p[i],:],28,28))
end
savefig("tmp.png")

ec = run(sess, encoded, feed_dict=Dict(x=>x_train[p,:]))
ec2 = (ec[1:3,:]+ec[4:6,:])/2

y = placeholder(Float64, shape=[nothing, 32])
variable_scope("decoder", reuse=AUTO_REUSE) do
    net = dense(y, 200, activation="relu")
    net = dense(net, 200, activation="relu")
    global imgs = dense(net, 784, activation="sigmoid")
end
IMGS = run(sess, imgs, feed_dict=Dict(y=>ec2))
# IMGS = IMGS .> 0.5
figure()
for i = 1:3
    subplot(1,3,i)
    imshow(reshape(IMGS[i,:],28,28))
end
savefig("tmp2.png")