text:
using ApproximateDPs
const ADP = ApproximateDPs
using FlightSims
const FS = FlightSims
using Transducers
using Plots
using Random, LinearAlgebra, ComponentArrays
using DynamicPolynomials, UnPack
using DifferentialEquations, DataFrames
using FileIO

function initialise(; exact_integration=false)
    # setting
    env = TwoDimensionalNonlinearPolynomialSystem()
    u_explorer = (x, p, t) -> sin(5*t)
    __x = State(env)()
    __t = 0.0
    n = length(__x)
    __u = u_explorer(__x, (), __t)
    m = length(__u)
    u_norm_max = 5.0  # input constraint for numerical optimisation; TODO: can be deprecated
    adp = CTValueIterationADP(n, m, FSimZoo.RunningCost(env), u_norm_max, exact_integration)
    env, u_explorer, adp
end

function explore(dir_log, env, adp, u_explorer;
        file_name="exploration.jld2",
        Δt=0.003, tf=1.8,  # the ettings in the paper and implemented code by T. Bian
    )
    file_path = joinpath(dir_log, file_name)
    prob, sol = nothing, nothing
    df = nothing
    x0 = State(env)(-2.9, -2.9)  # the setting in the paper
    if adp.exact_integration
        X0 = ComponentArray(x=x0, ∫Ψ=zeros(1, adp.N_Ψ))  # append system, 1 x N
        appended_dynamics! = function (env::TwoDimensionalNonlinearPolynomialSystem)
            @Loggable function dynamics!(dX, X, p, t; u)
                @unpack x, ∫Ψ = X
                @log state = x
                @log ∫Ψ = ∫Ψ
                @log input = u
                Dynamics!(env)(dX.x, X.x, (), t; u=u)
                dX.∫Ψ = ADP.Ψ(adp)(dX.x, u)
            end
        end
        prob, df = sim(
                       X0,
                       apply_inputs(appended_dynamics!(env); u=u_explorer);
                       tf=tf,
                       savestep=Δt
                      )
    else
        prob, df = sim(
                       x0,
                       apply_inputs(Dynamics!(env); u=u_explorer);
                       tf=tf,
                       savestep=Δt
                      )
    end
    FileIO.save(file_path, Dict("env" => env, "prob" => prob, "sol" => df.sol))
    if adp.exact_integration
        ∫Ψs = df.sol |> Map(datum -> datum.∫Ψ) |> collect
    end
    ts = df.time
    xs = df.sol |> Map(datum -> datum.state) |> collect
    us = df.sol |> Map(datum -> datum.input) |> collect
    df
end

function train!(adp; max_iter=100, w_tol=0.01)
    @show adp.V̂.param
    i = 0
    w_prev = deepcopy(adp.V̂.param)
    stop_conds = function(i, w_diff_norm)
        stop_conds_dict = Dict(
                              :max_iter => i >= max_iter,
                              :w_tol => w_diff_norm < w_tol,
                             )
    end
    display_res = function (adp::CTValueIterationADP)
        @unpack n = adp
        @polyvar x[1:n]
        @show adp.V̂(x)
    end
    while true
        i += 1
        @show i
        lr = 1 / (i+10)
        ADP.update!(adp, lr)
        @show adp.V̂.param
        display_res(adp)
        stop_conds_dict = stop_conds(i, norm(adp.V̂.param-w_prev))
        if any(values(stop_conds_dict))
            @show stop_conds_dict
            break
        else
            w_prev = deepcopy(adp.V̂.param)
        end
    end
end

function demonstrate(env, adp; Δt=0.01, tf=10.0)
    x0 = State(env)(-2.9, -2.9)
    prob, df = sim(
                   x0,
                   apply_inputs(Dynamics!(env); u=ADP.approximate_optimal_input(adp));
                   tf=tf,
                   savestep=Δt
                  )
    xs = df.sol |> Map(datum -> datum.state) |> collect
    plot(df.time, hcat(xs...)')
end

"""
Main codes for demonstration of continuous-time value-iteration adaptive dynamic programming (CT-VI-ADP) [1].
# References
[1] T. Bian and Z.-P. Jiang, “Value Iteration, Adaptive Dynamic Programming, and Optimal Control of Nonlinear Systems,” in 2016 IEEE 55th Conference on Decision and Control (CDC), Las Vegas, NV, USA, Dec. 2016, pp. 3375–3380. doi: 10.1109/CDC.2016.7798777.
[2] Forked repo of the code implementation by T. Bian, https://github.com/JinraeKim/nonlinear-VI/blob/main/nonlinear_sys.R
# Settings
`exact_integration`:
if `true`, `∫Ψ` is calculated with the dynamical system simultaneously (i.e., regarded as the true solution).
if `false`, it is calculated by numerical integration (default setting may be trapezoidal method.)
`Δt`: data sampling period
"""
function main(; seed=1, dir_log="data/main/continuous_time_vi_adp")
    Random.seed!(seed)
    env, u_explorer, adp = initialise(; exact_integration=true)
    df = explore(dir_log, env, adp, u_explorer;
                 Δt=0.003)
    # @show df.sol
    ADP.set_data!(adp, df)
    train!(adp;
           max_iter=100, w_tol=1e-2)
    demonstrate(env, adp)
end


text:
import java.io.IOException;
import java.io.InputStreamReader;
import java.io.PrintStream;
import java.util.ArrayList;

public class Screen {
    static final String ANSI_CLS = "\u001b[2J";
    static final String ANSI_HOME = "\u001b[H";

    private final SimpleIRCClient parent;

    // ANNOTATION: field erasure type
    private ArrayList{L /parent.clearHistory T} buffer;
    private final PrintStream out;
    private final InputStreamReader isr;
    private final ArrayList inputBuffer;
    private String incomplete;

    public Screen(SimpleIRCClient parent) {
    	this.parent = parent;
    	this.buffer = new ArrayList();
    	this.inputBuffer = new ArrayList();
	this.out = System.out;
	this.isr = new InputStreamReader(System.in);
    	this.incomplete = "";
    }

    public String getInput() {
	this.bufferInput();
	if (!inputBuffer.isEmpty()) {
	    String input = (String) inputBuffer.get(0);
	    inputBuffer.remove(0);
	    return input;
    	} else {
	    return null;
    	}
    }

    public void bufferInput() {
	try {
	    int c = -1;
	    while (isr.ready() && (c = isr.read()) != -1) {
		if (c == '\n') {
		    inputBuffer.add(incomplete);
		    incomplete = "";
		} else {
		    incomplete += (char) c;
		}
	    }
	} catch (Throwable e) {}
    }

    public void flush() {
    	this.bufferInput();
    	this.clearScreen();
	int drawn = 0;
	try {
	    for (int i = 0; i < buffer.size(); i++) {
		// Declassify output to the screen
		// ANNOTATION: declassification
    		out.println(({L})buffer.get(i));
		drawn++;
	    }
	} catch (Throwable e) {}
	for (int i = drawn; drawn < 20; drawn++) {
	    out.println();
    	}
	out.print("> " + incomplete);
    	out.flush();
    }

    public void writeLine(String line) {
	try {
	    buffer.add(line);
	    while (buffer.size() > 20) {
		buffer.remove(0);
	    }
	} catch (Throwable e) {}
    }

    public void clearScreen() {
    	out.print(ANSI_CLS + ANSI_HOME);
    	out.flush();
    }
}

ext: jl

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
"""
    ψ_Stommel(x, y; τ₀, β, r)

Stommel gyre streamfunction taken from equation (19.39) of Vallis (2nd edition) wher L is the
length of each side of the square domain, τ₀ is the magnitude of the wind stress, β is the
Rossby parameter, and r is the linear drag, or Rayleigh friction, acting on the vertically
integrated velocity.
"""
@inline ψ_Stommel(x, y; L, τ₀, β, r) = τ₀*π/β * (1 - x/L - exp(-x/L/ϵ(r, β, L))) * sin(π*y/L)

@inline u_Stommel(x, y; L, τ₀, β, r) = τ₀*π^2/β/L * (1 - x/L - exp(-x/L/ϵ(r, β, L))) * cos(π*y/L)
@inline v_Stommel(x, y; L, τ₀, β, r) = τ₀*π/β/L * (1 - 1/ϵ(r, β, L) * exp(-x/L/ϵ(r, β, L))) * sin(π*y/L)

@inline u_Stommel_max(; L, τ₀, β) = τ₀*π^2/β/L
@inline v_Stommel_max(; L, τ₀, β) = τ₀*π/β/L

@inline ϵ(r, β, L) = r / (β*L)


text:
module BinningAnalysisPlots

using BinningAnalysis, RecipesBase, StatsBase

include("recipes.jl")

end # module


text:
include("../src/ImportFolders.jl")
using Manager
using Helper

orgs = Creature.Organism[]
len  = 500

for i=1:len
  push!(orgs, Creature.create())
end

function run()
  local i::Int

  while true
    for i=1:len
      Helper.save(orgs[i], "code-before.jevo")
      Mutator.mutate(orgs[i])
      orgs[i].codeFn = eval(orgs[i].code)
      try
        Helper.save(orgs[i], "code-after.jevo")
        orgs[i].codeFn(orgs[i])
      catch e
      end
    end
  end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
#
# Result
# can have either of the following:
# - recordset
# - number of rows affected
#
# pending results can be
# - checked for readiness
# - cancelled
#
# results must be closed

# number of records to fetch with every fetchnext
const DEFAULT_FETCH_SIZE = 1024

type PendingResult
    session::HiveSession
    handle::TOperationHandle
    status::Nullable{TGetOperationStatusResp}
end

typealias RowCount Float64

type ResultSet
    session::HiveSession
    handle::TOperationHandle
    schema::Nullable{TTableSchema}
    rowset::Nullable{TRowSet}
    position::Int
    fetchsize::Int
    eof::Bool

    function ResultSet(session::HiveSession, handle::TOperationHandle)
        new(session, handle, Nullable{TOperationHandle}(), Nullable{TRowSet}(), 0, DEFAULT_FETCH_SIZE, false)
    end
end

typealias Result Union{ResultSet, RowCount, PendingResult}

result(pending::PendingResult) = isready(pending) ? result(pending.session, true, pending.handle) : pending
result(rs::ResultSet) = rs
function result(session::HiveSession, ready::Bool, handle::TOperationHandle)
    ready || (return PendingResult(session, handle, Nullable{TGetOperationStatusReq}()))

    # hasResultSet == true => a result set (possibly empty) can be fetched, modifiedRowCount is not set.
    # hasResultSet == false => modifiedRowCount is set
    handle.hasResultSet && (return ResultSet(session, handle))

    # modifiedRowCount >= 0 => number of rows affected is known
    # modifiedRowCount < 0 => operation is capable fo modifying rows, but the count is unknown. e.g. LOAD DATA
    return handle.modifiedRowCount
end

eof(r::ResultSet) = r.eof

close(pending::PendingResult) = close(pending.session, pending.handle)
close(resultset::ResultSet) = close(resultset.session, resultset.handle)
function close(session::HiveSession, handle::TOperationHandle)
    conn = session.conn
    request = thriftbuild(TCloseOperationReq, Dict(:operationHandle => handle))
    response = CloseOperation(conn.client, request)
    check_status(response.status)
end

cancel(pending::PendingResult) = cancel(pending.session, pending.handle)
cancel(resultset::ResultSet) = cancel(resultset.session, resultset.handle)
function cancel(session::HiveSession, handle::TOperationHandle)
    conn = session.conn
    request = thriftbuild(TCancelOperationReq, Dict(:operationHandle => handle))
    response = CancelOperation(conn.client, request)
    check_status(response.status)
end

function status(pending::PendingResult)
    if !isnull(pending.status)
        response = get(pending.status)
    else
        conn = pending.session.conn
        request = thriftbuild(TGetOperationStatusReq, Dict(:operationHandle => pending.handle))
        response = GetOperationStatus(conn.client, request)
        check_status(response.status)
        if response.operationState in READY_STATUS
            pending.status = Nullable(response)
        end
    end
    response
end

const READY_STATUS = (TOperationState.FINISHED_STATE, TOperationState.CANCELED_STATE, TOperationState.CLOSED_STATE, TOperationState.ERROR_STATE)
const RESULT_STATUS = (TOperationState.FINISHED_STATE)
isready(::ResultSet) = true
isready(::RowCount) = true
isready(pending::PendingResult) = status(pending).operationState in READY_STATUS
hasresult(pending::PendingResult) = status(pending).operationState in RESULT_STATUS
function sqlerror(pending::PendingResult)
    isready(pending) || (return utf8(""))
    hasresult(pending) && (return utf8(""))
    response = get(pending.status)
    sqlstate = isfilled(response, :sqlState) ? getfield(response, :sqlState) : utf8("")
    errormessage = isfilled(response, :errorMessage) ? getfield(response, :errorMessage) : utf8("")
    errorcode = isfilled(response, :errorCode) ? getfield(response, :errorCode) : Int32(0)
    utf8("Error. $errormessage ($sqlstate, $errorcode)")
end

const SCHEMA_CACHE = Dict{Union{AbstractString,Symbol}, TTableSchema}()

function schema(rs::ResultSet; cached::Union{AbstractString,Symbol,Void}=nothing)
    if !isnull(rs.schema)
        sch = get(rs.schema)
    else
        if (cached !== nothing) && haskey(SCHEMA_CACHE, cached)
            sch = SCHEMA_CACHE[cached]
        else
            conn = rs.session.conn
            request = thriftbuild(TGetResultSetMetadataReq, Dict(:operationHandle => rs.handle))
            response = GetResultSetMetadata(conn.client, request)
            check_status(response.status)
            sch = response.schema
            (cached !== nothing) && (SCHEMA_CACHE[cached] = sch)
        end
        rs.schema = Nullable(sch)
    end
    sch
end

# Fetch rows from the server corresponding to a particular OperationHandle.
fetchsize(rs::ResultSet) = rs.fetchsize
fetchsize!(rs::ResultSet, sz::Integer) = (rs.fetchsize = sz)

function fetchfirst(rs::ResultSet, nrows::Integer=rs.fetchsize)
    rs.position = 0
    fetchnext(rs, nrows)
end

function fetchnext(rs::ResultSet, nrows::Integer=rs.fetchsize)
    rs.eof && return nothing
    conn = rs.session.conn
    orientation = rs.position > 0 ? TFetchOrientation.FETCH_NEXT : TFetchOrientation.FETCH_FIRST
    request = thriftbuild(TFetchResultsReq, Dict(
        :operationHandle => rs.handle,
        :orientation => orientation,
        :maxRows => nrows))
    response = FetchResults(conn.client, request)
    check_status(response.status)

    rs.eof = !response.hasMoreRows
    rowset = response.results
    rs.rowset = Nullable(rowset)
    rs.position = rowset.startRowOffset + length(rowset.rows)
    rs
end

##
# Iterators
# - records, record iterator: fetches one record at a time as a tuple
# - dataframes, dataframe iterator: fetches one batch of records at a time, returns a dataframe
# - dataframe(dataframe iterator): fetches and returns the first batch

type DataFrameIterator
    rs::ResultSet

    function DataFrameIterator(rs::ResultSet, fetchsz::Integer=DEFAULT_FETCH_SIZE)
        fetchsize!(rs, fetchsz)
        new(rs)
    end
end

dataframes(rs::ResultSet, fetchsz::Integer=DEFAULT_FETCH_SIZE) = DataFrameIterator(rs, fetchsz)
dataframe(rs::ResultSet, fetchsz::Integer=DEFAULT_FETCH_SIZE) = dataframe(dataframes(rs, fetchsz))
function dataframe(iter::DataFrameIterator)
    df = reduce(vcat, iter)
    close(iter.rs)
    df
end

start(iter::DataFrameIterator) = iter.rs.eof
done(iter::DataFrameIterator, state) = state
function next(iter::DataFrameIterator, state)
    rs = iter.rs
    fetchnext(rs)

    rowset = get(rs.rowset)
    sch = schema(rs)
    ncols = length(sch.columns)
    colvecs = DataArray[]

    if isfilled(rowset, :columns)
        @logmsg("reading columns")
        # ref: https://issues.apache.org/jira/browse/HIVE-3746
        for col in rowset.columns
            push!(colvecs, julia_type(col))
        end
    else
        for T in coltypes(sch)
            push!(colvecs, DataArray(T[]))
        end
        @logmsg("reading rows for coltypes: $(coltypes(sch))")
        for row in rowset.rows
            colidx = 1
            for col in row.colVals
                push!(colvecs[colidx], julia_type(col))
                colidx += 1
            end
        end
    end

    df = DataFrame()
    for colidx in 1:ncols
        colname = sch.columns[colidx].columnName
        @logmsg("$colname: $(colvecs[colidx])")
        df[Symbol(colname)] = colvecs[colidx]
    end
    df, rs.eof
end

type RecordIterator
    dfiter::DataFrameIterator
    dfpos::Int
    df::Nullable{DataFrame}

    RecordIterator(rs::ResultSet) = new(DataFrameIterator(rs), 1, Nullable{DataFrame}())
end

records(rs::ResultSet) = RecordIterator(rs)

start(iter::RecordIterator) = start(iter.dfiter)
function done(iter::RecordIterator, state)
    done(iter.dfiter, state) || (return false)
    isnull(iter.df) && return false
    df = get(iter.df)
    iter.dfpos > nrow(df)
end
function next(iter::RecordIterator, state)
    if isnull(iter.df) || (iter.dfpos > nrow(get(iter.df)))
        df, state = next(iter.dfiter, state)
        iter.df = Nullable(df)
        iter.dfpos = 1
    else
        df = get(iter.df)
    end
    recs = tuple(convert(Array, df[iter.dfpos,:])...)
    iter.dfpos += 1
    recs, state
end


text:
function TIER_td_loss(;loss=Flux.mse, name=:Qavg, s_key=:s, a_key=:a, weight=nothing)
    (π, 𝒫, 𝒟, y; info=Dict(), z) -> begin
        Q = value(π, vcat(𝒟[s_key], z), 𝒟[a_key])

        # Store useful information
        Zygote.ignore() do
            info[name] = mean(Q)
        end

        loss(Q, y, agg = isnothing(weight) ? mean : weighted_mean(𝒟[weight]))
    end
end

function TIER_double_Q_loss(;name1=:Q1avg, name2=:Q2avg, kwargs...)
    l1 = TIER_td_loss(;name=name1, kwargs...)
    l2 = TIER_td_loss(;name=name2, kwargs...)

    (π, 𝒫, 𝒟, y; info=Dict(), z=𝒟[:z]) -> begin
        .5f0*(l1(critic(π).N1, 𝒫, 𝒟, y, info=info, z=z) + l2(critic(π).N2, 𝒫, 𝒟, y, info=info,z=z))
    end
end

TIER_TD3_actor_loss(π, 𝒫, 𝒟; info = Dict()) = -mean(value(critic(π).N1, vcat(𝒟[:s], 𝒟[:z]), action(actor(π), vcat(𝒟[:s], 𝒟[:z]))))

function TIER_TD3_target(π, 𝒫, 𝒟, γ::Float32; i, z=𝒟[:z])
    ap, _ = exploration(𝒫[:π_smooth], vcat(𝒟[:sp], z), π_on=actor(π), i=i)
    y = 𝒟[:r] .+ γ .* (1.f0 .- 𝒟[:done]) .* min.(value(critic(π), vcat(𝒟[:sp], z), ap)...)
end

TIER_action_regularization(π, 𝒟) = Flux.mse(action(actor(π), vcat(𝒟[:s], 𝒟[:z])), 𝒟[:a])
TIER_action_value_regularization(π, 𝒟) = begin
    v = value(critic(π), vcat(𝒟[:s], 𝒟[:z]), 𝒟[:a])
    v isa Tuple && (v = v[1])
    Flux.mse(v, 𝒟[:value])
end

function TIER(;π,
               observation_model,
               S,
               ΔN,
               𝒫=(;),
               A=action_space(π),
               buffer_size=1000,
               latent_dim,
               N_experience_replay,
               N_experience_obs,
               replay_store_weight = (D)->1f0,
               bayesian_inference,
               solver=TD3,
               required_columns=Symbol[],
               c_opt=(;),
               a_opt=(;),
               obs_opt=(;),
               a_loss=TIER_TD3_actor_loss,
               c_loss=TIER_double_Q_loss(),
               target_fn=TIER_TD3_target,
               zprior=MvNormal(zeros(latent_dim), I),
               kwargs...)

    required_columns = unique([:episode_end, :weight, required_columns...])
    buffer = ExperienceBuffer(S, A, buffer_size, required_columns)
    buffer.data[:z] = zeros(Float32, latent_dim, capacity(buffer))

    # This experience buffer is for recalling behavior
    buffer_er = ExperienceBuffer(S, A, N_experience_replay, [required_columns..., :value])
    buffer_er.data[:z] = zeros(Float32, latent_dim, capacity(buffer_er))

    # this experience buffer is for learning the latent embedding
    buffer_obs = ExperienceBuffer(S, A, N_experience_obs, [required_columns..., :value])
    buffer_obs.data[:z] = zeros(Float32, latent_dim, capacity(buffer_obs))

    # Buffer used to train for latent
    obs_opt = TrainingParams(;obs_opt...)
    𝒟obs = buffer_like(buffer_obs, capacity=obs_opt.batch_size)

    𝒫 = (;buffer_er, buffer_obs, obs_opt, z_dist=Any[zprior], zs=Any[zprior], observation_model, 𝒫...)

    # Define the solver
    𝒮 = solver(;π=π,
              S=S,
              𝒫=𝒫,
              buffer=buffer,
              c_opt=(;regularizer=BatchRegularizer(buffers=[buffer_er], batch_size=128, λ=0.5f0, loss=TIER_action_value_regularization, c_opt...)),
              a_opt=(;regularizer=BatchRegularizer(buffers=[buffer_er], batch_size=128, λ=0.5f0, loss=TIER_action_regularization, a_opt...)),
              extra_buffers = [buffer_er],
              buffer_fractions = [0.5, 0.5],
              a_loss=a_loss,
              c_loss=c_loss,
              target_fn=target_fn,
              kwargs...
              )

    function TIER_cb(D; 𝒮, info=Dict())
        # update the z distribution
        𝒮.𝒫.z_dist[1] = bayesian_inference(observation_model, D, 𝒮.𝒫.z_dist[1], info=info)

        # Set the agent's best estimate and record
        zbest = mean(𝒮.𝒫.z_dist[1])
        𝒮.agent.π.z = zbest
        push!(𝒮.𝒫[:zs], deepcopy(𝒮.𝒫.z_dist[1]))

        # Fill the buffer with latent estimate, value and computed weight
        D[:z] = repeat(zbest, 1, length(D[:r]))
        D[:value] = mean(value(critic(𝒮.agent.π), vcat(D[:s], D[:z]), D[:a]))
        D[:weight] .= replay_store_weight(D)

        # Add this buffer to our experience replay and observation buffers
        push_reservoir!(buffer_er, D, weighted=true)
        push_reservoir!(buffer_obs, D)

        # Train the obs model
        for j=1:obs_opt.epochs
                rand!(𝒟obs, buffer_obs, buffer, fracs=[0.5, 0.5])
                train!(Flux.params(observation_model), (;kwargs...) -> obs_opt.loss(observation_model, 𝒟obs; kwargs...), obs_opt, info=info)
        end
    end

    𝒮.post_sample_callback = TIER_cb
    𝒮
end


text:
using Dash
using DashHtmlComponents
using DashCoreComponents

app = dash(external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"])

app.layout = html_div() do
        html_h1("Hello Dash"),
        html_div("Dash.jl: Julia interface for Dash"),
        dcc_graph(
            id = "example-graph",
            figure = (
                data = [
                    (x = [1, 2, 3], y = [100, 1, 2], type = "bar", name = "SF"),
                    (x = [1, 2, 3], y = [2, 4, 5], type = "bar", name = "Montréal"),
                ],
                layout = (title = "Dash Data Visualization 2",)
            )
        )
    end

run_server(app, "0.0.0.0", 8080)

max_stars_repo_name: krystophny/profit

max_issues_repo_name: krystophny/profit

max_forks_repo_name: krystophny/profit


text:
export
    Serial

"""
    struct Serial{I,O} <: SFunc{I,O}

A sequence of sfuncs in series.

Although this could be implemented as a special case of NetworkSFunc,
the serial composition allows an easier and more efficient implementation of operations.
All but the first sfunc in the sequence will have a single input; the output of each
sfunc feeds into the input of the next sfunc.

To work properly, most of the operations on `Serial` need the support of the intermediate
sfuncs, given an input range. Rather than compute this each time, and to avoid having
the non-support operations take a size argument, support is memoized, and must be called
before other operations like logcpdf are called. The `support_memo` is a dictionary whose
keys are tuples of parent ranges and whose values are the support computed for those
parent ranges, along with the target size for which they were computed.
Storing the target size enables refinement algorithms that increase the size and improve
the support.

# Additional supported operators
- `support`
- `support_quality`
- `sample`
- `cpdf`
- `bounded_probs`
- `make_factors`
- `compute_pi`
- `send_lambda`

# Type parameters
- I the input type of the first sfunc
- O the output type of the last sfunc
"""
struct Serial{I,O} <: SFunc{I,O}
    components :: NTuple{N,SFunc} where N
    support_memo :: Dict
    Serial(I,O,sfuncs) = new{I,O}(sfuncs, Dict())
end

@impl begin
    struct SerialSupport end
    function support(sf::Serial{I,O}, parranges, size, curr) where {I,O}
        if parranges in keys(sf.support_memo)
            (sup,sz) = sf.support_memo[parranges]
            if sz >= size
                return last(sup)
            end
        end

        # Adequate support has not been found - compute it now
        prs = parranges
        compsups = Vector[]
        for component in sf.components
            sup = support(component, prs, size, output_type(component)[])
            push!(compsups, sup)
            prs = (sup,)
        end
        for c in curr
            if !(c in sup)
                push!(sup, c)
            end
        end
        sf.support_memo[parranges] = (tuple(compsups...), size)
        return sup
    end
end

@impl begin
    struct SerialSupportQuality end

    function support_quality(sf::Serial{I,O}, parranges) where {I,O}
        if !(parranges in keys(sf.support_memo))
            error("Support must be called before support_quality for Serial")
        end
        compsups = sf.support_memo[parranges][1]
        rank = support_quality_rank(:CompleteSupport)
        prs = parranges
        for (comp,sup) in zip(sf.components, compsups)
            q = support_quality_rank(support_quality(comp, prs))
            rank = min(rank, q)
            prs = (sup,)
        end
        return support_quality_from_rank(rank)
    end
end

@impl begin
    struct SerialSample end
    function sample(sf::Serial{I,O}, i::I) where {I,O}
        x = i
        for component in sf.components
            x = (sample(component, x),)
        end
        return x[1]
    end
end

function _checksup(sf, i)
    for k in keys(sf.support_memo)
        if all([i[j] in k[j] for j in 1:length(k)])
            return sf.support_memo[k][1]
        end
    end
    return nothing
end

@impl begin
    struct SerialCpdf end
    function cpdf(sf::Serial{I,O}, i::I, o::O) where {I,O}
        compsups = _checksup(sf, i)
        if isnothing(compsups)
            error("No support found for parent values in logcpdf for Serial")
        end
        us = [i]
        pis = [1.0]
        n = length(sf.components)
        for j in 1:n-1
            comp = sf.components[j]
            sup = compsups[j]
            newpis = zeros(Float64, length(sup))
            for (u,pi) in zip(us,pis)
                for l in 1:length(sup)
                    x = sup[l]
                    newpis[l] += cpdf(comp, u, x) * pi
                end
            end
            us = [(x,) for x in sup]
            pis = newpis
        end
        result = 0.0
        finalcomp = sf.components[n]
        for (u, pi) in zip(us,pis)
            result += cpdf(finalcomp, u, o) * pi
        end
        return result
    end
end

@impl begin
    struct SerialBoundedProbs end

    function bounded_probs(sf::Serial{I,O}, range::Vector{<:O},
            parranges) where {I,O}
        probs = Float64[]
        combos = cartesian_product(parranges)
        for i in combos
            for o in range
                push!(probs, cpdf(sf, tuple(i...), o))
            end
        end
        return (probs, probs)
    end
end

@impl begin
    struct SerialMakeFactors end

    function make_factors(sf::Serial{I,O}, range::Vector{<:O}, parranges,
            id, parids) where {I,O}
        if !(parranges in keys(sf.support_memo))
            error("Support must be called before make_factors for Serial")
        end
        compsups = sf.support_memo[parranges][1]
        inids = parids
        outid = nextkey()
        prs = parranges
        lowers = Factor[]
        uppers = Factor[]
        for i in 1:length(sf.components)
            outid = i == length(sf.components) ? id : nextkey()
            (ls, us) = make_factors(sf.components[i], compsups[i], prs, outid, inids)
            inid = (outid,)
            prs = (compsups[i],)
            append!(lowers, ls)
            append!(uppers, us)
        end
        return(lowers, uppers)
    end
end

@impl begin
    struct SerialComputePi end

    function compute_pi(sf::Serial{I,O}, range::Vector{<:O}, parranges,
            incoming_pis) where {I,O}
        ps = zeros(Float64, length(range))
        combos = cartesian_product(parranges)
        m = length(parranges)
        for combo in combos
            parpis = [cpdf(incoming_pis[i], (), combo[i]) for i in 1:m]
            ipi = reduce(*, parpis)
            for (j,o) in enumerate(range)
                ps[j] += cpdf(sf, tuple(combo...), o) * ipi
            end
        end
        return Cat(range, ps)
    end
end

function _incpis(components, supports, parranges, incoming_pis)
    incpis = Tuple[incoming_pis]
    prs = parranges
    for i in 1:length(components)-1
        comp = components[i]
        rng = supports[i]
        pi = compute_pi(comp, rng, prs, incpis[i])
        push!(incpis, (pi,))
        prs = (rng,)
    end
    return incpis
end

function _lambdas(components, supports, lambda, incpis)
    lambdas = Score[lambda]
    lam = lambda
    for i = length(components):-1:2
        comp = components[i]
        rng = supports[i]
        prs = (supports[i-1],)
        ipis = incpis[i]
        lam = send_lambda(comp, lam, rng, prs, ipis, 1)
        pushfirst!(lambdas, lam)
    end
    return lambdas
end

@impl begin
    struct SerialSendLambda end

    function send_lambda(sf::Serial{I,O}, lambda, range, parranges, incoming_pis,
            parent_idx) where {I,O}
        if !(parranges in keys(sf.support_memo))
            error("Support must be called before send_lambda for Serial")
        end
        compsups = sf.support_memo[parranges][1]
        incpis = _incpis(sf.components, compsups, parranges, incoming_pis)
        lambdas = _lambdas(sf.components, compsups, lambda, incpis)
        return send_lambda(sf.components[1], lambdas[1], compsups[1],
            parranges, incoming_pis, parent_idx)
    end
end
#=
@impl begin
    struct SerialInitialStats end
    initial_stats(sf::Serial) = map(initial_stats, sf.components)
end

@impl begin
    struct SerialAccumulateStats end
    accumulate_stats(sf::Serial, existing_stats, new_stats) =
        [accumulate_stats(sf.components[i], existing_stats[i], new_stats[i])
            for i in 1:length(sf.components)]
end

@impl begin
    struct SerialExpectedStats end
    function expected_stats(sf::Serial, range, parranges, incoming_pis, child_lambda)
        if !(parranges in keys(sf.support_memo))
            error("Support must be called before expected_stats for Serial")
        end
        compsups = sf.support_memo[parranges][1]
        incpis = _incpis(sf.components, compsups, parranges, incoming_pis)
        lambdas = _lambdas(sf.components, compsups, child_lambda, incpis)
        stats = Any[]
        prs = parranges
        for i in 1:length(sf.components)
            comp = sf.components[i]
            sup = compsups[i]
            ipi = incpis[i]
            lam = lambdas[i]
            push!(stats, expected_stats(comp, sup, prs, ipi, lam))
            prs = (sup,)
        end
        return tuple(stats...)
    end
end

@impl begin
    struct SerialMaximizeStats end
    function maximize_stats(sf::Serial, stats)
        ps = [maximize_stats(sf.components[i], stats[i]) for i in 1:length(sf.components)]
        return tuple(ps...)
    end
end
=#


text:
# We keep generating taskgraphs super simple - just generate a random graph.
function make_taskgraph(ntasks, nedges)
    tasks = [TaskgraphNode(string(i)) for i in 1:ntasks]

    # Randomly create edges.
    #
    # Since the Mapper supports multiple edges from a source to destination, we
    # don't have to worry about checking if an edge already exists.
    edges = map(1:nedges) do _
        source = rand(1:ntasks)
        dest = rand(1:ntasks)

        return TaskgraphEdge(string(source), string(dest))
    end

    return Taskgraph("taskgraph", tasks, edges)
end


text:
using MultivariatePolynomials
using JuMP

using RecipesBase

export Ellipsoid

struct Ellipsoid{T}
    Q::Matrix{T}
    c::Vector{T}
end

@recipe function f(ell::Ellipsoid)
    @assert LinearAlgebra.checksquare(ell.Q) == 2
    αs = range(0, stop=2π, length=1024)
    ps = [[cos(α), sin(α)] for α in αs]
    r = [sqrt(dot(p, ell.Q * p)) for p in ps]
    seriestype --> :shape
    legend --> false
    ell.c[1] .+ cos.(αs) ./ r, ell.c[2] .+ sin.(αs) ./ r
end

struct LiftedEllipsoid{T}
    P::Matrix{T}
end

function LiftedEllipsoid(ell::Ellipsoid)
    md = ell.Q*ell.c
    δ = ell.c'*md-1
    d = -md
    D = ell.Q
    P = [δ d'
         d D]
    LiftedEllipsoid(P)
end

Base.convert(::Type{Ellipsoid{T}}, ell::LiftedEllipsoid) where T = convert(Ellipsoid{T}, Ellipsoid(ell))
function Bbβλ(P)
    n = LinearAlgebra.checksquare(P) - 1
    ix = 1 .+ (1:n)
    β = P[1, 1]
    b = P[1, ix]
    B = P[ix, ix]
    λ = dot(b, B \ b) - β
    @assert λ >= 0
    B, b, β, λ
end
function Ellipsoid(ell::LiftedEllipsoid)
    # P is
    # λ * [c'Qc-1  -c'Q
    #         -Qc   Q]
    # Let P be [β b'; b B]
    # We have
    # β = λ c'Qc - λ
    # b = -λ Qc <=> Q^{-1/2}b = -λ Q^{1/2} c
    # hence
    # λ c'Qc = β + λ
    # λ^2 c'Qc = b'Q^{-1}b = λ b'B^{-1}b <=> λ c'Qc = b'B^{-1}b
    # Hence λ = b'B^{-1}b - β
    B, b, β, λ = Bbβλ(ell.P)
    c = -(B \ b)
    Q = B / λ
    Ellipsoid(Q, c)
end

@recipe function f(ell::LiftedEllipsoid)
    Ellipsoid(ell)
end

abstract type QuadCone{T, P, S} end

# It will not really be the center and the center for z = 0 is the the same as for <h, x> = 0
struct CenterPoint{T}
    h::Vector{T}
end
struct CenterQuadCone{T, P<:AbstractPolynomial{T}, S} <: QuadCone{T, P, S}
    p::P
    Q::Symmetric{S, Matrix{S}}
    h::Vector{Float64} # h is the center
    H::Matrix{Float64}
    vol::S
end
CenterQuadCone(p::P, Q::Symmetric{S, Matrix{S}}, c, H, vol::S) where {T, P<:AbstractPolynomial{T}, S} = CenterQuadCone{T, P, S}(p, Q, c, H, vol)
JuMP.value(p::CenterQuadCone) = CenterQuadCone(JuMP.value(p.p), Symmetric(JuMP.value.(p.Q)), p.h, p.H, JuMP.value(p.vol))
_β(m, h::CenterPoint{T}) where T = -one(T)
_b(m, h::CenterPoint{T}) where T = zeros(T, length(h.h))
QuadCone(p, Q, b, β, h::CenterPoint, H, vol) = CenterQuadCone(p, Q, h.h, H, vol)

samecenter(l1::CenterQuadCone, l2::CenterQuadCone) = l1.h == l2.h

struct InteriorPoint{T}
    h::Vector{T}
end
struct InteriorQuadCone{T, P<:AbstractPolynomial{T}, S} <: QuadCone{T, P, S}
    p::P
    Q::Symmetric{S, Matrix{S}}
    b::Vector{S}
    β::S
    h::Vector{Float64} # h is an interior point
    H::Matrix{Float64}
    vol::S
end
InteriorQuadCone(p::P, Q::Symmetric{S, Matrix{S}}, b::Vector{S}, β::S, c, H, vol::S) where {T, P<:AbstractPolynomial{T}, S} = InteriorQuadCone{T, P, S}(p, Q, b, β, c, H, vol)
JuMP.value(p::InteriorQuadCone) = InteriorQuadCone(JuMP.value(p.p), Symmetric(JuMP.value.(p.Q)), JuMP.value.(p.b), JuMP.value(p.β), p.h, p.H, JuMP.value(p.vol))
_β(m, h::InteriorPoint) = @variable m
_b(m, h::InteriorPoint) = @variable m [1:length(h.h)]
QuadCone(p, Q, b, β, h::InteriorPoint, H, vol) = InteriorQuadCone(p, Q, b, β, h.h, H, vol)

samecenter(l1, l2) = false

_householder(h) = householder([1; h.h]) # We add 1, for z

function lyapconstraint(_p::Function, N, s, l, y, t, m, cone, λuser)
    u = source(s, t)
    v = target(s, t)
    σ = symbol(s.automaton, t)
    startp = lhs(l[u].p, y, s.resetmaps[σ])
    E = s.resetmaps[σ].E
    newp = ATrp(_p(v), y, E)
    if v in N && samecenter(l[u], l[v])
        expr = newp - startp
        @constraint m expr in cone
        1.
    else
        if λuser === nothing
            λ = @variable m lower_bound=0.0
        else
            λ = λuser
        end
        expr = λ * newp - startp
        @constraint m expr in cone
        λ
    end
end

ellipsoid(p::QuadCone{T, P, JuMP.VariableRef}) where {T, P<:AbstractPolynomial{T}} = ellipsoid(JuMP.value(p))

function _HPH(D, d, δ, H)
    P = [δ d'
         d D]
    HPH = H * P * H
end
_HPH(p::CenterQuadCone) = _HPH(p.Q, zeros(size(p.Q, 1)), -1., p.H)
_HPH(p::InteriorQuadCone) = _HPH(p.Q, p.b, p.β, p.H)
function ellipsoid(p::QuadCone)
    LiftedEllipsoid(inv(_HPH(p)))
end


text:
my_fun(x, y) = 2*x+ 3*y

der_my_fun(x, y) = ForwardDiff.derivative(x -> my_fun(x, y), x )


text:
"""
```
    sample_points_on_ellipse(A::Real, B::Real, H::Real, K::Real, τ::Real, N::Integer, α₁::Real, α₂::Real)
```
    Samples N data points in the angle range [α₁, α₂] for an ellipse specified by  semi-major (A) semi-minor (B) axes,
    centroid (H,K) and orientation (τ). All angles are assumed to be specified in radians. The results are returned
    as `[Observations](@ref)`.
"""
function sample_points_on_ellipse(A::Real, B::Real, H::Real, K::Real, τ::Real, N::Integer, α₁::Real, α₂::Real)
    ℳ = [SVector(0.0,0.0) for n = 1:N]
    for (n,α) in enumerate(range(α₁, stop = α₂, length = N))
        x = H + A*cos(α)*cos(τ) - B*sin(α)*sin(τ)
        y = K + A*cos(α)*sin(τ) + B*sin(α)*cos(τ)
        𝐦 = SVector(x,y)
        ℳ[n] = 𝐦
    end
    data = tuple(ℳ)
    return Observations(data)
end


text:
"""
    imshow(arg1; kw...)

Is a simple front end to the [`grdimage`](@ref) program that accepts GMTgrid, GMTimage, 2D array
of floats or strings with file names of grids or images. The normal options of the *grdimage* program
also apply here but some clever guessing of suitable necessary parameters is done if they are not
provided. Contrary to other image producing modules the "show' keyword is not necessary to display the
image. Here it is set by default. If user wants to use *imshow* to create layers of a more complex fig
he can use *show=false* for the intermediate layers.

# Examples
```julia-repl
# Plot vertical shaded illuminated view of the Mexican hat
julia> G = gmt("grdmath -R-15/15/-15/15 -I0.3 X Y HYPOT DUP 2 MUL PI MUL 8 DIV COS EXCH NEG 10 DIV EXP MUL =");
julia> imshow(G, frame="a", shade="+a45")

# Plot a random heat map
julia> imshow(rand(128,128))

# Display a web downloaded jpeg image wrapped into a sinusoidal projection
julia> imshow("http://larryfire.files.wordpress.com/2009/07/untooned_jessicarabbit.jpg", region="d", frame="g", proj="I15", img_in="r", fmt=:jpg)
```
See also: [`grdimage`](@ref)
"""
function imshow(arg1; first=true, kw...)
	# Take a 2D array of floats and turn it into a GMTgrid or if input is a string assume it's a file name
	# In this later case try to figure if it's a grid or an image and act accordingly.
	is_image = false
	if (isa(arg1, String))		# If it's string it has to be a file name. Check extension to see if is an image
		ffname, ext = splitext(arg1)
		ext = lowercase(ext)
		if (ext == ".jpg" || ext == ".tif" || ext == ".tiff" || ext == ".png" || ext == ".bmp" || ext == ".gif")
			is_image = true
		end
		G = arg1
	elseif (isa(arg1, Array{UInt8}))
		G = mat2img(arg1; kw...)
	else
		G = mat2grid(arg1)
	end

	d = KW(kw)
	see = (!haskey(d, :show)) ? true : see = d[:show]	# No explicit 'show' keyword means show=true

	if (is_image)
		if (haskey(d, :D) || haskey(d, :img_in) || haskey(d, :image_in))	# OK, user set -D so don't repeat
			grdimage(G; first=first, show=see, kw...)
		else
			grdimage(G; first=first, D=1, show=see, kw...)
		end
	else
		imshow(G; first=first, kw...)
	end
end

function imshow(arg1::GMTgrid; first=true, kw...)
	# Here the default is to show, but if a 'show' was used let it rule
	d = KW(kw)
	see = (!haskey(d, :show)) ? true : see = d[:show]	# No explicit 'show' keyword means show=true
	opt_p = parse_common_opts(d, "", [:p], first)
	if (opt_p == "")
		grdimage("", arg1; first=first, show=see, kw...)
	else
		zsize = ((val = find_in_dict(d, [:JZ :Jz :zscale :zsize])[1]) !== nothing) ? val : 5
		srf = ((val = find_in_dict(d, [:Q :surf :surftype])[1]) !== nothing) ? val : "i100"
		grdview("", arg1; first=first, show=see, p=opt_p[4:end], JZ=zsize, Q=srf, kw...)
	end
end

function imshow(arg1::GMTimage; first=true, kw...)
	# Here the default is to show, but if a 'show' was used let it rule
	d = KW(kw)
	see = (!haskey(d, :show)) ? true : see = d[:show]	# No explicit 'show' keyword means show=true
	grdimage("", arg1; first=first, D=true, show=see, kw...)
end


text:
export neighborhood,parametric,cartesian,jacobian,normal,shapefuns,shapeders,get_shape_curl
mutable struct subd_point{T}
    cart::Vector{T}
    paracoords::Vector{T}
    J::Array{T,2}
    detJ::T
    normal::SVector{3,T}
    shapes::Matrix{T}
    shapeders::Matrix{T}
end
Base.getindex(p::subd_point{T}, i::Int) where {T} = p.cart[i]
function neighborhood(chart::subd_chart{T}, u) where {T}
    vertices = chart.vertices
    Nv = chart.N
    val = Nv - 6
    shapefunc = shape_function(u[1],u[2],val)
    shape_der = shape_function_der(u[1],u[2],val)
    coords = zeros(3)
    t1 = zeros(3)
    t2 = zeros(3)
    for id = 1:3
        for iv = 1:Nv
            coords[id] += vertices[iv][id]*shapefunc[iv]
            t1[id] += vertices[iv][id]*shape_der[iv,1]
            t2[id] += vertices[iv][id]*shape_der[iv,2]
        end
    end
    J = [t1 t2]'
    norm = [t1[2] * t2[3] - t1[3] * t2[2];
            t1[3] * t2[1] - t1[1] * t2[3];
            t1[1] * t2[2] - t1[2] * t2[1]]
    detJ = sqrt(norm[1]*norm[1] + norm[2]*norm[2] + norm[3]*norm[3])
    norm = norm./detJ
    return subd_point{T}(coords,u,J,detJ,norm,shapefunc,shape_der)
end

function parametric(d::subd_point)
    return d.paracoords
end

function cartesian(d::subd_point)
    v = d.cart
    return v
end

function jacobian(d::subd_point)
    return d.detJ
end

function normal(d::subd_point)
    return d.normal
end

function shapefuns(d::subd_point)
    return d.shapes
end

function shapeder(d::subd_point)
    return d.shapeders
end

function JacobMatrix33(d::subd_point)
    Jacob = d.J
    norm = d.normal
    return [Jacob;norm']
end

function Jacob_inv(d::subd_point)
    J = JacobMatrix33(d)
    invJ = inv(J)
    return invJ[:,1:2]
end

function get_shape_curl(d::subd_point)
    norm = d.normal
    sders = d.shapeders
    invJ = Jacob_inv(d)
    grad = invJ*sders'
    ncol = size(grad,2)
    curl = zeros(size(grad))
    for i = 1 : ncol
        curl[:,i] = cross(norm,grad[:,i])
    end
    return curl'
end

max_stars_repo_name: krcools/CompScienceMeshes

max_issues_repo_name: krcools/CompScienceMeshes

max_forks_repo_name: krcools/CompScienceMeshes


text:
println("""
If this is the first time running Ahorn this might take a while.
The window might stay blank for a long time, this is normal as packages are precompiling in the background.
""")

using Ahorn

Ahorn.displayMainWindow()


text:
      module BHDL
      using Reexport

include("ROOT/placer/main.jl")

@reexport using .placer
end

max_stars_repo_name: bhdl/bhdl

max_issues_repo_name: bhdl/bhdl

max_forks_repo_name: bhdl/bhdl


text:
# Test build_optimizer_model!
@testset "build_optimizer_model!" begin
    # initialize model
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1],
                        derivative_method = OrthogonalCollocation(3))
    @infinite_parameter(m, 0 <= pars[1:2] <= 1, supports = [0, 1])
    @infinite_variable(m, 1 >= x(par) >= 0, Int)
    @infinite_variable(m, y(par, pars) == 2, Bin, start = 0)
    @point_variable(m, x(0), x0)
    @point_variable(m, y(0, [0, 0]), 0 <= y0 <= 1, Int)
    @hold_variable(m, 0 <= z <= 1, Bin)
    @hold_variable(m, w == 1, Int, start = 1)
    @deriv(x, par)
    data1 = DiscreteMeasureData(par, [1, 1], [0, 1])
    meas1 = measure(x - w, data1)
    meas2 = measure(y, data1)
    @constraint(m, c1, x + par - z == 0)
    @constraint(m, c2, z + x0 >= -3)
    @constraint(m, c3, meas1 + z == 0)
    @BDconstraint(m, c4(par in [0.5, 1]), meas2 - 2y0 + x <= 1)
    @constraint(m, c5, meas2 == 0)
    @objective(m, Min, x0 + meas1)
    # test normal usage
    @test isa(build_optimizer_model!(m), Nothing)
    @test optimizer_model_ready(m)
    @test num_variables(optimizer_model(m)) == 14
    # test repeated build
    @test isa(build_optimizer_model!(m), Nothing)
    @test optimizer_model_ready(m)
    @test num_variables(optimizer_model(m)) == 14
end

# Test optimizer model querying methods
@testset "Optimizer Model Queries" begin
    # initialize model
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1],
                        derivative_method = OrthogonalCollocation(3))
    @infinite_variable(m, x(par))
    @point_variable(m, x(0), x0)
    @hold_variable(m, z)
    d1 = @deriv(x, par)
    data1 = DiscreteMeasureData(par, [1, 1], [0, 1])
    meas1 = measure(x - z, data1)
    @constraint(m, c1, x + par - z == 0)
    @constraint(m, c2, z + x0 >= -3)
    @constraint(m, c3, meas1 + z == 0)
    f = parameter_function(sin, par)
    build_optimizer_model!(m)
    tm = optimizer_model(m)
    tdata = transcription_data(tm)
    # Test optimizer_model_variable
    @testset "optimizer_model_variable" begin
        # test normal usage
        @test optimizer_model_variable(x, label = All) == transcription_variable(x, label = All)
        @test optimizer_model_variable(x, label = All, ndarray = true) == transcription_variable(x, label = All, ndarray = true)
        @test optimizer_model_variable(x0) == transcription_variable(x0)
        @test optimizer_model_variable(z) == transcription_variable(z)
        @test optimizer_model_variable(d1, label = InternalLabel) == transcription_variable(d1, label = InternalLabel)
        @test optimizer_model_variable(f) == [0, sin(1)]
        # test fallback
        @test_throws ErrorException optimizer_model_variable(x, Val(:Bad), my_key = true)
    end
    # Test variable_supports
    @testset "variable_supports" begin
        # test finite variable fallback
        @test InfiniteOpt.variable_supports(tm, dispatch_variable_ref(z), :key) == ()
        @test InfiniteOpt.variable_supports(tm, dispatch_variable_ref(x0), :key, label = All) == ()
        # test fallback
        @test_throws ErrorException InfiniteOpt.variable_supports(tm, x, Val(:Bad))
    end
    # Test supports (variables)
    @testset "supports (Variables)" begin
        # test normal usage
        @test supports(x) == [(0.,), (1.,)]
        @test supports(x, ndarray = true) == [(0.,), (1.,)]
        @test supports(x, label = All) == [(0.,), (0.5,), (1.,)]
        @test supports(meas1) == ()
        @test supports(d1, label = InternalLabel) == [(0.5,)]
        @test supports(f, label = All) == [(0.,), (0.5,), (1.,)]
    end
    # Test optimizer_model_expression
    @testset "optimizer_model_expression" begin
        # test variable references
        @test optimizer_model_expression(x, label = All) == transcription_variable(x, label = All)
        @test optimizer_model_expression(z) == transcription_variable(z)
        @test optimizer_model_expression(x0) == transcription_variable(x0)
        @test optimizer_model_expression(x0, ndarray = true) == transcription_variable(x0)
        # test expression without variables
        expr = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}) + 42
        @test optimizer_model_expression(expr) == zero(AffExpr) + 42
        # test normal expressions
        xt = transcription_variable(x, label = All)
        zt = transcription_variable(z)
        @test optimizer_model_expression(x^2 + z) == [xt[1]^2 + zt, xt[3]^2 + zt]
        @test optimizer_model_expression(x^2 + z, ndarray = true) == [xt[1]^2 + zt, xt[3]^2 + zt]
        @test optimizer_model_expression(x^2 + z, label = All) == [xt[1]^2 + zt, xt[2]^2 + zt, xt[3]^2 + zt]
        @test optimizer_model_expression(2z - 3) == 2zt - 3
        @test optimizer_model_expression(2 * f) == [zero(AffExpr), zero(AffExpr) + sin(1) * 2]
        # test fallback
        @test_throws ErrorException optimizer_model_expression(c1, Val(:Bad), my_key = true)
    end
    # Test expression_supports
    @testset "expression_supports" begin
        # test normal usage
        @test InfiniteOpt.expression_supports(tm, x, Val(:TransData)) == [(0.,), (1.,)]
        @test InfiniteOpt.expression_supports(tm, x^2 - x0, Val(:TransData)) == [(0.,), (1.,)]
        @test InfiniteOpt.expression_supports(tm, x^2 - x0, Val(:TransData), label = All) == [(0.,), (0.5,), (1.,)]
        # test fallback
        @test_throws ErrorException InfiniteOpt.expression_supports(tm, z^2, Val(:Bad))
    end
    # Test supports (expressions)
    @testset "supports (Expressions)" begin
        # test normal usage
        @test supports(x) == [(0.,), (1.,)]
        @test supports(2x - z + x0 + 43) == [(0.,), (1.,)]
        @test supports(2x - z + x0 + 43, ndarray = true) == [(0.,), (1.,)]
        expr = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}) + 42
        @test supports(expr, label = All) == ()
    end
    # Test optimizer_model_constraint
    @testset "optimizer_model_constraint" begin
        # test normal usage
        @test optimizer_model_constraint(c1) == transcription_constraint(c1)
        @test optimizer_model_constraint(c2, label = All) == transcription_constraint(c2, label = All)
        @test optimizer_model_constraint(c2, label = All, ndarray = true) == transcription_constraint(c2, label = All, ndarray = true)
        @test optimizer_model_constraint(c3) == transcription_constraint(c3)
        # test fallback
        @test_throws ErrorException optimizer_model_constraint(c1, Val(:Bad), my_key = true)
    end
    # Test constraint_supports
    @testset "constraint_supports" begin
        # test normal usage
        @test InfiniteOpt.constraint_supports(tm, c1) == [(0.,), (1.,)]
        @test InfiniteOpt.constraint_supports(tm, c1, label = All) == [(0.,), (0.5,), (1.,)]
        @test InfiniteOpt.constraint_supports(tm, c1, label = All, ndarray = true) == [(0.,), (0.5,), (1.,)]
        # test fallback
        @test_throws ErrorException InfiniteOpt.constraint_supports(tm, c1, Val(:Bad))
    end
    # Test supports (constraints)
    @testset "supports (Constraints)" begin
        # test normal usage
        @test supports(c1, label = UserDefined) == [(0.,), (1.,)]
        # test fallback
        @test supports(c2) == ()
    end
end

# Test optimize!
@testset "JuMP.optimize!" begin
    # initialize model
    mockoptimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                             eval_objective_value=false)
    m = InfiniteModel(mockoptimizer)
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_parameter(m, 0 <= pars[1:2] <= 1, supports = [0, 1])
    @infinite_variable(m, 1 >= x(par) >= 0, Int)
    @infinite_variable(m, y(par, pars) == 2, Bin, start = 0)
    @point_variable(m, x(0), x0)
    @point_variable(m, y(0, [0, 0]), 0 <= y0 <= 1, Int)
    @hold_variable(m, 0 <= z <= 1, Bin)
    @hold_variable(m, w == 1, Int, start = 1)
    data1 = DiscreteMeasureData(par, [1, 1], [0, 1])
    meas1 = measure(x - w, data1)
    meas2 = measure(y, data1)
    @constraint(m, c1, x + par - z == 0)
    @constraint(m, c2, z + x0 >= -3)
    @constraint(m, c3, meas1 + z == 0)
    @BDconstraint(m, c4(par in [0.5, 1]), meas2 - 2y0 + x <= 1)
    @constraint(m, c5, meas2 == 0)
    @objective(m, Min, x0 + meas1)
    # test normal usage
    @test isa(optimize!(m), Nothing)
    @test optimizer_model_ready(m)
    @test num_variables(optimizer_model(m)) == 8
end

# Test JuMP.result_count
@testset "JuMP.result_count" begin
    # Setup the infinite model
    optimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                         eval_objective_value=false)
    m = InfiniteModel(optimizer)
    tm = optimizer_model(m)
    model = MOIU.Model{Float64}()
    JuMP.optimize!(tm)
    mockoptimizer = JuMP.backend(tm).optimizer.model
    MOI.set(mockoptimizer, MOI.ResultCount(), 2)
    MOI.set(mockoptimizer, MOI.TerminationStatus(), MOI.OPTIMAL)
    @test result_count(m) == 2
end


text:
using GeoJSON, FactCheck, Turf

# Point
pt = GeoJSON.Point(Float64[1,0])
flipped = Turf.flip(pt)
@fact flipped.coordinates --> roughly([0,1])
@fact pt.coordinates --> roughly([1,0]) # not mutated
pt = GeoJSON.Feature(GeoJSON.Point(Float64[1,0]))
flipped = Turf.flip(pt)
@fact flipped.geometry.coordinates --> roughly([0,1])
@fact pt.geometry.coordinates --> roughly([1,0]) # not mutated
line = GeoJSON.LineString(Vector{Float64}[[1,0],[1,0]])
flipped = Turf.flip(line)
@fact flipped.coordinates[1] --> roughly([0,1])
@fact flipped.coordinates[2] --> roughly([0,1])
polygon = GeoJSON.Polygon(Vector{Vector{Float64}}[Vector{Float64}[[1,0],[1,0],[1,12]],
                                                  Vector{Float64}[[.2,.2],[.3,.3],[.1,.2]]])
flipped = Turf.flip(polygon)
@fact flipped.coordinates[1][1] --> roughly([0,1])
@fact flipped.coordinates[1][2] --> roughly([0,1])
@fact flipped.coordinates[1][3] --> roughly([2,1])
@fact flipped.coordinates[2][3] --> roughly([.2,.1])
pt1 = GeoJSON.Feature(GeoJSON.Point(Float64[1,0]))
pt2 = GeoJSON.Feature(GeoJSON.Point(Float64[1,0]))
fc = GeoJSON.FeatureCollection(GeoJSON.Feature[pt1,pt2])
flipped = Turf.flip(fc)
@fact flipped.features[1].geometry.coordinates --> roughly([0,1])
@fact flipped.features[2].geometry.coordinates --> roughly([0,1])

#explode a polygon
polygon = GeoJSON.Polygon(Vector{Vector{Float64}}[Vector{Float64}[[0,0],[10,0],[10,10],[0,10],[0,0]]])
pt1 = GeoJSON.Feature(GeoJSON.Point(Float64[0,0]))
pt2 = GeoJSON.Feature(GeoJSON.Point(Float64[0,10]))
pt3 = GeoJSON.Feature(GeoJSON.Point(Float64[10,10]))
pt4 = GeoJSON.Feature(GeoJSON.Point(Float64[10,0]))
fc = GeoJSON.FeatureCollection(GeoJSON.Feature[pt1,pt2,pt3,pt4])
exploded = Turf.explode(polygon)
@fact typeof(exploded) --> GeoJSON.FeatureCollection
for i=1:length(fc.features)
    @fact typeof(exploded.features[i].geometry) --> GeoJSON.Point
    @fact exploded.features[i].geometry.coordinates --> roughly(fc.features[i].geometry.coordinates)
end
# explode a featurecollection
polygonFC = GeoJSON.FeatureCollection(GeoJSON.Feature[GeoJSON.Feature(polygon)])
exploded = Turf.explode(polygon)
@fact typeof(exploded) --> GeoJSON.FeatureCollection
for i=1:length(fc.features)
    @fact typeof(exploded.features[i].geometry) --> GeoJSON.Point
    @fact exploded.features[i].geometry.coordinates --> roughly(fc.features[i].geometry.coordinates)
end
# explode a single point
point = GeoJSON.Point(Float64[0,0])
fc = GeoJSON.FeatureCollection(GeoJSON.Feature[GeoJSON.Feature(GeoJSON.Point(Float64[0,0]))])
exploded = Turf.explode(point)
@fact typeof(exploded) --> GeoJSON.FeatureCollection
for i=1:length(fc.features)
    @fact typeof(exploded.features[i].geometry) --> GeoJSON.Point
    @fact exploded.features[i].geometry.coordinates --> roughly(fc.features[i].geometry.coordinates)
end
# explode a linestring
line = GeoJSON.LineString(Vector{Float64}[[0,0],[1,1],[0,1],[0,0]])
pt1 = GeoJSON.Feature(GeoJSON.Point(Float64[0,0]))
pt2 = GeoJSON.Feature(GeoJSON.Point(Float64[1,1]))
pt3 = GeoJSON.Feature(GeoJSON.Point(Float64[0,1]))
pt4 = GeoJSON.Feature(GeoJSON.Point(Float64[0,0]))
fc = GeoJSON.FeatureCollection(GeoJSON.Feature[pt1,pt2,pt3,pt4])
exploded = Turf.explode(line)
@fact typeof(exploded) --> GeoJSON.FeatureCollection
for i=1:length(fc.features)
    @fact typeof(exploded.features[i].geometry) --> GeoJSON.Point
    @fact exploded.features[i].geometry.coordinates --> roughly(fc.features[i].geometry.coordinates)
end

pt1 = GeoJSON.Feature(GeoJSON.Point(Float64[50,51]))
pt2 = GeoJSON.Feature(GeoJSON.Point(Float64[100,101]))
multipoint = Turf.combine(GeoJSON.FeatureCollection(GeoJSON.Feature[pt1,pt2]))
@fact typeof(multipoint) --> GeoJSON.MultiPoint
@fact multipoint.coordinates[1] --> roughly([50,51])
@fact multipoint.coordinates[2] --> roughly([100,101])

line1 = GeoJSON.Feature(GeoJSON.LineString(Vector{Float64}[[102,-10],[130,4]]))
line2 = GeoJSON.Feature(GeoJSON.LineString(Vector{Float64}[[40,-20],[150,18]]))
multiline = Turf.combine(GeoJSON.FeatureCollection(GeoJSON.Feature[line1,line2]))
@fact typeof(multiline) --> GeoJSON.MultiLineString
@fact multipoint.coordinates[1][1] --> roughly([102,-10])
@fact multipoint.coordinates[1][2] --> roughly([130,4])
@fact multipoint.coordinates[2][1] --> roughly([40,-20])
@fact multipoint.coordinates[2][2] --> roughly([150,18])

poly1 = GeoJSON.Feature(GeoJSON.Polygon(Vector{Vector{Float64}}[Vector{Float64}[[20,0],[101,0],[101,1],[100,1],[100,0]]]))
poly2 = GeoJSON.Feature(GeoJSON.Polygon(Vector{Vector{Float64}}[Vector{Float64}[[30,0],[102,0],[103,1]]]))
multipoly = Turf.combine(GeoJSON.FeatureCollection(GeoJSON.Feature[poly1,poly2]))
@fact typeof(multipoly) --> GeoJSON.MultiPolygon
@fact multipoint.coordinates[1][1][1] --> roughly([20,0])
@fact multipoint.coordinates[1][1][2] --> roughly([101,0])
@fact multipoint.coordinates[1][1][3] --> roughly([101,1])
@fact multipoint.coordinates[1][1][4] --> roughly([100,1])
@fact multipoint.coordinates[1][1][5] --> roughly([100,0])
@fact multipoint.coordinates[2][1][1] --> roughly([30,0])
@fact multipoint.coordinates[2][1][2] --> roughly([102,0])
@fact multipoint.coordinates[2][1][3] --> roughly([103,1])

@fact Turf.isclockwise(Vector{Float64}[[0,0],[1,1],[1,0],[0,0]]) --> true
@fact Turf.isclockwise(Vector{Float64}[[0,0],[1,0],[1,1],[0,0]]) --> false

hourglass = GeoJSON.parse("""{
        "type": "Polygon",
        "coordinates": [
          [
            [
              -12.034835815429688,
              8.901183448260598
            ],
            [
              -12.060413360595701,
              8.899826693726117
            ],
            [
              -12.036380767822266,
              8.873199368734273
            ],
            [
              -12.059383392333983,
              8.871418491385919
            ],
            [
              -12.034835815429688,
              8.901183448260598
            ]
          ]
        ]
      }""")
hourglasskinks = Turf.kinks(hourglass)
@fact typeof(hourglasskinks) --> GeoJSON.FeatureCollection
@fact length(hourglasskinks.features) --> 2
@fact typeof(hourglasskinks.features[1]) --> GeoJSON.Point

triple = GeoJSON.parse("""{
        "type": "Polygon",
        "coordinates": [
          [
            [
              -44.384765625,
              1.0106897682409128
            ],
            [
              -53.4375,
              0.4833927027896987
            ],
            [
              -43.154296875,
              -6.402648405963884
            ],
            [
              -53.173828125,
              -6.708253968671543
            ],
            [
              -43.857421875,
              -13.752724664396975
            ],
            [
              -54.09667968749999,
              -14.944784875088372
            ],
            [
              -53.3935546875,
              -11.867350911459308
            ],
            [
              -44.384765625,
              1.0106897682409128
            ]
          ]
        ]
      }""")
triplekinks = Turf.kinks(triple)
@fact length(triplekinks.features) --> 6

feature = GeoJSON.parse("""{
    "type": "Feature",
    "geometry": {
      "type": "Polygon",
      "coordinates": [
        [
          [
            -12.034835815429688,
            8.901183448260598
          ],
          [
            -12.060413360595701,
            8.899826693726117
          ],
          [
            -12.036380767822266,
            8.873199368734273
          ],
          [
            -12.059383392333983,
            8.871418491385919
          ],
          [
            -12.034835815429688,
            8.901183448260598
          ]
        ]
      ]
    },
    "properties": null
  }""")

featurekinks = Turf.kinks(feature)
@fact length(featurekinks.features) --> 2


text:
export
    NGSIMEnv,
    reset,
    step,
    observation_space_spec,
    action_space_spec,
    obs_names,
    render

#=
Description:
    NGSIM env that plays NGSIM trajectories, allowing the agent to take the place
    of one of the vehicles in the trajectory
=#
type NGSIMEnv <: Env
    trajdatas::Vector{ListRecord}
    trajinfos::Vector{Dict}
    roadways::Vector{Roadway}
    roadway::Union{Void, Roadway} # current roadway
    scene::Scene
    rec::SceneRecord
    ext::MultiFeatureExtractor
    egoid::Int # current id of relevant ego vehicle
    ego_veh::Union{Void, Vehicle} # the ego vehicle
    traj_idx::Int # current index into trajdatas
    t::Int # current timestep in the trajdata
    h::Int # current maximum horizon for egoid
    H::Int # maximum horizon
    primesteps::Int # timesteps to prime the scene
    Δt::Float64

    # settings
    terminate_on_collision::Bool
    terminate_on_off_road::Bool

    # metadata
    epid::Int # episode id
    render_params::Dict # rendering options
    infos_cache::Dict # cache for infos intermediate results
    function NGSIMEnv(
            params::Dict;
            trajdatas::Union{Void, Vector{ListRecord}} = nothing,
            trajinfos::Union{Void, Vector{Dict}} = nothing,
            roadways::Union{Void, Vector{Roadway}} = nothing,
            reclength::Int = 5,
            Δt::Float64 = .1,
            primesteps::Int = 50,
            H::Int = 50,
            terminate_on_collision::Bool = true,
            terminate_on_off_road::Bool = true,
            render_params::Dict = Dict("zoom"=>5., "viz_dir"=>"/tmp"))
        param_keys = keys(params)
        @assert in("trajectory_filepaths", param_keys)

        # optionally overwrite defaults
        reclength = get(params, "reclength", reclength)
        primesteps = get(params, "primesteps", primesteps)
        H = get(params, "H", H)
        for (k,v) in get(params, "render_params", render_params)
            render_params[k] = v
        end
        terminate_on_collision = get(params, "terminate_on_collision", terminate_on_collision)
        terminate_on_off_road = get(params, "terminate_on_off_road", terminate_on_off_road)

        # load trajdatas if not provided
        if trajdatas == nothing || trajinfos == nothing || roadways == nothing
            trajdatas, trajinfos, roadways = load_ngsim_trajdatas(
                params["trajectory_filepaths"],
                minlength=primesteps + H
            )
        end

        # build components
        scene_length = max_n_objects(trajdatas)
        scene = Scene(scene_length)
        rec = SceneRecord(reclength, Δt, scene_length)
        ext = build_feature_extractor(params)
        infos_cache = fill_infos_cache(ext)
        return new(
            trajdatas,
            trajinfos,
            roadways,
            nothing,
            scene,
            rec,
            ext,
            0, nothing, 0, 0, 0, H, primesteps, Δt,
            terminate_on_collision, terminate_on_off_road,
            0, render_params, infos_cache
        )
    end
end
function reset(
        env::NGSIMEnv;
        offset::Int=env.H + env.primesteps,
        egoid::Union{Void,Int}=nothing,
        start::Union{Void,Int}=nothing,
        traj_idx::Int=1)
    # sample the trajectory, ego vehicle
    env.traj_idx, env.egoid, env.t, env.h = sample_trajdata_vehicle(
        env.trajinfos,
        offset,
        traj_idx,
        egoid,
        start
    )

    # update / reset containers
    env.epid += 1
    empty!(env.rec)
    empty!(env.scene)

    # prime
    for t in env.t:(env.t + env.primesteps)
        update!(env.rec, get!(env.scene, env.trajdatas[env.traj_idx], t))
    end
    # set the ego vehicle
    env.ego_veh = env.scene[findfirst(env.scene, env.egoid)]
    # set the roadway
    env.roadway = env.roadways[env.traj_idx]
    # env.t is the next timestep to load
    env.t += env.primesteps + 1
    # enforce a maximum horizon
    env.h = min(env.h, env.t + env.H)
    return get_features(env)
end

#=
Description:
    Propagate a single vehicle through an otherwise predeterined trajdata

Args:
    - env: environment to be stepped forward
    - action: array of floats that can be converted into an AccelTurnrate
=#
function _step!(env::NGSIMEnv, action::Array{Float64})
    # convert action into form
    ego_action = AccelTurnrate(action...)
    # propagate the ego vehicle
    ego_state = propagate(
        env.ego_veh,
        ego_action,
        env.roadway,
        env.Δt
    )
    # update the ego_veh
    env.ego_veh = Entity(env.ego_veh, ego_state)

    # load the actual scene, and insert the vehicle into it
    get!(env.scene, env.trajdatas[env.traj_idx], env.t)
    vehidx = findfirst(env.scene, env.egoid)
    orig_veh = env.scene[vehidx] # for infos purposes
    env.scene[vehidx] = env.ego_veh

    # update rec with current scene
    update!(env.rec, env.scene)

    # compute info about the step
    step_infos = Dict{String, Float64}()
    step_infos["rmse_pos"] = sqrt(abs2((orig_veh.state.posG - env.ego_veh.state.posG)))
    step_infos["rmse_vel"] = sqrt(abs2((orig_veh.state.v - env.ego_veh.state.v)))
    step_infos["rmse_t"] = sqrt(abs2((orig_veh.state.posF.t - env.ego_veh.state.posF.t)))
    step_infos["x"] = env.ego_veh.state.posG.x
    step_infos["y"] = env.ego_veh.state.posG.y
    step_infos["s"] = env.ego_veh.state.posF.s
    step_infos["phi"] = env.ego_veh.state.posF.ϕ
    return step_infos
end

function _extract_rewards(env::NGSIMEnv, infos::Dict{String, Float64})
    r = 0
    if infos["is_colliding"] == 1
        r -= 1
    end
    if infos["is_offroad"] == 1
        r -= 1
    end
    return r
end

function Base.step(env::NGSIMEnv, action::Array{Float64})
    step_infos = _step!(env, action)
    # compute features and feature_infos
    features = get_features(env)
    feature_infos = _compute_feature_infos(env, features)
    # combine infos
    infos = merge(step_infos, feature_infos)
    # update env timestep to be the next scene to load
    env.t += 1
    # compute terminal
    if env.t >= env.h
        terminal = true
    elseif env.terminate_on_collision && infos["is_colliding"] == 1
        terminal = true
    elseif env.terminate_on_off_road && (abs(infos["markerdist_left"]) > 3 && abs(infos["markerdist_right"]) > 3)
        terminal = true
    else
        terminal = false
    end
    reward = _extract_rewards(env, infos)
    return features, reward, terminal, infos
end
function _compute_feature_infos(env::NGSIMEnv, features::Array{Float64})
    is_colliding = features[env.infos_cache["is_colliding_idx"]]
    markerdist_left = features[env.infos_cache["markerdist_left_idx"]]
    markerdist_right = features[env.infos_cache["markerdist_right_idx"]]
    is_offroad = features[env.infos_cache["out_of_lane_idx"]]
    return Dict{String, Float64}(
        "is_colliding"=>is_colliding,
        "markerdist_left"=>markerdist_left,
        "markerdist_right"=>markerdist_right,
        "is_offroad"=>is_offroad
    )
end
function AutoRisk.get_features(env::NGSIMEnv)
    veh_idx = findfirst(env.scene, env.egoid)
    pull_features!(env.ext, env.rec, env.roadway, veh_idx)
    return deepcopy(env.ext.features)
end
function observation_space_spec(env::NGSIMEnv)
    low = zeros(length(env.ext))
    high = zeros(length(env.ext))
    feature_infos = feature_info(env.ext)
    for (i, fn) in enumerate(feature_names(env.ext))
        low[i] = feature_infos[fn]["low"]
        high[i] = feature_infos[fn]["high"]
    end
    infos = Dict("high"=>high, "low"=>low)
    return (length(env.ext),), "Box", infos
end
action_space_spec(env::NGSIMEnv) = (2,), "Box", Dict("high"=>[4.,.15], "low"=>[-4.,-.15])
obs_names(env::NGSIMEnv) = feature_names(env.ext)

#=
Description:
    Render the scene

Args:
    - env: environment to render

Returns:
    - img: returns a (height, width, channel) image to display
=#
function render(
        env::NGSIMEnv;
        egocolor::Vector{Float64}=[1.,0.,0.],
        camtype::String="follow",
        static_camera_pos::Vector{Float64}=[0.,0.],
        camera_rotation::Float64=0.,
        canvas_height::Int=800,
        canvas_width::Int=800)
    # define colors for all the vehicles
    carcolors = Dict{Int,Colorant}()
    egocolor = ColorTypes.RGB(egocolor...)
    for veh in env.scene
        carcolors[veh.id] = veh.id == env.egoid ? egocolor : colorant"green"
    end

    # define a camera following the ego vehicle
    if camtype == "follow"
        cam = AutoViz.CarFollowCamera{Int}(env.egoid, env.render_params["zoom"])
    elseif camtype == "static"
        cam = AutoViz.StaticCamera(VecE2(static_camera_pos...), env.render_params["zoom"])
    else
        error("invalid camera type $(camtype)")
    end
    stats = [
        CarFollowingStatsOverlay(env.egoid, 2),
        NeighborsOverlay(env.egoid, textparams = TextParams(x = 600, y_start=300))
    ]

    # rendermodel for optional rotation
    # note that for this to work, you have to comment out a line in AutoViz
    # src/overlays.jl:27 `clear_setup!(rendermodel)` in render
    rendermodel = RenderModel()
    camera_rotate!(rendermodel, deg2rad(camera_rotation))

    # render the frame
    frame = render(
        env.scene,
        env.roadway,
        stats,
        rendermodel = rendermodel,
        cam = cam,
        car_colors = carcolors,
        canvas_height=canvas_height,
        canvas_width=canvas_width
    )

    # save the frame
    if !isdir(env.render_params["viz_dir"])
        mkdir(env.render_params["viz_dir"])
    end
    ep_dir = joinpath(env.render_params["viz_dir"], "episode_$(env.epid)")
    if !isdir(ep_dir)
        mkdir(ep_dir)
    end
    filepath = joinpath(ep_dir, "step_$(env.t).png")
    write_to_png(frame, filepath)

    # load and return the frame as an rgb array
    img = PyPlot.imread(filepath)
    return img
end


text:
using SmallInts
using Base.Test

for (n,I) in enumerate((Int1, Int2, Int4))
    nbits = 1 << (n-1)
    minval = (Int8(-1) << (nbits-1))
    maxval = (Int8(1) << (nbits-1)) - Int8(1)
    for x in minval:maxval
        @test Int8(I(x)) === x
        @test I(x).val === x
        @test rem(x, I) === I(x)
        @test convert(I, x) === I(x)
        @test typemin(I).val == ~typemax(I).val

        @test rem(I(x), I) === I(x)
        @test rem(I(x), Bool) === rem(x, Bool)
        @test rem(rem(x, Bool), I) === I(rem(x, Bool) ? 1 : 0)
        @test convert(I, I(x)) === I(x)
        @test convert(I, rem(x, Bool)) === rem(rem(x, Bool), I)
        if x == rem(x, Bool)
            @test convert(Bool, I(x)) === rem(I(x), Bool)
        else
            @test_throws InexactError convert(Bool, I(x))
        end

        if x < 0
            @test leading_zeros(I(x)) == 0
            @test leading_ones(I(x)) == leading_zeros(~I(x))
        else
            @test leading_zeros(I(x)) ==
                leading_zeros(x) - leading_zeros(maxval) + 1
            @test leading_ones(I(x)) == 0
        end
        @test trailing_zeros(I(x)) ==
            min(trailing_zeros(x), trailing_ones(maxval) + 1)
        if x == 0
            @test trailing_ones(I(x)) == 0
        elseif x == -1
            @test trailing_ones(I(x)) == nbits
        else
            @test trailing_ones(I(x)) == trailing_ones(x)
        end

        @test ~I(x) === rem(~x, I)
        @test +I(x) === rem(+x, I)
        @test -I(x) === rem(-x, I)
        @test abs(I(x)) === rem(abs(x), I)

        for y in minval:maxval
            @test I(x) & I(y) === I(x & y)
            @test I(x) | I(y) === I(x | y)
            @test I(x) $ I(y) === I(x $ y)

            if y < 0
                @test I(x) << Int(y) === I(0)
                @test I(x) >> Int(y) === I(x < 0 ? -1 : 0)
                @test I(x) >>> Int(y) === I(0)
            else
                @test I(x) << Int(y) === rem(x << y, I)
                @test I(x) >> Int(y) === I(x >> y)
                @test I(x) >>> Int(y) ===
                    rem((x >>> y) & (Int8(-1) >>> (8 - nbits + y)), I)
            end

            @test (I(x) == I(y)) === (x == y)
            @test (I(x) != I(y)) === (x != y)
            @test (I(x) <= I(y)) === (x <= y)
            @test (I(x) < I(y)) === (x < y)
            @test (I(x) >= I(y)) === (x >= y)
            @test (I(x) > I(y)) === (x > y)

            @test I(x) + I(y) === rem(x + y, I)
            @test I(x) - I(y) === rem(x - y, I)
            @test I(x) * I(y) === rem(x * y, I)
            if y == 0 || (x == minval && y == -1)
                @test_throws DivideError div(I(x), I(y))
                @test_throws DivideError fld(I(x), I(y))
                @test_throws DivideError cld(I(x), I(y))
            else
                @test div(I(x), I(y)) === I(div(x, y))
                @test fld(I(x), I(y)) === I(fld(x, y))
                @test cld(I(x), I(y)) === I(cld(x, y))
            end
            if y == 0
                @test_throws DivideError rem(I(x), I(y))
                @test_throws DivideError mod(I(x), I(y))
            else
                @test rem(I(x), I(y)) === I(rem(x, y))
                @test mod(I(x), I(y)) === I(mod(x, y))
            end
        end
    end

    @test_throws InexactError convert(I, -9)
    @test_throws InexactError convert(I, 8)
end

for (n,U) in enumerate((UInt1, UInt2, UInt4))
    nbits = 1 << (n-1)
    maxval = (0x01 << nbits) - 0x01
    for x in 0x00:maxval
        @test UInt8(U(x)) === x
        @test U(x).val === x
        @test rem(x, U) === U(x)
        @test convert(U, x) === U(x)

        @test rem(U(x), U) === U(x)
        @test rem(U(x), Bool) === rem(x, Bool)
        @test rem(rem(x, Bool), U) === U(rem(x, Bool) ? 1 : 0)
        @test convert(U, U(x)) === U(x)
        @test convert(U, rem(x, Bool)) === rem(rem(x, Bool), U)
        if x == rem(x, Bool)
            @test convert(Bool, U(x)) === rem(U(x), Bool)
        else
            @test_throws InexactError convert(Bool, U(x))
        end

        @test leading_zeros(U(x)) ==
            leading_zeros(x) - leading_zeros(maxval)
        @test leading_ones(U(x)) == leading_zeros(~U(x))
        @test trailing_zeros(U(x)) ==
            min(trailing_zeros(x), trailing_ones(maxval))
        @test trailing_ones(U(x)) == trailing_ones(x)

        @test ~U(x) === rem(~x, U)
        @test +U(x) === rem(+x, U)
        @test -U(x) === rem(-x, U)
        @test abs(U(x)) === rem(abs(x), U)

        for y in 0x00:maxval
            @test U(x) & U(y) === U(x & y)
            @test U(x) | U(y) === U(x | y)
            @test U(x) $ U(y) === U(x $ y)

            @test U(x) << Int(y) === rem(x << y, U)
            @test U(x) >> Int(y) === U(x >> y)
            @test U(x) >>> Int(y) === U(x >>> y)

            @test (U(x) == U(y)) === (x == y)
            @test (U(x) != U(y)) === (x != y)
            @test (U(x) <= U(y)) === (x <= y)
            @test (U(x) < U(y)) === (x < y)
            @test (U(x) >= U(y)) === (x >= y)
            @test (U(x) > U(y)) === (x > y)

            @test U(x) + U(y) === rem(x + y, U)
            @test U(x) - U(y) === rem(x - y, U)
            @test U(x) * U(y) === rem(x * y, U)
            if y == 0
                @test_throws DivideError div(U(x), U(y))
                @test_throws DivideError rem(U(x), U(y))
                @test_throws DivideError fld(U(x), U(y))
                @test_throws DivideError mod(U(x), U(y))
                @test_throws DivideError cld(U(x), U(y))
            else
                @test div(U(x), U(y)) === U(div(x, y))
                @test rem(U(x), U(y)) === U(rem(x, y))
                @test fld(U(x), U(y)) === U(fld(x, y))
                @test mod(U(x), U(y)) === U(mod(x, y))
                @test cld(U(x), U(y)) === U(cld(x, y))
            end
        end
    end

    @test_throws InexactError convert(U, -1)
    @test_throws InexactError convert(U, 16)
end


text:
using DistStat, Test, LinearAlgebra, Pkg
type=[Float64,Float32]

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    using CUDA
    ArrayType = CuArray
else
    ArrayType = Array
end

for T in type
    test1=[1.0 1.0 1.0; -1.0 0.0 2.0; 0.0 1.0 3.0]
    test2=[1.0 2.0 3.0; -1.0 0.0 1.0; 2.0 0.0 1.0]
    test1_conv=ArrayType{T}(test1); test2_conv=ArrayType{T}(test2)
    answer1=Array{T}(undef,3,3); answer2=Array{T}(undef,3,3)
    r1=ArrayType{T}(undef, 3, 3); r2=ArrayType{T}(undef,3,3)

    for i in 1:3
        for j in 1:3
            answer1[i,j]=sqrt(sum((test1[:,i]-test1[:,j]).^2))
            answer2[i,j]=sqrt(sum((test1[:,i]-test2[:,j]).^2))
        end
    end
    answer1=ArrayType{T}(answer1)
    answer2=ArrayType{T}(answer2)

    result1=DistStat.euclidean_distance!(r1,test1_conv,test1_conv)
    result2=DistStat.euclidean_distance!(r2,test1_conv,test2_conv)

    @test isapprox(answer1,result1)
    @test isapprox(answer2,result2)

    zero1=[0.0, 0.0, 0.0]

    @test isapprox(result1[diagind(result1)],zero1)
    @test all(result1 .>= 0)
    @test all(result2 .>= 0)

    test=reshape(collect(-24:24),7,7)
    test_conv=ArrayType{T}(test)
    test_dist=distribute(test_conv)
    cols=test_dist.partitioning[DistStat.Rank()+1][2]

    r=MPIArray{T,2,ArrayType}(undef,7,7)
    answer=Array{T}(undef,7,7)
    for i in 1:7
        for j in 1:7
            answer[i,j]=sqrt(sum((test[:,i]-test[:,j]).^2))
        end
    end
    answer=ArrayType{T}(answer)

    result=DistStat.euclidean_distance!(r,test_dist)

    @test isapprox(answer[:,cols],result.localarray)

    zero2=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

    @test isapprox(answer[diagind(answer)],zero2)

    @test all(result.localarray .>= 0.0)

end


text:
## Dynamics
abstract type AbstractModel <: AbstractFunction end
output_dim(model::AbstractModel) = state_dim(model)

abstract type ContinuousDynamics <: AbstractModel end
@inline dynamics(model::ContinuousDynamics, z::AbstractKnotPoint) =
    dynamics(model, state(z), control(z), time(z))
@inline dynamics(model::ContinuousDynamics, xdot, z::AbstractKnotPoint) =
    dynamics!(model, xdot, state(z), control(z), time(z))

@inline evaluate(model::ContinuousDynamics, x, u, p) =
    dynamics(model, x, u, p.t)
@inline evaluate!(model::ContinuousDynamics, ẋ, x, u, p) =
    dynamics!(model, ẋ, x, u, p.t)

@inline dynamics(model::ContinuousDynamics, x, u, t) = dynamics(model, x, u)
@inline dynamics!(model::ContinuousDynamics, ẋ, x, u, t) = dynamics!(model, ẋ, x, u)
@inline jacobian!(::FunctionSignature, ::UserDefined, model::ContinuousDynamics, J, ẋ, z) =
    jacobian!(model, J, ẋ, state(z), control(z), time(z))


text:
using Distributed, AzManagers, Random, TOML, Test, HTTP, AzSessions, JSON, Pkg

include(joinpath(homedir(), "azmanagers-setup.jl"))

session = AzSession(;protocal=AzClientCredentials, client_id=client_id, client_secret=client_secret)

azmanagers_pinfo = Pkg.project()
pkgs=TOML.parse(read(joinpath(dirname(azmanagers_pinfo.path),"Manifest.toml"), String))
pkg=pkgs["AzManagers"][1]
azmanagers_rev=get(pkg, "repo-rev", "")

@testset "AzManagers, addprocs, ppi=$ppi" for ppi in (1,)
    ninstances = 1
    group = "test$(randstring('a':'z',4))"

    # Set up iteration vars
    url = "https://management.azure.com/subscriptions/$subscriptionid/resourceGroups/$resourcegroup/providers/Microsoft.Compute/virtualMachineScaleSets/$group?api-version=2019-12-01"
    tppi = ppi*ninstances                       # Total number of Julia processes in the entire scale set

    #
    # Unit Test 1 - Create scale set and start Julia processes
    #
    addprocs(templatename, ninstances;
        waitfor = true,
        subscriptionid = subscriptionid,
        resourcegroup = resourcegroup,
        ppi = ppi,
        group = group,
        session = session)

    # Verify that the scale set is present
    _r = HTTP.request("GET", url, Dict("Authorization"=>"Bearer $(token(session))"); verbose=0)
    @test _r.status == 200

    #
    # Unit Test 2 - Verify that (Total # of Julia processes specified) == (Total # of Julia processes actual)
    #
    @test nworkers() == tppi

    #
    # Unit Test 3 - Verify that there is a healthy connection to each node
    #
    myworkers = [remotecall_fetch(gethostname, workers()[i]) for i=1:tppi]

    #
    # Unit Test 4 - Verify that none of the created Julia processes are the master Julia process
    #
    master = gethostname()
    unique_workers = unique(myworkers)

    @test length(unique_workers) == ninstances
    for worker in myworkers
        @test master != worker
    end

    #
    # Unit Test 5 - Verify that the cloud-init startup script ran successfully
    #
    for i = 1:tppi
        @test_broken remotecall_fetch(isfile, workers()[i], ".git-credentials")
    end

    #
    # Unit Test 6 - Delete the Julia processes, scale set instances and the scale set itself
    #

    # First, verify that the scale set is present
    _r = HTTP.request("GET", url, Dict("Authorization"=>"Bearer $(token(session))"); verbose=0)
    @test _r.status == 200

    @info "Deleting cluster..."
    rmprocs(workers())

    # Last, verify that the scale set has been deleted
    itry = 0
    while true
        itry += 1
        try
            HTTP.request("GET", url, Dict("Authorization"=>"Bearer $(token(session))"); verbose=0)
        catch _e
            e = JSON.parse(String(_e.response.body))
            if _e.status == 404 && e["error"]["code"] == "ResourceNotFound"
                @info "Cluster deleted!"
                break
            end
            if itry == 10
                @warn "cluster not deleted"
                break
            end
        end
        sleep(5)
    end
end

@testset "environment, addproc" begin
    mkpath("myproject")
    cd("myproject")
    Pkg.activate(".")
    Pkg.add("AzSessions")
    Pkg.add("Distributed")
    Pkg.add("JSON")
    Pkg.add("HTTP")

    Pkg.add(PackageSpec(name="AzManagers", rev=azmanagers_rev))

    r = randstring('a':'z',4)
    bname = "test$r"

    testvm = addproc(templatename; basename=bname, session=session, customenv=true)
    testjob = @detachat testvm begin
        using Pkg
        pinfo = Pkg.project()
        write(stdout, "project path is $(pinfo.path)\n")
    end
    wait(testjob)
    @test contains(read(testjob), "myproject")
    rmproc(testvm; session=session)
end

@testset "environment, addprocs" begin
    mkpath("myproject")
    cd("myproject")
    Pkg.activate(".")
    Pkg.add("AzSessions")
    Pkg.add("Distributed")
    Pkg.add("JSON")
    Pkg.add("HTTP")

    Pkg.add(PackageSpec(name="AzManagers", rev=azmanagers_rev))

    group = "test$(randstring('a':'z',4))"

    addprocs(templatename, 1; waitfor=true, group=group, session=session, customenv=true)
    @everywhere using Pkg
    pinfo = remotecall_fetch(Pkg.project, workers()[1])
    rmprocs(workers())

    @test contains(pinfo.path, "myproject")
end

@testset "AzManagers, addproc" begin
    r = randstring('a':'z',4)
    basename = "test$r"
    testvm = addproc(templatename; basename=basename, session=session)
    testjob = @detachat testvm begin
        write(stdout, "write to stdout\n")
        write(stderr, "write to stderr\n")
    end
    wait(testjob)
    @test read(testjob) == "write to stdout\n"
    @test read(testjob; stdio=stderr) == "write to stderr\n"
    rmproc(testvm; session=session)

    testvm = addproc(templatename, name=basename, session=session)
    testjob = @detachat testvm begin
        write(stdout, "write to stdout\n")
        write(stderr, "write to stderr\n")
    end
    wait(testjob)
    @test read(testjob) == "write to stdout\n"
    @test read(testjob; stdio=stderr) == "write to stderr\n"
    rmproc(testvm; session=session)
end

@testset "AzManagers, detach" for kwargs in ( (dummy="dummy"), )

    #
    # Unit Test 1 - Create a detached job and persist the server
    #
    job1 = @detach vm(;vm_template=templatename, session=session, persist=true) begin
        write(stdout, "job1 - stdout string")
        write(stderr, "job1 - stderr string")
    end

    #
    # Unit Test 2 - Send a new job to the server started above
    #
    job2 = @detachat job1.vm begin
        write(stdout, "job2 - stdout string")
        write(stderr, "job2 - stderr string")
    end

    # Wait for jobs to finish
    wait(job1)
    wait(job2)

    @test status(job1) == "done"
    @test read(job1;stdio=stdout) == "job1 - stdout string"
    @test read(job1;stdio=stderr) == "job1 - stderr string"

    @test status(job2) == "done"
    @test read(job2;stdio=stdout) == "job2 - stdout string"
    @test read(job2;stdio=stderr) == "job2 - stderr string"

    #
    # Unit Test 3 - shut-down the detached server
    #
    rmproc(job1.vm; session=session)

    #
    # Unit Test 4 - create a new job on a new detached server that auto-destructs upon completion of its work
    #
    job3 = @detach vm(;vm_template=templatename, session=session, persist=false) begin
    end
end

@testset "AzManagers, detach, variablebundle" begin
    r = randstring('a':'z',4)
    basename = "test$r"
    testvm = addproc(templatename; basename=basename, session=session)
    variablebundle!(a=1.0,b=3.14)
    testjob = @detachat testvm begin
        if variablebundle(:a) ≈ 1.0 && variablebundle(:b) ≈ 3.14
            write(stdout, "passed")
        else
            write(stdout, "failed")
        end
    end
    wait(testjob)
    @test contains(read(testjob), "passed")
end


text:
@testset "Conjugate Gradient" begin
    # TODO: Investigate the exceptions (could be they just need more iterations?)
    # Test Optim.cg for all differentiable functions in MultivariateProblems.UnconstrainedProblems.examples
    # Trigonometric gets stuck in a local minimum?
    skip = ("Trigonometric", "Extended Powell")
    run_optim_tests(ConjugateGradient(),
                    skip=skip,
                    convergence_exceptions = (("Powell", 1), ("Powell", 2), ("Polynomial", 1),
                                              ("Extended Rosenbrock", 1),
                                              ("Extended Powell", 1),
                                              ("Extended Powell", 2)),
                    minimum_exceptions = (("Paraboloid Diagonal", 1)),
                    minimizer_exceptions = (("Paraboloid Diagonal", 1),
                                            ("Extended Powell", 1),
                                            ("Extended Powell", 2)),
                    f_increase_exceptions = (("Hosaki"),),
                    iteration_exceptions = (("Paraboloid Diagonal", 10000),),
                    show_name = debug_printing)

    @testset "matrix input" begin
        cg_objective(X, B) = sum(abs2, X .- B)/2

        function cg_objective_gradient!(G, X, B)
            for i = 1:length(G)
                G[i] = X[i]-B[i]
            end
        end

        Random.seed!(1)
        B = rand(2,2)
        results = Optim.optimize(X -> cg_objective(X, B), (G, X) -> cg_objective_gradient!(G, X, B), rand(2,2), ConjugateGradient())
        @test Optim.converged(results)
        @test Optim.minimum(results) < 1e-8
    end

    @testset "Undefined beta_k behaviour" begin
        # Ref #669
        # If y_k is zero, then betak is undefined.
        # Check that we don't produce NaNs

        f(x) = 100 - x[1] + exp(x[1] - 100)
        g!(grad, x) = grad[1] = -1 + exp(x[1] - 100)
        res = optimize(f, g!, [0.0], ConjugateGradient(alphaguess=LineSearches.InitialStatic(alpha=1.0), linesearch=LineSearches.BackTracking()))
        @test Optim.converged(res)
        @test Optim.minimum(res) ≈ 1.0
    end
end


text:
using Mill, JSON, BSON, Flux, JsonGrinder, Test
using HierarchicalUtils
import HierarchicalUtils: NodeType, children, InnerNode, LeafNode, printtree
using JsonGrinder: DictEntry, suggestextractor, schema
using Mill: reflectinmodel

j1 = JSON.parse("""{"a": 4, "b": {"a":[1,2,3],    "b": 1},"c": { "a": {"a":[1,2,3],"b":[4,5,6]}}}""",inttype=Float64)
j2 = JSON.parse("""{"a": 4, "c": {"a":{"a":[2,3], "b":[5,6]}}}""")
j3 = JSON.parse("""{"a": 4, "b": {"a":[1,2,3,4],  "b": 1}}""")
j4 = JSON.parse("""{"a": 4, "b": {}}""")
j5 = JSON.parse("""{"b": {}}""")
j6 = JSON.parse("""{}""")

sch = schema([j1,j2,j3,j4,j5,j6])
ext = suggestextractor(sch)

@testset "printtree" begin
    @test buf_printtree(ext, trav=true) ==
    """
    Dict [""]
      ├── a: Float32 ["E"]
      ├── b: Dict ["U"]
      │        ├── a: Array of ["Y"]
      │        │        └── Float32 ["a"]
      │        └── b: Float32 ["c"]
      └── c: Dict ["k"]
               └── a: Dict ["s"]
                        ├── a: Array of ["u"]
                        │        └── Float32 ["v"]
                        └── b: Array of ["w"]
                                 └── Float32 ["x"]"""
end

@testset "nnodes" begin
    @test nnodes(ext) == 12
    @test nnodes(ext[:a]) == 1
    @test nnodes(ext[:b]) == 4
    @test nnodes(ext[:c]) == 6
end

@testset "nleafs" begin
    @test nleafs(ext[:a]) + nleafs(ext[:b]) + nleafs(ext[:c]) == nleafs(ext)
end

@testset "children" begin
    @test children(ext) == (; :a=>ext[:a], :b=>ext[:b], :c=>ext[:c])
    @test children(ext[:a]) == ()
    @test children(ext[:b]) == (; :a=>ext[:b][:a], :b=>ext[:b][:b])
    @test children(ext[:b][:a]) == (ext[:b][:a].item,)
    @test children(ext[:b][:b]) == ()
    @test children(ext[:c]) == (; :a=>ext[:c][:a])
    @test children(ext[:c][:a]) == (; :a=>ext[:c][:a][:a], :b=>ext[:c][:a][:b])
    @test children(ext[:c][:a][:a]) == (ext[:c][:a][:a].item,)
    @test children(ext[:c][:a][:b]) == (ext[:c][:a][:b].item,)
end

@testset "nchildren" begin
    @test nchildren(ext) == 3
    @test nchildren(ext[:a]) == 0
    @test nchildren(ext[:b]) == 2
    @test nchildren(ext[:b][:a]) == 1
    @test nchildren(ext[:b][:b]) == 0
    @test nchildren(ext[:c]) == 1
    @test nchildren(ext[:c][:a]) == 2
    @test nchildren(ext[:c][:a][:a]) == 1
    @test nchildren(ext[:c][:a][:b]) == 1
end

@testset "getindex on strings" begin
    @test ext[""] == ext
    @test ext["E"] == ext[:a]
    @test ext["U"] == ext[:b]
    @test ext["Y"] == ext[:b][:a]
    @test ext["a"] == ext[:b][:a].item
    @test ext["c"] == ext[:b][:b]
    @test ext["k"] == ext[:c]
    @test ext["s"] == ext[:c][:a]
    @test ext["u"] == ext[:c][:a][:a]
    @test ext["v"] == ext[:c][:a][:a].item
    @test ext["w"] == ext[:c][:a][:b]
    @test ext["x"] == ext[:c][:a][:b].item
end

@testset "NodeIterator" begin
    @test collect(NodeIterator(ext)) == [ext[""], ext["E"], ext["U"], ext["Y"], ext["a"], ext["c"],
        ext["k"], ext["s"], ext["u"], ext["v"], ext["w"], ext["x"]]
end

@testset "LeafIterator" begin
    @test collect(LeafIterator(ext)) == [ext["E"], ext["a"], ext["c"], ext["v"], ext["x"]]
end

@testset "TypeIterator" begin
    @test collect(TypeIterator(ExtractArray, ext)) == [ext["Y"], ext["u"], ext["w"]]
end

@testset "show" begin
    e = ExtractCategorical(["a","b"])
    @test buf_printtree(e) == """Categorical d = 3"""

    e = ExtractOneHot(["a","b"], "name", nothing)
    @test buf_printtree(e) == """OneHot d = 3"""

    dict = Dict("a" => ExtractArray(ExtractScalar(Float64,2,3)),"b" => ExtractArray(ExtractScalar(Float64,2,3)))
    br = ExtractDict(dict)
    @test buf_printtree(br, trav=true) ==
    """
    Dict [""]
      ├── a: Array of ["E"]
      │        └── Float64 ["M"]
      └── b: Array of ["U"]
               └── Float64 ["c"]"""

    dict = Dict("a" => ExtractScalar(Float64,2,3),"b" => ExtractScalar(Float64), "c" => ExtractArray(ExtractScalar(Float64,2,3)))
    br = ExtractDict(dict)
    @test buf_printtree(br, trav=true) ==
    """
    Dict [""]
      ├── a: Float64 ["E"]
      ├── b: Float64 ["U"]
      └── c: Array of ["k"]
               └── Float64 ["s"]"""

    other1 = Dict("a" => ExtractArray(ExtractScalar(Float64,2,3)),"b" => ExtractArray(ExtractScalar(Float64,2,3)))
    br1 = ExtractDict(other1)
    other = Dict("a" => ExtractArray(br1), "b" => ExtractScalar(Float64,2,3))
    br = ExtractDict(other)
    @test buf_printtree(br, trav=true) ==
    """
    Dict [""]
      ├── a: Array of ["E"]
      │        └── Dict ["M"]
      │              ├── a: Array of ["O"]
      │              │        └── Float64 ["P"]
      │              └── b: Array of ["Q"]
      │                       └── Float64 ["R"]
      └── b: Float64 ["U"]"""
end


text:
#=
  Bayesian Data Analysis of Wine Quality dataset

  I'm continue to testingMvNormal models to see how well it does on
  various datasets.

  From http://archive.ics.uci.edu/ml/datasets/Wine+Quality
  """
  Data Set Information:

  The two datasets are related to red and white variants of the Portuguese "Vinho Verde" wine. For more details, consult: [Web Link] or the reference [Cortez et al., 2009]. Due to privacy and logistic issues, only physicochemical (inputs) and sensory (the output) variables are available (e.g. there is no data about grape types, wine brand, wine selling price, etc.).

  These datasets can be viewed as classification or regression tasks. The classes are ordered and not balanced (e.g. there are many more normal wines than excellent or poor ones). Outlier detection algorithms could be used to detect the few excellent or poor wines. Also, we are not sure if all input variables are relevant. So it could be interesting to test feature selection methods.

  Attribute Information:

  For more information, read [Cortez et al., 2009].
  Input variables (based on physicochemical tests):
  1 - fixed acidity
  2 - volatile acidity
  3 - citric acid
  4 - residual sugar
  5 - chlorides
  6 - free sulfur dioxide
  7 - total sulfur dioxide
  8 - density
  9 - pH
  10 - sulphates
  11 - alcohol
  Output variable (based on sensory data):
  12 - quality (score between 0 and 10)
  """

  Dataset: http://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/

  The datasets are in the files wine_quality-white.csv and winequality-red.csv.
  Here I use the dataset for white wines, and will see the problem as a
  classification problem.

  Summary for white wines:

  Num correct: 2622.0
  Num incorrect: 2276.0
  Accuracy: 0.5353205389955084
  Confusion matrix:
  7×7 Matrix{Float64}:
  0.0  0.0    8.0    11.0    1.0  0.0  0.0
  0.0  0.0    0.0     0.0    0.0  0.0  0.0
  0.0  0.0  673.0   935.0   12.0  0.0  0.0
  0.0  0.0  277.0  1757.0  164.0  0.0  0.0
  0.0  0.0   19.0   669.0  192.0  0.0  0.0
  0.0  0.0    1.0   115.0   59.0  0.0  0.0
  0.0  0.0    0.0     2.0    3.0  0.0  0.0
  num bad: 2364

  Real bad: when the difference between true and predicted >= 2:
  Number of real bad predictions (diff >= 2): 245  (0.0500204164965292%)

  Summary for red wines:
Predictions:
Num correct: 982.0
Num incorrect: 617.0
Accuracy: 0.6141338336460288
Confusion matrix:
7×7 Matrix{Float64}:
 0.0  0.0    9.0    1.0   0.0  0.0  0.0
 0.0  0.0    0.0    0.0   0.0  0.0  0.0
 0.0  0.0  518.0  214.0   2.0  0.0  0.0
 0.0  0.0  183.0  433.0  22.0  0.0  0.0
 0.0  0.0    6.0  162.0  31.0  0.0  0.0
 0.0  0.0    0.0   11.0   7.0  0.0  0.0
 0.0  0.0    0.0    0.0   0.0  0.0  0.0
num bad: 653

Real bad: when the difference between true and predicted >= 2:
Number of real bad predictions (diff >= 2): 47  (0.029393370856785492%)

=#
using Turing, StatsPlots, Distributions
using DataFrames, StatsBase
using CSV
include("jl_utils.jl")

#
# Read dataset and preprocessing
#
df = CSV.read("winequality-white.csv",DataFrame,header=true)
# df = CSV.read("winequality-red.csv",DataFrame,header=true)

display(describe(df))
num_rows,num_cols = size(df)
data = Matrix(df)
x = data[:,1:end-1]
y = data[:,end]

categories = unique(y) |> sort
println("Mean value for each categories: ", categories)
display(categories .|> s -> (s,mean(x[y.==s,:], dims=1)))
println("\n")

# Standardize the dataset (better and faster)
# x = StatsBase.standardize(UnitRangeTransform,x) # This give very bad result
means = mean(x,dims=1)
stds = std(x,dims=1)
maxes = [maximum(x[:,i]) for i in 1:size(x)[2]]
x = (x .- means) ./ stds
# x = (x .- means) ./ maxes'

println("mean(x): ", mean(x,dims=1))
println("std(x): ", std(x,dims=1))

@model function wine_quality_model(x,y=missing, ::Type{T} = Float64) where {T}
    n,p = size(x)

    if y === missing
      y = Vector{T}(undef, n)
    end

    a ~ Normal(0,1)
    params ~ filldist(Normal(0,1),p)
    sigma ~ Uniform(0,1)

    # for i in 1:n
    #   mu = a .+ dot(x[i,:], params)
    #   y[i] ~ Normal(mu,sigma)
    # end

    mu = a .+ x*params
    y ~ MvNormal(mu, sigma)

    return y

end

model = wine_quality_model(x,y)

# println("Prior:")
# chns = sample(model, Prior(), 10_000)
# display(chns)

println("\nModel:")
# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 1_000)
# chns = sample(model, PG(5), 1_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, HMC(0.1,6), 1_000)
chns = sample(model, NUTS(), 1_000)

chns_orig = chns
display(chns)

# display(plot(chns))
println("\nPredictions:")

# mpred = iris_model(x[50:51,:],missing)
mpred = wine_quality_model(x,missing)
pred = predict(mpred,chns)
# display(pred)

cats = Dict(3=>1,4=>3,5=>3,6=>4,7=>5,8=>6,9=>7)
confusion, bads = confusion_matrix(y,pred,cats,false)
println("Confusion matrix:")
display(confusion)
# println("\nBad classfied:")
# display(bads)
println("num bad: $(length(bads))")
println("\nReal bad: when the difference between true and predicted >= 2:")
num_real_bad = 0
for (b1,b2) in bads
  if abs(b2[1]-b2[2]) >= 2
    # println("b21:$(b2[1]) b22:$(b2[2])")
    global num_real_bad += 1
  end
end
println("Number of real bad predictions (diff >= 2): $num_real_bad  ($(num_real_bad / length(y))%)")

# println("\nbads: credible_intervals")
# for b in bads
#   println("b:",b)
#   ix = b[1]
#   println("x[$ix]: $(x[ix,:])  true cat (y[$ix]): $(y[ix])  (pred mean: ", mean(pred["y[$ix]"]))
#   # credible_interval(pred,"y[$(b[1])]",0.3)
# end

# println("\nMean value for each category: ", 1:3)
# display(1:3 .|> s -> (s,mean(x[y.==s,:], dims=1)))
# println("\n")

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank