text:
import DrWatson: quickactivate
quickactivate(@__DIR__, "Chemostat_Folsom2014")

@time begin
    import SparseArrays
    import Base.Threads: @threads, threadid

    # -------------------------------------------------------------------
    import Chemostat_Folsom2014
    const ChF = Chemostat_Folsom2014

    const iJR = ChF.iJR904
    const Fd = ChF.FolsomData # experimental data
    const Bd = ChF.BegData    # cost data

    # -------------------------------------------------------------------
    # run add "https://github.com/josePereiro/Chemostat" in the
    # julia Pkg REPL for installing the package
    import Chemostat
    import Chemostat.LP.MathProgBase

    const Ch = Chemostat
    const ChU = Ch.Utils
    const ChSS = Ch.SteadyState
    const ChLP = Ch.LP
    const ChEP = Ch.MaxEntEP
    const ChSU = Ch.SimulationUtils

    import ChemostatPlots
    const ChP = ChemostatPlots

    import UtilsJL
    const UJL = UtilsJL

    using Statistics
    using Base.Threads
    using ProgressMeter
    using Serialization

    # -------------------------------------------------------------------
    using Plots, FileIO
    import GR
    GR.inline("png")

end

## -------------------------------------------------------------------
INDEX = ChU.load_data(iJR.MAXENT_VARIANTS_INDEX_FILE; verbose = false);

# -------------------------------------------------------------------
const ME_Z_OPEN_G_OPEN        = :ME_Z_OPEN_G_OPEN
const ME_MAX_POL              = :ME_MAX_POL
const ME_Z_EXPECTED_G_MOVING  = :ME_Z_EXPECTED_G_MOVING
const ME_Z_FIXXED_G_BOUNDED   = :ME_Z_FIXXED_G_BOUNDED
const ME_Z_EXPECTED_G_BOUNDED = :ME_Z_EXPECTED_G_BOUNDED

ALL_METHODS = [
    # ME_Z_OPEN_G_OPEN,
    ME_MAX_POL,
    # ME_Z_FIXXED_G_BOUNDED,
    # ME_Z_EXPECTED_G_BOUNDED,
    # ME_Z_EXPECTED_G_MOVING
]

# -------------------------------------------------------------------
fileid = "2.1"
function mysavefig(p, pname; params...)
    fname = UJL.mysavefig(p, string(fileid, "_", pname), iJR.MODEL_FIGURES_DIR; params...)
    @info "Plotting" fname
end
myminmax(a::Vector) = (minimum(a), maximum(a))
EXCH_FLX_IDERS = ["GLC", "PYR", "SUCC", "LAC", "FORM", "AC", "O2", "CO2"]

Fd_rxns_map = iJR.load_rxns_map()
Fd_mets_map = iJR.load_mets_map()

## -------------------------------------------------------------------
# collect DAT
include("2.2.1_collect_DAT.jl")

## -------------------------------------------------------------------
EXPS = DAT[:EXPS]

exp_colors = let
    colors = Plots.distinguishable_colors(length(EXPS))
    Dict(exp => color for (exp, color) in zip(EXPS, colors))
end

ider_colors = let
    iders = [EXCH_FLX_IDERS; "D"]
    colors = Plots.distinguishable_colors(length(iders))
    Dict(ider => color for (ider, color) in zip(iders, colors))
end

method_colors = Dict(
    ME_Z_OPEN_G_OPEN => :red,
    ME_MAX_POL => :orange,
    ME_Z_FIXXED_G_BOUNDED => :green,
    ME_Z_EXPECTED_G_MOVING => :purple,
    ME_Z_EXPECTED_G_BOUNDED => :blue,
)

## -------------------------------------------------------------------
# PLOTS
# -------------------------------------------------------------------

## -------------------------------------------------------------------
# polytope box volume
include("2.2.2_pol_vox_volume.jl")

## -------------------------------------------------------------------
# Var study
include("2.2.3_var_study.jl")

## -------------------------------------------------------------------
# MSE study
include("2.2.4_MSE_study.jl")

## -------------------------------------------------------------------
# proj 2D
include("2.2.5_proj2d_study.jl")

## -------------------------------------------------------------------
# Beta tendencies
include("2.2.6_Beta_tendencies.jl")

## -------------------------------------------------------------------
# flx correlations
include("2.2.7_flx_correlations.jl")

## -------------------------------------------------------------------
# inbounds study
include("2.2.8_inbound_study.jl")

## -------------------------------------------------------------------
# marginals
include("2.2.9_marginals.jl")

## -------------------------------------------------------------------
# legends
include("2.2.10_leyends.jl")


text:
export FiniteDiff

#TODO add boundary conditions

"""
`FiniteDiff([domainType=Float64::Type,] dim_in::Tuple, direction = 1)`

`FiniteDiff(x::AbstractArray, direction = 1)`

Creates a `LinearOperator` which, when multiplied with an array `x::AbstractArray{N}`, returns the discretized gradient over the specified `direction` obtained using forward finite differences.

```julia
julia> FiniteDiff(Float64,(3,))
δx  ℝ^3 -> ℝ^2

julia> FiniteDiff((3,4),2)
δy  ℝ^(3, 4) -> ℝ^(3, 3)

julia> all(FiniteDiff(ones(2,2,2,3),1)*ones(2,2,2,3) .== 0)
true

```

"""
struct FiniteDiff{T,N,D,C <: CartesianIndices{N}} <: LinearOperator
	dim_in::NTuple{N,Int}
    idx::C
	function FiniteDiff{T,N,D}(dim_in) where {T,N,D}
		D > N && error("direction is bigger the number of dimension $N")
        idx = CartesianIndices(([i == D ? (2:d) : (1:d) for (i,d) in enumerate(dim_in)]...,))
        new{T,N,D,typeof(idx)}(dim_in,idx)
	end
end

# Constructors
#default constructor
FiniteDiff(domainType::Type, dim_in::NTuple{N,Int}, dir::Int64 = 1) where {N} =
FiniteDiff{domainType,N,dir}(dim_in)

FiniteDiff(dim_in::NTuple{N,Int}, dir::Int64 = 1) where {N} =
FiniteDiff(Float64, dim_in, dir)

FiniteDiff(x::AbstractArray{T,N}, dir::Int64 = 1) where {T,N}  = FiniteDiff(eltype(x), size(x), dir)

# Mappings

@generated function mul!(y::AbstractArray{T,N},
                         L::FiniteDiff{T,N,D},
                         b::AbstractArray{T,N}) where {T,N,D}
	z = zeros(Int,N)
	z[D] = 1
	idx = CartesianIndex(z...)
	ex = quote
        for I in L.idx
			y[I-$idx] = b[I]-b[I-$idx]
		end
		return y
	end
end

@generated function mul!(y::AbstractArray{T,N},
                         L::AdjointOperator{FiniteDiff{T,N,D,C}},
                         b::AbstractArray{T,N}) where {T,N,D,C}
	z = zeros(Int,N)
	z[D] = 1
	idx = CartesianIndex(z...)
	ex = quote
		for I in CartesianIndices(size(y))
			y[I] =
			I[$D] == 1 ? -b[I]  :
			I[$D] == size(y,$D) ?   b[I-$idx]  : -b[I]+b[I-$idx]
		end
		return y
	end
end

# Properties

domainType(L::FiniteDiff{T, N}) where {T, N} = T
codomainType(L::FiniteDiff{T, N}) where {T, N} = T

function size(L::FiniteDiff{T,N,D}) where {T,N,D}
	dim_out = [L.dim_in...]
	dim_out[D] = dim_out[D]-1
	return ((dim_out...,), L.dim_in)
end

fun_name(L::FiniteDiff{T,N,1}) where  {T,N} = "δx"
fun_name(L::FiniteDiff{T,N,2}) where  {T,N} = "δy"
fun_name(L::FiniteDiff{T,N,3}) where  {T,N} = "δz"
fun_name(L::FiniteDiff{T,N,D}) where {T,N,D}  = "δx$D"

is_full_row_rank(L::FiniteDiff) = true

ext: jl


text:
#module map_covid19
using Dates, DataFrames, ColorSchemes, PlotlyJS

#include("load_covid19_data.jl")

reasonable_resolution() = (1440, 900)

# time index
#for ti = 8:length(covid19us[1].confirmed)
#for ti = length(covid19us[1].confirmed)-length(covid19us[1].confirmed)+8 : length(covid19us[1].confirmed)
for ti = length(covid19us[1].confirmed)-10 : length(covid19us[1].confirmed)
#ti = 200
#ti = length(covid19us[1].confirmed)

display(covid19us[1].time[ti])

local cconfirmed = [covid19us[i].confirmed[ti] for i in 1:length(covid19us)];
local dcconfirmed = [mean(covid19us[i].confirmed[ti-7+1:ti] - covid19us[i].confirmed[ti-8+1:ti-1]) for i in 1:length(covid19us)];
local dcconfirmedinst = [mean(covid19us[i].confirmed[ti] - covid19us[i].confirmed[ti-1]) for i in 1:length(covid19us)];
local cdeath = [covid19us[i].death[ti] for i in 1:length(covid19us)];
local dcdeath = [mean(covid19us[i].death[ti-7+1:ti] - covid19us[i].death[ti-8+1:ti-1]) for i in 1:length(covid19us)];
local dcdeathinst = [mean(covid19us[i].death[ti] - covid19us[i].death[ti-1]) for i in 1:length(covid19us)];

# extracting county indices with land area and population
local cind = findall((countyarea .> 0) .& (cpop .> 0));

local tind = 1:length(covid19us[1].confirmed);
local t = covid19us[1].time[tind];

# plot a map of confirmed cases
#cval = cconfirmed[cind] ./ (countyarea[cind] .* cpop[cind]);
local cval = cconfirmed[cind] ./ cpop[cind];
local log10cval = log10.(cval);
#sval = countyarea[cind] .* cpop[cind];
local sval = cpop[cind];
local log10sval = log10.(sval);

# daily confirmed cases per 100k averaged over 7 days
local dCONFIRMEDpc = dcconfirmed[cind] ./ cpop[cind] .* 1e5;
local bind = findall(dCONFIRMEDpc .<= 0);
local dCONFIRMEDpc[bind] .= NaN;

# daily instant confirmed cases per 100k
local dCONFIRMEDpcinst = dcconfirmedinst[cind] ./ cpop[cind] .* 1e5;
local bind = findall(dCONFIRMEDpcinst .<= 0);
local dCONFIRMEDpcinst[bind] .= NaN;

# daily death per 100k averaged over 7 days
local dDEATHpc = dcdeath[cind] ./ cpop[cind] .* 1e5;
local bind = findall(dDEATHpc .<= 0);
local dDEATHpc[bind] .= NaN;

# daily instant death per 100k
local dDEATHpcinst = dcdeathinst[cind] ./ cpop[cind] .* 1e5;
local bind = findall(dDEATHpcinst .<= 0);
local dDEATHpcinst[bind] .= NaN;

local df = DataFrame(LON = clon[cind],
    LAT = clat[cind],
    POPULATION = cpop[cind],
    log10POPULATION = log10.(cpop[cind]),
    logPOPULATION = log.(cpop[cind]),
    nrootPOPULATION = (cpop[cind]) .^ (1.0/3)/5,

    CONFIRMEDpc = cconfirmed[cind] ./ cpop[cind] .* 1e5,
    CONFIRMED = cconfirmed[cind],
    log10CONFIRMEDpc = log10.(cconfirmed[cind] ./ cpop[cind] .* 1e5),

    dCONFIRMEDpc = dCONFIRMEDpc,
    dCONFIRMED = dcconfirmed[cind],
    log10dCONFIRMEDpc = log10.(dCONFIRMEDpc),

    dCONFIRMEDpcinst = dCONFIRMEDpcinst,
    dCONFIRMEDinst = dcconfirmedinst[cind],
    log10dCONFIRMEDpcinst = log10.(dCONFIRMEDpcinst),

    DEATHpc = cdeath[cind] ./ cpop[cind] .* 1e5,
    DEATH = cdeath[cind],
    log10DEATHpc = log10.(cdeath[cind] ./ cpop[cind] .* 1e5),

    dDEATHpc = dDEATHpc,
    dDEATH = dcdeath[cind],
    log10dDEATHpc = log10.(dDEATHpc),

    dDEATHpcinst = dDEATHpcinst,
    dDEATHinst = dcdeathinst[cind],
    log10dDEATHpcinst = log10.(dDEATHpcinst),
    );

function pmapconfirmed()
    pmapsg = PlotlyJS.scattergeo(; locationmode="USA-states",
        lat = df[!, :LAT],
        lon = df[!, :LON],
        marker_size=df[!, :nrootPOPULATION],
        marker_color=df[!, :log10CONFIRMEDpc],
        marker_colorscale="Jet",
        marker_showscale = true,
        marker_cmin = 2,
        marker_cmax = 4,

        hoverinfo="text",
        hovertext=
            "County: " .* county[cind] .* ", " .* state[cind] .* "
" .*
            "Population: " .* string.(df[!, :POPULATION]) .* "
" .*
            "Confirmed Cases (actual): " .* string.(df[!, :CONFIRMED]) .* "
" .*
            "Confirmed Cases per 100k: " .* string.(round.(df[!, :CONFIRMEDpc], digits=2)) .* "
" .*
            "Cumulative Deaths (actual): " .* string.(df[!, :DEATH]) .* "
" .*
            "Cumulative Deaths per 100k: " .* string.(round.(df[!, :DEATHpc], digits=2)) .* "
" .*
            "% infected with COVID19: " .* string.(round.(df[!, :CONFIRMED] ./ df[!, :POPULATION] .* 100, digits=2)) .* "
" .*
            "% death rate: " .* string.(round.(df[!, :DEATH] ./ df[!, :CONFIRMED] .*100, digits=2))

        );

    geo = attr(scope="usa",
        projection_type="ablbers usa",
        showland=true,
        landcolor="rgb(220, 220, 220)",
        subunitwidth=1,
        countrywidth=1,
        subunitcolor="rgb(255,255,255)",
        countrycolor="rgb(255,255,255)",
        );

    title = attr(text = "Cumulative Confirmed Cases per 100k " * string(t[ti])[1:10] * " (log10)
https://truedichotomy.github.io/covid19_public/
", yref = "paper", y=1.0)

    layout = Layout(; title=title, showlegend=false, geo=geo)
    PlotlyJS.plot(pmapsg, layout)
end

function pmapdconfirmed()
    pmapsg = PlotlyJS.scattergeo(; locationmode="USA-states",
        lat = df[!, :LAT],
        lon = df[!, :LON],
        #marker_size=df[!, :logPOPULATION]/1.3,
        marker_size=df[!, :nrootPOPULATION],
        marker_color=df[!, :log10dCONFIRMEDpc],
        marker_colorscale="Jet",
        marker_showscale = true,
        marker_cmin = -0.5,
        marker_cmax = 2.5,
        hoverinfo="text",
        #hovertext=string.(df[!, :dCONFIRMED]) .* " " .* county[cind],
        hovertext=
            "County: " .* county[cind] .* ", " .* state[cind] .* "
" .*
            "Population: " .* string.(df[!, :POPULATION]) .* "
" .*
            "Confirmed Cases (actual): " .* string.(df[!, :CONFIRMED]) .* "
" .*
            "Confirmed Cases per 100k: " .* string.(round.(df[!, :CONFIRMEDpc], digits=2)) .* "
" .*
            "Daily New Cases (actual): " .* string.(df[!, :dCONFIRMED]) .* "
" .*
            "Daily New Cases per 100k: " .* string.(round.(df[!, :dCONFIRMEDpc], digits=2)) .* "
" .*
            "% infected with COVID19: " .* string.(round.(df[!, :CONFIRMED] ./ df[!, :POPULATION] .* 100, digits=2)) .* "
" .*
            "% death rate: " .* string.(round.(df[!, :DEATH] ./ df[!, :CONFIRMED] .*100, digits=2))
        );

    geo = attr(scope="usa",
        projection_type="ablbers usa",
        showland=true,
        landcolor="rgb(220, 220, 220)",
        subunitwidth=1,
        countrywidth=1,
        subunitcolor="rgb(255,255,255)",
        countrycolor="rgb(255,255,255)",
        );

    title = attr(text = "Daily New Cases per 100k " * string(t[ti])[1:10] * " (log10, 7 day avg.)
https://truedichotomy.github.io/covid19_public/
", yref = "paper", y=1.0)

    layout = Layout(; title=title, showlegend=false, geo=geo)
    PlotlyJS.plot(pmapsg, layout)
end

function pmapdconfirmedinst()
    pmapsg = PlotlyJS.scattergeo(; locationmode="USA-states",
        lat = df[!, :LAT],
        lon = df[!, :LON],
        marker_size=df[!, :nrootPOPULATION],
        #marker_size=df[!, :logPOPULATION]/1.3,
        marker_color=df[!, :log10dCONFIRMEDpcinst],
        marker_colorscale="Jet",
        marker_showscale = true,
        marker_cmin = -0.5,
        marker_cmax = 2.5,
        hoverinfo="text",
        #hovertext=string.(df[!, :dCONFIRMED]) .* " " .* county[cind],
        hovertext=
            "County: " .* county[cind] .* ", " .* state[cind] .* "
" .*
            "Population: " .* string.(df[!, :POPULATION]) .* "
" .*
            "Confirmed Cases (actual): " .* string.(df[!, :CONFIRMED]) .* "
" .*
            "Confirmed Cases per 100k: " .* string.(round.(df[!, :CONFIRMEDpc], digits=2)) .* "
" .*
            "Daily New Cases (actual): " .* string.(df[!, :dCONFIRMEDinst]) .* "
" .*
            "Daily New Cases per 100k: " .* string.(round.(df[!, :dCONFIRMEDpcinst], digits=2)) .* "
" .*
            "% infected with COVID19: " .* string.(round.(df[!, :CONFIRMED] ./ df[!, :POPULATION] .* 100, digits=2)) .* "
" .*
            "% death rate: " .* string.(round.(df[!, :DEATH] ./ df[!, :CONFIRMED] .*100, digits=2))
        );

    geo = attr(scope="usa",
        projection_type="ablbers usa",
        showland=true,
        landcolor="rgb(220, 220, 220)",
        subunitwidth=1,
        countrywidth=1,
        subunitcolor="rgb(255,255,255)",
        countrycolor="rgb(255,255,255)",
        );

    title = attr(text = "Daily New Cases per 100k " * string(t[ti])[1:10] * " (log10)
https://truedichotomy.github.io/covid19_public/
", yref = "paper", y=1.0)

    layout = Layout(; title=title, showlegend=false, geo=geo)
    PlotlyJS.plot(pmapsg, layout)
end

function pmapdeath()
    pmapsg = PlotlyJS.scattergeo(; locationmode="USA-states",
        lat = df[!, :LAT],
        lon = df[!, :LON],
        marker_size=df[!, :nrootPOPULATION],
        #marker_size=df[!, :logPOPULATION]/1.3,
        marker_color=df[!, :log10DEATHpc],
        marker_colorscale="Jet",
        marker_showscale = true,
        marker_cmin = 0,
        marker_cmax = 2.5,
        hoverinfo="text",
        #hovertext=string.(df[!, :DEATH]) .* " " .* county[cind],
        hovertext=
            "County: " .* county[cind] .* ", " .* state[cind] .* "
" .*
            "Population: " .* string.(df[!, :POPULATION]) .* "
" .*
            "Confirmed Cases (actual): " .* string.(df[!, :CONFIRMED]) .* "
" .*
            "Confirmed Cases per 100k: " .* string.(round.(df[!, :CONFIRMEDpc], digits=2)) .* "
" .*
            "Cumulative Deaths (actual): " .* string.(df[!, :DEATH]) .* "
" .*
            "Cumulative Deaths per 100k: " .* string.(round.(df[!, :DEATHpc], digits=2)) .* "
" .*
            "% infected with COVID19: " .* string.(round.(df[!, :CONFIRMED] ./ df[!, :POPULATION] .* 100, digits=2)) .* "
" .*
            "% death rate: " .* string.(round.(df[!, :DEATH] ./ df[!, :CONFIRMED] .*100, digits=2))
        );

    geo = attr(scope="usa",
        projection_type="ablbers usa",
        showland=true,
        landcolor="rgb(220, 220, 220)",
        subunitwidth=1,
        countrywidth=1,
        subunitcolor="rgb(255,255,255)",
        countrycolor="rgb(255,255,255)",
        );

    title = attr(text = "Cumulative Death per 100k " * string(t[ti])[1:10] * " (log10)
https://truedichotomy.github.io/covid19_public/
", yref = "paper", y=1.0)

    layout = Layout(; title=title, showlegend=false, geo=geo)
    PlotlyJS.plot(pmapsg, layout)
end

function pmapddeath()
    pmapsg = PlotlyJS.scattergeo(; locationmode="USA-states",
        lat = df[!, :LAT],
        lon = df[!, :LON],
        #marker_size=df[!, :logPOPULATION]/1.3,
        marker_size=df[!, :nrootPOPULATION],
        marker_color=df[!, :log10dDEATHpc],
        marker_colorscale="Jet",
        marker_showscale = true,
        marker_cmin = -1.5,
        marker_cmax = 1,
        hoverinfo="text",
        #hovertext=string.(df[!, :dDEATH]) .* " " .* county[cind],
        hovertext=
            "County: " .* county[cind] .* ", " .* state[cind] .* "
" .*
            "Population: " .* string.(df[!, :POPULATION]) .* "
" .*
            "Confirmed Cases (actual): " .* string.(df[!, :DEATH]) .* "
" .*
            "Confirmed Cases per 100k: " .* string.(round.(df[!, :DEATHpc], digits=2)) .* "
" .*
            "Daily New Deaths (actual): " .* string.(df[!, :dDEATH]) .* "
" .*
            "Daily New Deaths per 100k: " .* string.(round.(df[!, :dDEATHpc], digits=2)) .* "
" .*
            "% infected with COVID19: " .* string.(round.(df[!, :CONFIRMED] ./ df[!, :POPULATION] .* 100, digits=2)) .* "
" .*
            "% death rate: " .* string.(round.(df[!, :DEATH] ./ df[!, :CONFIRMED] .*100, digits=2))
        );

    geo = attr(scope="usa",
        projection_type="ablbers usa",
        showland=true,
        landcolor="rgb(220, 220, 220)",
        subunitwidth=1,
        countrywidth=1,
        subunitcolor="rgb(255,255,255)",
        countrycolor="rgb(255,255,255)",
        );

    title = attr(text = "Daily New Deaths per 100k " * string(t[ti])[1:10] * " (log10, 7 day avg.)
https://truedichotomy.github.io/covid19_public/
", yref = "paper", y=1.0)

    layout = Layout(; title=title, showlegend=false, geo=geo)
    PlotlyJS.plot(pmapsg, layout)
end

function pmapddeathinst()
    pmapsg = PlotlyJS.scattergeo(; locationmode="USA-states",
        lat = df[!, :LAT],
        lon = df[!, :LON],
        #marker_size=df[!, :logPOPULATION]/1.3,
        marker_size=df[!, :nrootPOPULATION],

        marker_color=df[!, :log10dDEATHpcinst],
        marker_colorscale="Jet",
        marker_showscale = true,
        marker_cmin = -1.5,
        marker_cmax = 1,
        hoverinfo="text",
        #hovertext=string.(df[!, :dDEATH]) .* " " .* county[cind],
        hovertext=
            "County: " .* county[cind] .* ", " .* state[cind] .* "
" .*
            "Population: " .* string.(df[!, :POPULATION]) .* "
" .*
            "Confirmed Cases (actual): " .* string.(df[!, :DEATH]) .* "
" .*
            "Confirmed Cases per 100k: " .* string.(round.(df[!, :DEATHpc], digits=2)) .* "
" .*
            "Daily New Deaths (actual): " .* string.(df[!, :dDEATHinst]) .* "
" .*
            "Daily New Deaths per 100k: " .* string.(round.(df[!, :dDEATHpcinst], digits=2)) .* "
" .*
            "% infected with COVID19: " .* string.(round.(df[!, :CONFIRMED] ./ df[!, :POPULATION] .* 100, digits=2)) .* "
" .*
            "% death rate: " .* string.(round.(df[!, :DEATH] ./ df[!, :CONFIRMED] .*100, digits=2))
        );

    geo = attr(scope="usa",
        projection_type="ablbers usa",
        showland=true,
        landcolor="rgb(220, 220, 220)",
        subunitwidth=1,
        countrywidth=1,
        subunitcolor="rgb(255,255,255)",
        countrycolor="rgb(255,255,255)",
        );

    title = attr(text = "Daily New Deaths per 100k " * string(t[ti])[1:10] * " (log10)
https://truedichotomy.github.io/covid19_public/
", yref = "paper", y=1.0)

    layout = Layout(; title=title, showlegend=false, geo=geo)
    PlotlyJS.plot(pmapsg, layout)
end

#figoutdir = "/Volumes/GoogleDrive/My Drive/COVID19/";
local figoutdir = "../covid19_public/ts_maps/";

#PlotlyJS._js_path = "https://cdn.plot.ly/plotly-latest.min.js";

#PlotlyJS.savehtml(pmapconfirmed(), figoutdir * "covid19map_confirmed.html", :remote);
#PlotlyJS.savehtml(pmapdconfirmed(), figoutdir * "covid19map_delta_confirmed_7days.html", :remote);
#PlotlyJS.savehtml(pmapdconfirmedinst(), figoutdir * "covid19map_delta_confirmed_latest.html", :remote);
#PlotlyJS.savehtml(pmapdeath(), figoutdir * "covid19map_death.html", :remote);
#PlotlyJS.savehtml(pmapddeath(), figoutdir * "covid19map_delta_death_7days.html", :remote);
#PlotlyJS.savehtml(pmapddeathinst(), figoutdir * "covid19map_delta_death_latest.html", :remote);
#end

#p = pmapdconfirmed()
local mio = open(figoutdir * "covid19map_delta_confirmed_7days_" * string(t[ti]) * ".jpg", "w")
PlotlyJS.savefig(mio, pmapdconfirmed(); format="jpeg", width=1280, height=720, scale = 3);
close(mio)

local mio = open(figoutdir * "covid19map_confirmed_7days_" * string(t[ti]) * ".jpg", "w")
PlotlyJS.savefig(mio, pmapconfirmed(); format="jpeg", width=1280, height=720, scale = 3);
close(mio)

local mio = open(figoutdir * "covid19map_death_7days_" * string(t[ti]) * ".jpg", "w")
PlotlyJS.savefig(mio, pmapdeath(); format="jpeg", width=1280, height=720, scale = 3);
close(mio)

end


text:
export euler, is_cut_edge

function euler(G::SimpleDigraph{T}, u::T, v::T) where {T}
    notrail = T[]
    #check in_degrees and out_degrees of start and end vertex first
    if u == v
        if in_deg(G,u) != out_deg(G,u)
            return notrail
        end
    else
        if out_deg(G,u) - in_deg(G,u) != 1 ||
            in_deg(G,v) - out_deg(G,v) != 1
            return notrail
        end
    end

    #check if the undirected graph has an euler path
    simpleG = simplify(G)
    if length(euler(simpleG,u,v)) == 0
        return notrail
    end
    #check if connected
    #check if in_deg == out_deg for all other vertices

    GG = deepcopy(G)
    return euler_work!(GG, u)
end

euler(G::SimpleDigraph,u) = euler(G,u,u)
euler(G::SimpleDigraph) = euler(G,first(G.V))

# determine if an edge in a directed graph is a cut edge
function is_cut_edge(G::SimpleDigraph{T}, u::T, v::T) where {T}
    if !has(G,u,v)
        error("No such edge in this graph")
    end

    SimpleGraphs.delete!(G,u,v)
    P = find_path(G,v,u)
    if (length(P) == 0)
        add!(G,u,v)
        return true
    else
        add!(G,u,v)
        return false
    end
end

# helper function to determine if there is euler path
# function euler_work!(G::SimpleDigraph{T}, u::T) where {T}
function euler_work!(G::SimpleDigraph{T}, u::T) where {T}
    trail = T[]
    while true
        Nu = out_neighbors(G,u)
        if (length(Nu)) == 0
            push!(trail, u)
            return trail
        end
        if length(Nu) == 1
            v = Nu[1]
            push!(trail,u)
            SimpleGraphs.delete!(G,u,v)
            u = v
        else
            for w in Nu
                if !is_cut_edge(G,u,w)
                    SimpleGraphs.delete!(G,u,w)
                    push!(trail,u)
                    u = w
                    break
                end
            end
        end
    end
    error("Not gonna happen")
end


text:
function create_export_context(vendor::AbstractVendor, product::AbstractProduct, data::AbstractData, keystr::String)
    # The controller is responsible for instantiating the dialog GUI,
    # handling user input and tiggering the file importer once the user
    # confirms the selection of a particular file.
    control = FileDialogControl(true)
    # The model store the path that the user has confirmed, as well as the
    # path that the user has selected but not necessarily confirmed.
    model = FileDialogModel(Path(pwd(),""), Path(pwd(),""))
    # The display properties include the title of the file dialog and what
    # action string is to be displayed on the button.
    properties = FileDialogDisplayProperties(caption = "Save File###"*keystr, action = "Save###"*keystr, width = Cfloat(640), height = Cfloat(395))
    # The importer is triggered once the user has confirmed the selection of a
    # particular path.
    store_csv = CSVExporter(false, CSVSchema(vendor, product, data))
    # Construct a context and assign roles to complete the requisite export interactions.
    ExportContext(control, model, properties, store_csv)
end

function (store::CSVExporter{<:AbstractSchema})(path::Path, eda::ElectrodermalData, li::LabelledIntervals)
    @show "Storing Electrodermal Data"
    df = DataFrame(Data=String[])
    timestamp = get_timestamp(eda)
    frequency = get_sampling_frequency(eda)
    raw_eda = get_unprocessed_eda(eda)
    eda = get_eda(eda)
    # TODO The exporting of filtered vs unfiltered data needs to be refactored.
    if norm(raw_eda - eda) < 1e-15
        df = DataFrame(Raw = vcat(timestamp, frequency, raw_eda), Unfiltered = vcat(timestamp, frequency, eda))
    else
        df = DataFrame(Raw = vcat(timestamp, frequency, raw_eda), Filtered = vcat(timestamp, frequency, eda))
    end
    #df = DataFrame(Raw = vcat(timestamp, frequency, raw_eda))
    labelled_intervals = get_labelled_intervals(li)
    for (key, labelled_interval) in pairs(labelled_intervals)
        nested_interval = labelled_interval.nested_interval
        start = get_start(nested_interval)
        stop = get_stop(nested_interval)
        indicator = [ start <= i <= stop ? 1.0 : 0.0 for i = 1:length(raw_eda)]
        df = hcat(df, DataFrame([vcat(timestamp, frequency, indicator)], [Symbol(key)]))
    end
    store(path, df)
end


text:
#==============================================================================#
# Options to give in parameter file
#==============================================================================#

#------------------------------------------------------------------------------#
# Files
#------------------------------------------------------------------------------#

# set genotype, map, and pedigree file names
#genotype_file = "ALGP2_Cycle_1_5_Geno_BLUPF90.txt_clean"
genotype_file = "genotypes.txt"
map_file      = "ChrInfoBLUPF90.txt"
ped_file      = "dam_1_5.ped"

#------------------------------------------------------------------------------#
# General Options
#------------------------------------------------------------------------------#

#------------------------------------------------------------------------------#
# Map File Options
#------------------------------------------------------------------------------#

# Map file options
map_file_delimiter = ' '      # [char] needs to be a character with ' not "
map_file_header    = false    # [logical] is there a header on the map file
map_file_start_row = 1        # [Int] Line to start reading from

#------------------------------------------------------------------------------#
# Genotype File Options
#------------------------------------------------------------------------------#

# SNP format:
# "dense" =      ID1     0501102020102052
# "space" =      ID1     0 5 0 1 1 0 2 0 2 0 1 0 2 0 5 2
SNP_format = "dense"

# headers on map file or genotype file
genotype_file_delimiter   = ' '        # [char] set what separates in the genotype file
genotype_file_start_row   = 1          # [Int] row to start reading
genotype_missing_value    = 5          # [Int] set missing value for genotypes (will take anything not 0,1,2)

# NOT USED
n_skip_genotype_file      = 0          # [Int] number of lines to skip in the genotype file
genotype_header           = false      # [logical] ask if there is a header
genotype_start_col        = 2          # [Int]

# number of SNPs to sample
#n_sample_M = 0              # Select all SNPs
n_sample_M = 50000                  # [Int] number to sample from genotype file

#------------------------------------------------------------------------------#
# QC Options
#------------------------------------------------------------------------------#

# Call Rates (minimum proportion of non-missing genotypes for each)
call_rate_animals  = 0.90
call_rate_SNPs = 0.90

# MAF threshold
MAF_threshold = 0.05

# Exclude IDs
exclude_animals = []
exclude_SNPs    = []
exclude_Chr     = []
exclude_SNP     = []

#----------------------------------------#
# PCA Options
#----------------------------------------#

# set number of Principal components to use and plot bivariate scatterplots
n_PCs = 3

#----------------------------------------#
# Save Options
#----------------------------------------#

# save clean SNPs
save_clean = "yes"

# Save Raw Relationship Matrices
saveG = "no"
saveA = "no"
saveH = "no"

# Save Inverses of relationship matrices
saveGinv = "no"
saveAinv = "no"
saveHinv = "no"


text:
const test_bounds = [Bounds.Ellipsoid, Bounds.MultiEllipsoid]
const test_props = [Proposals.Uniform(), Proposals.RWalk(ratio=0.9, walks=75), Proposals.RStagger(ratio=0.9, walks=75, scale=0.5), Proposals.Slice(slices=10), Proposals.RSlice()]

@testset "$(nameof(bound)), $(nameof(typeof(proposal)))" for bound in test_bounds, proposal in test_props
    @testset "Correlated Gaussian Conjugate Prior - ndims=$D" for D in [2, 4, 8]
        if D == 8 && (proposal isa Proposals.RWalk || proposal isa Proposals.RStagger)
            # TODO evidence estimates are terrible for D=8
            continue
        end
        model, logz = Models.CorrelatedGaussian(D)

        sampler = Nested(D, 50D; bounds=bound, proposal=proposal)

        chain, state = sample(rng, model, sampler; dlogz=0.01)
        chain_res = sample(chain, Weights(vec(chain[:weights])), length(chain))
        # test posteriors
        vals = Array(chain_res)
        means = mean(vals, dims=1)
        tols = 3 .* std(vals, mean=means, dims=1) # 3-sigma
        μ = fill(2.0, D)
        Σ = fill(0.95, D, D)
        Σ[diagind(Σ)] .= 1
        expected = Σ * ((Σ + I) \ μ)
        @test all(@.(abs(means - expected) < tols))

        # logz (5-sigma)
        @test state.logz ≈ logz atol = 5state.logzerr
    end

    @testset "Gaussian Shells" begin
        model, logz = Models.GaussianShells()

        sampler = Nested(2, 1000; bounds=bound, proposal=proposal)

        chain, state = sample(rng, model, sampler; dlogz=0.01)

        # logz
        @test state.logz ≈ logz atol = 5state.logzerr
    end

    @testset "Gaussian Mixture Model" begin
        σ = 0.1
        μ1 = ones(2)
        μ2 = -ones(2)
        inv_σ = diagm(0 => fill(1 / σ^2, 2))

        function logl(x)
        dx1 = x .- μ1
        dx2 = x .- μ2
        f1 = -dx1' * (inv_σ * dx1) / 2
        f2 = -dx2' * (inv_σ * dx2) / 2
        return logaddexp(f1, f2)
    end

        prior(X) = 10 .* X .- 5
        model = NestedModel(logl, prior)

        analytic_logz = log(4π * σ^2 / 100)

        spl = Nested(2, 1000, bounds=bound, proposal=proposal)
        chain, state = sample(rng, model, spl; dlogz=0.2)
        chain_res = sample(chain, Weights(vec(chain[:weights])), length(chain))

        diff = state.logz - analytic_logz
        atol = 5state.logzerr
        if diff > atol
        @warn "logz estimate is poor" bound proposal error = diff tolerance = atol
    end

        @test state.logz ≈ analytic_logz atol = atol # within 3σ
        xmodes = sort!(findpeaks(chain_res[:, 1, 1])[1:2])
        @test xmodes[1] ≈ -1 atol = 2σ
        @test xmodes[2] ≈ 1 atol = 2σ
        ymodes = sort!(findpeaks(chain_res[:, 2, 1])[1:2])
        @test ymodes[1] ≈ -1 atol = 2σ
        @test ymodes[2] ≈ 1 atol = 2σ
    end
end


text:
import Base: getindex, length, eltype
export primetables, primelookup, primetableinfo, primetablefilename
export primetabletype, primesievetype

# A single table of π(x) with constant increment between values of x
immutable PrimeTable
    data::Array{Int128,1}  # values of π(x)
    incr::Int128           # increment to x for successive elements
    maxn::Int128           # largest value in data. (The last element)
    expn::Int              # log10(incr), an integer
end

const Zero = zero(Int128)

length(t::PrimeTable) = length(t.data)
getindex(t::PrimeTable,i) = (t.data)[i]
eltype(t::PrimeTable) = eltype(t.data)
getindex(t::Array{PrimeTable}, i::Int,j::Int) = t[i][j]
primetabletype() = eltype(primetables[1])
primesievetype() = UInt64

#  Return a list (pi-tab, min, rem), where `pi-tab' is
#  the value of prime pi function at argument `min',
#  and `rem'= `x'-`min'"
function piandrem(t::PrimeTable, x)
    (q,rem) = divrem(x,t.incr)
   return ((t.data)[q], t.incr * q, rem)
end

# Find the table with finest increments such that
# x falls on or between incrments and look up value for x
function piandrem{T<:Real}(x::T)
    j = 0
    @inbounds for i in 1:length(primetables) x <= primetables[i].maxn && (j = i; break) end
    j == 0 && error("x is too large!")
    piandrem(primetables[j],x)
end

function primelookup(x)
    x = conv128(x)
    j = 0
    @inbounds for i in 1:length(primetables) x < primetables[i].maxn && (j = i ; break) end
    j == 0 && error("x is too large!")
   (j,piandrem(primetables[j],x))
end

# Look up prime pi in table, compute remaining primes
function _countprimes(stop)
    stop < 10 && return  convert(Int128,ntcountprimes(stop))
    (count,i,rem) = piandrem(Int128(stop))
    return rem == Zero ? convert(Int128,count) :
    convert(Int128,count+ntcountprimes(i,i+rem))
end

function _countprimes(start,stop)
    local n1,n2,count1,count2
    if start < 10
        n1 = ntcountprimes(start)
        count1 = zero(n1)
    else
        (count1,i1,rem1) = piandrem(Int128(start))
        n1 = rem1 == Zero ? Zero : ntcountprimes(i1,i1+rem1)
    end
    if stop < 10
        n2 = ntcountprimes(stop)
        count2 = zero(n2)
    else
        (count2,i2,rem2) = piandrem(Int128(stop))
        n2 = rem2 == Zero ? Zero : ntcountprimes(i2,i2+rem2)
    end
    convert(Int128, count2 - count1 + n2 - n1)
end

const _primetablefilename = Pkg.dir("PrimeSieve") * "/data/primetables128bin.dat"
primetablefilename() = _primetablefilename

# Read the tables from a binary data file.  First Int is number of
# tables.  Second Int is number of elements in first table.  Next come
# all the elements in the table, then the number of elements for the
# next table, etc.
function _readbintables()
    fn = primetablefilename()
    if stat(fn).inode == 0
        error("Can't find file containing prime tables, $fn\n"
              * "Maybe your package installation is corrupt.")
    end
    mystream = open(fn)
    buf = zeros(Int,1)
    read!(mystream,buf)
    numtables = buf[1]
    bintables = Array{Array{Int128,1}}(0)
    for itab in 1:numtables
        read!(mystream,buf)
        numprimes = buf[1]
        a = Array{Int128}(numprimes)
        read!(mystream,a)
        push!(bintables,a)
    end
    close(mystream)
    return bintables
end

# Read binary tables and make table data structures.
function loadprimetables()
    bintables = _readbintables()
    tables = Array{PrimeTable}(length(bintables))
    base = Int128(10)
    for i in 1:length(bintables)
        data = bintables[i]
        expn = i
        incr = base^expn
        maxn = incr * length(data)
        pt = PrimeTable(data,incr,maxn,expn)
        tables[i] = pt
    end
    tables
end

function primetableinfo()
    dtype = eltype(primetables[1])
    println("Tables of π(x). element type: $dtype. Listed are: table number, increment in x (and first value of x),")
    println("number of entries in the table, largest x in table.\n")
    println("table  incr    tab len  max x")
    for i in 1:length(primetables)
        t = primetables[i]
        l = length(t)
        ip = rpad("$i",6)
        incr = rpad("10^$i",7)
        ll = int(log10(l))
        len = rpad("10^$ll",8)
        maxn = "10^$(ll+i)"
        println("$ip $incr $len $maxn")
    end
end

# Read the tables now.
const primetables = loadprimetables();


text:
using NLPModels
using CUTEst

function cutest2spg(probname)
    prob = CUTEstModel(probname)

    f(x) = obj(prob, x)
    function g!(g, x)
        g .= grad(prob, x)
    end
    x0 = copy(prob.meta.x0)
    lower = copy(prob.meta.lvar)
    upper = copy(prob.meta.uvar)
    return f, g!, x0, lower, upper, prob
end


text:
using Printf
using Plots

n =100
del, E, F ,B = init(n)

x = collect(0:n+1)*pi/(n+1)
U0 = 2.0*sin.(x)
x =x[2:n+1]
U0 = U0[2:n+1]

Ut = zeros(n)
Ux= 2.0*cos.(x)

w1 = del
w2 = del.*(E.+1.0)
w0 = hcat(w1,w2)

function runSim(w0,U0,Ut,Ux,tEnd)

tStart = 0;
t1 =tStart;
deltaT = 0.005;
#tEnd = 1.5;
w = zeros(n,2)

while t1 < tEnd

w, dt = FVMIBL(w0,U0,Ut,Ux);
ww = w
del = ww[:,1]
E = (ww[:,2]./ww[:,1]) .- 1.0
ncell = length(del)
Csep = zeros(ncell-1)
#xtrunc = x[2:end-1]
#for i = 2:ncell-1
    #Csep[i] = (del[i] - del[i-1])/(del[i+1] - del[i])
    #if Csep[i] > 5.
    #    println("singularity (separation) detected at t=$t1, x=$(x[i]./pi), Csep=$(Csep[i])")
    #end
#end
display(plot(x/pi,[del, E], xticks = 0:0.1:1, layout=(2,1), legend = false))
sleep(0.05)
w0 = w;

@printf("Time :%1.10f , Time step size %1.10f \n", t1, dt);
#map(x -> @sprintf("Seperation occure at :%1.10f  \n",x), xtrunc[seperation]./pi);

t1 = t1+dt
end

return del, E
end

max_stars_repo_name: desanga/IBMFVM

max_issues_repo_name: desanga/IBMFVM

max_forks_repo_name: desanga/IBMFVM


text:
workspace()

include("file_io/read_key.jl")
include("file_io/write_key.jl")

include("file_io/read_input.jl")
include("file_io/read_output_generation.jl")
include("file_io/write_output.jl")

include("characteristic_matrix_ellipse.jl")
include("characteristic_matrix_ellipsoid.jl")
include("inverse_characteristic_matrix_ellipsoid.jl")
include("rotation_matrix.jl")
include("inverse_rotation_matrix.jl")

include("generate_cell_lists.jl")
include("overlap_cuboid_binary.jl")
include("axis_aligned_bounding_box.jl")
include("intersect_box_box.jl")

foo = @__FILE__
@eval @everywhere f = $foo
#@everywhere println(f)
@everywhere (program_file_dir, program_file_name) = splitdir(f)
@everywhere include(joinpath(program_file_dir, "characterize_S2.jl"))
@everywhere include(joinpath(program_file_dir, "characterize_S3.jl"))
@everywhere include(joinpath(program_file_dir, "characterize_zeta2.jl"))
@everywhere include(joinpath(program_file_dir, "signed_distance_mod.jl"))
@everywhere include(joinpath(program_file_dir, "position_mod.jl"))
@everywhere include(joinpath(program_file_dir, "rotation_matrix.jl"))
@everywhere include(joinpath(program_file_dir, "generate_random_unit_quaternion.jl"))
@everywhere include(joinpath(program_file_dir, "quaternion_mult.jl"))
@everywhere include(joinpath(program_file_dir, "rotate.jl"))

function run_characterization()
	# Inititalization of random number generation device.
	random_seed::Int64 = convert(Int64, time_ns())
	srand(random_seed)

	# Start time.
	t_start_ns::Int64 = convert(Int64, time_ns())

	# Change current folder to the folder where this script lies.
	(program_file_dir::String, program_file_name::String) = splitdir(PROGRAM_FILE)
	program_file_dir = abspath(program_file_dir)
	cd(program_file_dir)

	# Assert that input is file and store path.
	input_file_path::String = ""
	if isfile(ARGS[1])
		input_file_path = ARGS[1]
	else
		println("No input file specified or specified input file does not exist. Aborting.")
		return nothing
	end

	# Read characerization input from file.
	println(join(("Reading characerization input from file ", input_file_path, "...")))
	(	output_generation_path::String,
		number_of_samples::Int64,
		d2::Array{Float64, 1},
		d3::Array{Float64, 1},
		theta3::Array{Float64, 1},
		number_of_cells_x::Int64,
		number_of_cells_y::Int64,
		number_of_cells_z::Int64,
		output_file_path::String) = read_input(input_file_path)

	# Read generation output from file.
	println(join(("Reading generation output from file ", output_generation_path, "...")))
	(	particle_type::String,
		R::Array{Float64, 2},
		Lx::Float64,
		Ly::Float64,
		Lz::Float64,
		phi::Float64,
		X::Array{Float64, 1},
		Y::Array{Float64, 1},
		Z::Array{Float64, 1},
		Q0::Array{Float64, 1},
		Q1::Array{Float64, 1},
		Q2::Array{Float64, 1},
		Q3::Array{Float64, 1},
		execution_time_generation::Float64) = read_output_generation(output_generation_path)
	number_of_particles::Int64 = length(X)

	# Characteristic/rotation matrix entries.
	A11::Array{Float64, 1} = zeros(number_of_particles)
	A12::Array{Float64, 1} = zeros(number_of_particles)
	A13::Array{Float64, 1} = zeros(number_of_particles)
	A21::Array{Float64, 1} = zeros(number_of_particles)
	A22::Array{Float64, 1} = zeros(number_of_particles)
	A23::Array{Float64, 1} = zeros(number_of_particles)
	A31::Array{Float64, 1} = zeros(number_of_particles)
	A32::Array{Float64, 1} = zeros(number_of_particles)
	A33::Array{Float64, 1} = zeros(number_of_particles)

	a11::Float64 = 0.0
	a12::Float64 = 0.0
	a13::Float64 = 0.0
	a21::Float64 = 0.0
	a22::Float64 = 0.0
	a23::Float64 = 0.0
	a31::Float64 = 0.0
	a32::Float64 = 0.0
	a33::Float64 = 0.0
	if particle_type != "sphere"
		for current_particle = 1:number_of_particles
			(a11, a12, a13, a21, a22, a23, a31, a32, a33) = inverse_rotation_matrix(Q0[current_particle], Q1[current_particle], Q2[current_particle], Q3[current_particle])
			A11[current_particle] = a11
			A12[current_particle] = a12
			A13[current_particle] = a13
			A21[current_particle] = a21
			A22[current_particle] = a22
			A23[current_particle] = a23
			A31[current_particle] = a31
			A32[current_particle] = a32
			A33[current_particle] = a33
		end
	end

	# Create cell lists.
	cell_overlap::Float64 = 0.0
	cell_lists::Array{Array{Int64, 1}, 3} = generate_cell_lists(particle_type,
																R,
																Lx,
																Ly,
																Lz,
																X,
																Y,
																Z,
																Q0,
																Q1,
																Q2,
																Q3,
																number_of_cells_x,
																number_of_cells_y,
																number_of_cells_z,
																cell_overlap)

	# Run characterization.
	number_of_workers::Int64 = nworkers() # This is determined by the the '-p' input flag to Julia.
	number_of_samples_per_worker::Array{Int64, 1} = convert(Array{Int64, 1}, floor(number_of_samples / number_of_workers) * ones(number_of_workers))
	number_of_samples_remaining::Int64 = number_of_samples - sum(number_of_samples_per_worker)
	number_of_samples_per_worker[1:number_of_samples_remaining] += 1
	#println(number_of_samples_per_worker)

	# Compute S2.
	S2::Array{Float64, 1} = zeros(size(d2))
#	S2 = @parallel (+) for current_worker = 1:number_of_workers
#		characterize_S2(
#			particle_type,
#			R,
#			Lx,
#			Ly,
#			Lz,
#			X,
#			Y,
#			Z,
#			Q0,
#			Q1,
#			Q2,
#			Q3,
#			A11,
#			A12,
#			A13,
#			A21,
#			A22,
#			A23,
#			A31,
#			A32,
#			A33,
#			number_of_samples_per_worker[current_worker],
#			d2,
#			cell_lists)
#	end
#	S2 /= convert(Float64, number_of_samples)

	# Compute S3.
	number_of_distances::Int64 = length(d3)
	number_of_angles::Int64 = length(theta3)
	S3::Array{Float64, 3} = zeros(number_of_distances, number_of_distances, convert(Int64, number_of_angles/2 + 1))
	S3 = @parallel (+) for current_worker = 1:number_of_workers
		characterize_S3( particle_type, R, Lx, Ly, Lz, X, Y, Z, Q0, Q1, Q2, Q3,
			A11, A12, A13, A21, A22, A23, A31, A32, A33,
			number_of_samples_per_worker[current_worker], d3, theta3, cell_lists)
	end
	S3 /= convert(Float64, number_of_samples * number_of_angles)

	# Compute zeta2.
	S1::Float64 = 1.0 - phi
	integrand::Array{Float64, 3} = zeros(number_of_distances, number_of_distances, convert(Int64, number_of_angles/2 + 1))
    for current_r1 = 1:number_of_distances
        for current_r2 = 1:number_of_distances
            for current_angle = 1:convert(Int64, number_of_angles/2 + 1)
                integrand[current_r1, current_r2, current_angle] = 1.0 / d3[current_r1] * 1.0 / d3[current_r2] *
					0.5 * sin(theta3[current_angle]) * (3.0 * cos(theta3[current_angle])^2 - 1.0) *
                    ( S3[current_r1, current_r2, current_angle] - S3[current_r1, current_r1, 1] * S3[current_r2, current_r2, 1] / S1 )
            end
        end
    end
	integral::Float64 = 0.0
	dr1::Float64 = 0.0
	dr2::Float64 = 0.0
	dtheta::Float64 = 0.0
	for current_r1 = 1:number_of_distances
		if current_r1 == 1
			dr1 = 0.5 * (d3[2] - d3[1])
		elseif current_r1 == number_of_distances
			dr1 = 0.5 * (d3[number_of_distances] - d3[number_of_distances - 1])
		else
			dr1 = 0.5 * (d3[current_r1 + 1] - d3[current_r1 - 1])
		end

        for current_r2 = 1:number_of_distances
			if current_r2 == 1
				dr2 = 0.5 * (d3[2] - d3[1])
			elseif current_r2 == number_of_distances
				dr2 = 0.5 * (d3[number_of_distances] - d3[number_of_distances - 1])
			else
				dr2 = 0.5 * (d3[current_r2 + 1] - d3[current_r2 - 1])
			end

            for current_angle = 1:convert(Int64, number_of_angles/2 + 1)
				if current_angle == 1
					dtheta = 0.5 * (theta3[2] - theta3[1])
				elseif current_angle == number_of_distances
					dtheta = 0.5 * (theta3[convert(Int64, number_of_angles/2 + 1)] - theta3[convert(Int64, number_of_angles/2 + 1) - 1])
				else
					dtheta = 0.5 * (theta3[current_angle + 1] - theta3[current_angle - 1])
				end

				integral += integrand[current_r1, current_r2, current_angle] * dr1 * dr2 * dtheta
			end
		end
	end
    zeta2::Float64 = 1.0 - 9.0/(2.0 * phi * (1.0 - phi)) * integral
	println(zeta2)

	# Kill all workers.
	#rmprocs(workers(); waitfor = typemax(Int))

	t_finish_ns::Int64 = convert(Int64, time_ns())
	t_exec::Float64 = convert(Float64, t_finish_ns - t_start_ns) / 1e9

	# Write output.
	write_output(
		output_file_path,
		number_of_samples,
		d2,
		S2,
		d3,
		theta3,
		S3,
		zeta2,
		t_exec)
	println(join(("Output written to ", output_file_path, ".")))
	println("Finished.")

	nothing
end

run_characterization()


text:
using Push
using Base.Test

cfg_path = joinpath(dirname(@__FILE__), "configuration/code.cfg")
cfg = Push.load_configuration(cfg_path)

# CODE.QUOTE
s = Push.run("(CODE.QUOTE X)", cfg)
@test s.code == {:X}
s = Push.run("(CODE.QUOTE -12)", cfg)
@test s.code == {-12}
s = Push.run("(CODE.QUOTE -5.1)", cfg)
@test s.code == {convert(Float32, -5.1)}
s = Push.run("(CODE.QUOTE TRUE)", cfg)
@test s.code == {true}
s = Push.run("(CODE.QUOTE CODE.QUOTE)", cfg)
@test s.code == {convert(Symbol, "CODE.QUOTE")}
s = Push.run("(CODE.QUOTE A CODE.QUOTE (X Y Z))", cfg)
@test s.code == {:A, {:X, :Y, :Z}}

# CODE.=
s = Push.run("(CODE.QUOTE X CODE.QUOTE X CODE.=)", cfg)
@test s.boolean == [true] && isempty(s.code)
s = Push.run("(CODE.QUOTE X CODE.QUOTE Y CODE.=)", cfg)
@test s.boolean == [false] && isempty(s.code)
s = Push.run("(CODE.QUOTE CODE.QUOTE CODE.QUOTE CODE.QUOTE CODE.=)", cfg)
@test s.boolean == [true] && isempty(s.code)

# CODE.APPEND
s = Push.run("(CODE.QUOTE (X Y) CODE.QUOTE (A B) CODE.APPEND)", cfg)
@test s.code == {{:A, :B, :X, :Y}} && isempty(s.name)
s = Push.run("(CODE.QUOTE (X Y) CODE.QUOTE Z CODE.APPEND)", cfg)
@test s.code == {{:Z, :X, :Y}} && isempty(s.name)
s = Push.run("(CODE.QUOTE A CODE.QUOTE B CODE.APPEND)", cfg)
@test s.code == {{:B, :A}} && isempty(s.name)
s = Push.run("(CODE.QUOTE A CODE.APPEND)", cfg)
@test s.code == {:A} && isempty(s.name)

# CODE.ATOM
s = Push.run("(CODE.ATOM)", cfg)
@test isempty(s.code) && isempty(s.boolean)
s = Push.run("(CODE.QUOTE 32 CODE.ATOM)", cfg)
@test isempty(s.code) && s.boolean == [true]
s = Push.run("(CODE.QUOTE 32.0 CODE.ATOM)", cfg)
@test isempty(s.code) && s.boolean == [true]
s = Push.run("(CODE.QUOTE X CODE.ATOM)", cfg)
@test isempty(s.code) && s.boolean == [true]
s = Push.run("(CODE.QUOTE FALSE CODE.ATOM)", cfg)
@test isempty(s.code) && s.boolean == [true]
s = Push.run("(CODE.QUOTE (FALSE) CODE.ATOM)", cfg)
@test isempty(s.code) && s.boolean == [false]

# CODE.CAR
#
#
# WHAT IS CAR OF AN EMPTY LIST?
#
#
s = Push.run("(CODE.QUOTE X CODE.CAR)", cfg)
@test s.code == {:X}
s = Push.run("(CODE.CAR)", cfg)
@test isempty(s.code)
s = Push.run("(CODE.QUOTE X CODE.QUOTE (1 2 3 4 5) CODE.CAR)", cfg)
@test s.code == {:X, 1}
s = Push.run("(CODE.QUOTE () CODE.CAR)", cfg)
@test s.code == {{}}

# CODE.CDR
s = Push.run("(CODE.CDR)", cfg)
@test isempty(s.code)
s = Push.run("(CODE.QUOTE X CODE.CDR)", cfg)
@test s.code == {{}}
s = Push.run("(CODE.QUOTE (X) CODE.CDR)", cfg)
@test s.code == {{}}
s = Push.run("(CODE.QUOTE (X Y Z) CODE.CDR)", cfg)
@test s.code == {{:Y, :Z}}

# CODE.CONS
s = Push.run("(CODE.QUOTE Z CODE.QUOTE (X Y) CODE.CONS)", cfg)
@test s.code == {{:Z, :X, :Y}}
s = Push.run("(CODE.QUOTE Z CODE.CONS)", cfg)
@test s.code == {:Z}
s = Push.run("(CODE.QUOTE (Z) CODE.QUOTE (X Y) CODE.CONS)", cfg)
@test s.code == {{{:Z}, :X, :Y}}
s = Push.run("(CODE.QUOTE B CODE.QUOTE A CODE.CONS)", cfg)
@test s.code == {{:B, :A}}
s = Push.run("(CODE.QUOTE (B) CODE.QUOTE A CODE.CONS)", cfg)
@test s.code == {{{:B}, :A}}

# CODE.DUP
s = Push.run("(CODE.DUP)", cfg)
@test isempty(s.code)
s = Push.run("(CODE.QUOTE X CODE.DUP)", cfg)
@test s.code == {:X, :X}

# CODE.EXTRACT
s = Push.run("(CODE.EXTRACT)", cfg)
@test isempty(s.code)
s = Push.run("(CODE.QUOTE (A B C D) CODE.EXTRACT)", cfg)
@test s.code == {{:A, :B, :C, :D}}
s = Push.run("(0 CODE.QUOTE (A B C D) CODE.EXTRACT)", cfg)
@test s.code == {{:A, :B, :C, :D}}
s = Push.run("(1 CODE.QUOTE (A B C D) CODE.EXTRACT)", cfg)
@test s.code == {:A}
s = Push.run("(2 CODE.QUOTE (A B C D) CODE.EXTRACT)", cfg)
@test s.code == {:B}
s = Push.run("(3 CODE.QUOTE (A B C D) CODE.EXTRACT)", cfg)
@test s.code == {:C}
s = Push.run("(4 CODE.QUOTE (A B C D) CODE.EXTRACT)", cfg)
@test s.code == {:D}
s = Push.run("(5 CODE.QUOTE (A B C D) CODE.EXTRACT)", cfg)
@test s.code == {{:A, :B, :C, :D}}
s = Push.run("(6 CODE.QUOTE (A B C D) CODE.EXTRACT)", cfg)
@test s.code == {:A}
s = Push.run("(-1 CODE.QUOTE (A B C D) CODE.EXTRACT)", cfg)
@test s.code == {:A}
s = Push.run("(-2 CODE.QUOTE (A B C D) CODE.EXTRACT)", cfg)
@test s.code == {:B}
s = Push.run("(-5 CODE.QUOTE (A B C D) CODE.EXTRACT)", cfg)
@test s.code == {{:A, :B, :C, :D}}
s = Push.run("(-6 CODE.QUOTE (A B C D) CODE.EXTRACT)", cfg)
@test s.code == {:A}
s = Push.run("(0 CODE.QUOTE A CODE.EXTRACT)", cfg)
@test s.code == {:A}
s = Push.run("(1 CODE.QUOTE ((A B C) (D E) (F)) CODE.EXTRACT)", cfg)
@test s.code == {{:A, :B, :C}}
s = Push.run("(2 CODE.QUOTE ((A B C) (D E) (F)) CODE.EXTRACT)", cfg)
@test s.code == {:A}
s = Push.run("(6 CODE.QUOTE ((A B C) (D E) (F)) CODE.EXTRACT)", cfg)
@test s.code == {:D}

# CODE.FLUSH
s = Push.run("(CODE.QUOTE 76 CODE.QUOTE (A B C D) CODE.FLUSH)", cfg)
@test s.code == {}
s = Push.run("(CODE.FLUSH)", cfg)
@test s.code == {}

# CODE.FROMBOOLEAN
s = Push.run("(A CODE.FROMBOOLEAN)", cfg)
@test s.name == {:A} && isempty(s.code)
s = Push.run("(FALSE FALSE TRUE CODE.FROMBOOLEAN)", cfg)
@test s.code == {true} && s.boolean == {false, false}
s = Push.run("(TRUE TRUE FALSE CODE.FROMBOOLEAN)", cfg)
@test s.code == {false} && s.boolean == {true, true}

# CODE.FROMINTEGER
s = Push.run("(7 CODE.FROMINTEGER)", cfg)
@test s.code == {7} && isempty(s.integer)
s = Push.run("(1 2 3 CODE.FROMINTEGER)", cfg)
@test s.integer == {1, 2} && s.code == {3}
s = Push.run("(3 2 1 CODE.FROMINTEGER)", cfg)
@test s.integer == {3, 2} && s.code == {1}

# CODE.FROMNAME
s = Push.run("(CODE.FROMNAME)", cfg)
@test isempty(s.code)
s = Push.run("(NAME.QUOTE ADD CODE.FROMNAME)", cfg)
@test s.code == {:ADD} && isempty(s.name)

# CODE.IF
s = Push.run("(CODE.QUOTE X CODE.QUOTE Y CODE.IF)", cfg)
@test s.code == {:X, :Y}
s = Push.run("(TRUE CODE.QUOTE X CODE.IF)", cfg)
@test s.code == {:X} && s.boolean == {true}
s = Push.run("(TRUE CODE.QUOTE (A B C) CODE.QUOTE (X Y Z) CODE.IF)", cfg)
@test s.name == {:A, :B, :C} && isempty(s.code) && isempty(s.boolean)
s = Push.run("(FALSE CODE.QUOTE (A B C) CODE.QUOTE (X Y Z) CODE.IF)", cfg)
@test s.name == {:X, :Y, :Z} && isempty(s.code) && isempty(s.boolean)

# CODE.LENGTH
s = Push.run("(CODE.LENGTH)", cfg)
@test isempty(s.integer)
s = Push.run("(CODE.QUOTE X CODE.LENGTH)", cfg)
@test s.integer == {1}
s = Push.run("(CODE.QUOTE (X) CODE.LENGTH)", cfg)
@test s.integer == {1}
s = Push.run("(CODE.QUOTE (X Y Z) CODE.LENGTH)", cfg)
@test s.integer == {3}
s = Push.run("(CODE.QUOTE ((A B C) D E) CODE.LENGTH)", cfg)
@test s.integer == {3}

# CODE.LIST
s = Push.run("(CODE.QUOTE A CODE.LIST)", cfg)
@test s.code == {:A}
s = Push.run("(CODE.QUOTE Z CODE.QUOTE A CODE.QUOTE B CODE.LIST)", cfg)
@test s.code == {:Z, {:B, :A}}
s = Push.run("(CODE.QUOTE Z CODE.QUOTE (1 2 3 4) CODE.QUOTE (A B C) CODE.LIST)", cfg)
@test s.code == {:Z, {{:A, :B, :C}, {1, 2, 3, 4}}}
s = Push.run("(CODE.QUOTE Z CODE.QUOTE (1 2 3 4) CODE.QUOTE A CODE.LIST)", cfg)
@test s.code == {:Z, {:A, {1, 2, 3, 4}}}

# CODE.NOOP
s = Push.run("(CODE.NOOP)", cfg)
@test isempty(s.code) && isempty(s.boolean) && isempty(s.integer) && isempty(s.float) && isempty(s.name) && isempty(s.exec)
s = Push.run("(1 CODE.NOOP 5)", cfg)
@test isempty(s.code) && isempty(s.boolean) && s.integer == {1, 5} && isempty(s.float) && isempty(s.name) && isempty(s.exec)

# CODE.NTH
s = Push.run("(CODE.NTH)", cfg)
@test isempty(s.code)
s = Push.run("(0 CODE.QUOTE X CODE.NTH)", cfg)
@test s.code == {:X}
s = Push.run("(919 CODE.QUOTE X CODE.NTH)", cfg)
@test s.code == {:X}
s = Push.run("(98 CODE.QUOTE () CODE.NTH)", cfg)
@test s.code == {{}}
s = Push.run("(CODE.QUOTE (1 2 3) CODE.NTH)", cfg)
@test s.code == {{1, 2, 3}}
s = Push.run("(0 CODE.QUOTE (1 2 3) CODE.NTH)", cfg)
@test s.code == {1}
s = Push.run("(1 CODE.QUOTE (1 2 3) CODE.NTH)", cfg)
@test s.code == {2}
s = Push.run("(2 CODE.QUOTE (1 2 3) CODE.NTH)", cfg)
@test s.code == {3}
s = Push.run("(3 CODE.QUOTE (1 2 3) CODE.NTH)", cfg)
@test s.code == {1}
s = Push.run("(-1 CODE.QUOTE (1 2 3) CODE.NTH)", cfg)
@test s.code == {2}

# CODE.NULL
s = Push.run("(CODE.NULL)", cfg)
@test isempty(s.code) && isempty(s.boolean)
s = Push.run("(CODE.QUOTE () CODE.NULL)", cfg)
@test isempty(s.code) && s.boolean == [true]
s = Push.run("(CODE.QUOTE (1 2 3) CODE.NULL)", cfg)
@test isempty(s.code) && s.boolean == [false]
s = Push.run("(CODE.QUOTE 1 CODE.NULL)", cfg)
@test isempty(s.code) && s.boolean == [false]

# CODE.POP
s = Push.run("(CODE.QUOTE X CODE.QUOTE Y CODE.POP)", cfg)
@test s.code == {:X}

# CODE.ROT
s = Push.run("(CODE.QUOTE A CODE.QUOTE B CODE.ROT)", cfg)
@test s.code == {:A, :B}
s = Push.run("(CODE.QUOTE Z CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.ROT)", cfg)
@test s.code == {:Z, :C, :A, :B}

# CODE.SIZE
s = Push.run("(CODE.QUOTE X CODE.SIZE)", cfg)
@test isempty(s.code) && s.integer == {1}
s = Push.run("(CODE.QUOTE (X) CODE.SIZE)", cfg)
@test isempty(s.code) && s.integer == {2}
s = Push.run("(CODE.QUOTE (X Y Z) CODE.SIZE)", cfg)
@test isempty(s.code) && s.integer == {4}
s = Push.run("(CODE.QUOTE ((X) (Y Z)) CODE.SIZE)", cfg)
@test isempty(s.code) && s.integer == {6}
s = Push.run("(CODE.SIZE)", cfg)
@test isempty(s.code) && isempty(s.integer)

# CODE.STACKDEPTH
s = Push.run("(CODE.STACKDEPTH)", cfg)
@test s.integer == {0}
s = Push.run("(CODE.QUOTE X CODE.STACKDEPTH)", cfg)
@test s.integer == {1}
s = Push.run("(CODE.QUOTE X CODE.QUOTE 3 CODE.STACKDEPTH)", cfg)
@test s.integer == {2}

# CODE.SWAP
s = Push.run("(CODE.QUOTE A CODE.SWAP)", cfg)
@test s.code == {:A}
s = Push.run("(CODE.QUOTE Z CODE.QUOTE A CODE.QUOTE B CODE.SWAP)", cfg)
@test s.code == {:Z, :B, :A}

# CODE.DO
s = Push.run("((CODE.QUOTE CODE.= CODE.DUP) CODE.DO)", cfg)
@test isempty(s.code) && s.boolean == {true}
s = Push.run("(CODE.QUOTE (X Y Z) CODE.DO)", cfg)
@test isempty(s.code) && s.name == {:X, :Y, :Z}
s = Push.run("(CODE.QUOTE (1 2 3 INTEGER.+ INTEGER.+) CODE.DO)", cfg)
@test isempty(s.code) && s.integer == {6}

# CODE.DO*
s = Push.run("((CODE.DUP CODE.QUOTE CODE.=) CODE.DO*)", cfg)
@test isempty(s.code) && isempty(s.boolean)
s = Push.run("(CODE.QUOTE (X Y Z) CODE.DO*)", cfg)
@test isempty(s.code) && s.name == {:X, :Y, :Z}
s = Push.run("(CODE.QUOTE (1 2 3 INTEGER.+ INTEGER.+) CODE.DO*)", cfg)
@test isempty(s.code) && s.integer == {6}
s = Push.run("(CODE.QUOTE (CODE.STACKDEPTH) CODE.DO*)", cfg)
@test isempty(s.code) && s.integer == {0}

# CODE.DO*RANGE
s = Push.run("(CODE.QUOTE (X) CODE.DO*RANGE)", cfg)
@test s.code == {{:X}} && isempty(s.name)
s = Push.run("(0 CODE.QUOTE (X) CODE.DO*RANGE)", cfg)
@test s.code == {{:X}} && isempty(s.name) && s.integer == {0}
s = Push.run("(0 4 CODE.QUOTE (X) CODE.DO*RANGE)", cfg)
@test s.name == {:X, :X, :X, :X, :X}

# CODE.DO*COUNT
s = Push.run("(CODE.QUOTE (2 INTEGER.*) CODE.DO*COUNT)", cfg)
@test s.code == {{2, convert(Symbol, "INTEGER.*")}} && isempty(s.integer)
s = Push.run("(6 CODE.QUOTE (2 INTEGER.*) CODE.DO*COUNT)", cfg)
@test s.integer == {0, 2, 4, 6, 8, 10}

# CODE.DO*TIMES
s = Push.run("(CODE.QUOTE (INTEGER.DUP) CODE.DO*TIMES)", cfg)
@test s.code == {{convert(Symbol, "INTEGER.DUP")}} && isempty(s.integer)
s = Push.run("(6 CODE.QUOTE (INTEGER.DUP) CODE.DO*TIMES)", cfg)
@test isempty(s.code) && isempty(s.integer)
s = Push.run("(5 CODE.QUOTE (X) CODE.DO*TIMES)", cfg)
@test s.name == {:X, :X, :X, :X, :X}
s = Push.run("(0 5 CODE.QUOTE (1 INTEGER.+) CODE.DO*TIMES)", cfg)
@test isempty(s.code) && s.integer == {5}
s = Push.run("(0 5 CODE.QUOTE (INTEGER.DUP 1 INTEGER.+) CODE.DO*TIMES)", cfg)
@test isempty(s.code) && s.integer == {0, 1, 2, 3, 4, 5}

# CODE.NTHCDR
s = Push.run("(CODE.QUOTE (A B C) CODE.NTHCDR)", cfg)
@test s.code == {{:A, :B, :C}}
s = Push.run("(0 CODE.QUOTE () CODE.NTHCDR)", cfg)
@test s.code == {{}} && isempty(s.integer)
s = Push.run("(10 CODE.QUOTE X CODE.NTHCDR)", cfg)
@test s.code == {{}} && isempty(s.integer)
s = Push.run("(0 CODE.QUOTE (A B C D E) CODE.NTHCDR)", cfg)
@test s.code == {{:B, :C, :D, :E}} && isempty(s.integer)
s = Push.run("(-1 CODE.QUOTE (A B C D E) CODE.NTHCDR)", cfg)
@test s.code == {{:C, :D, :E}} && isempty(s.integer)
s = Push.run("(3 CODE.QUOTE (A B C D E) CODE.NTHCDR)", cfg)
@test s.code == {{:E}} && isempty(s.integer)
s = Push.run("(-3 CODE.QUOTE (A B C D E) CODE.NTHCDR)", cfg)
@test s.code == {{:E}} && isempty(s.integer)
s = Push.run("(4 CODE.QUOTE (A B C D E) CODE.NTHCDR)", cfg)
@test s.code == {{}} && isempty(s.integer)
s = Push.run("(-4 CODE.QUOTE (A B C D E) CODE.NTHCDR)", cfg)
@test s.code == {{}} && isempty(s.integer)
s = Push.run("(5 CODE.QUOTE (A B C D E) CODE.NTHCDR)", cfg)
@test s.code == {{:B, :C, :D, :E}} && isempty(s.integer)
s = Push.run("(-5 CODE.QUOTE (A B C D E) CODE.NTHCDR)", cfg)
@test s.code == {{:B, :C, :D, :E}} && isempty(s.integer)

# CODE.SHOVE
s = Push.run("(CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.QUOTE D CODE.SHOVE)", cfg)
@test s.code == {:A, :B, :C, :D}
s = Push.run("(-98 CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.QUOTE D CODE.SHOVE)", cfg)
@test s.code == {:A, :B, :C, :D}
s = Push.run("(0 CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.QUOTE D CODE.SHOVE)", cfg)
@test s.code == {:A, :B, :C, :D}
s = Push.run("(1 CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.QUOTE D CODE.SHOVE)", cfg)
@test s.code == {:A, :B, :D, :C}
s = Push.run("(2 CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.QUOTE D CODE.SHOVE)", cfg)
@test s.code == {:A, :D, :B, :C}
s = Push.run("(3 CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.QUOTE D CODE.SHOVE)", cfg)
@test s.code == {:D, :A, :B, :C}
s = Push.run("(78 CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.QUOTE D CODE.SHOVE)", cfg)
@test s.code == {:D, :A, :B, :C}

# CODE.YANK
s = Push.run("(CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.QUOTE D CODE.YANK)", cfg)
@test s.code == {:A, :B, :C, :D}
s = Push.run("(0 CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.QUOTE D CODE.YANK)", cfg)
@test s.code == {:A, :B, :C, :D}
s = Push.run("(-871 CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.QUOTE D CODE.YANK)", cfg)
@test s.code == {:A, :B, :C, :D}
s = Push.run("(1 CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.QUOTE D CODE.YANK)", cfg)
@test s.code == {:A, :B, :D, :C}
s = Push.run("(2 CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.QUOTE D CODE.YANK)", cfg)
@test s.code == {:A, :C, :D, :B}
s = Push.run("(3 CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.QUOTE D CODE.YANK)", cfg)
@test s.code == {:B, :C, :D, :A}
s = Push.run("(78 CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.QUOTE D CODE.YANK)", cfg)
@test s.code == {:B, :C, :D, :A}

# CODE.YANKDUP
s = Push.run("(CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.QUOTE D CODE.YANKDUP)", cfg)
@test s.code == {:A, :B, :C, :D}
s = Push.run("(0 CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.QUOTE D CODE.YANKDUP)", cfg)
@test s.code == {:A, :B, :C, :D, :D}
s = Push.run("(-871 CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.QUOTE D CODE.YANKDUP)", cfg)
@test s.code == {:A, :B, :C, :D, :D}
s = Push.run("(1 CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.QUOTE D CODE.YANKDUP)", cfg)
@test s.code == {:A, :B, :C, :D, :C}
s = Push.run("(2 CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.QUOTE D CODE.YANKDUP)", cfg)
@test s.code == {:A, :B, :C, :D, :B}
s = Push.run("(3 CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.QUOTE D CODE.YANKDUP)", cfg)
@test s.code == {:A, :B, :C, :D, :A}
s = Push.run("(78 CODE.QUOTE A CODE.QUOTE B CODE.QUOTE C CODE.QUOTE D CODE.YANKDUP)", cfg)
@test s.code == {:A, :B, :C, :D, :A}

# CODE.INSERT
# Pushes the result of inserting the second item of the CODE stack into the
# first item, at the position indexed by the item at the top of the INTEGER
# stack. Indexing is computed as in CODE.EXTRACT.
s = Push.run("(CODE.QUOTE D CODE.QUOTE (A B C) CODE.INSERT)", cfg)
@test s.code == {:D, {:A, :B, :C}}
s = Push.run("(0 CODE.QUOTE D CODE.QUOTE (A B C) CODE.INSERT)", cfg)
@test s.code == {:D} && isempty(s.integer)
s = Push.run("(1 CODE.QUOTE X CODE.QUOTE ((A B C) (D E) (F)) CODE.INSERT)", cfg)
@test s.code == {{:X, {:D, :E}, {:F}}} && isempty(s.integer)
s = Push.run("(-1 CODE.QUOTE X CODE.QUOTE ((A B C) (D E) (F)) CODE.INSERT)", cfg)
@test s.code == {{:X, {:D, :E}, {:F}}} && isempty(s.integer)
s = Push.run("(2 CODE.QUOTE X CODE.QUOTE ((A B C) (D E) (F)) CODE.INSERT)", cfg)
@test s.code == {{{:X, :B, :C}, {:D, :E}, {:F}}} && isempty(s.integer)
s = Push.run("(-2 CODE.QUOTE X CODE.QUOTE ((A B C) (D E) (F)) CODE.INSERT)", cfg)
@test s.code == {{{:X, :B, :C}, {:D, :E}, {:F}}} && isempty(s.integer)
s = Push.run("(5 CODE.QUOTE X CODE.QUOTE ((A B C) (D E) (F)) CODE.INSERT)", cfg)
@test s.code == {{{:A, :B, :C}, :X, {:F}}} && isempty(s.integer)
s = Push.run("(-5 CODE.QUOTE X CODE.QUOTE ((A B C) (D E) (F)) CODE.INSERT)", cfg)
@test s.code == {{{:A, :B, :C}, :X, {:F}}} && isempty(s.integer)
s = Push.run("(4 CODE.QUOTE X CODE.QUOTE (A B C) CODE.INSERT)", cfg)
@test s.code == {:X} && isempty(s.integer)
s = Push.run("(6 CODE.QUOTE X CODE.QUOTE (A B C) CODE.INSERT)", cfg)
@test s.code == {{:A, :X, :C}} && isempty(s.integer)
s = Push.run("(-6 CODE.QUOTE X CODE.QUOTE (A B C) CODE.INSERT)", cfg)
@test s.code == {{:A, :X, :C}} && isempty(s.integer)

# CODE.CONTAINER
# Pushes the "container" of the 2nd CODE stack item within the first CODE
# stack item onto the CODE stack. If second item contains the first
# anywhere (i.e. in any nested list) then the container is the smallest
# sub-list that contains but is not equal to the first instance.
s = Push.run("(CODE.QUOTE X CODE.QUOTE (A B C) CODE.CONTAINER)", cfg)
@test s.code == {{}}
s = Push.run("(CODE.QUOTE X CODE.QUOTE (X A B) CODE.CONTAINER)", cfg)
@test s.code == {{:X, :A, :B}}
s = Push.run("(CODE.QUOTE X CODE.QUOTE (A ((B C) X D) E) CODE.CONTAINER)", cfg)
@test s.code == {{{:B, :C}, :X, :D}}
s = Push.run("(CODE.QUOTE X CODE.QUOTE ((A X) B C X) CODE.CONTAINER)", cfg)
@test s.code == {{:A, :X}}

# CODE.DISCREPANCY
s = Push.run("(CODE.QUOTE X CODE.QUOTE X CODE.DISCREPANCY)", cfg)
@test s.integer == {0}
s = Push.run("(CODE.QUOTE X CODE.QUOTE Y CODE.DISCREPANCY)", cfg)
@test s.integer == {2}
s = Push.run("(CODE.QUOTE (X Y Z) CODE.QUOTE (X Y A B) CODE.DISCREPANCY)", cfg)
@test s.integer == {3}

# CODE.MEMBER
s = Push.run("(CODE.QUOTE (A B C) CODE.MEMBER)", cfg)
@test s.code == {{:A, :B, :C}} && isempty(s.boolean)
s = Push.run("(CODE.QUOTE A CODE.QUOTE (A B C) CODE.MEMBER)", cfg)
@test isempty(s.code) && s.boolean == [true]
s = Push.run("(CODE.QUOTE (A) CODE.QUOTE (A B C) CODE.MEMBER)", cfg)
@test isempty(s.code) && s.boolean == [false]
s = Push.run("(CODE.QUOTE A CODE.QUOTE B CODE.MEMBER)", cfg)
@test isempty(s.code) && s.boolean == [false]
s = Push.run("(CODE.QUOTE A CODE.QUOTE A CODE.MEMBER)", cfg)
@test isempty(s.code) && s.boolean == [true]
s = Push.run("(CODE.QUOTE A CODE.QUOTE ((A) B C) CODE.MEMBER)", cfg)
@test isempty(s.code) && s.boolean == [false]

# CODE.CONTAINS
s = Push.run("(CODE.QUOTE (1 2 3) CODE.QUOTE 1 CODE.CONTAINS)", cfg)
@test isempty(s.code) && s.boolean == [true]
s = Push.run("(CODE.QUOTE (1 2 3) CODE.QUOTE 4 CODE.CONTAINS)", cfg)
@test isempty(s.code) && s.boolean == [false]
s = Push.run("(CODE.QUOTE 1 CODE.QUOTE 1 CODE.CONTAINS)", cfg)
@test isempty(s.code) && s.boolean == [false]
s = Push.run("(CODE.QUOTE () CODE.QUOTE 1 CODE.CONTAINS)", cfg)
@test isempty(s.code) && s.boolean == [false]
s = Push.run("(CODE.QUOTE () CODE.QUOTE () CODE.CONTAINS)", cfg)
@test isempty(s.code) && s.boolean == [false]
s = Push.run("(CODE.QUOTE 1 CODE.QUOTE () CODE.CONTAINS)", cfg)
@test isempty(s.code) && s.boolean == [false]
s = Push.run("(CODE.QUOTE () CODE.QUOTE () CODE.CONTAINS)", cfg)
@test isempty(s.code) && s.boolean == [false]
s = Push.run("(CODE.QUOTE ((1 2) 3) CODE.QUOTE (1 2) CODE.CONTAINS)", cfg)
@test isempty(s.code) && s.boolean == [true]
s = Push.run("(CODE.QUOTE (((2 1) 9 8) 2 3) CODE.QUOTE 1 CODE.CONTAINS)", cfg)
@test isempty(s.code) && s.boolean == [true]

# CODE.SUBST
# -- LEFT OUT FOR NOW

# CODE.RAND

# CODE.POSITION
# Does this operate like EXTRACT, or like MEMBER?

# CODE.INSTRUCTIONS

# CODE.DEFINE
s = Push.run("(ADD CODE.QUOTE (5 3 INTEGER.+) CODE.DEFINE ADD)", cfg)
@test isempty(s.code) && isempty(s.name) && s.integer == {8}

# CODE.DEFINITION


text:
"""
*Gamble*

`Gamble` constructs a gamble object with probability vector `p` and outcome vector `v`.

- `p`: probability vector
- `v`: outcome vector

Constructor
````julia
Gamble(;p=[.5,.5], v=[10.0,0.0])
````
"""
mutable struct Gamble{T1,T2}
    p::T1
    v::T2
end

function Gamble(;p=[.5,.5], v=[10.0,0.0])
    return Gamble(p, v)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 61d79b04-7db7-11eb-1f8f-cbd858393030
md"""
Engle-Granger Two-Step Error Correction Model

$y_{t} = log(Y_{t}), x_{t} = log(X_{t})$

Order 1

$y_{t} = c + \gamma x_{t} + z_{t}$

$\delta y_{t} = y_{t} - y_{t - 1}  = \alpha z_{t - 1} + \epsilon_{t}$

Test for Cointegration
+ Ad-Fuller Test on residuals $z_{t}$
"""

# ╔═╡ 3f05c2d8-7dba-11eb-3c50-f169294281fc
md"""

https://hudsonthames.org/correct-backtest-methodology-pairs-trading/

https://docs.google.com/presentation/d/1oFpv7OUi3W9F2D30rEQRy_Ra5ejLtlCMnk9nTJnMDr4/edit#slide=id.gb70b23aaa3_1_0

https://drive.google.com/file/d/1ZeJD81OrKln8QDxm1sU63ivRgqXCcQbS/view

https://docs.google.com/presentation/d/19xUNdddCsdXg2KizcUMVvdjhTnA8Cuq293brXB8rSWM/edit#slide=id.p

https://docs.google.com/presentation/d/1Cb9NizOkvCx6-p6X3AVZyxjviub7kF2sCPy2gBiT9nk/edit#slide=id.p

https://hudsonthames.org/copula-for-pairs-trading-overview-of-common-strategies/

https://hudsonthames.org/optimal-stopping-in-pairs-trading-ornstein-uhlenbeck-model/

Tim Leung, Optimal Mean Reversion Trading with Transaction Costs and Stop-Loss Exit (2015)
https://arxiv.org/abs/1411.5062

Statistical Arbitrage Strategy in Multi-Asset Market Using Time Series Analysis

https://www.scirp.org/journal/paperinformation.aspx?paperid=100361

"""

# ╔═╡ Cell order:
# ╠═61d79b04-7db7-11eb-1f8f-cbd858393030
# ╠═3f05c2d8-7dba-11eb-3c50-f169294281fc


text:
module Endowment

if  pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

if "$(pwd())\\data" ∉ LOAD_PATH
  push!(LOAD_PATH,"$(pwd())\\data")
end

#Use this to print
#=
using Weave
codeName = "train"
weave(Pkg.dir("$(pwd())\\$(codeName).jl"),
  informat="script",
  out_path = "$(pwd())\\output\\$(codeName)_Appendix.html",
  doctype = "md2html")
=#
##############0.7/1.0 compatability hacks
#NOTE: Delete these upon 0.7

#############Dependencies####################
using DataFrames, Distributions, CSV, GZip, ForwardDiff,
  Gadfly, NLopt, Measures, Formatting, HTTP,
  CodecZlib

using CT

##################Exports########################
export initializePar, ###Acquire NCCSData
cleanup,
getNCCSFiles,
getNCCSFilesIndex,
preProcessData,
validateSingleTypeCSV,
readSingleTypeCSV,
readMetaIndex,
formMetaData,
MetaTuple,
MetaDictionary,
computeEndowmentReturns

#############################Shared Constants###########
const OUTPUT_PATH = pwd() * "\\output"
const WORKING_PATH = pwd() * "\\working"
const DATA_PATH = pwd() * "\\data"
const NCCS_PATH = DATA_PATH * "\\NCCS"
const DOCUMENTATION_PATH = DATA_PATH * "\\documentation"

const METADATA_NAME = "NCCSMetadata"
const META_INDEX_NAME = "NCCSMetaIndex"
const META_PATH = DOCUMENTATION_PATH * "\\raw meta"
const SUB_TABLES_PATH = META_PATH * "\\subtables"
const FILES_INDEX_NAME = "NCCSFiles"
const KEEPER_NAME = "Keeper List"
const Nothing = Void #NOTE: 0.7 compatabiltiy hack, delete on ugrade
const VAR_TYPES = Dict(:Numeric=>Float64, :Character=>String)

const Conformed = Symbol
const MConformed = Union{Symbol, Missing}

#default compressor algorithms
I_CSV_STREAMER(p::String) = GzipDecompressorStream(open(p))
I_JLS_STREAMER(p::String) = GzipDecompressorStream(open(p))
O_JLS_STREAMER(p::String) = GzipCompressorStream(open(p, "w"))

include("EndowmentCode\\EndowmentPar.jl")
include("EndowmentCode\\MetaTuples.jl")
include("EndowmentCode\\MetaDictionaries.jl")
include("EndowmentCode\\EndowmentGetNCCS.jl")
include("EndowmentCode\\EndowmentData.jl")

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

using MXNet

#--------------------------------------------------------------------------------
# define lenet

# input
data = mx.Variable(:data)

# first conv
conv1 = @mx.chain mx.Convolution(data, kernel=(5,5), num_filter=20)  =>
                  mx.Activation(act_type=:tanh) =>
                  mx.Pooling(pool_type=:max, kernel=(2,2), stride=(2,2))

# second conv
conv2 = @mx.chain mx.Convolution(conv1, kernel=(5,5), num_filter=50) =>
                  mx.Activation(act_type=:tanh) =>
                  mx.Pooling(pool_type=:max, kernel=(2,2), stride=(2,2))

# first fully-connected
fc1   = @mx.chain mx.Flatten(conv2) =>
                  mx.FullyConnected(num_hidden=500) =>
                  mx.Activation(act_type=:tanh)

# second fully-connected
fc2   = mx.FullyConnected(fc1, num_hidden=10)

# softmax loss
lenet = mx.SoftmaxOutput(fc2, name=:softmax)

#--------------------------------------------------------------------------------
# load data
batch_size = 100
include("mnist-data.jl")
train_provider, eval_provider = get_mnist_providers(batch_size; flat=false)

#--------------------------------------------------------------------------------
# fit model
model = mx.FeedForward(lenet, context=mx.gpu())

# optimizer
optimizer = mx.SGD(η=0.05, μ=0.9, λ=0.00001)

# fit parameters
mx.fit(model, optimizer, train_provider, n_epoch=20, eval_data=eval_provider)

max_stars_repo_name: Vikas-kum/incubator-mxnet

max_issues_repo_name: Vikas-kum/incubator-mxnet

max_forks_repo_name: Vikas-kum/incubator-mxnet


text:
using ImageFiltering
using LAP_julia: prepare_gaussian_filters

filter_speed_df = DataFrame(
    index = Int[],
    filter_alg = Symbol[],
    img_size = Int[],
    khs = Int[],
    bench = BenchmarkTools.Trial[])

kernel_half_sizes = cat(collect(1:51), collect(61:10:101), collect(141:40:501), collect(501:98:795), dims=1);
img_sizes = [100, 200, 400, 800, 1600, 3200];
filter_algs = [Algorithm.FIR, Algorithm.FFT]

##
df = deepcopy(filter_speed_df)
let index = 0
    for img_size in img_sizes

        rand_img = rand(img_size, img_size)

        # limit the max kernel size to 1/4 the image
        kernel_limit = img_size/4
        khs_modified = filter(x -> x <= kernel_limit, kernel_half_sizes)

        for khs in khs_modified

            forward_ker, backward_ker = prepare_gaussian_filters(khs)

            for alg in filter_algs
                filt_result = similar([1.0], (img_size, img_size))

                # do the filtering, timed
                bench = @benchmark imfilter!($alg, $filt_result, $rand_img, $forward_ker[1], "symmetric")

                index = index + 1
                println("at index: ", index, " img_size:", img_size, " khs: ", khs)
                push!(df, Dict(:index => index,
                               :alg => Symbol(alg),
                               :img_size => img_size,
                               :khs => khs,
                               :bench => bench))
            end
        end
    end
end

filt_result = similar([1.0], (img_size, img_size))
img = gen_chess(25,2)
alg = Algorithm.FIR
res = AbstractResource{alg}
forward_ker, backward_ker = prepare_gaussian_filters(13)
# do the filtering, timed
imfilter(Float64, rand_img, forward_ker[1], "symmetric", Algorithm.FIR)
bench = @benchmark imfilter!($res, $filt_result, $rand_img, $forward_ker[1], "symmetric")


text:
## Julia program for Metropolization of the Gibbs Sampler
## author: weiya <szcfweiya@gmail.com>
## date: 2018-08-25

## example 7.3.8

## data
a = [0.06, 0.14, 0.11, 0.09];
b = [0.17, 0.24, 0.19, 0.20];
x = [9, 15, 12, 7, 8];

function ex738(T)
    z = ones(T+1, 4)
    for t = 1:T
        for i = 1:2
            atrial = 0
            for k = 1:x[i]
                u = rand()
                if u <= a[i]*mu # ? sample mu first
            end
        end
    end
    println(1)
end

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
module State

import Base: convert, copy, getindex, setindex!, length,
             show, start, done, next, endof, eltype, isapprox

export QuantumState, QUBIT0, QUBIT1, BELL_STATE,
       frombloch, fromstates, fromvector, randomstate

# A n-qubit quantum state is represented by the Kronecker-product of their
# qubits, stored in `vector`, and the number of bits that make up that state,
# stored in `bits`
type QuantumState
    vector::Vector{Complex{Float64}}
    bits::Int
end

# The two constants representing the |0> and |1> states (respectively)
const QUBIT0 = QuantumState([1,0], 1)
const QUBIT1 = QuantumState([0,1], 1)
const BELL_STATE  = QuantumState(1 / sqrt(2) .* [1,0,0,1], 2)

# Constructs a new quantum state from the given state vector. It automatically
# calculates the number of bits
function fromvector{T<:Number}(state::Vector{T})
    numbits = floor(Int, log2(length(state)))
    return QuantumState(state, numbits)
end

# Constructs a new quantum state by taking the Kronecker product of any number
# of other quantum states
function fromstates(state0::QuantumState, states::QuantumState...)
    state::Vector{Complex{Float64}} = mapfoldl(x -> x.vector, kron, state0.vector, states)
    return fromvector(state)
end

# Returns the 1-qubit quantum state that represents the (r, theta, phi) location
# on the Bloch-sphere
function frombloch(theta::Float64, phi::Float64)
    state = [cos(theta / 2), exp(im * phi) * sin(theta / 2)]
    return QuantumState(state, 1)
end

function isapprox(q1::QuantumState, q2::QuantumState, atol=0.0001)
    for (x, y) in zip(q1, q2)
        if !isapprox(x, y, atol = atol)
            return false
        end
    end
    return true
end

function randomstate(nbits::Int)
    vec = [rand() + rand() * im for _ in 1:2^nbits]
    n = sqrt(sumabs2(vec))
    vec = vec .* (1.0 / n)
    return QuantumState(vec, nbits)
end

convert(state::QuantumState) = state.vector
copy(state::QuantumState) = QuantumState(copy(state.vector), copy(state.bits))

length(state::QuantumState) = length(state.vector)
getindex(state::QuantumState, i::Int) = state.vector[i]
setindex!(state::QuantumState, c::Complex{Float64}, i::Int) = state.vector[i] = c

start(state::QuantumState) = start(state.vector)
done(state::QuantumState, current) = done(state.vector, current)
next(state::QuantumState, current) = next(state.vector, current)
endof(state::QuantumState) = endof(state.vector)
eltype(::Type{QuantumState}) = Complex{Float64}

function show(io::IO, state::QuantumState)
    for i = 1:length(state)
        basis = i - 1

        real = Base.real(state[i])
        resign = real < 0.0 ? '-' : ' '
        real = abs(real)

        imag = Base.imag(state[i])
        imsign = imag < 0.0 ? '-' : '+'
        imag = abs(imag)

        @printf("%c%04.4f %c %04.4fi    |%s>    P: %04.4f\n",
                resign, real, imsign, imag, bin(basis, state.bits), abs2(state[i]))
    end
end

end

max_stars_repo_name: leios/quantum

max_issues_repo_name: leios/quantum

max_forks_repo_name: leios/quantum


text:
include("Grp/Agnostic.jl")
include("Grp/GenGrp.jl")
include("Grp/SmallGroups.jl")
include("Grp/Morphisms.jl")


text:
type BatesModel{CalibratedModelType} <: Model{CalibratedModelType}
  modT::CalibratedModelType
  hestonModel::HestonModel
  nu::ConstantParameter
  delta::ConstantParameter
  lambda::ConstantParameter
  process::BatesProcess
  common::ModelCommon
end

function BatesModel(process::BatesProcess)
  hm = HestonModel(process.hestonProcess)
  nu = ConstantParameter([process.nu], NoConstraint())
  delta = ConstantParameter([process.delta], PositiveConstraint())
  lambda = ConstantParameter([process.lambda], PositiveConstraint())

  bm = BatesModel(CalibratedModelType(), hm, nu, delta, lambda, process, ModelCommon())

  generate_arguments!(bm)

  return bm
end

get_theta(bm::BatesModel) = get_theta(bm.hestonModel)
get_kappa(bm::BatesModel) = get_kappa(bm.hestonModel)
get_sigma(bm::BatesModel) = get_sigma(bm.hestonModel)
get_rho(bm::BatesModel) = get_rho(bm.hestonModel)
get_v0(bm::BatesModel) = get_v0(bm.hestonModel)
get_nu(bm::BatesModel) = bm.nu.data[1]
get_delta(bm::BatesModel) = bm.delta.data[1]
get_lambda(bm::BatesModel) = bm.lambda.data[1]

generate_arguments!(bm::BatesModel) =
  bm.process = BatesProcess(get_risk_free_rate(bm.process), get_dividend_yield(bm.process), get_s0(bm.process),
              get_v0(bm), get_kappa(bm), get_theta(bm), get_sigma(bm), get_rho(bm),
              get_nu(bm), get_delta(bm), get_lambda(bm))


text:
function read_trees(filename)

    open(filename) do f
        data = readlines(f)
        return transpose(map(x->(x=='#'), hcat(collect.(data)...)))
    end

end

function hit_count(trees, down, right)
    nrow, ncol = size(trees)
    hit = 0
    j = 1
    for i in 1:down:nrow
        hit += trees[i,j]
        j = mod1(j+right,ncol)
    end
    return hit
end

trees = read_trees("input03.txt")
println(hit_count(trees, 1, 3))
println(hit_count(trees, 1, 1) *
        hit_count(trees, 1, 3) *
        hit_count(trees, 1, 5) *
        hit_count(trees, 1, 7) *
        hit_count(trees, 2, 1) )

max_stars_repo_name: pnavaro/AdventOfCode

max_issues_repo_name: pnavaro/AdventOfCode

max_forks_repo_name: pnavaro/AdventOfCode


text:
module ExtensibleUnions

export addtounion!, # extensible_unions.jl
       extensiblefunction!, # extensible_functions.jl
       extensibleunion!, # extensible_unions.jl
       @ExtensibleUnion, # extensible_unions.jl
       ExtensibleUnion,  # extensible_unions.jl
       members,          # extensible_unions.jl
       isextensiblefunction, # extensible_functions.jl
       @extensible,      # extensible_functions.jl
       isextensibleunion # extensible_unions.jl

const _registry_extensibleunion_to_genericfunctions = Dict{Any, Any}()
const _registry_extensibleunion_to_members = Dict{Type, Set{Type}}()
const _registry_genericfunctions_to_extensibleunions = Dict{Any, Any}()

include("code_transformation.jl")
include("extensible_functions.jl")
include("extensible_unions.jl")
include("init.jl")
include("reset.jl")
include("update_methods.jl")

end # module


text:
### Data analysis ###

# this helper function works like
# aggregate(df, by_sym, mean)
# but uses column "weight_sym" as weights
function aggregate_wmean(df,by_sym,weight_sym)
  function wmean(df, by_sym, weight_sym)
    ndf=DataFrame()
    for sym in names(df)
      if sym==by_sym || sym==weight_sym continue end
      ndf[Symbol(sym,"_mean")]=mean(df[sym],WeightVec(df[weight_sym]))
    end
    return ndf
  end
  by(df,by_sym,d->wmean(d,by_sym,weight_sym))
end

function aggregate_by_age(data)
  workminutes_cutoff = 60
  # this would work in principle but cannot take the weights into account
  # work_df = aggregate(data[data[:work].>= workminutes_cutoff,:], :age, mean)
  # nwork_df = aggregate(data[data[:work].< workminutes_cutoff,:], :age, mean)

  work_df = aggregate_wmean(data[(data[:work].>= workminutes_cutoff) & (data[:employment] .== 1),:], :age, :weight)
  nwork_df = aggregate_wmean(data[(data[:work].< workminutes_cutoff) & (data[:employment] .== 1),:], :age, :weight)

  # save data for reproducibility
  mkpath("../data/processed")
  writetable("../data/processed/_work.csv",work_df)
  writetable("../data/processed/_nwork.csv",nwork_df)

  #remove "_mean" suffix from column names
  names_df = names(work_df)
  for sym in names_df
    if sym==:age continue end
    nsym = Symbol(replace(string(sym), "_mean", ""))
    rename!(work_df, sym, nsym)
    rename!(nwork_df, sym, nsym)
  end
  return work_df, nwork_df
end


text:
println("Starting to generate data and PDF. This may take a while.")

import Pkg

Pkg.activate(".")
Pkg.update()
Pkg.instantiate()

using BayesianAssignment

build()


text:
try
    using Revise
    ENV["JULIA_REVISE"] = "auto"
catch e
    @warn "Error initializing Revise" exception=(e, catch_backtrace())
end


text:
@test_throws ErrorException eccentricity(g4)
z = eccentricity(a1, distmx1)
@test z == [6.2, 4.2, 6.2]
@test diameter(z) == diameter(a1, distmx1) == 6.2
@test periphery(z) == periphery(a1, distmx1) == [1,3]
@test radius(z) == radius(a1, distmx1) == 4.2
@test center(z) == center(a1, distmx1) == [2]

z = eccentricity(a2, distmx2)
@test z == [6.2, 4.2, 6.1]
@test diameter(z) == diameter(a2, distmx2) == 6.2
@test periphery(z) == periphery(a2, distmx2) == [1]
@test radius(z) == radius(a2, distmx2) == 4.2
@test center(z) == center(a2, distmx2) == [2]

@test size(LightGraphs.DefaultDistance()) == (typemax(Int), typemax(Int))
d = LightGraphs.DefaultDistance(3)
@test size(d) == (3, 3)
@test d[1,1] == 1
@test d[1:2, 1:2] == LightGraphs.DefaultDistance(2)
@test d == transpose(d) == ctranspose(d)


text:
using Oceananigans
using CairoMakie
# using LazyGrids
# using UnicodePlots

grid = RegularRectilinearGrid(size=(128, 128), x=(-5, 5), y=(-5, 5),
topology=(Periodic, Periodic, Flat))
model = NonhydrostaticModel(grid=grid, tracers=:c, buoyancy=nothing)

initial_c(x, y, z) = exp.(-x.^2 - y.^2)
set!(model, u=1, c=initial_c)

xs =  LinRange(-5, 5, 128)
ys =  LinRange(-5, 5, 128)

#  i = 1

# while i <= 3

#     i += 1

# end
simulation = Simulation(model, Δt=1e-2, stop_iteration = 10)
run!(simulation)

function plot_tracer(simulation)
    c = simulation.model.tracers.c
    fig = Figure(resolution=(700, 450), fontsize=18, font="sans")
    ax = fig[1, 1] = Axis(fig, xlabel="x", ylabel="y")
    cMatrixTrimmed = interior(model.tracers.c)[:, :, 1]'
    CairoMakie.heatmap!(ax, xs, ys, cMatrixTrimmed, colormap = :deep)
    display(fig)
end

simulation.callbacks[:plotter] = Callback(plot_tracer, schedule=IterationInterval(10))

    #fig = Figure(resolution=(700, 450), fontsize=18, font="sans")
    #ax = fig[1, 1] = Axis(fig, xlabel="x", ylabel="y")
    #cMatrixTrimmed = interior(model.tracers.c)[:, :, 1]'
    #CairoMakie.heatmap!(xs, ys, cMatrixTrimmed, colormap = :deep)

    #fig

# 0:130 auf 0:130 Datensatz wird auf 1:128 auf 1:128 getrimmt
# matrixTrimmed = model.tracers.c.data[1:128,1:128, 1]
#CairoMakie.heatmap(interior(model.tracers.c)[:, :, 1]')
# returns a view of f that excludes halo points

# (xx, yy) = ndgrid_array(xs, ys)

# us = [-x for x in xs]
# vs = [2y for y in ys]
# u2(x, y) = -x
# v2(x, y) = 2y

# # wird jetzt durch model.tracers.c.data ersetzt
# # zz = exp.(-xx.^2 - yy.^2)

# vel1  =-xx
# vel2 = 2*yy

# xsToPlot = xs[begin:2:end]

# #xsToPlot1 = xs[1:3:128]
# ysToPlot = ys[begin:2:end]
# #ysToPlot1 = ys[1:3:128]
# #(xx, yy) = ndgrid_array(xsToPlot, ysToPlot)
# #a[(1:length(a)) .% 3 .!= 0]

# u = ones(2)
# v = ones(2)

# arrows(xsToPlot, ysToPlot, interior(model.tracers.c), v, arrowsize=5, lengthscale=0.3,
# normalize = false)


text:
#=
Project:      MIRTdev
File:         image_geom_mri.jl
Description:  Image geom function
Author:       Rajas Gupta
Date Created: October 1st 2019
Date Update:  October 7th 2019

2007-02-22, Jeff Fessler, University of Michigan (image_geom_mri.m)
=#

export image_geom_mri

#using MIRT: image_geom

"""
    ig = image_geom_mri(varargin)

Same as `image_geom()` but default offsets are 0.5 so that
image pixel indices go from -N/2 to N/2-1.

"""
image_geom_mri( ; kwargs...) =
    image_geom( ; offset_x = .5, offset_y = .5, offset_z = .5, kwargs...)


text:
using Test

using LinearAlgebra

using CLIMA.TurbulenceConvection.TriDiagSolvers
TDMA = TriDiagSolvers

@testset "TriDiagSolvers" begin
  N = 2:10
  for n in N
    dl = rand(n-1)
    du = rand(n-1)
    d = rand(n);

    A = Array(Tridiagonal(dl, d, du))
    b = rand(length(d))

    x_correct = inv(A)*b

    xtemp = zeros(n)
    gamma = zeros(n-1)
    beta = zeros(n)
    x_TDMA = zeros(n)

    solve_tridiag!(x_TDMA, b, dl, d, du, n, xtemp, gamma, beta)
    tol = eps(Float32)

    err = [abs(x-y) for (x, y) in zip(x_correct, x_TDMA)]
    # @show eps(Float64)
    # @show "Full", err
    if !all([x<tol for x in err])
        @show tol
        @show "Full", err
    end
    @assert all([x<tol for x in err])

    init_β_γ!(beta, gamma, dl, d, du, n)
    solve_tridiag_stored!(x_TDMA, b, dl, beta, gamma, n, xtemp)

    err = [abs(x-y) for (x, y) in zip(x_correct, x_TDMA)]
    # @show "Stored", err
    @assert all([x<tol for x in err])

    dl_mod = zeros(n)
    dl_mod[2:end] = dl
    du_mod = zeros(n)
    du_mod[1:end-1] = du
    TDMA.solve_tridiag_old(n, b, dl_mod, d, du_mod)
    x_TDMA = b
    err = [abs(x-y) for (x, y) in zip(x_correct, x_TDMA)]
    # @show "old", err
    @assert all([x<tol for x in err])

  end
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

"""
Exception thrown when an error occurred calling MXNet API.
"""
abstract type AbstractMXError <: Exception end

"General MXNet API error"
struct MXError <: AbstractMXError
  msg::String

  MXError(s::AbstractString) = new(string(s))
end

Base.show(io::IO, e::AbstractMXError) = print(io, e.msg)


text:
hook(f, x) = x

@adjoint! hook(f, x) = x, x̄ -> (nothing, f(x̄),)

macro showgrad(x)
  :(hook($(esc(x))) do x̄
      println($"D($x) = ", repr(x̄))
      x̄
    end)
end


text:
"""
$(TYPEDEF)

The abstract representation of the cost for a player in a game.
"""
abstract type PlayerCost end

"""
    $(FUNCTIONNAME)(x::SVector{nx}, u::SVector{nu}, t::AbstractFloat)

A functor to evaluate the cost representation.
"""
function (pc::PlayerCost) end

"""
    $FUNCTIONNAME(pc::PlayerCost, x::SVector{nx}, u::SVector{nu}, t::AbstractFloat)

Computes the quadratic approximation of the cost.

We provide a convencience default below.
"""
function quadraticize end

"""
    $(FUNCTIONNAME)(pc::PlayerCost)

Returns the index of the player this cost is associated with.
"""
function player_id end

"""
$(TYPEDEF)

Represents a parameter free player cost that is defined by a funciton f.
"""
struct FunctionPlayerCost{F} <: PlayerCost
    f::F
end
(pc::FunctionPlayerCost)(args...) = pc.f(args...)

"""
$(TYPEDEF)

Represents the quadratic players costs for a singel player in a game.

# Fields

$(TYPEDFIELDS)
"""
struct QuadraticPlayerCost{nx, nu, TL<:SVector{nx}, TQ<:SMatrix{nx, nx},
                           TRL<:SVector{nu}, TRQ<:SMatrix{nu, nu}} <: PlayerCost
    "The linear state cost"
    l::TL
    "The qudratic state cost matrix"
    Q::TQ
    "The linear term of the input cost."
    r::TRL
    "A square matrix to represent the quadratic control cost for this player"
    R::TRQ
end

(pc::QuadraticPlayerCost)(::Any, x::SVector, u::SVector) = pc(x, u)

function (pc::QuadraticPlayerCost)(x::SVector, u::SVector)
    return 1//2 * x'*pc.Q*x + pc.l'*x + 1//2 * u'*pc.R*u + pc.r'*u
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

"""

Parameters
----------
distval
    a charasteristic measure to skip element pair, 0..5 => near, 10+ => far
    5 means that distance of slave element midpoint and point to project
    is 5 times larger than length of element
"""
mutable struct Contact <: BoundaryProblem
    dimension :: Int
    rotate_normals :: Bool
    finite_sliding :: Bool
    friction :: Bool
    dual_basis :: Bool
    use_forwarddiff :: Bool
    forwarddiff_assemble_in_pieces :: Bool
    minimum_active_set_size :: Int
    distval :: Float64
    remove_from_set :: Bool # allow removal of non-potential contact pairs
    allow_quads :: Bool
    remove_nodes :: Vector{Int}
    always_in_contact :: Bool
    update_contact_pairing :: Bool
    iteration :: Int
    contact_state_in_first_iteration :: Symbol
    alpha :: Float64
    drop_tolerance :: Float64
    store_fields :: Vector{AbstractString}
end

function Contact()
    default_fields = ["element area", "contact area", "weighted gap",
        "contact pressure", "active nodes", "inactive nodes", "stick nodes",
        "slip nodes", "complementarity condition", "contact error"]
    return Contact(
        -1,      # dimension
        false,   # rotate_normals
        false,   # finite_sliding
        false,   # friciton
        true,    # dual basis
        false,   # use forwarddiff
        false,   # when using forwarddiff, assemble interface in pieces
        0,       # minimum active set size
        5.0,     # distance value for contact detection
        false,   # allow removal of non-potential contact pairs
        false,   # allow quadrangles in contact discretization
        [],      # remove these nodes always from set
        false,   # mainly for debugging, do not remove inactive nodes
        true,    # update contact pairing on each loop
        1,       # iteration counter
        :AUTO,   # contact state in first iteration, AUTO, INACTIVE, ACTIVE
        0.0,     # alpha basis transform parameter
        1.0e-9,  # drop tolerance
        default_fields)
end

function assemble!(problem::Problem{Contact}, time::Real)
    if problem.properties.dimension == -1
        problem.properties.dimension = dim = size(first(problem.elements), 1)
    end
    dimension = Val{problem.properties.dimension}
    finite_sliding = Val{problem.properties.finite_sliding}
    friction = Val{problem.properties.friction}
    use_forwarddiff = Val{problem.properties.use_forwarddiff}
    assemble!(problem, time, dimension, finite_sliding, friction, use_forwarddiff)
    problem.properties.iteration += 1
end

const ContactElements2D = Union{Seg2}