text:
# # DataFrameMacros.jl

# DataFrameMacros.jl is an opinionated take on DataFrame manipulation in Julia with a syntax geared towards clarity, brevity and convenience.
# It offers macros that translate expressions into [DataFrames.jl](https://github.com/JuliaData/DataFrames.jl) function calls.

# Here is a simple example:

using DataFrameMacros #hide
using DataFrames #hide
df = DataFrame(name = ["Mary Louise Parker", "Thomas John Fisher"])
@transform(df, :middle_initial = split(:name)[2][1] * ".")

# Unlike DataFrames.jl, most operations are **row-wise** by default.
# This often results in cleaner code that's easier to understand and reason about, especially when string or object manipulation is involved.
# Such operations often don't have a clean broadcasting syntax, for example, `somestring[2]` is easier to read than `getindex.(somestrings, 2)`.
# The same is true for `someobject.property` and `getproperty.(someobjects, :property)`.

# The following macros are currently available:
# - `@transform` / `@transform!`
# - `@select` / `@select!`
# - `@groupby`
# - `@combine`
# - `@subset` / `@subset!`
# - `@sort` / `@sort!`
# - `@unique`

# Together with [Chain.jl](https://github.com/jkrumbiegel/Chain.jl), you get a convient syntax for chains of transformations:

using DataFrameMacros
using DataFrames
using Chain
using Random
using Statistics
Random.seed!(123)

df = DataFrame(
    id = shuffle(1:5),
    group = rand('a':'b', 5),
    weight_kg = randn(5) .* 5 .+ 60,
    height_cm = randn(5) .* 10 .+ 170)

@chain df begin
    @subset(:weight_kg > 50)
    @transform(:BMI = :weight_kg / (:height_cm / 100) ^ 2)
    @groupby(iseven(:id), :group)
    @combine(:mean_BMI = mean(:BMI))
    @sort(sqrt(:mean_BMI))
end

# ## Design choices

# These are the most important aspects that differ from other packages ([DataFramesMeta.jl](https://github.com/JuliaData/DataFramesMeta.jl) in particular):

# - All macros except `@combine` work **row-wise** by default. This reduces syntax complexity in most cases because no broadcasting is necessary. A flag macro (`@c` or `@r`) can be used to switch between row/column-based mode when needed.
# - `@groupby` and `@sort` allow using arbitrary expressions including multiple columns, without having to `@transform` first and repeat the new column names.
# - Column expressions are interpolated into the macro with `$`.
# - Keyword arguments to the macro-underlying functions work by separating them from column expressions with the `;` character.
# - Target column names are written with `:` symbols to avoid visual ambiguity (`:newcol = ...`). This also allows to use `AsTable` as a target like in DataFrames.jl.
# - The flag macro can also include the character `m` to switch on automatic `passmissing` in row-wise mode.
# - There is also a `@t` flag macro, which extracts every `:sym = expression` expression and collects the new symbols in a named tuple, while setting the target to `AsTable`.

# ## Examples

# TOC

using DataFrameMacros
using DataFrames
using Random
using Statistics
Random.seed!(123)

df = DataFrame(
    id = shuffle(1:5),
    group = rand('a':'b', 5),
    weight_kg = randn(5) .* 5 .+ 60,
    height_cm = randn(5) .* 10 .+ 170)

# ### @select
@select(df, :height_m = :height_cm / 100)
#-
@select(df, AsTable = (w = :weight_kg, h = :height_cm))

# ### @transform
@transform(df, :weight_g = :weight_kg * 1000)
#-
@transform(df, :BMI = :weight_kg / (:height_cm / 100) ^ 2)
# #### column flag @c
@transform(df, :weight_z = @c (:weight_kg .- mean(:weight_kg)) / std(:weight_kg))

# ### @groupby & @combine
g = @groupby(df, iseven(:id))
#-
@combine(g, :total_weight_kg = sum(:weight_kg))

# ### @sort

@sort(df, sqrt(:height_cm) / :weight_kg; rev = true)

# ### @unique

namedf = DataFrame(name = ["Joe Smith", "Eric Miller", "Frank Smith"])
@unique(namedf, last(split(:name)))

# ### interpolating column expressions

# If you have a variable or expression that you want to use as a column identifier, interpolate it into the macro with `$`.

the_column = :weight_kg
@combine(df, :total_weight = sum($the_column))
#-

a_string = "weight"
@combine(df, :total_weight = sum($(a_string * "_kg")))

# You can use strings and integers directly with `$`.
@combine(df, :sum = sum($"weight_kg" .* $4))

# ### passmissing flag @m

df = DataFrame(name = ["joe", "jim", missing, "james"])

@transform(df, :cap_name = @m uppercasefirst(:name))

# ### escaping symbols

# The `$` symbol usually signals that an expression is to be used as a column identifier.
# The only exception is `$` in front of a bare symbol.
# In that case, it signals that the symbol should be left as it is.

df = DataFrame(color = [:red, :green, :blue])
@transform(df, :is_red = :color == $:red)

# ### `@t` flag macro for automatic `AsTable`

# To use `AsTable` as a target, you usually have to construct a NamedTuple in the passed function.
# You can avoid both passing `AsTable` explicitly and constructing the NamedTuple by using the `@t` flag macro.
# All expressions of the type `:symbol = expression` are collected, the `:symbol`s are replaced with anonymous variables, and these variables are collected in a NamedTuple as the return value automatically.

df = DataFrame(a = 1:3, b = 4:6)
df2 = @transform df @t begin
    x = :a + :b
    :y = x * 2
    :z = x + 4
end

# ### block syntax

# You can pass a begin/end block to every macro instead of multiple separate arguments.

df = DataFrame(
    id = shuffle(1:5),
    group = rand('a':'b', 5),
    weight_kg = randn(5) .* 5 .+ 60,
    height_cm = randn(5) .* 10 .+ 170)

@transform df begin
    :weight_g = :weight_kg / 1000
    :BMI = :weight_kg / (:height_cm / 100) ^ 2
    :weight_z = @c (:weight_kg .- mean(:weight_kg)) / std(:weight_kg)
end


text:
 #mtdna_model = @reaction_network mtdna begin
 # b, wildtype --> 2wildtype
 # d, wildtype  --> null
 # bmut, mut --> 2mut
 # dmut, mut --> null
 # m, wildtype --> wildtype + mut
 # # can we have wildtype --> 2mut
#end b d bmut dmut m

using Random
using Colors
using Plots
theme(:ggplot2)

function update(t, vals)
  wt = vals[1];
  mut = vals[2];
  haz = [ifelse(sum(vals)<target,b*wt,0.0), d*wt, ifelse(sum(vals)<target,bmut*mut,0.0), dmut*mut, m*wt];
  # ^ specifying the reactions. 1st = wt replication,  b, wildtype --> 2wildtype
  # 2nd = wt degredation, d, wildtype  --> null
  # 3rd = mut replication, bmut, mut --> 2mut
  # 4th = mut degredation, dmut, mut --> null
  # 5th = mutaion, m, wildtype --> wildtype + mut
  hazfrac = cumsum(haz/sum(haz));
  index = minimum(findall(hazfrac .- Random.rand() .> 0));
  tnew = t + (1.0/sum(haz))*(Random.randexp());
  valsnew = [wt + delta_wt[index], mut + delta_mut[index]]
  return((tnew,valsnew));
end

b = 0.001875*365
d = 0.001875*365*0.86 # Note, slightly lower than b, so that when total is low, population grows
bmut = 0.001875*365
dmut = 0.001875*365*0.86 # Note, slightly lower than b, so that when total is low, population grows
m = (5e-6)*365
target = 100.0

vals0 = [100, 0]
tmax = 100.0

delta_wt = [1,-1,0,0,0];
delta_mut = [0,0,1,-1,1];

simres = [vals0]
times = [0.0]

while (times[end] < tmax) & (sum(simres[end]) > 0)
 global simres, times
 (tnew, valsnew) = update(times[end],simres[end])
 push!(simres,valsnew)
 push!(times,tnew)
end

sol = (t = times, u = simres)

plot(sol.t, [sol.u[i][1] for i in 1:length(sol.t)], labels="Wildtype", color= "blue", lw=1.0)
plot!(sol.t, [sol.u[i][2] for i in 1:length(sol.t)], labels="Mutant", color="red", lw=1.0)
plot!(sol.t, [sol.u[i][1]+sol.u[i][2] for i in 1:length(sol.t)], legend=:left, labels="Total", color="white", lw=0.2, title="Low Initial Mutation Load Where b<bmut", xaxis="Time (Years)", yaxis="Population Size")

plot!(sol.t, [sol.u[i][1] for i in 1:length(sol.t)], labels="", color= "blue",lw=0.4)
plot!(sol.t, [sol.u[i][2] for i in 1:length(sol.t)], labels="", color="red", lw=0.4)
plot!(sol.t, [sol.u[i][1]+sol.u[i][2] for i in 1:length(sol.t)], legend=:left,layout=4, labels="", color="white", lw=0.2)

function sim_gillespie(;b,d,bmut,dmut,m,target,vals0,tmax)
    function update(t, vals)
      wt = vals[1];
      mut = vals[2];
      haz = [ifelse(sum(vals)<target,b*wt,0.0), d*wt, ifelse(sum(vals)<target,bmut*mut,0.0), dmut*mut, m*wt];
      # ^ specifying the reactions. 1st = wt replication,  b, wildtype --> 2wildtype
      # 2nd = wt degredation, d, wildtype  --> null
      # 3rd = mut replication, bmut, mut --> 2mut
      # 4th = mut degredation, dmut, mut --> null
      # 5th = mutaion, m, wildtype --> wildtype + mut
      hazfrac = cumsum(haz/sum(haz));
      index = minimum(findall(hazfrac .- Random.rand() .> 0));
      tnew = t + (1.0/sum(haz))*(Random.randexp());
      valsnew = [wt + delta_wt[index], mut + delta_mut[index]]
      return((tnew,valsnew));
    end
    delta_wt = [1,-1,0,0,0];
    delta_mut = [0,0,1,-1,1];

    simres = [vals0]
    times = [0.0]

    while (times[end] < tmax) & (sum(simres[end]) > 0)
     (tnew, valsnew) = update(times[end],simres[end])
     push!(simres,valsnew)
     push!(times,tnew)
    end
    sol = (t = times, u = simres)
    return(sol)
end

sol = sim_gillespie(;b = 0.001875*365,d = 0.001875*365*0.86,bmut = 0.001875*365,dmut = 0.001875*365*0.86,m = (5e-6)*365, target = 100.0,vals0 = [100, 0],tmax = 100.0)

nsims=150
transparency = 0.4
thick = 0.2
solutions =[sim_gillespie(;b = 0.001875*365,d = 0.001875*365*0.86,bmut = 0.001875*365*1.0,dmut = 0.001875*365*0.86,m = (5e-7)*365, target = 100.0,vals0 = [100, 0],tmax = 100.0) for i in 1:nsims]
panelA = plot(solutions[1].t,[solutions[1].u[i][1]+solutions[1].u[i][2] for i in 1:length(solutions[1].t)], labels="", color=RGBA(0,0,0,transparency), lw=thick)
for j in 2:nsims
    plot!(solutions[j].t,[solutions[j].u[i][1]+solutions[j].u[i][2] for i in 1:length(solutions[j].t)], labels="", color=RGBA(0,0,0,transparency), lw=thick, legend=false, xaxis="Time (Years)", yaxis="mtDNA copy number", title="High initial pathogenic variant load where b=bmut", title_location=:left)
end
panelB = plot(solutions[1].t,[solutions[1].u[i][2]/(solutions[1].u[i][1]+solutions[1].u[i][2]) for i in 1:length(solutions[1].t)], labels="", color=RGBA(1,0,0,transparency), lw=thick)
for j in 2:nsims
    plot!(solutions[j].t,[solutions[j].u[i][2]/(solutions[j].u[i][1]+solutions[j].u[i][2]) for i in 1:length(solutions[j].t)], labels="", color=RGBA(1,0,0,transparency), lw=thick, legend=false, xaxis="Time (Years)", yaxis="Pathogenic variant load")
    #plot!(rand(100), lw=0)
end

plot(panelA,panelB,layout=(1,2), fmt=:pdf)

savefig("D:\\himl.pdf")

nsims=1
transparency = 1.0
thick = 0.7
solutions =[sim_gillespie(;b = 0.001875*365,d = 0.001875*365*0.86,bmut = 0.001875*365*1.1,dmut = 0.001875*365*0.86,m = (5e-7)*365, target = 100.0,vals0 = [100, 0],tmax = 100.0) for i in 1:nsims]
panelC = plot(solutions[1].t,[solutions[1].u[i][1]+solutions[1].u[i][2] for i in 1:length(solutions[1].t)], labels="mtDNA copy number", color=RGBA(0,0,0,transparency), lw=thick, legend=false,  xaxis="Time (Years)", yaxis="Population Size", title="Copy Number of Cell")
panelB = plot(solutions[1].t,[solutions[1].u[i][2] for i in 1:length(solutions[1].t)], color=RGBA(0.5,0,0.5,transparency), lw=thick, legend=false,  xaxis="Time (Years)", yaxis="Population Size", title="Mutant mtDNA")
panelA = plot(solutions[1].t,[solutions[1].u[i][1] for i in 1:length(solutions[1].t)], color=RGBA(0,0,1,transparency), lw=thick, legend=false,  xaxis="Time (Years)", yaxis="Population Size", title="Wild-type mtDNA")
panelD = plot(solutions[1].t,[solutions[1].u[i][2]/(solutions[1].u[i][1]+solutions[1].u[i][2]) for i in 1:length(solutions[1].t)], title="Pathogenic variant load", color=RGBA(1,0,0,transparency), lw=thick, legend=false, xaxis="Time (Years)", yaxis="Fraction Value")

plot(panelA,panelB,panelC,panelD,layout=(2,2), fmt=:pdf)

savefig("D:\\4panel.pdf")

# What happens if a cells is born with no mutations?

# What happens if a cell inherits a lot of mutations?  70%?

# What happens in both of the scenarios above if bmut > b?

ext: jl


text:
function compute()
	N=100
	s1=sum(i for i=1:N)
	s2=sum(i*i for i=1:N)

	return s1^2-s2
end

print(compute())


text:
using JLD2

@doc """
    runMC(param::Parameter)
    runMC(params::AbstractArray{Parameter}
          ;
          parallel::Bool=false,
          autoID::Bool=true)

Runs Monte Carlo simulation(s) and returns calculated observables.

If a checkpoint file named `"\$(param["Checkpoint Filename Prefix"])_\$(param["ID"]).jld2"` exists and
`param["Checkpoint Interval"] > 0.0`, `runMC` loads this file and restarts the pending simulation.

# Keyward aruguments
- `autoID`: If true, `"ID"`s will be set (overwritten) as `params[i]["ID"] = i`.
- `parallel`: If true, runs simulations in parallel (uses `pmap` instead of `map`).

# Required keys in `param`
- "Model"
- "Update Method"

# Optional keys in `param`
- "MCS": The number of Monte Carlo steps after thermalization
    - Default: `8192`
- "Thermalization": The number of Monte Carlo steps for thermalization
    - Default: `MCS>>3`
- "Seed": The initial seed of the random number generator, `MersenneTwister`
    - Default: determined randomly (see `Random.srand`)
- "Checkpoint Filename Prefix": See above document.
    - Default: `"cp"`
- "ID": See above document.
    - Default: `0`
- "Checkpoint Interval": Time interval between writing checkpoint file in seconds.
    - Default: `0.0`, this means that NO checkpoint file will be loaded and saved.
"""
function runMC(params::AbstractArray{T}; parallel::Bool=false, autoID::Bool=true) where T<:Dict
    map_fn = ifelse(parallel, pmap, map)
    return map_fn(enumerate(params)) do idp
        id,p = idp
        if autoID
            p["ID"] = id
        end
        return runMC(p)
    end
end

function runMC(param::Parameter)
    model = param["Model"](param)
    if "Seed" in keys(param)
        srand(model, param["Seed"])
    end
    ret = runMC(model, param)
    return ret
end

function runMC(model, param::Parameter)
    MODEL = typeof(model)
    verbose = get(param, "Verbose", false) :: Bool
    if verbose
        println("Start: ", param)
    end
    cp_filename = @sprintf("%s_%d.jld2",
                           get(param, "Checkpoint Filename Prefix", "cp")::String,
                           get(param, "ID", 0)::Int)
    cp_interval = get(param, "Checkpoint Interval", 0.0) :: Float64
    if VERSION > v"0.6.4" && cp_interval != 0.0
        Compat.@info """"Checkpoint Interval" is set to 0.0 automatically since current JLD2.jl cannot save/load Random.MersenneTwister object properly in julia-0.7"""
        cp_interval = 0.0
    end
    tm = time()

    MCS = get(param, "MCS", 8192) :: Int
    Therm = get(param, "Thermalization", MCS>>3) :: Int

    mcs = 0
    MCS += Therm
    obs = BinningObservableSet()
    makeMCObservable!(obs, "Time per MCS")
    makeMCObservable!(obs, "MCS per Second")

    if cp_interval > 0.0 && ispath(cp_filename)
        jld = jldopen(cp_filename)
        model = jld["model"] :: MODEL
        obs = jld["obs"] :: BinningObservableSet
        mcs = convert(Int, jld["mcs"]) :: Int
        close(jld)
    end

    if "UpdateMethod" in keys(param)
        warn("\"UpdateMethod\" is deprecated. Use instead \"Update Method\".")
        param["Update Method"] = param["UpdateMethod"] :: Function
    end
    update! = param["Update Method"] :: Function
    estimator = get(param, "Estimator", default_estimator(model, update!)) :: Function
    p = convert_parameter(model, param)

    while mcs < MCS
        if mcs < Therm
            update!(model, p...)
        else
            t = @elapsed begin
                st = update!(model,p...)
                localobs = estimator(model, p..., st)
            end
            obs["Time per MCS"] << t
            obs["MCS per Second"] << 1.0/t
            accumulateObservables!(model, obs, localobs)
        end
        mcs += 1
        if cp_interval > 0.0 && time() - tm > cp_interval
            jld = jldopen(cp_filename, "w")
            jld["model"] = model
            jld["obs"] = obs
            jld["mcs"] = mcs
            close(jld)
            tm += cp_interval
         end
    end

    if cp_interval > 0.0
        jld = jldopen(cp_filename, "w")
        jld["model"] = model
        jld["obs"] = obs
        jld["mcs"] = mcs
        close(jld)
    end

    jk = postproc(model, param, obs)

    if verbose
        println("Finish: ", param)
    end
    return jk
end

@doc """
    accumulateObservables!(model, obs::MCObservableSet, localobs::Dict)

Accumulates `localobs` into `obs`. For example, `obs["Energy"] << localobs["Energy"]`.
"""
function accumulateObservables!(::Model, obs::MCObservableSet, localobs::Measurement)
    if length(obs) < 3
        @inbounds for key in keys(localobs)
            makeMCObservable!(obs, key)
            obs[key] << localobs[key]
        end
    else
        @inbounds for key in keys(localobs)
            obs[key] << localobs[key]
        end
    end
    return obs
end

@doc """
    postproc(model::Model, param::Parameter, obs::MCObservableSet)

Post process of observables. For example, Specific heat will be calculated from energy, energy^2, and temperature.
"""
function postproc end

@doc doc"""
    postproc(model::Union{Ising, Potts}, param::Parameter, obs::MCObservableSet)

# Observables to be calculated
In the following, $m$ is total magnetization per site and $\epsilon$ is total energy per site.

- `"Binder Ratio"`
    - $R := \frac{\left \langle m^4 \right \rangle}{\left \langle m^2 \right\rangle^2}$
- `"Susceptibility"`
    - $\chi := \frac{N}{T}\left(\left\langle m^2\right\rangle\right)$
- `"Connected Susceptibility"`
    - $\frac{N}{T}\left(\left\langle m^2\right\rangle - \left\langle |m| \right\rangle^2\right)$
- `"Specific Heat"`
    - $\frac{N}{T^2}\left(\left\langle \epsilon^2\right\rangle - \left\langle \epsilon \right\rangle^2\right)$
"""
function postproc(model::Union{Ising, Potts}, param::Parameter, obs::MCObservableSet)
    nsites = numsites(model)
    T = param["T"] :: Float64
    beta = 1.0/T

    jk = jackknife(obs)
    jk["Binder Ratio"] = jk["Magnetization^4"] / (jk["Magnetization^2"]^2)
    jk["Susceptibility"] = (nsites*beta)*jk["Magnetization^2"]
    jk["Connected Susceptibility"] = (nsites*beta)*(jk["Magnetization^2"] - jk["|Magnetization|"]^2)
    jk["Specific Heat"] = (nsites*beta*beta)*(jk["Energy^2"] - jk["Energy"]^2)
    return jk
end

@doc doc"""
    postproc(model::Union{Clock, XY}, param::Parameter, obs::MCObservableSet)

# Observables to be calculated
In the following, $m$ is total magnetization per site and $\epsilon$ is total energy per site.

- `"Binder Ratio x"`
    - $\frac{\left \langle m_x^4 \right \rangle}{\left \langle m_x^2 \right\rangle^2}$
- `"Binder Ratio y"`
    - $\frac{\left \langle m_y^4 \right \rangle}{\left \langle m_y^2 \right\rangle^2}$
- `"Binder Ratio"`
    - $\frac{\left \langle |m|^4 \right \rangle}{\left \langle |m|^2 \right\rangle^2}$
- `"Susceptibility x"`
    - $\frac{N}{T}\left(\left\langle m_x^2\right\rangle\right)$
- `"Susceptibility y"`
    - $\frac{N}{T}\left(\left\langle m_y^2\right\rangle\right)$
- `"Susceptibility y"`
    - $\frac{N}{T}\left(\left\langle |m|^2\right\rangle\right)$
- `"Connected Susceptibility x"`
    - $\frac{N}{T}\left(\left\langle m_x^2\right\rangle - \left\langle |m_x| \right\rangle^2\right)$
- `"Connected Susceptibility y"`
    - $\frac{N}{T}\left(\left\langle m_y^2\right\rangle - \left\langle |m_y| \right\rangle^2\right)$
- `"Connected Susceptibility"`
    - $\frac{N}{T}\left(\left\langle |m|^2\right\rangle - \left\langle |m| \right\rangle^2\right)$
- `"Specific Heat"`
    - $\frac{N}{T^2}\left(\left\langle \epsilon^2\right\rangle - \left\langle \epsilon \right\rangle^2\right)$
"""
function postproc(model::Union{Clock, XY}, param::Parameter, obs::MCObservableSet)
    nsites = numsites(model)
    T = param["T"]
    beta = 1.0/T

    jk = jackknife(obs)
    jk["Binder Ratio x"] = jk["Magnetization x^4"] / (jk["Magnetization x^2"]^2)
    jk["Binder Ratio y"] = jk["Magnetization y^4"] / (jk["Magnetization y^2"]^2)
    jk["Binder Ratio"] = jk["|Magnetization|^4"] / (jk["|Magnetization|^2"]^2)
    jk["Susceptibility x"] = (nsites*beta)*jk["Magnetization x^2"]
    jk["Susceptibility y"] = (nsites*beta)*jk["Magnetization y^2"]
    jk["Susceptibility"] = (nsites*beta)*jk["|Magnetization|^2"]
    jk["Connected Susceptibility x"] = (nsites*beta)*(jk["Magnetization x^2"] - jk["|Magnetization x|"]^2)
    jk["Connected Susceptibility y"] = (nsites*beta)*(jk["Magnetization y^2"] - jk["|Magnetization y|"]^2)
    jk["Connected Susceptibility"] = (nsites*beta)*(jk["|Magnetization|^2"] - jk["|Magnetization|"]^2)
    jk["Specific Heat"] = (nsites*beta*beta)*(jk["Energy^2"] - jk["Energy"]^2)
    return jk
end

@doc doc"""
    postproc(model::QuantumXXZ, param::Parameter, obs::MCObservableSet)

# Observables to be calculated
In the following, $s$ is sign of weight, $m$ is total magnetization per site,
and $\epsilon$ is total energy per site.

- `"Magnetization"`
    - $\left\langle m s\right\rangle\Big/\left\langle s \right\rangle$
- `"|Magnetization|"`
    - $\left\langle |m| s\right\rangle\Big/\left\langle s \right\rangle$
- `"Magnetization^2"`
    - $\left\langle m^2 s\right\rangle\Big/\left\langle s \right\rangle$
- `"Magnetization^4"`
    - $\left\langle m^4 s\right\rangle\Big/\left\langle s \right\rangle$
- `"Energy"`
    - $\left\langle \epsilon s\right\rangle\Big/\left\langle s \right\rangle$
- `"Energy^2"`
    - $\left\langle \epsilon^2 s\right\rangle\Big/\left\langle s \right\rangle$
- `"Binder Ratio"`
    - $\frac{\left \langle m^4 \right \rangle}{\left \langle m^2 \right\rangle^2}$
- `"Susceptibility"`
    - $\frac{N}{T}\left(\left\langle m^2\right\rangle\right)$
- `"Connected Susceptibility"`
    - $\frac{N}{T}\left(\left\langle m^2\right\rangle - \left\langle |m| \right\rangle^2\right)$
- `"Specific Heat"`
    - $\frac{N}{T^2}\left(\left\langle \epsilon^2\right\rangle - \left\langle \epsilon \right\rangle^2\right)$
"""
function postproc(model::QuantumXXZ, param::Parameter, obs::MCObservableSet)
    nsites = numsites(model)
    T = param["T"]
    beta = 1.0/T

    jk = jackknife(obs)

    for oname in ("Magnetization", "|Magnetization|",
                  "Magnetization^2", "Magnetization^4",
                  "Energy", "Energy^2",
                 )
        jk[oname] = jk["Sign * $oname"] / jk["Sign"]
    end

    jk["Binder Ratio"] = jk["Magnetization^4"] / (jk["Magnetization^2"]^2)
    jk["Susceptibility"] = (nsites*beta)*jk["Magnetization^2"]
    jk["Connected Susceptibility"] = (nsites*beta)*(jk["Magnetization^2"] - jk["|Magnetization|"]^2)
    jk["Specific Heat"] = (nsites*beta*beta)*(jk["Energy^2"] - jk["Energy"]^2)
    return jk
end


text:
"""
    add_backward_links!(g, w)

Add a backward link to each directed link in `g` with weight `w`.
"""
function add_backward_links!(g::SimpleWeightedDiGraph, w::Real)
    for e in collect(edges(g))
        (u, v) = src(e), dst(e)
        if !has_edge(g, v, u)
            add_edge!(g, v, u, w)
        end
    end
    return g
end

add_backward_links(g, w) = add_backward_links!(copy(g), w)

add_backward_links!(w::Real) = g -> add_backward_links!(g,w)
add_backward_links(w::Real) = g -> add_backward_links(g,w)


text:
include("..\\libs\\FunctionsRobot.jl")

function countPartition!(r)
    back_path = BackPath(r, (Sud, West))
    print(counterPartition!(r))
    BackPath(r, (Sud, West))
    back!(r, back_path)
end


text:
using Plots
using ImageMagick
using FileIO
using Images
using ImageView
using FFTW

function convert_to_float(a::RGB)
    a = float(a)
    return (a.r + a.g + a.b)/3.0
end

function convert_to_float(a::RGBA)
    a = float(a)
    return (a.r + a.g + a.b + a.alpha)*0.25
end

function convert_to_rgba(r, i)
    return RGBA(r, 0, i, 1)
end

function normalize!(a)
    factor = 1/(sum(abs2.(a)))
    a .*= factor
end

function square_animation()
    x = zeros(100)
    x[1:5] .= 1
    x[end-5:end] .= 1
    for i = 1:length(x)-11
        x[i+5] = 1
        x[i] = 0
        if i < 6
            x[end-6+i] = 0
        end
        filename = "output/time"*lpad(i, 5, "0")*".png"
        time_plot = plot(x)
        savefig(time_plot, filename)

        y = fftshift(fft(x))
        filename = "output/frequency"*lpad(i, 5, "0")*".png"
        freq_plot = plot(real(y))
        freq_plot = plot!(freq_plot, imag(y))
        savefig(freq_plot, filename)

    end
end

function fourier_challenge(input_filename::String, output_filename::String)
    img = load(input_filename)
    img_real = float(red.(img))
    img_imag = float(blue.(img))

    array = img_real + im*img_imag
    array = fftshift(ifft(array))
    #array = ifft(array)
    array = abs2.(array)
    normalize!(array)
    img_final = convert_to_rgba.(real(array), imag(array))
    save(output_filename, map(clamp01nan, img_final))
end


text:
mutable struct WT_Image
    kind::Int32 #bytes per pixel
    width::Int32
    height::Int32
    text::Cstring #NULL for TIFF
    array::Ptr{UInt8} #data
end

mutable struct Line_Params
    offset::Float32
    angle::Float32
    width::Float32
    score::Float32
end

mutable struct Seed
    xpnt::Int32
    ypnt::Int32
    xdir::Int32
    ydir::Int32
end

mutable struct Whisker2
    id::Int32
    time::Int32
    len::Int32
    x::Ptr{Float32}
    y::Ptr{Float32}
    thick::Ptr{Float32}
    scores::Ptr{Float32}
end

function Whisker2(w::Whisker1)
    id=w.id
    time=w.time
    len=w.len
    x=pointer(w.x)
    y=pointer(w.y)
    thick=pointer(w.thick)
    scores=pointer(w.scores)
    Whisker2(id,time,len,x,y,thick,scores)
end

struct JT_Measurements
    row::Int32
    fid::Int32
    wid::Int32
    state::Int32 #1 if Whisker, 0 if not. Defaults to 0 when measurements are first made
    face_x::Int32
    face_y::Int32
    col_follicle_x::Int32
    col_follicle_y::Int32
    valid_velocity::Int32
    n::Int32
    face_axis::Cuchar
    data::Ptr{Float64}
    velocity::Ptr{Float64}
end

mutable struct JT_Params
    paramMIN_LENPRJ::Int32
    paramMIN_LENSQR::Int32
    paramMIN_LENGTH::Int32
    paramDUPLICATE_THRESHOLD::Float32
    paramFRAME_DELTA::Int32
    paramHALF_SPACE_TUNNELING_MAX_MOVES::Int32
    paramHALF_SPACE_ASSYMETRY_THRESH::Float32
    paramMAX_DELTA_OFFSET::Float32
    paramMAX_DELTA_WIDTH::Float32
    paramMAX_DELTA_ANGLE::Float32
    paramMIN_SIGNAL::Float32
    paramWIDTH_MAX::Float32
    paramWIDTH_MIN::Float32
    paramWIDTH_STEP::Float32
    paramANGLE_STEP::Float32
    paramOFFSET_STEP::Float32
    paramTLEN::Int32
    paramMIN_SIZE::Int32
    paramMIN_LEVEL::Int32
    paramHAT_RADIUS::Float32
    paramSEED_THRESH::Float32
    paramSEED_ACCUM_THRESH::Float32
    paramSEED_ITERATION_THRESH::Float32
    paramSEED_ITERATIONS::Int32
    paramSEED_SIZE_PX::Int32
    paramSEED_ON_GRID_LATTICE_SPACING::Int32
    paramSEED_METHOD::Int32
    paramIDENTITY_SOLVER_SHAPE_NBINS::Int32
    paramIDENTITY_SOLVER_VELOCITY_NBINS::Int32
    paramCOMPARE_IDENTITIES_DISTS_NBINS::Int32
    paramHMM_RECLASSIFY_BASELINE_LOG2::Float32
    paramHMM_RECLASSIFY_VEL_DISTS_NBINS::Int32
    paramHMM_RECLASSIFY_SHP_DISTS_NBINS::Int32
    paramSHOW_PROGRESS_MESSAGES::Cchar
    paramSHOW_DEBUG_MESSAGES::Cchar
end


text:
using DataDeps

register(DataDep(
  "auto-mpg",
  "https://archive.ics.uci.edu/ml/datasets/Auto+MPG",
  "https://archive.ics.uci.edu/ml/machine-learning-databases/auto-mpg/auto-mpg.data",
  "48b830e11feee5572525f8f1691ddb9d38d3d7b7063edcd8fca672c2a5e17d8d",
  post_fetch_method=(path -> begin
    # Inconsistent use of tabs and whitespace delimiting
    # Convert to whitespace
    write(path, replace(read(path, String), '\t' => ' '))

    UCIData.process_dataset(path,
        target_index=1,
        feature_indices=2:8,
        id_indices=9,
        categoric_indices=8,
        delim=' ',
        ignorerepeated=true,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
# ------------ Static map ------------ #

struct StaticMap{T, N, K, B} <: AddressMap{K}
    tree::NamedTuple{T, <: NTuple{N, <: AddressMap{K}}}
    isempty::B
    function StaticMap{K}() where K
        new{(), 0, K, true}(NamedTuple(), Val(true))
    end
    function StaticMap{K}(nt::NamedTuple) where K
        new{keys(nt), length(nt), K, Val{isempty(nt)}}(nt, Val(isempty(nt)))
    end
    function StaticMap{K}(tree::Dict{Symbol, AddressMap{K}}) where K
        nt = (; tree...)
        new{keys(nt), length(tree), K, Val{isempty(tree)}}(nt, Val(isempty(true)))
    end
end
StaticMap(tree::Dict{Any, AddressMap{K}}) where K = StaticMap{K}(tree)
function StaticMap(k::A, v) where A <: Address
    nt = NamedTuple{(k, )}((Value(v), ))
    StaticMap{Value}(nt)
end
function StaticMap(k::A, v::S) where {A <: Address, S <: StaticMap}
    nt = NamedTuple{(k, )}((v, ))
    StaticMap{Value}(nt)
end

@inline get_address_schema(sm::StaticMap{T, N, K, B}) where {T, N, K, B} = T
@inline get_address_schema(sm::Type{StaticMap{T, N, K, B}}) where {T, N, K, B} = T

@inline shallow_iterator(sm::StaticMap) = zip(keys(sm.tree), values(sm.tree))

@inline function get(sm::StaticMap{T, Value}, addr, fallback) where T
    haskey(sm, addr) || return fallback
    return getindex(sm, addr)
end

@inline function get_sub(sm::StaticMap, addr::A) where A <: Address
    haskey(sm.tree, addr) && return getindex(sm.tree, addr)
    Empty()
end
@inline get_sub(sm::StaticMap, addr::Tuple{}) = Empty()

@inline getindex(sm::StaticMap, addrs...) = get_value(get_sub(sm, addrs))

@inline Base.isempty(sm::StaticMap{T, K, Val{false}}) where {T, K, B} = false
@inline Base.isempty(sm::StaticMap{T, K, Val{true}}) where {T, K, B} = true

@inline haskey(sm::StaticMap, addr::A) where A <: Address = haskey(sm.tree, addr) && has_value(get_sub(sm, addr))

@inline has_sub(sm::StaticMap, addr::A) where A <: Address = haskey(sm.tree, addr)

function set_sub(sm::StaticMap{T, K, B}, addr::A, v::AddressMap{<: K}) where {T, K, B, A <: Address}
    d = Dict{Symbol, AddressMap{K}}(map(keys(sm.tree)) do k
                                        k => getindex(sm, k)
                                    end)
    d[addr] = v
    StaticMap(d)
end
@inline set_sub(sm::StaticMap{T, K, B}, addr::Tuple{A}, v::AddressMap{<: K}) where {T, K, B, A <: Address} = set_sub(sm, addr[1], v)
function set_sub(sm::StaticMap{T, K, B}, addr::Tuple, v::AddressMap{<: K}) where {T, K, B}
    hd, tl = addr[1], addr[2 : end]
    d = Dict{Symbol, AddressMap{K}}()
    if !haskey(sm.tree, hd)
        d[hd] = StaticMap{K}()
    end
    set_sub(d[hd], tl, v)
    StaticMap(d)
end

function merge(sel1::StaticMap{T, K, B1},
               sel2::StaticMap{L, K, B2}) where {T, K, B1, L, B2}
    tree = Dict{Any, AddressMap{K}}()
    inter = intersect(keys(sel1.tree), keys(sel2.tree))
    check = false
    for k in setdiff(keys(sel2.tree), keys(sel1.tree))
        tree[k] = sel2.tree[k]
    end
    for k in setdiff(keys(sel1.tree), keys(sel2.tree))
        tree[k] = sel1.tree[k]
    end
    for k in inter
        tree[k], b = merge(sel1.tree[k], sel2.tree[k])
        check = check || b
    end
    return StaticMap{K}((; tree...)), !isempty(inter) || check
end

merge(sm::StaticMap, ::Empty) = deepcopy(sm), false
merge(::Empty, sm::StaticMap) = deepcopy(sm), false

function collect!(par::T, addrs::Vector{Any}, chd::Dict{Any, Any}, tr::S, meta) where {T <: Tuple, S <: StaticMap}
    iterate(tr) do (k, v)
        collect!((par..., k), addrs, chd, v, meta)
    end
end
function collect!(addrs::Vector{Any}, chd::Dict{Any, Any}, tr::S, meta) where S <: StaticMap
    iterate(tr) do (k, v)
        collect!((k, ), addrs, chd, v, meta)
    end
end

function static(k::A, v) where A <: Address
    StaticMap(k, v)
end
function static(k::Tuple{A}, v) where A <: Address
    StaticMap(k[1], v)
end
function static(k::Tuple, v) where A <: Address
    static(k[1], static(k[2 : end], v))
end
function static(v::Vector{Pair{Tuple{Symbol}, K}}) where K
    ks = map(v) do k
        k[1][1]
    end
    vs = map(v) do k
        Value(k[2])
    end
    nt = NamedTuple{tuple(ks...)}(vs)
    StaticMap{Value}(nt)
end
function static(v::Vector{Pair{T, K}}) where {T <: Tuple, K}
    if !isempty(v)
        k, l = v[1]
        sm = static(k, l)
        for (k, v) in v[2 : end]
            sm, _ = merge(sm, static(k, v))
        end
        return sm
    end
end

@inline static(dm::DynamicMap{Value}) = StaticMap(dm.tree)


text:
## model
zfixed   = [-1e6,   0, 40, 60]
rho      = [1e12,   20, 1, 100]
# operator
ntimesperdecade = 10
nfreqsperdecade = 5
nkᵣeval = 60
zTx = -120
zRx = -80
x_rx = -115.
y_rx = 0.
rx_roll = 0.
rx_pitch = 0.
rx_yaw = 0.
tx_roll = 10. #ACHTUNG!
tx_pitch = 10. #ACHTUNG!
tx_yaw = 0.
addprimary = true
tempest = transD_GP.TEMPEST1DInversion.Bfield(
  				     ntimesperdecade = ntimesperdecade,
                 nfreqsperdecade = nfreqsperdecade,
					  nkᵣeval = nkᵣeval,
                 times  = times,
                 ramp   = ramp,
					  zTx    = zTx,
					  zRx    = zRx,
					  x_rx   = x_rx,
					  y_rx   = y_rx,
					  rx_roll = rx_roll,
					  rx_pitch = rx_pitch,
					  rx_yaw = rx_yaw,
					  tx_roll = tx_roll,
					  tx_pitch = tx_pitch,
					  tx_yaw = tx_yaw,
					  addprimary = addprimary)
# X and Z components in fT from GA-AEM
ross_primary = [33.0699      -13.689]
ross_secondary = [
6.148033      -8.96926
   4.254406     -7.237281
   3.515803     -6.429036
   3.037861     -5.838994
   2.680435      -5.36655
   2.369648     -4.949131
   2.011034      -4.45433
    1.59146     -3.831225
   1.147514     -3.094577
  0.7279429     -2.289014
  0.3895116      -1.50605
  0.1698442    -0.8633888
0.05873083    -0.4281239
0.01569159    -0.1862741
0.002626163     -0.069275]

if addprimary
	ross_secondary .+= ross_primary
end

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
import AbstractTrees
import AbstractTrees.children
import AbstractTrees.printnode
using AbstractTrees: print_tree
using TSML

# From  https://github.com/JuliaOpt/Convex.jl/blob/master/src/utilities/tree_print.jl
# Plot it with: https://github.com/JuliaPlots/GraphRecipes.jl#abstracttrees-trees

# extennd abstract trees
children(d::Dict) = [p for p in d]
children(p::Pair) = AbstractTrees.children(p[2])

function printnode(io::IO, p::Pair)
  v = AbstractTrees.children(p[2])
  str = isempty(v)  ? string(p[1], ": ", p[2]) : string(p[1], ": ")
  print(io, str)
end

function showtree(p::Transformer)
  if isnothing(p.model)
    return
  elseif :transformers in keys(p.model) # pipeline
    println(:pipeline)
    print_tree(p.model[:transformers])
  else
    println(:transfomer)
    print_tree(p.model) # any transformer except pipeline
  end
end

function test_print_tree()
  there = Dict(:ok=>:yes)
  bb=Dict(:hello=>there)
  a=Dict(:b => bb, :d=>:e)
  d = Dict(:a => 2,:d => Dict(:b => 4,:c => "Hello"),:e => 5.0)
  print_tree(a)
  print_tree(d)
end


text:
function forloop(startNumber,endNumber,increment)
    i = startNumber;
    while (i < endNumber)
        println("i = ", i);
        i= i+increment;
    end;
end;

forloop(0,5,1);

max_stars_repo_name: fahimfarhan/legendary-coding-odyssey

max_issues_repo_name: fahimfarhan/legendary-coding-odyssey

max_forks_repo_name: fahimfarhan/legendary-coding-odyssey


text:
#==========================================================================================#

# TYPE

mutable struct Poisson <: MLE

    sample::Microdata
    β::Vector{Float64}
    V::AbstractMatrix{Float64}

    Poisson() = new()
end

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

# CONSTRUCTOR

function Poisson(MD::Microdata)
    obj        = Poisson()
    obj.sample = MD
    return obj
end

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

# ESTIMATION

function _fit!(obj::Poisson, ::UnitWeights)

    has_offset = haskey(obj.sample.mapping, :offset)

    y  = getvector(obj, :response)
    x  = getmatrix(obj, :control)
    μ  = Array{Float64}(undef, length(y))
    r  = Array{Float64}(undef, length(y))
    xx = Array{Float64}(undef, size(x)...)

    β₀ = vcat(fill(0.0, size(x, 2) - 1), log(mean(y)))

    has_offset && (xo = getmatrix(obj, :offset, :control))

    function L(β::Vector)
        has_offset ? mul!(μ, xo, vcat(1.0, β)) : mul!(μ, x, β)
        return sum(exp.(μ) .- y .* μ)
    end

    function G!(g::Vector, β::Vector)
        has_offset ? mul!(μ, xo, vcat(1.0, β)) : mul!(μ, x, β)
        r .= exp.(μ) .- y
        mul!(g, transpose(x), r)
    end

    function LG!(g::Vector, β::Vector)

        has_offset ? mul!(μ, xo, vcat(1.0, β)) : mul!(μ, x, β)
        ll = 0.0

        @inbounds for (i, (yi, μi)) in enumerate(zip(y, μ))
            ηi   = exp(μi)
            r[i] = ηi - yi
            ll  += ηi - yi * μi
        end

        mul!(g, transpose(x), r)

        return ll
    end

    function H!(h::AbstractMatrix, β::Vector)
        has_offset ? mul!(μ, xo, vcat(1.0, β)) : mul!(μ, x, β)
        xx .= x .* exp.(μ)
        mul!(h, transpose(x), xx)
    end

    res = optimize(TwiceDifferentiable(L, G!, LG!, H!, β₀), β₀, Newton())

    if Optim.converged(res)
        obj.β = Optim.minimizer(res)
    else
        throw("likelihood maximization did not converge")
    end
end

function _fit!(obj::Poisson, w::AbstractWeights)

    has_offset = haskey(obj.sample.mapping, :offset)

    y  = getvector(obj, :response)
    x  = getmatrix(obj, :control)
    μ  = Array{Float64}(undef, length(y))
    r  = Array{Float64}(undef, length(y))
    xx = Array{Float64}(undef, size(x)...)

    β₀ = vcat(fill(0.0, size(x, 2) - 1), log(mean(y, w)))

    has_offset && (xo = getmatrix(obj, :offset, :control))

    function L(β::Vector)
        has_offset ? mul!(μ, xo, vcat(1.0, β)) : mul!(μ, x, β)
        return sum(exp.(μ) .- y .* μ, w)
    end

    function G!(g::Vector, β::Vector)
        has_offset ? mul!(μ, xo, vcat(1.0, β)) : mul!(μ, x, β)
        r .= w .* (exp.(μ) .- y)
        mul!(g, transpose(x), r)
    end

    function LG!(g::Vector, β::Vector)

        has_offset ? mul!(μ, xo, vcat(1.0, β)) : mul!(μ, x, β)
        ll = 0.0

        @inbounds for (i, (yi, μi, wi)) in enumerate(zip(y, μ, w))
            ηi   = exp(μi)
            r[i] = wi * (ηi - yi)
            ll  += wi * (ηi - yi * μi)
        end

        mul!(g, transpose(x), r)

        return ll
    end

    function H!(h::AbstractMatrix, β::Vector)
        has_offset ? mul!(μ, xo, vcat(1.0, β)) : mul!(μ, x, β)
        r  .= w .* exp.(μi)
        xx .= x .* r
        mul!(h, transpose(x), xx)
    end

    res = optimize(TwiceDifferentiable(L, G!, LG!, H!, β₀), β₀, Newton())

    if Optim.converged(res)
        obj.β = Optim.minimizer(res)
    else
        throw("likelihood maximization did not converge")
    end
end

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

# SCORE (DERIVATIVE OF THE LIKELIHOOD FUNCTION)

score(obj::Poisson) = Diagonal(residuals(obj)) * getmatrix(obj, :control)

# EXPECTED JACOBIAN OF SCORE × NUMBER OF OBSERVATIONS

function jacobian(obj::Poisson, ::UnitWeights)

    x = getmatrix(obj, :control)

    if haskey(obj.sample.mapping, :offset)
        xo = getmatrix(obj, :offset, :control)
        v  = xo * vcat(1.0, obj.β)
    else
        v  = x * obj.β
    end

    v .= exp.(v)

    return - crossprod(x, v)
end

function jacobian(obj::Poisson, w::AbstractWeights)

    x = getmatrix(obj, :control)

    if haskey(obj.sample.mapping, :offset)
        xo = getmatrix(obj, :offset, :control)
        v  = xo * vcat(1.0, obj.β)
    else
        v  = x * obj.β
    end

    v .= w .* exp.(v)

    return - crossprod(x, v)
end

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

# LINEAR PREDICTOR

function linear_predictor(obj::Poisson, MD::Microdata)

    (getnames(obj, :control) != getnames(MD, :control)) && throw("missing variables")

    if haskey(obj.sample.mapping, :offset)
        return getmatrix(MD, :offset, :control) * vcat(1.0, obj.β)
    else
        return getmatrix(MD, :control) * obj.β
    end
end

# FITTED VALUES

predict(obj::Poisson, MD::Microdata) = exp.(linear_predictor(obj, MD))

# DERIVATIVE OF FITTED VALUES

function jacobexp(obj::Poisson)

    if haskey(obj.sample.mapping, :offset)
        v = getmatrix(obj, :offset, :control) * vcat(1.0, obj.β)
    else
        v = getmatrix(obj, :control) * obj.β
    end

    v .= exp.(v)

    return Diag(v) * getmatrix(obj, :control)
end

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

# UTILITIES

coefnames(obj::Poisson) = getnames(obj, :control)

# LIKELIHOOD FUNCTION

function _loglikelihood(obj::Poisson, w::AbstractWeights)
    y = getvector(obj, :response)
    μ = linear_predictor(obj)
    return sum(y .* μ .- exp.(μ) .- loggamma.(1.0 .+ y), w)
end

# LIKELIHOOD FUNCTION UNDER NULL MODEL

function _nullloglikelihood(obj::Poisson, w::AbstractWeights)
    if haskey(obj.sample.mapping, :offset)
        y = getvector(obj, :response)
        o = getmatrix(obj, :offset)
        β = log(mean(y, w) / mean(exp.(o), w))
        return sum(y .* (o .+ β) .- exp.(o .+ β) .- loggamma.(1.0 .+ y), w)
    else
        y = getvector(obj, :response)
        η = mean(y, w)
        μ = log(η)
        return sum(y .* μ .- η .- loggamma.(1.0 .+ y), w)
    end
end

# DEVIANCE

function _deviance(obj::Poisson, w::AbstractWeights)
    y = getvector(obj, :response)
    μ = linear_predictor(obj)
    return 2 * sum(xlogx.(y) .- y .* (1.0 .+ μ) .+ exp.(μ), w)
end

# DEVIANCE UNDER NULL MODEL

function _nulldeviance(obj::Poisson, w::AbstractWeights)
    y = getvector(obj, :response)
    μ = mean(y, w)
    return 2 * sum(xlogx.(y) .- y .* (1.0 .+ log(μ)) .+ μ, w)
end

max_stars_repo_name: lbittarello/Microeconometrics

max_issues_repo_name: lbittarello/Microeconometrics

max_forks_repo_name: lbittarello/Microeconometrics


text:
@testset "mono_mode" begin
  stats1 = ripqp(QuadraticModel(qps1), iconf = InputConfig(presolve = false))
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  stats2 = ripqp(QuadraticModel(qps2))
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats1.status == :acceptable

  stats3 = ripqp(QuadraticModel(qps3))
  @test isapprox(stats3.objective, 5.32664756, atol = 1e-2)
  @test stats3.status == :acceptable
end

@testset "K2_5" begin
  stats1 = ripqp(QuadraticModel(qps1), display = false, iconf = InputConfig(sp = K2_5LDLParams()))
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  stats2 = ripqp(
    QuadraticModel(qps2),
    display = false,
    iconf = InputConfig(sp = K2_5LDLParams(), mode = :multi),
  )
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable

  stats3 = ripqp(
    QuadraticModel(qps3),
    display = false,
    iconf = InputConfig(sp = K2_5LDLParams(regul = :dynamic)),
  )
  @test isapprox(stats3.objective, 5.32664756, atol = 1e-2)
  @test stats3.status == :acceptable
end

@testset "KrylovK2" begin
  for precond in [:Identity, :Jacobi]
    for kmethod in [:minres, :minres_qlp, :dqgmres]
      stats2 = ripqp(
        QuadraticModel(qps2),
        display = false,
        iconf = InputConfig(
          sp = K2KrylovParams(uplo = :L, kmethod = kmethod, preconditioner = precond),
          history = true,
        ),
        itol = InputTol(
          max_iter = 50,
          max_time = 20.0,
          ϵ_rc = 1.0e-2,
          ϵ_rb = 1.0e-2,
          ϵ_pdd = 1.0e-2,
        ),
      )
      @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-1)
      @test stats2.status == :acceptable

      stats3 = ripqp(
        QuadraticModel(qps3),
        display = false,
        iconf = InputConfig(
          sp = K2KrylovParams(uplo = :U, kmethod = kmethod, preconditioner = precond),
        ),
        itol = InputTol(
          max_iter = 50,
          max_time = 20.0,
          ϵ_rc = 1.0e-2,
          ϵ_rb = 1.0e-2,
          ϵ_pdd = 1.0e-2,
        ),
      )
      @test isapprox(stats3.objective, 5.32664756, atol = 1e-1)
      @test stats3.status == :acceptable
    end
  end
end

@testset "KrylovK2_5" begin
  for kmethod in [:minres, :minres_qlp]
    stats1 = ripqp(
      QuadraticModel(qps1),
      display = true,
      iconf = InputConfig(
        sp = K2_5KrylovParams(kmethod = kmethod, preconditioner = :Identity),
        history = true,
      ),
      itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-2, ϵ_rb = 1.0e-2, ϵ_pdd = 1.0e-2),
    )
    @test isapprox(stats1.objective, -1.59078179, atol = 1e-1)
    @test stats1.status == :acceptable

    stats2 = ripqp(
      QuadraticModel(qps2),
      display = false,
      iconf = InputConfig(
        sp = K2_5KrylovParams(uplo = :U, kmethod = kmethod, preconditioner = :Jacobi),
        solve_method = :IPF,
      ),
      itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-2, ϵ_rb = 1.0e-2, ϵ_pdd = 1.0e-2),
    )
    @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-1)
    @test stats2.status == :acceptable

    stats3 = ripqp(
      QuadraticModel(qps3),
      display = false,
      iconf = InputConfig(sp = K2_5KrylovParams(kmethod = kmethod, preconditioner = :Jacobi)),
      itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-2, ϵ_rb = 1.0e-2, ϵ_pdd = 1.0e-2),
    )
    @test isapprox(stats3.objective, 5.32664756, atol = 1e-1)
    @test stats3.status == :acceptable
  end
end

@testset "K2 structured LP" begin
  for kmethod in [:tricg, :trimr, :gpmr]
    stats4 = ripqp(
      QuadraticModel(qps4),
      display = false,
      iconf = InputConfig(sp = K2StructuredParams(kmethod = kmethod), solve_method = :IPF),
    )
    @test isapprox(stats4.objective, -4.6475314286e02, atol = 1e-2)
    @test stats4.status == :acceptable
  end
end

@testset "K2.5 structured LP" begin
  for kmethod in [:tricg, :trimr, :gpmr]
    stats4 = ripqp(
      QuadraticModel(qps4),
      display = false,
      iconf = InputConfig(sp = K2_5StructuredParams(kmethod = kmethod), solve_method = :IPF),
    )
    @test isapprox(stats4.objective, -4.6475314286e02, atol = 1e-2)
    @test stats4.status == :acceptable
  end
end


text:
export MtlFunction

const MTLFunction = Ptr{MtFunction}

mutable struct MtlFunction
    handle::MTLFunction
    lib::MtlLibrary

    "Get a handle to a kernel function in a Metal Library."
    function MtlFunction(lib::MtlLibrary, name::String)
        handle = mtNewFunctionWithName(lib, name)
        handle == C_NULL && throw(KeyError(name))
        obj = new(handle, lib)
        finalizer(unsafe_destroy!, obj)
        return obj
    end
end

function unsafe_destroy!(fun::MtlFunction)
    mtRelease(fun.handle)
end

Base.unsafe_convert(::Type{MTLFunction}, fun::MtlFunction) = fun.handle

Base.:(==)(a::MtlFunction, b::MtlFunction) = a.handle == b.handle
Base.hash(fun::MtlFunction, h::UInt) = hash(mod.handle, h)

## properties

Base.propertynames(::MtlFunction) = (:device, :label, :name, :functionType)

function Base.getproperty(fun::MtlFunction, f::Symbol)
    if f === :device
        MtlDevice(mtFunctionDevice(fun))
    elseif f == :label
        ptr = mtFunctionLabel(fun)
        ptr == C_NULL ? nothing : unsafe_string(ptr)
    elseif f == :name
        unsafe_string(mtFunctionName(fun))
    elseif f == :functionType
        mtFunctionType(fun)
    else
        getfield(fun, f)
    end
end

## display

function Base.show(io::IO, ::MIME"text/plain", fun::MtlFunction)
    println(io, "MtlFunction:")
    println(io, " name:   ", fun.name)
    println(io, " type:   ", fun.functionType)
    println(io, " device: ", fun.device)
    print(io,   " label:  ", fun.label)
end

function Base.show(io::IO, fun::MtlFunction)
    print(io, "MtlFunction($(fun.name))")
end


text:
@doc raw"""
    HomogeneousPoissonPointProcess{T<:Vector{Float64}} <: AbstractSpatialPointProcess{T}

Homegeneous [Poisson point process](https://en.wikipedia.org/wiki/Poisson_point_process) with intensity ``\beta > 0``, denoted ``\operatorname{Poisson}(\beta)``.

``\operatorname{Poisson}(\beta)`` has density (w.r.t. the homogenous Poisson point process with unit intensity ``\operatorname{Poisson}(1)``) proportional to

```math
    \prod_{x \in X}
        \beta.
```
"""
struct HomogeneousPoissonPointProcess{T<:Vector{Float64}} <: AbstractSpatialPointProcess{T}
    β::Float64
    window::AbstractSpatialWindow{Float64}
end

function Base.show(io::IO, pp::HomogeneousPoissonPointProcess{T}) where {T}
    print(io, "HomogeneousPoissonPointProcess{$T}\n- β = $(pp.β)\n- window = $(pp.window)")
end

@doc raw"""
    HomogeneousPoissonPointProcess(
        β::Real,
        window::Union{Nothin,AbstractSpatialWindow}=nothing
    )

Construct a [`PRS.HomogeneousPoissonPointProcess`](@ref) with intensity `β` restricted to `window`.

Default window (`window=nothing`) is [`PRS.SquareWindow`](@ref)`()`.

```jldoctest; output = true
using PartialRejectionSampling

β = 40
win = PRS.SquareWindow(zeros(2), 1)
PRS.HomogeneousPoissonPointProcess(β, win)

# output

HomogeneousPoissonPointProcess{Array{Float64,1}}
- β = 40.0
- window = SquareWindow [0.0, 1.0]^2
```
"""
function HomogeneousPoissonPointProcess(
    β::Real,
    window::Union{Nothing,AbstractSpatialWindow}=nothing
)
    @assert β > 0
    win = isnothing(window) ? SquareWindow() : window
    return HomogeneousPoissonPointProcess{Vector{Float64}}(β, win)
end

"""
    intensity(pp::HomogeneousPoissonPointProcess) = pp.β
"""
intensity(pp::HomogeneousPoissonPointProcess) = pp.β

## Sampling

"""
    generate_sample(
        [rng::Random.AbstractRNG,]
        pp::HomogeneousPoissonPointProcess{Vector{T}},
        win::Union{Nothing,AbstractWindow}=nothing
    )::Matrix{T} where {T<:Float64}

Generate an exact sample from [`PRS.HomogeneousPoissonPointProcess`](@ref) on window `win`.
Sampled points are stored as columns of the output matrix.

Default window (`win=nothing`) is `window(pp)`.
"""
function generate_sample(
    rng::Random.AbstractRNG,
    pp::HomogeneousPoissonPointProcess{Vector{T}},
    win::Union{Nothing,AbstractWindow}=nothing
)::Matrix{T} where {T<:Float64}
    window_ = isnothing(win) ? window(pp) : win
    n = rand(rng, Distributions.Poisson(pp.β * volume(window_)))
    return rand(rng, window_, n)
end

# function generate_sample(
#     pp::HomogeneousPoissonPointProcess,
#     win::Union{Nothing,AbstractWindow}=nothing
# )
#     return generate_sample(Random.default_rng(), pp, win=win)
# end

@doc raw"""
    generate_sample_poisson_union_balls(
        [rng::Random.AbstractRNG,]
        β::Real,
        centers::Matrix,
        radius::Real,
        win::Union{Nothing,AbstractWindow}=nothing
    )::Matrix{Float64}

Generate an exact sample from a homogenous [`PRS.HomogeneousPoissonPointProcess`](@ref) Poisson(β) on ``\bigcup_{i} B(c_i, r)`` (union of balls centered at ``c_i`` with the same radius ``r``).

If `win ≂̸ nothing`, returns the points falling in `win`.

!!! hint

    Use the independence property of the Poisson point process on disjoint subsets in order to

    - Sample from Poisson(β) on ``B(c_1, r)``,
    - Sample from Poisson(β) on ``B(c_2, r) \setminus B(c_1, r)``,
    - Sample from Poisson(β) on ``B(c_j, r) \setminus \bigcup_{i<j} B(c_i, r)``,
    - ...
"""
function generate_sample_poisson_union_balls(
    rng::Random.AbstractRNG,
    β::Real,
    centers::Matrix,
    radius::Real,
    win::Union{Nothing,AbstractWindow}=nothing
)::Matrix{Float64}
    d = size(centers, 1)
    !(isnothing(win)) && @assert dimension(win) == d

    𝒫 = Distributions.Poisson(β * volume(BallWindow(zeros(d), radius)))
    points = Matrix{Float64}(undef, d, 0)
    for (i, c) in enumerate(eachcol(centers))
        n = rand(rng, 𝒫)
        n == 0 && continue
        proposed = rand(rng, BallWindow(c, radius), n)
        if i == 1
            points = hcat(points, proposed)
            continue
        end
        centers_ = @view centers[:, 1:i-1]
        accept = vec(all(pairwise_distances(centers_, proposed) .> radius, dims=1))
        points = hcat(points, proposed[:, accept])
    end

    if isnothing(win) || isempty(points)
        return points
    end
    return points[:, vec(mapslices(x -> x in win, points; dims=1))]
end

function generate_sample_poisson_union_balls(
    β::Real,
    centers::Matrix,
    radius::Real,
    win::Union{Nothing,AbstractWindow}=nothing
)
    generate_sample_poisson_union_balls(Random.default_rng(), β, centers, radius, win)
end


text:
using ElectronDisplay
ElectronDisplay.CONFIG.showable = showable
ElectronDisplay.CONFIG.single_window = true
ElectronDisplay.CONFIG.focus = false
using ImageMagick, FileIO
using WGLMakie, Makie, Test
using Pkg
path = normpath(joinpath(dirname(pathof(Makie)), "..", "ReferenceTests"))
Pkg.develop(PackageSpec(path = path))
WGLMakie.activate!()
using ReferenceTests
using ReferenceTests: database_filtered

excludes = Set([
    "Streamplot animation",
    "Transforming lines",
    "image scatter",
    "Line GIF",
    "surface + contour3d",
    # Hm weird, looks like some internal JSServe error missing an Observable:
    "Errorbars x y low high",
    "Rangebars x y low high",
    # These are a bit sad, since it's just missing interpolations
    "FEM mesh 2D",
    "FEM polygon 2D",
    # missing transparency & image
    "Wireframe of a Surface",
    "Image on Surface Sphere",
    "Surface with image",
    # Marker size seems wrong in some occasions:
    "Hbox",
    "UnicodeMarker",
    # Not sure, looks pretty similar to me! Maybe blend mode?
    "Test heatmap + image overlap",
    "Stars",
    "heatmaps & surface",
    "OldAxis + Surface",
    "Order Independent Transparency",
    "Record Video",
    "fast pixel marker",
    "Animated surface and wireframe"
])

database = database_filtered(excludes)

basefolder = joinpath(@__DIR__, "reference_test_output")
rm(basefolder; force=true, recursive=true)
mkdir(basefolder)

refimage_set = "refimages"
tests_root_folder = joinpath(basefolder, refimage_set)
mkdir(tests_root_folder)

tests_record_folder = joinpath(tests_root_folder, "recorded")
mkdir(tests_record_folder)

ReferenceTests.record_tests(database, recording_dir = tests_record_folder)

tests_refimages_download_folder = ReferenceTests.download_refimages(; name=refimage_set)
tests_refimages_folder = joinpath(tests_root_folder, "reference")
cp(tests_refimages_download_folder, tests_refimages_folder)

missing_refimages, scores = ReferenceTests.record_comparison(tests_root_folder)
ReferenceTests.test_comparison(missing_refimages, scores; threshold = 0.032)


text:
using ScriptIt

#using Khepri
#backend(autocad)
#traceability(true)
#excluded_modules([Base, Khepri, Atom, Atom.CodeTools, Base.CoreLogging,ScriptIt])
#ScriptIt.startScript()

max_stars_repo_name: Daniel-Fermoselle/Script-It-Julia

max_issues_repo_name: Daniel-Fermoselle/Script-It-Julia

max_forks_repo_name: Daniel-Fermoselle/Script-It-Julia


text:
using ImageBase
using ImageBase: varmult
using Statistics
using Test

@testset "Reductions" begin
    _abs(x::Colorant) = mapreducec(abs, +, 0, x)

    @testset "sumfinite, meanfinite, varfinite" begin
        for T in generate_test_types([N0f8, Float32], [Gray, RGB])
            A = rand(T, 5, 5)
            s12 = sum(A, dims=(1,2))
            @test eltype(s12) <: Union{T, float(T), float64(T)}

            @test sumfinite(A) ≈ sum(A)
            @test sumfinite(A, dims=1) ≈ sum(A, dims=1)
            @test sumfinite(A, dims=(1, 2)) ≈ sum(A, dims=(1, 2))

            @test meanfinite(A) ≈ mean(A)
            @test meanfinite(A, dims=1) ≈ mean(A, dims=1)
            @test meanfinite(A, dims=(1, 2)) ≈ mean(A, dims=(1, 2))

            @test varfinite(A) ≈ varmult(⋅, A)
            @test varfinite(A, dims=1) ≈ varmult(⋅, A, dims=1)
            @test varfinite(A, dims=(1, 2)) ≈ varmult(⋅, A, dims=(1, 2))

            # test NaN/Inf
            if eltype(T) != N0f8
                A = rand(T, 5, 5) .- 0.5 .* oneunit(T)
                A[1] = Inf
                @test sum(A) ≈ A[1]
                @test sum(abs, A) ≈ A[1]
                @test sumfinite(A) ≈ sum(A[2:end])
                @test sumfinite(abs, A) ≈ sum(abs, A[2:end])
                A[1] = NaN
                @test isnan(sum(A))
                @test isnan(sum(abs, A))
                @test sumfinite(A) ≈ sum(A[2:end])
                @test sumfinite(abs, A) ≈ sum(abs, A[2:end])

                A = rand(T, 5, 5) .- 0.5 .* oneunit(T)
                A[1] = Inf
                @test mean(A) ≈ A[1]
                @test mean(abs, A) ≈ A[1]
                @test meanfinite(A) ≈ mean(A[2:end])
                @test meanfinite(abs, A) ≈ mean(abs, A[2:end])
                A[1] = NaN
                @test isnan(mean(A))
                @test isnan(mean(abs, A))
                @test meanfinite(A) ≈ mean(A[2:end])
                @test meanfinite(abs, A) ≈ mean(abs, A[2:end])

                A = rand(T, 5, 5)
                A[1] = Inf
                @test isnan(varmult(⋅, A))
                @test varfinite(A) ≈ varmult(⋅, A[2:end])
                A[1] = NaN
                @test isnan(varmult(⋅, A))
                @test varfinite(A) ≈ varmult(⋅, A[2:end])
            end
        end

        A = [NaN, 1, 2, 3]
        @test meanfinite(A, dims=1) ≈ [2]
        @test varfinite(A, dims=1) ≈ [1]

        A = [NaN NaN 1;
            1 2 3]
        vf = varfinite(A, dims=2)
        @test isnan(vf[1])

        A = [NaN 1 2 3;
            NaN 6 5 4]
        mf = meanfinite(A, dims=1)
        vf = varfinite(A, dims=1)
        @test isnan(mf[1])
        @test mf[2:end] ≈ [3.5,3.5,3.5]
        @test isnan(vf[1])
        @test vf[2:end] ≈ [12.5,4.5,0.5]

        @test meanfinite(A, dims=2) ≈ reshape([2, 5], 2, 1)
        @test varfinite(A, dims=2) ≈ reshape([1, 1], 2, 1)

        @test meanfinite(A, dims=(1,2)) ≈ [3.5]
        @test varfinite(A, dims=(1,2)) ≈ [3.5]

        # Ensure we're consistant with our decision to `abs2` in ColorVectorSpace
        # See also: https://github.com/JuliaGraphics/ColorVectorSpace.jl/blob/master/README.md#abs-and-abs2
        A = rand(Gray, 4, 4)
        @test varfinite(A) ≈ varfinite(RGB.(A))
        A[1] = Inf
        @test varfinite(A) ≈ varfinite(RGB.(A))
        A[1] = NaN
        @test varfinite(A) ≈ varfinite(RGB.(A))
    end

    @testset "minfinite, maxfinite, maxabsfinite" begin
        for T in generate_test_types([N0f8, Float32], [Gray, ])
            A = rand(T, 5, 5)
            @test @inferred(minimum_finite(A)) == minimum(A)
            @test @inferred(maximum_finite(A)) == maximum(A)
            @test minimum_finite(A; dims=1) == minimum(A; dims=1)
            @test maximum_finite(A; dims=1) == maximum(A; dims=1)
            @test_broken @inferred maximum_finite(A; dims=1)
            @test_broken @inferred minimum_finite(A; dims=1)

            @test maximum_finite(abs2, A) == maximum(abs2, A)
            @test_broken @inferred maximum(abs2, A)
            @test_broken @inferred minimum(abs2, A)

            if eltype(T) != N0f8
                A = rand(T, 5, 5) .- 0.5 * rand(T, 5, 5)
                A[1] = Inf

                @test @inferred(minimum_finite(A)) == minimum(A[2:end])
                @test @inferred(maximum_finite(A)) == maximum(A[2:end])
                @test minimum_finite(abs, A) == minimum(abs, A[2:end])
                @test maximum_finite(abs2, A) == maximum(abs2, A[2:end])
            end
        end

        # minimum_finite and maximum_finite for RGB are processed per channel
        A = rand(RGB{Float32}, 5, 5)
        @test minimum_finite(A) == RGB(minimum(channelview(A), dims=(2, 3))...)
        @test maximum_finite(A) == RGB(maximum(channelview(A), dims=(2, 3))...)

        # Container of abstract type
        A = Any[1, 2, Gray(0.3)]
        @test minimum_finite(A) == 0.3
        @test maximum_finite(A) == 2.0
        @test Base.return_types(minimum_finite, (typeof(A),)) == Base.return_types(minimum, (typeof(A), ))
    end

end


text:
"""
There are also identification robust versions of likelihood ratio and
lagrange multiplier test. Moreire (2003)[@moreira2003] proposed a conditional
likelihood ratio (CLR) test for weakly identified linear IV
models. Kleibergen (2005)[@kleibergen2005] developed a Lagrange multiplier (often called
the KLM) test and extended Moreira's CLR test to weakly identified GMM
models.  More recently, Andrews and Guggenberge (2015) [@andrews2015]
and Andrews and Guggenberge (2017) [@andrews2017] showed the
validity of these tests under more general conditions. These tests are
somewhat more complicated than the AR test, but they have the
advantage that they are often more powerful. The AR test statistic has
a \\chi^2_{m} distribution, where m is the number of moments. The
CLR and KLM statistics under strong identification have \\chi^2_k
distributions (as does the Wald statistic), where k is the number of
parameters. Consequently, when the model is overidentified, the CLR
and LM tests are more powerful than the AR test.

Here is an implementation of the KLM and CLR statistics. The names of
variables roughly follows the notation of Andrews and
Guggenberge(2017

"""

"""
    statparts(gi::Function)

Return a function.

compute common components of klm, rk, & clr stats
follows notation of Andrews & Guggenberger 2017, section 3.1
"""
function statparts(gi::Function)
  # compute common components of klm, rk, & clr stats
  # follows notation of Andrews & Guggenberger 2017, section 3.1
    function P(A::AbstractMatrix) # projection matrix
    A*pinv(A'*A)*A'
    end
    function(θ)
        giθ = gi(θ)
        p = length(θ)
        (n, k) = size(giθ)
        Ω = cov(giθ)
        gn=mean(gi(θ), dims=1)'
        #G = ForwardDiff.jacobian(θ->mean(gi(θ),dims=1),θ)
        Gi= ForwardDiff.jacobian(gi,θ)
        Gi = reshape(Gi, n , k, p)
        G = mean(Gi, dims=1)
        Γ = zeros(eltype(G),p,k,k)
        D = zeros(eltype(G),k, p)
        for j in 1:p
        for i in 1:n
          Γ[j,:,:] += (Gi[i,:,j] .- G[1,:,j]) * giθ[i,:]'
        end
        Γ[j,:,:] ./= n
        D[:,j] = G[1,:,j] - Γ[j,:,:]*inv(Ω)*gn
        end
        return(n,k,p,gn, Ω, D, P)
    end
end

"""
    klm(gi::Function)

Return test statistic for KLM method.
"""

function klm(gi::Function)
    SP = statparts(gi)
    function(θ)
        (n,k,p,gn, Ω, D, P) = SP(θ)
        return n*(gn'*Ω^(-1/2)*P(Ω^(-1/2)*D)*Ω^(-1/2)*gn)[1]
    end
end

"""
    clr(gi::Function)

Return test statistic for CLR method.
"""
function clr(gi::Function)
  SP = statparts(gi)
  function(θ)
        (n,k,p,gn, Ω, D, P) = SP(θ)
        rk = eigmin(n*D'*inv(Ω)*D)
        AR  = (n*gn'*inv(Ω)*gn)[1]
        lm = (n*(gn'*Ω^(-1/2)*P(Ω^(-1/2)*D)*Ω^(-1/2)*gn))[1]
        lr = 1/2*(AR - rk + sqrt( (AR-rk)^2 + 4*lm*rk))

        # simulate to find p-value
        S = 5000
        function randc(k,p,r,S)
        χp = rand(Chisq(p),S)
        χkp = rand(Chisq(k-p),S)
        0.5.*(χp .+ χkp .- r .+
              sqrt.((χp .+ χkp .- r).^2 .+ 4 .* χp.*r))
        end
        csim = randc(k,p,rk,S)
        return mean(csim.<=lr)
    end
end


text:
abstract type LinkedList{T} end

Base.eltype(::Type{<:LinkedList{T}}) where T = T

mutable struct Cons{T} <: LinkedList{T}
    head::T
    tail::Cons{T}
    Cons{T}(h, t) where T = new{T}(h, t)
    Cons{T}() where T = new{T}()
end

cons(h, t::LinkedList{T}) where {T} = Cons{T}(h, t)

nil(T) = Cons{T}()
nil() = nil(Any)
isnil(x::Cons) = !isdefined(x, :tail)

head(x::Cons) = x.head
tail(x::Cons) = x.tail

Base.:(==)(x::Cons, y::Cons) = begin
    isnil(x) && isnil(y) && return true
    return !isnil(x) && !isnil(y) && (x.head == y.head) && (x.tail == y.tail)
end

function Base.show(io::IO, l::LinkedList{T}) where T
    if isnil(l)
        if T === Any
            print(io, "nil()")
        else
            print(io, "nil(", T, ")")
        end
    else
        print(io, "list(")
        show(io, head(l))
        for t in tail(l)
            print(io, ", ")
            show(io, t)
        end
        print(io, ")")
    end
end

list() = nil()

function list(elts...)
    l = nil()
    for i=length(elts):-1:1
        l = cons(elts[i],l)
    end
    return l
end

function list(elts::T...) where T
    l = nil(T)
    for i=length(elts):-1:1
        l = cons(elts[i],l)
    end
    return l
end

function Base.length(l::Cons)
    n = 0
    for i in l
        n += 1
    end
    return n
end

function Base.map(f::Base.Callable, l::Cons{T}) where T
    isnil(l) && return l
    first = f(l.head)
    l2 = cons(first, nil(typeof(first) <: T ? T : typeof(first)))
    for h in l.tail
        l2 = cons(f(h), l2)
    end
    reverse(l2)
end

function Base.filter(f::Function, l::LinkedList{T}) where T
    l2 = nil(T)
    for h in l
        if f(h)
            l2 = cons(h, l2)
        end
    end
    reverse(l2)
end

function Base.reverse(l::LinkedList{T}) where T
    l2 = nil(T)
    for h in l
        l2 = cons(h, l2)
    end
    return l2
end

function Base.copy(l::Cons)
    isnil(l) && return l
    l2 = reverse(reverse(l))
end

Base.cat(lst::LinkedList) = lst

function Base.cat(lst::LinkedList, lsts::LinkedList...)
    T = typeof(lst).parameters[1]
    n = length(lsts)
    for i = 1:n
        T2 = typeof(lsts[i]).parameters[1]
        T = typejoin(T, T2)
    end

    l2 = nil(T)
    for h in lst
        l2 = cons(h, l2)
    end

    for i = 1:n
        for h in lsts[i]
            l2 = cons(h, l2)
        end
    end

    reverse(l2)
end

function Base.iterate(l::LinkedList, state::Cons = l)
    isnil(state) && return nothing
    state.head, state.tail
end


text:
using POMDPs
using MCTS
using POMDPModels

function POMDPs.simulate(mdp::POMDP,
                  policy::Policy,
                  initialstate::Any,
                  rng::AbstractRNG=MersenneTwister(rand(Uint32)),
                  eps::Float64=0.0)

    disc = 1.0
    r = 0.0
    s = deepcopy(initialstate)

    trans_dist = create_transition_distribution(mdp)

    while disc > eps && !isterminal(mdp, s)

        a = action(policy, s)
        r += disc*reward(mdp, s, a)

        transition!(trans_dist, mdp, s, a)
        rand!(rng, s, trans_dist)

        disc*=discount(mdp)
    end

    return r
end

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

function run_batch(n::Int64,
                   mdp::POMDP,
                   policy::Policy,
                   initialstate::Any;
                   rng=MersenneTwister(rand(Uint32)),
                   eps=0.0)
    rewards = zeros(n)
    space = states(mdp)
    s = create_state(mdp)
    rand!(s, space)
    while s.bumped && s.done && inreward(s, mdp.reward_states)
        rand!(s, space)
    end
    for i = 1:n
        rewards[i] = simulate(mdp, policy, initialstate, rng, eps)
    end
    return mean(rewards)
end

function inreward(s::Any, rstates::Any)
    for rs in rstates
        if s == rs
            return true
        end
    end
    return false
end

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

mdp = LegacyGridWorld(10,10)
rewards = zeros(20, 9)
initialstate = GridWorldState(1,1)
n=300

for (i,d) in enumerate(1:20), (j,ec) in enumerate(0.0:0.5:4.0)
  println("On: $d, $ec, $i, $j")
  mcts = MCTSSolver(depth=d, exploration_constant=ec)
  policy = MCTSPlanner(mcts, mdp)
  rewards[i,j] = run_batch(n, mdp, policy, initialstate,eps=0.5)
end


text:
module TestVariance
    using StreamStats
    using Distributions
    using Base.Test

    # Means of uniform draws
    for n in rand(1:1_000_000, 100)
        xs = rand(n)
        stat = StreamStats.Variance()
        for x in xs
            update!(stat, x)
        end
        online_v = var(stat)
        online_s = std(stat)
        online_vs = state(stat)
        online_n = nobs(stat)
        batch_v = var(xs)
        batch_s = std(xs)
        @test_approx_eq(online_v, batch_v)
        @test_approx_eq(online_vs, batch_v)
        @test_approx_eq(online_s, batch_s)
        @test online_n == n
    end
end


text:
abstract LaneFollowingDriver <: DriverModel{LaneFollowingAccel}
track_longitudinal!(model::LaneFollowingDriver, v_ego::Float64, v_oth::Float64, headway::Float64) = model # do nothing by default

function observe!(model::LaneFollowingDriver, scene::Scene1D, roadway::StraightRoadway, egoid::Int)

    vehicle_index = findfirst(scene, egoid)

    fore_res = get_neighbor_fore(scene, vehicle_index, roadway)

    v_ego = scene[vehicle_index].state.v
    v_oth = scene[fore_res.ind].state.v
    headway = fore_res.Δs

    track_longitudinal!(model, v_ego, v_oth, headway)

    return model
end


text:
using DataFrames, Gadfly

set_default_plot_size(21cm, 8cm)

D = DataFrame(x=1:6, y=[0.39, 0.26, 0.31, 0.43, 0.1, 0.78], Shape=repeat(["a","b","c"], outer=2))
coord = Coord.cartesian(xmin=0, xmax=7, ymin=0, ymax=1.0)
p1 = plot(D, x=:x, y=:y, color=:x,  coord,
    Scale.color_discrete, Geom.point, Geom.hair,
    Guide.title("Scale.color_discrete, Theme(alphas=[0.5])"),
    Theme(alphas=[0.5], discrete_highlight_color=identity,
        point_size=2mm)
)
p2 = plot(D, x=:x, y=:y, alpha=:x, shape=:Shape, coord,
    Scale.alpha_discrete, Geom.point, Geom.hair,
    Guide.title("Scale.alpha_discrete, Theme(default_color=\"green\")"),
    Theme(default_color="green",   discrete_highlight_color=c->"gray",
        point_size=2mm, alphas=[0.0,0.2,0.4,0.6,0.8,1.0])
)
hstack(p1,p2)


text:
using MultivariateSeries

f = (u,v) -> 0.5*cos(0.7*pi*(u+v))+0.4*sin(4*pi*u)-0.2*cos(pi*v)
x = @ring x1 x2
L = monoms(x,5)
T = 10
mnt = (V->f(V[1]/10,V[2]/10))
sigma = series(mnt, L)

w, Xi = decompose(sigma)


text:
function potential_matrix(p::Param)

    #Variables
    xmin        = p.xmin
    xmax        = p.xmax
    Nx          = p.xnodes

    ymin        = p.ymin
    ymax        = p.ymax
    Ny          = p.ynodes

    #Discretize space
    dx      = (xmax - xmin)/(Nx-1)
    x       = [xmin + (i-1)*dx for i in 1:Nx]

    dy      = (ymax - ymin)/(Ny-1)
    y       = [ymin + (i-1)*dy for i in 1:Ny]

    if(p.Boundaries == :periodic )
        dx      = (xmax - xmin)/(Nx)
        x       = [xmin + (i-1)*dx for i in 1:Nx]

        dy      = (ymax - ymin)/(Ny)
        y       = [ymin + (i-1)*dy for i in 1:Ny]
    end

    #Empty out argument matrix
    A = spzeros(Nx , Ny)

    for (ix , posx) in enumerate(x) , (iy , posy) in enumerate(y)

        A[ix,iy] = 0.0#1.0*((posx^2 + posy^2)/25^2)^10

    end
    #make diagonal matrice

    return sparse( Diagonal(reshape(A , Nx*Ny ) ) )

 end


text:
using ImageNoise
using Documenter, DemoCards

templates, theme = cardtheme()
examples, postprocess_cb = makedemos("examples", templates)

format = Documenter.HTML(edit_link = "master",
                         prettyurls = get(ENV, "CI", nothing) == "true",
                         assets = [theme])

makedocs(modules  = [ImageNoise],
         format   = format,
         sitename = "ImageNoise",
         pages    = [
             "Home" => "index.md",
             "Examples" => examples,
             "References" => "reference.md"
        ])

postprocess_cb()

deploydocs(repo = "github.com/johnnychen94/ImageNoise.jl.git")


text:
module DataMod

using Base.Iterators: Zip
using CSV
using DataFrames
using MLJ: macro_f1score
using Random
using Statistics
using StatsBase
using StatsModels
using ..Helpers

export Data, Reader
export read_proben, read_uci, read_toy, preader, ureader, readers_full, readers_preliminary, readers_multiclass
export prediction_accuracy, prediction_error, f1

struct Data
    df::DataFrame
    xdf::DataFrame
    xmat::Matrix{Float64}
    ydf::DataFrame
    ymat::VecOrMat{N} where N <: Number
    y1d::Vector{Int}
    xs::Vector{Symbol}
    ys::Vector{Symbol}
    nclass::Int
    function Data(df::DataFrame, xs::Vector{Symbol}, ys::Vector{Symbol})
        xdf = select(df, xs)
        ydf = select(df, ys)
        ymat = Matrix(ydf)
        new(df, xdf, Matrix(xdf), ydf, ymat, ys1d(ymat), xs, ys, length(ys))
    end
    function Data(data::Data, selxs::Vector{Symbol})
        selxdf = select(data.df, selxs)
        new(select(data.df, [selxs; data.ys]), selxdf, Matrix(selxdf), data.ydf, data.ymat, data.y1d, selxs, data.ys, length(data.ys))
    end
    function Data(data::Data, enhanced_df::DataFrame, enhanced_xs::Vector{Symbol})
        xdf = select(enhanced_df, enhanced_xs)
        new(enhanced_df, xdf, Matrix(xdf), data.ydf, data.ymat, data.y1d, enhanced_xs, data.ys, data.nclass)
    end
end
Base.show(io::IO, data::Data) = Base.show(io, data.df)
Base.iterate(d::Data) = (d, nothing)
Base.iterate(::Data, ::Any) = nothing
Base.length(d::Data) = 1

struct Reader
    name::String
    read::Function
    m::Int
    λ::Float64
    uci::Bool
    cv::Bool
    Reader(dataset::String, read::Function, uci::Bool) = new(dataset, read, 0, 0., uci, uci && !isfile("uci/$dataset" * "_fix_test.csv"))
    Reader(dataset::String, read::Function, m::Int, uci::Bool) = new(dataset, read, m, 0., uci, uci && !isfile("uci/$dataset" * "_fix_test.csv"))
    Reader(dataset::String, read::Function, m::Int, λ::Float64, uci::Bool) = new(dataset, read, m, λ, uci, uci && !isfile("uci/$dataset" * "_fix_test.csv"))
end
Base.show(io::IO, r::Reader) = begin
    cv = r.cv ? " - using " * string(getfield(r.read, 3)) * "-fold CV with RNG " * string(getfield(r.read, 2)) : ""
    repo = r.uci ? "UCI" * cv : "Proben1"
    print(io, "Reader of '$(r.name)' dataset (m = $(r.m), λ = $(r.λ)) of $repo")
end

function same_class(ŷ::Any, y::Any)::Bool
    if length(y) > 1
        argmax(ŷ) == argmax(y)
    else
        (ŷ[1] >= .5) == (y[1] >= .5)
    end
end
get_preds(models::Vector{<:RegressionModel}, data::Data)::Matrix{Float64} =
    hcat([predict(m, data.xmat) for m in models]...)
get_preds(model::RegressionModel, data::Data)::Matrix{Float64} =
    predict(model, data.xmat)
function prediction_accuracy(models::Union{<:RegressionModel, Vector{<:RegressionModel}}, data::Data)::Float64
    preds = get_preds(models, data)
    correct = sum(same_class.(eachrow(preds), eachrow(data.ymat)))::Int
    return correct / nrow(data.df)
end
function prediction_error(models::Union{<:RegressionModel, Vector{<:RegressionModel}}, data::Data)::Float64
    return 1 - prediction_accuracy(models, data)
end
function f1(models::Union{<:RegressionModel, Vector{<:RegressionModel}}, data::Data)
    predmat = get_preds(models, data)
    predy1d = ys1d(predmat)
    return macro_f1score(predy1d, data.y1d)
end

function kfold_cv_part(dat::DataFrame, rng::AbstractRNG, k::Int)::Vector{Tuple{DataFrame, DataFrame}}
    n = nrow(dat)
    idx = shuffle(rng, 1:n)
    rawsize = n / k
    nceil = k * (rawsize % 1)
    size = floor(Int, rawsize)
    sizes = [[size for _ in 1:(k - nceil)]; [size + 1 for _ in 1:nceil]]
    sum = 0
    idxs = []
    for s in sizes
        push!(idxs, idx[(sum + 1):(sum + s)])
        sum += s
    end
    return [(dat[setdiff(idx, i), :], dat[i, :]) for i in idxs]
end

function xnames(n::Integer; y::Bool = false)::Vector{Symbol}
    n == 1 && return y ? [:y] : [:x]
    return Symbol.(y ? 'y' : 'x', 1:n)
end
ynames(n::Integer)::Vector{Symbol} = xnames(n, y=true)

function read_proben(dataset::String, permutation::Int = 1, collapse_binary::Bool = true)::Tuple{Data, Data}#, Data}
    path = "proben1/$dataset/$dataset$permutation.dt"
    # 1 = bool_in; 2 = real_in; 3 = bool_out; 4 = bool_in; 5 = training_examples; 6 = validation_examples; 7 = test_examples
    heads = parse.(Int32, [x[2] for x in split.(Iterators.take(eachline(path), 7), '=')])
    nin = max(heads[1], heads[2])
    nout = max(heads[3], heads[4])
    ntrain = heads[5]
    nvalid = heads[6]
    ntest = heads[7]
    xs = xnames(nin)
    ys = ynames(nout)
    dat = CSV.read(path, DataFrame, header=[xs; ys], skipto=8, delim=' ', types=Float64, ignorerepeated=true)
    if collapse_binary && length(ys) == 2
        select!(dat, Not(ys[1]))
        deleteat!(ys, 1)
    end
    one_value = filter(x -> length(unique(dat[!, x])) == 1, xs)
    no_class = filter(y -> mean(dat[!, y]) == 0, ys)
    select!(dat, Not([one_value; no_class]))
    setdiff!(xs, one_value)
    setdiff!(ys, no_class)
    d = @inbounds dat[1:ntrain, :]
    # d_valid = @inbounds dat[ntrain+1:ntrain+nvalid, :]
    d_test = @inbounds dat[ntrain+nvalid+1:ntrain+nvalid+ntest, :]
    return (
        Data(d, xs, ys),
        # Data(d_valid, xs, ys),
        Data(d_test, xs, ys)
    )
end

function read_uci(dataset::String, rng::AbstractRNG = MersenneTwister(), k::Int = 4)::Union{Tuple{Data, Data}, Vector{Tuple{Data, Data}}}
    path = "uci/$dataset" * "_fix.csv"
    testpath = "uci/$dataset" * "_fix_test.csv"
    heads = parse.(Int32, [x[2] for x in split.(Iterators.take(eachline(path), 2), '=')])
    nin = heads[1]
    nout = heads[2]
    xs = xnames(nin)
    ys = ynames(nout)
    dat = CSV.read(path, DataFrame, header=[xs; ys], skipto=3, delim=',', types=Float64, ignorerepeated=true)
    if isfile(testpath)
        testdat = CSV.read(testpath, DataFrame, header=[xs; ys], skipto=3, delim=',', types=Float64, ignorerepeated=true)
        return (Data(dat, xs, ys), Data(testdat, xs, ys))
    else
        cv_dat = kfold_cv_part(dat, rng, k)
        return [(Data(dat, xs, ys), Data(testdat, xs, ys)) for (dat, testdat) in cv_dat]
    end
end

read_toy()::Data = read_proben("_toy")[1]

preader(dataset::String, m::Int = 0, λ::Float64 = 0.)::Reader = Reader(dataset, () -> read_proben(dataset), m, λ, false)
ureader(dataset::String, rng::AbstractRNG = MersenneTwister(), m::Int = 0, λ::Float64 = 0., k::Int = 4)::Reader = Reader(dataset, () -> read_uci(dataset, rng, k), m, λ, true)

function readers_full(start::Int = 1)::Vector{Reader}
    dp = ["cancer", "card", "gene", "glass", "heart", "horse", "soybean", "thyroid"]
    mp = [5, 5, 5, 5, 5, 5, 5, 145]
    lp = [.1, .5, 7.5, .1, .1, 1., .5, .1]
    du = ["abalone", "adult", "covtype", "madelon", "optdigits", "page-blocks", "pendigits",
          "poker", "satimage", "segmentation", "shuttle", "spect", "vehicle", "waveform"]
    mu = [110, 150, 180, 5, 175, 100, 185, 185, 180, 25, 70, 150, 115, 60]
    lu = [.1, .1, .1, 7., .5, .1, .1, .1, .1, .1, .1, .5, .1, 1.5]
    return [[preader(d, m, l) for (d, m, l) in zip(dp, mp, lp)]; [ureader(du[d], MersenneTwister(123d), mu[d], lu[d]) for d in eachindex(du)]][start:end]
end
function readers_preliminary(start::Int = 1)::Vector{Reader}
    dp = ["cancer", "card", "gene", "glass", "heart", "horse", "mushroom", "soybean", "thyroid"][start:end]
    return [preader(d) for d in dp]
end
function readers_multiclass()::Vector{Reader}
    dp = ["gene", "glass", "horse", "soybean", "thyroid"]
    du = ["abalone", "covtype", "optdigits", "page-blocks", "pendigits",
          "poker", "satimage", "segmentation", "shuttle", "vehicle", "waveform"]
    return [[preader(d) for d in dp]; [ureader(du[d], MersenneTwister(123d)) for d in eachindex(du)]]
end

end

using .DataMod

max_stars_repo_name: KasperNooteboom/thesis-rvfl-fs

max_issues_repo_name: KasperNooteboom/thesis-rvfl-fs

max_forks_repo_name: KasperNooteboom/thesis-rvfl-fs


text:
using .LazySets: LazySet, support_vector

function ConvexBodyProximityQueries.support(obj::LazySet{N}, dir::SVector{D, N}) where {D, N}
    return support_vector(dir, obj) |> SVector{D, N}
end


text:
using Roassal
using Test
using Gtk
# ------------------------------------
function testColor()
    aColor = RColor()
    @test aColor.r == aColor.g
    @test aColor.r == aColor.b
end

function testBox01()
    aColor = RColor()
    box1 = RBox(color=aColor)
    @test box1.color == aColor
    @test pos(box1) == (0,0)
    @test extent(box1) == (10,10)
    @test compute_encompassing_rectangle(box1) == (-5, -5, 10, 10)

    c = get_color(box1)
    @test c === aColor
    anotherColor = RColor(1.0, 1.0, 0.5)
    set_color!(box1, anotherColor)
    @test get_color(box1) === anotherColor
end

function testBox02()
    box = RBox()
    @test box.x == 0
    @test box.color != RColor()
end

function testCanvas()
    box1 = RBox(color=RColor(0.4, 0.1, 0.3))
    box = RBox()
    canvas = RCanvas()
    @test number_of_shapes(canvas) == 0
    add!(canvas, box1)
    @test number_of_shapes(canvas) == 1
    add!(canvas, box)
    @test number_of_shapes(canvas) == 2
end

function testCanvas02()
    canvas = RCanvas()
    @test isempty(get_shapes(canvas))
    box1 = RBox
    box2 = RBox
    add!(canvas, box1)
    add!(canvas, box2)
    @test length(shapeOf(canvas)) == 2
    @test get_shapes(canvas)[1] === box1
    @test get_shapes(canvas)[2] === box2
end

function testRemovingShape()
    canvas = RCanvas()
    box1 = RBox()
    box2 = RBox()
    box3 = RBox()
    add!(canvas, box1)
    add!(canvas, box2)
    add!(canvas, box3)

    removeShape!(canvas, box2)
    @test number_of_shapes(canvas) == 2
    @test get_shapes(canvas)[1] === box1
    @test get_shapes(canvas)[2] === box3

    removeShape!(canvas, box1)
    @test number_of_shapes(canvas) == 1
    @test get_shapes(canvas)[1] === box3

    # Not in the canvas!
    removeShape!(canvas, box1)
    @test number_of_shapes(canvas) == 1
    @test get_shapes(canvas)[1] === box3

    removeShape!(canvas, box3)
end

# ------------------------------------
function testFetchingShape()
    canvas = RCanvas()
    box = RBox()
    add!(canvas, box)
    translate_to!(box, 50, 50)
    @test pos(box) == (50, 50)

    #@test isPositionWithinEncompassingRectangleOf(box, 50, 50)
    #@test !isPositionWithinEncompassingRectangleOf(box, -50, 50)
    #@test !isPositionWithinEncompassingRectangleOf(box, -50, -50)

    @test get_shape_at_position(canvas, 50, 50) === box
    @test get_shape_at_position(canvas, 50, 50) == box
    @test get_shape_at_position(canvas, -50, -50) == canvas
    @test get_shape_at_position(canvas, -50, -50) === canvas
end

# ------------------------------------
function testCallbacks01()
    box = RBox()
    @test numberOfCallbacks(box) == 0
    add_callback!(box, Callback(:mouseMove, (event, shape) -> println("enter box")))
    @test numberOfCallbacks(box) == 1
end

function testCallbacks02()
    box = RBox()
    a = 0

    # Does not do much since there is no callback
    trigger_callback(box, :mouseMove, "not an event")
    trigger_callback(box, :mouseMove, "not an event")

    add_callback!(box, Callback(:mouseMove, () -> a = a + 1))
    @test numberOfCallbacks(box) == 1
    @test a == 0
    trigger_callback(box, :mouseMove, "not an event")
    @test a == 1
    trigger_callback(box, :mouseMove, "not an event")
    @test a == 2

    trigger_callback(box, :mouseMove, "not an event2")
    @test a == 3

    trigger_callback(box, :mouseMoveNot, "not an event")
    @test a == 3
end

function testCallbacks03()
    canvas = RCanvas()
    @test numberOfCallbacks(canvas) == 0
    a = 0

    # Does not do much since there is no callback
    trigger_callback(canvas, :mouseMove, "not an event")
    trigger_callback(canvas, :mouseMove, "not an event")

    add_callback!(canvas, Callback(:mouseMove, () -> a = a + 1))
    @test numberOfCallbacks(canvas) == 1
    @test a == 0
    trigger_callback(canvas, :mouseMove, "not an event")
    @test a == 1
    trigger_callback(canvas, :mouseMove, "not an event")
    @test a == 2

    trigger_callback(canvas, :mouseMove, "not an event2")
    @test a == 3

    trigger_callback(canvas, :mouseMoveNot, "not an event")
    @test a == 3
end

#= function testCallbacks04()
    canvas = RCanvas()
    t = nothing
    add_callback!(canvas, Callback(:mouseMove, () -> t = (event, shape)))
    @test isnothing(t)
    @test numberOfCallbacks(canvas) == 1

    trigger_callback(canvas, :mouseMove, "not an event")
    @test length(t) == 2
    @test t[1] == "not an event"
    @test t[2] == canvas
end =#

function testCallback05()
    canvas = RCanvas()
    box1 = RBox()
    box2 = RBox()
    add!(canvas, box1)
    add!(canvas, box2)
    add_callback!(get_shapes(canvas), Callback(:mouseMove, () -> 42))
    @test numberOfCallbacks(box1) == 1
    @test numberOfCallbacks(box2) == 1

    add_callback!(get_shapes(canvas)[1], Callback(:mouseMove, () -> 42))
    @test numberOfCallbacks(box1) == 2
    @test numberOfCallbacks(box2) == 1

    add_callback!(get_shapes(canvas), Callback(:mouseMove, () -> 42))
    @test numberOfCallbacks(box1) == 3
    @test numberOfCallbacks(box2) == 2
end

@testset "test - offset_from_canvas_to_screen" begin
    gtk = GtkCanvas()
    @test offset_from_canvas_to_screen(gtk) != (0, 0)
    @test offset_from_canvas_to_screen(gtk) != (0, 0)
    # @test (offset_from_canvas_to_screen(gtk) .+ offset_from_canvas_to_screen(gtk)) == (0, 0)
end

# ------------------------------------
@testset "Roassal" begin
    testColor()
    testBox01()
    testBox02()
    testCanvas()
    testFetchingShape()
    testCallbacks01()
    testCallbacks02()
    #testCallbacks04()
end

include("Interaction.jl");
include("Line.jl");


text:
module PolarimetryCalibration

include("mueller.jl")

export müller_rec, müller_sky, lin_pol_source

end # module

max_stars_repo_name: KajWiik/PolarimetryCalibration

max_issues_repo_name: KajWiik/PolarimetryCalibration

max_forks_repo_name: KajWiik/PolarimetryCalibration


text:
# %% markdown
# # ImageList -> Main
# %% codecell
#export
using FileIO
using Images
using Serialization
using Random
using CUDAapi
using Plots
import GR
using Images
using CuArrays
using ImageView
using Statistics, Random
import ProgressMeter
using Distributions
using Zygote
gr()
# %% markdown
# fromFolder
# %% codecell
#export
# Helper function to add path to array
function add_path(cat::String)
    temp_dir = readdir(joinpath(path, cat))
    return [joinpath(path, cat, x) for x in temp_dir],
    fill(cat, size(temp_dir, 1))
end
# %% markdown
# ## Class Distribution
# %% codecell
#export
function classDistribution(y)
    """
    Function to plot class distribution to see if balanced or not.
    """
    labels = unique(y)
    cnts = [sum(y .== i) for i in labels]
    display(plot(cnts, seriestype = [:bar]))
    return cnts, maximum(cnts)
end
# %% codecell
#export
path = "/media/subhaditya/DATA/COSMO/Datasets/catDog/"
# %% codecell
#export
# Define number of threads
Threads.nthreads() = length(Sys.cpu_info())
# %% codecell
collect(1:10)
# %% codecell
#export

"""
Function to create an array of images and labels -> when the directory structure is as follows
- main
    - category1
        - file1...
    -category2
        - file1...
    ...
"""
function fromFolder(path::String, imageSize = 64::Int64)
    @info path, imageSize
    categories = readdir(path)
    total_files =
        collect(Iterators.flatten([add_path(x)[1] for x in categories]))
    total_categories =
        collect(Iterators.flatten([add_path(x)[2] for x in categories]))
    distrib, max_dis = classDistribution(total_categories)

    indices_repeat = indexin(unique(total_categories), total_categories)
    # oversample
    total_add = max_dis .- distrib # get the differences to oversample
    oversample = false
    if sum(total_add) > 100
        @info "Oversampling"
        images = zeros((
            imageSize,
            imageSize,
            3,
            size(max_dis * length(unique(total_categories)), 1),
        ))
        oversample = true
        oversample_index = length(y) - sum(total_add)# keep a track of indices from the back
    else
        @info "No need to oversample"
        images = zeros((imageSize, imageSize, 3, size(total_categories, 1)))
        oversample = false
    end

    Threads.@threads for idx in collect(1:size(total_files, 1))
        img = channelview(imresize(
            load(total_files[idx]),
            (imageSize, imageSize),
        ))
        img = convert(Array{Float64}, img)
        images[:, :, :, idx] = permutedims(img, (2, 3, 1))
        #         @info oversample
        if oversample == true

            if idx in indices_repeat
                labelrep = findfirst(x -> x == idx, indices_repeat) # index in the repeated list
                to_repeat = total_add[labelrep] # no of times to repeat
                total_categories = vcat(
                    total_categories,
                    fill(total_categories[indices_repeat[labelrep]], to_repeat),
                )
                Threads.@threads for idx2 in collect(oversample_index:to_repeat)
                    images[:, :, :, idx2] =
                        images[:, :, :, indices_repeat[labelrep]]

                end

            end
        end
    end

    @info "Done loading images"

    return images, total_categories

end

# %% codecell
#export
X, y = fromFolder(path, 64);
# %% codecell
size(X), size(y)
# %% markdown
# ## Splitting
# %% codecell
at = 0.7
n = length(y)
idx = shuffle(1:n)
train_idx = view(idx, 1:floor(Int, at * n));
test_idx = view(idx, (floor(Int, at * n)+1):n);
# %% codecell
ytrain, ytest = y[train_idx, :], y[test_idx, :]
Xtrain, Xtest = X[:, :, :, train_idx], X[:, :, :, test_idx]
@info length(ytrain), length(ytest)
@info length(Xtrain), length(Xtest)
# %% codecell
#export
function splitter(pct_split = 0.7::Float16)
    """
    Splits into train/test by pct_split%
    """
    n = length(y)
    idx = shuffle(1:n)
    train_idx = view(idx, 1:floor(Int, pct_split * n))
    test_idx = view(idx, (floor(Int, pct_split * n)+1):n)
    ytrain, ytest = y[train_idx, :], y[test_idx, :]
    Xtrain, Xtest = X[:, :, :, train_idx], X[:, :, :, test_idx]
    return Xtrain, ytrain, Xtest, ytest
end
# %% codecell
#export
Xtrain, ytrain, Xtest, ytest = splitter(0.8);
# %% markdown
# ## Linear
# %% codecell
using Zygote
# %% codecell
W = rand(2, 5)
b = rand(2)
# %% codecell
Dense(x) = W * x .+ b
# %% codecell
function loss(x, y)
    ŷ = Dense(x)
    sum((y - ŷ) .^ 2)
end
# %% codecell
x, y = rand(5), rand(2)
# loss(x,y)
# %% codecell
α = 0.1
# %% codecell
x, y = rand(100), rand(2)
W = rand(2, 100)
b = rand(2)
# %% codecell
W = zeros(2, 100)
# %% codecell
for a in collect(1:50)
    gs = gradient(() -> loss(x, y), Params([W, b]))
    W̄ = gs[W]
    W .= α .* W̄
    ŷ = W * x .+ b
    @info sum((y - ŷ) .^ 2)
end
# %% markdown
# # Initialization
# - Zero Initialization: set all weights to 0
# - Normal Initialization: set all weights to random small numbers
# - Lecun Initialization: normalize variance
# - Xavier Intialization (glorot init)
# - Kaiming Initialization (he init)
# %% markdown
# ## Lecun
# It draws samples from a truncated normal distribution centered on 0 with stddev <- sqrt(1 / fan_in) where fan_in is the number of input units in the weight tensor..
# %% codecell
using Distributions
# %% codecell
#export
lecun_normal(fan_in) = return Distributions.Normal(0, sqrt(1 / fan_in))
# %% codecell
W = rand(lecun_normal(2), 2, 100)
b = rand(lecun_normal(2), 2)
# %% markdown
# ## Xavier Normal
# It draws samples from a truncated normal distribution centered on 0 with stddev = sqrt(2 / (fan_in + fan_out)) where fan_in is the number of input units in the weight tensor and fan_out is the number of output units in the weight tensor.
# %% codecell
#export
xavier_normal(fan_in, fan_out) =
    return Distributions.Normal(0, sqrt(2 / (fan_in + fan_out)))
# %% codecell
W = rand(xavier_normal(2, 100), 2, 100)
b = rand(xavier_normal(2, 2), 2)
# %% markdown
# # Xavier Uniform
# It draws samples from a uniform distribution within -limit, limit where limit is sqrt(6 / (fan_in + fan_out)) where fan_in is the number of input units in the weight tensor and fan_out is the number of output units in the weight tensor.
# %% codecell
#export
function xavier_uniform(fan_in, fan_out)
    limit = sqrt(6 / (fan_in + fan_out))
    return Distributions.Uniform(-limit, limit)
end
# %% codecell
W = rand(xavier_uniform(2, 100), 2, 100)
b = rand(xavier_uniform(2, 2), 2)
# %% markdown
# ## He Normal
# It draws samples from a truncated normal distribution centered on 0 with stddev = sqrt(2 / fan_in) where fan_in is the number of input units in the weight tensor.
# %% codecell
#export
he_normal(fan_in) = return Distributions.Normal(0, sqrt(2 / (fan_in)))
# %% codecell
W = rand(he_normal(2), 2, 100)
b = rand(he_normal(2), 2)
# %% markdown
# # He Uniform
# It draws samples from a uniform distribution within -limit, limit where limit is sqrt(6 / fan_in) where fan_in is the number of input units in the weight tensor.
# %% codecell
#export
function he_uniform(fan_in)
    limit = sqrt(6 / (fan_in))
    return Distributions.Uniform(-limit, limit)
end
# %% codecell
W = rand(he_uniform(2), 2, 100)
b = rand(he_uniform(2), 2)
# %% markdown
# # Batching
# %% markdown
# ## One hot
# %% codecell
ytrain
# %% codecell
labels = unique(ytrain);
encodedlabels = Dict(labels .=> collect(1:length(labels)))
# %% markdown
## One cold
# %% codecell
ytrain
# %% codecell
function onecold(y_enc)
    labels = unique(y_enc)
    encodedlabels = Dict(labels .=> collect(1:length(labels)))
    global ytrain
    for a in keys(encodedlabels)
        ytrain = replace(ytrain, a=>encodedlabels[a])
    end
end
# %% codecell
onecold(ytrain)

# %% codecell
ytrain
# %% codecell
# export
function testGen(c::Channel)
    for n=1:4
        put!(c,n);
    end
    put!(c,"stop");
end
# %% codecell
test_yi = Channel(testGen);

# %%
take!(test_yi)

## Batch generator
# %%
#export
function datagen(c::Channel)
    global rep_len
    for n=1:bs:rep_len
        put!(c,n);
    end
    put!(c,"stop");
end

#%%
bs = 64
rep_len = length(ytest)
bunchData = Channel(datagen);
#%% markdown
## Actual batches
#%%
for i in collect(1:round(rep_len/bs)+1)
    current_index = take!(bunchData)
    try
        x_batch,y_batch = Xtest[:, :, :, current_index:current_index+bs-1],ytest[current_index:current_index+bs-1]
        @info size(x_batch),size(y_batch)
    catch e
        x_batch,y_batch = Xtest[:, :, :, rep_len-bs:rep_len],ytest[rep_len-bs:rep_len]
        @info size(x_batch),size(y_batch)
    end
end

#%%
## Training loop with SGD
bs = 64
rep_len = length(ytest)
bunchData = Channel(datagen);
W = rand(he_uniform(rep_len), rep_len, 100)
b = rand(he_uniform(rep_len), rep_len)

for i in collect(1:round(rep_len/bs)+1)
    current_index = take!(bunchData)
    try
        x_batch,y_batch = Xtest[:, :, :, current_index:current_index+bs-1],ytest[current_index:current_index+bs-1]
        @info size(x_batch),size(y_batch)
    catch e
        x_batch,y_batch = Xtest[:, :, :, rep_len-bs:rep_len],ytest[rep_len-bs:rep_len]
        @info size(x_batch),size(y_batch)
    end
end