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

using Markdown
using InteractiveUtils

# ╔═╡ c41fb568-25bd-11ec-1f56-d1fbf2dd4e33
# Load up modules for working
begin
	using Pkg
	cd(dirname(@__FILE__))
	Pkg.activate(".")
	@info pwd()

	function ingredients(path::String)
		name = Symbol(basename(path))
		m = Module(name)
		Core.eval(m,
			Expr(:toplevel,
				 :(eval(x) = $(Expr(:core, :eval))($name, x)),
				 :(include(x) = $(Expr(:top, :include))($name, x)),
				 :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
				 :(include($path))))
		m
	end

	Analyse = ingredients("./src/Analyse.jl").Analyse
end

# ╔═╡ 770505bc-19cb-44c4-b8de-30e37234bc68
md"""
Notebook showing the entire pipeline for Twitter information operations report datasets

"""

# ╔═╡ d2c74016-f98a-4902-8c61-ba0a33195d47
md"""
we enlarge the default notebook width
"""

# ╔═╡ 7a60c92e-9d74-492b-970a-b271896c85d4
# Make cells wider
html"""<style>
/*              screen size more than:                     and  less than:                     */
@media screen and (max-width: 699px) { /* Tablet */
  /* Nest everything into here */
    main { /* Same as before */
        max-width: 1200px !important; /* Same as before */
        margin-right: 200px !important; /* Same as before */
    } /* Same as before*/

}

@media screen and (min-width: 700px) and (max-width: 1199px) { /* Laptop*/
  /* Nest everything into here */
    main { /* Same as before */
        max-width: 1200px !important; /* Same as before */
        margin-right: 200px !important; /* Same as before */
    } /* Same as before*/
}

@media screen and (min-width:1200px) and (max-width: 1920px) { /* Desktop */
  /* Nest everything into here */
    main { /* Same as before */
        max-width: 1200px !important; /* Same as before */
        margin-right: 200px !important; /* Same as before */
    } /* Same as before*/
}

@media screen and (min-width:1921px) { /* Stadium */
  /* Nest everything into here */
    main { /* Same as before */
        max-width: 1200px !important; /* Same as before */
        margin-right: 200px !important; /* Same as before */
    } /* Same as before*/
}
</style>
"""

# ╔═╡ 67adb4ee-a0dc-47ca-bd0d-67d3b2687ee4
md"""
We load up the analysis modules
"""

# ╔═╡ eebae507-271f-4935-8f13-ca25bd1907f6
datapath = "./bigdata"

# ╔═╡ b0525f84-427f-400c-8f42-189b84283d42
readdir("./bigdata/")

# ╔═╡ 07fa28e4-ea9f-4bd9-a0ad-38f3f8d94a47
md"""
## Get external tweet ids
the function below will identify the external tweets to recover and write them out to as csv file.
"""

# ╔═╡ f7549366-d538-4900-b4e2-85d28169925c
Analyse.get_externals(datapath)

# ╔═╡ 3036b5bd-77c5-441a-bc38-bf23bb56b12c
md"""## Download external tweet ids
cf. `./src/hydrate.py` & `./src/hydratorpipeline.py`"""

# ╔═╡ 3be71733-6e2f-4ec3-83ec-e46bcc872371
md"""
## Generate interaction networks

"""

# ╔═╡ 0ee608f7-5c64-4051-bc19-7f64e04305a5
begin
	msgfiles = "./bigdata/honduras_022020_tweets_csv_hashed.csv"
	usrfiles = "./bigdata/honduras_022020_users_csv_hashed.csv"
	G_rt, G_rp, G_rt_bp, G_rp_bp, G_rt_bp_proj, G_rp_bp_proj, Uinfo, RUinfo = Analyse.grapher(msgfiles, usrfiles)
end

# ╔═╡ 400c36ff-d96a-482a-8d92-48e5d6407bd1

# ╔═╡ Cell order:
# ╟─770505bc-19cb-44c4-b8de-30e37234bc68
# ╟─d2c74016-f98a-4902-8c61-ba0a33195d47
# ╟─7a60c92e-9d74-492b-970a-b271896c85d4
# ╟─67adb4ee-a0dc-47ca-bd0d-67d3b2687ee4
# ╠═c41fb568-25bd-11ec-1f56-d1fbf2dd4e33
# ╠═eebae507-271f-4935-8f13-ca25bd1907f6
# ╠═b0525f84-427f-400c-8f42-189b84283d42
# ╟─07fa28e4-ea9f-4bd9-a0ad-38f3f8d94a47
# ╠═f7549366-d538-4900-b4e2-85d28169925c
# ╟─3036b5bd-77c5-441a-bc38-bf23bb56b12c
# ╟─3be71733-6e2f-4ec3-83ec-e46bcc872371
# ╠═0ee608f7-5c64-4051-bc19-7f64e04305a5
# ╠═400c36ff-d96a-482a-8d92-48e5d6407bd1


text:
@testset "DataFrames" begin

    using DataFrames

    # DataFrame from Vector{Bedgraph.Record}.
    df = DataFrame(Bag.records)

    @test typeof(df) == DataFrame
    @test size(df) == (9,4)

    @test df[:chrom] == Bag.chroms
    @test df[:first] == Bag.firsts
    @test df[:last] == Bag.lasts
    @test df[:value] == Bag.values

    @test DataFrame(Bag.records) == Bag.records |> DataFrame

    # DataFrame from bedGraph file.
    df2 = DataFrame(load(Bag.file))

    @test typeof(df2) == DataFrame
    @test size(df2) == (9,4)

    @test df2[:chrom] == Bag.chroms
    @test df2[:first] == Bag.firsts
    @test df2[:last] == Bag.lasts
    @test df2[:value] == Bag.values

    @test DataFrame(load(Bag.file)) == load(Bag.file) |> DataFrame

    # DataFrame from headerless bedGraph file.
    df3 = DataFrame(load(Bag.file_headerless))
    @test typeof(df3) == DataFrame
    @test size(df3) == (9,4)

    @test df3[:chrom] == Bag.chroms
    @test df3[:first] == Bag.firsts
    @test df3[:last] == Bag.lasts
    @test df3[:value] == Bag.values

    @test DataFrame(load(Bag.file_headerless)) == load(Bag.file_headerless) |> DataFrame

    # Save and load from DataFrame.
    save(Bag.tmp_output_path, df)
    @test df == load(Bag.tmp_output_path) |> DataFrame

    df |> save(Bag.tmp_output_path)
    @test df == load(Bag.tmp_output_path) |> DataFrame

end # test DataFrames

ext: jl


text:
function getindexfromyear_rice_2010(year)
    baseyear = 2005

    if rem(year - baseyear, 10) != 0
        error("Invalid year")
    end

    return div(year - baseyear, 10) + 1
end

#Function to read a single parameter value from original RICE 2010 model.
function getparam_single(f, range::AbstractString, regions)
    vals= Array{Float64}(undef, length(regions))
    for (i,r) = enumerate(regions)
        data=readxl(f,"$(r)!$(range)")
        vals[i]=data[1]
    end
    return vals
end

#Function to read a time series of parameter values from original RICE 2010 model.
function getparam_timeseries(f, range::AbstractString, regions, T)
    vals= Array{Float64}(undef, T, length(regions))
    for (i,r) = enumerate(regions)
        data=readxl(f,"$(r)!$(range)")
        for n=1:T
            vals[n,i] = data[n]
        end
    end
    return vals
end


text:
const freudenstein_roth = let res_init=zeros(2), jac_init=zeros(2,2), x_init=[15.0, -2.0]

    function res(x, r)
        r[1] = -1.3e1 + x[1] - 2.0e0*x[2] + 5.0e0*x[2]^2 - x[2]^3
        r[2] = -2.9e1 + x[1] - 1.4e1*x[2] + x[2]^2 + x[2]^3
        return r
    end

    function jac(x, j)
        j[1,1] = 1.e0
        j[1,2] = -2.e0 + x[2]*(1.e1 - 3.e0*x[2])
        j[2,1] = 1.e0
        j[2,2] = -1.4e1 + x[2]*(2.e0 + 3.e0*x[2])
        return j
    end

    f(;scale=1, verbose=false, print_steps=false) = testone("freudenstein_roth", res, jac,
                                                   res_init, jac_init, x_init;
                                                   scale=scale, verbose=verbose,
                                                   print_steps=print_steps)
end


text:
# MIRT.jl

# note: using Reexport might be useful later

module MIRT

include("z-all.jl")

end # module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""ClusterRole is a cluster level, logical grouping of PolicyRules that can be referenced as a unit by a RoleBinding or ClusterRoleBinding. Deprecated in v1.17 in favor of rbac.authorization.k8s.io/v1 ClusterRole, and will no longer be served in v1.20.

    IoK8sApiRbacV1alpha1ClusterRole(;
        aggregationRule=nothing,
        apiVersion=nothing,
        kind=nothing,
        metadata=nothing,
        rules=nothing,
    )

    - aggregationRule::IoK8sApiRbacV1alpha1AggregationRule : AggregationRule is an optional field that describes how to build the Rules for this ClusterRole. If AggregationRule is set, then the Rules are controller managed and direct changes to Rules will be stomped by the controller.
    - apiVersion::String : APIVersion defines the versioned schema of this representation of an object. Servers should convert recognized schemas to the latest internal value, and may reject unrecognized values. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#resources
    - kind::String : Kind is a string value representing the REST resource this object represents. Servers may infer this from the endpoint the client submits requests to. Cannot be updated. In CamelCase. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - metadata::IoK8sApimachineryPkgApisMetaV1ObjectMeta : Standard object's metadata.
    - rules::Vector{IoK8sApiRbacV1alpha1PolicyRule} : Rules holds all the PolicyRules for this ClusterRole
"""
mutable struct IoK8sApiRbacV1alpha1ClusterRole <: SwaggerModel
    aggregationRule::Any # spec type: Union{ Nothing, IoK8sApiRbacV1alpha1AggregationRule } # spec name: aggregationRule
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    rules::Any # spec type: Union{ Nothing, Vector{IoK8sApiRbacV1alpha1PolicyRule} } # spec name: rules

    function IoK8sApiRbacV1alpha1ClusterRole(;aggregationRule=nothing, apiVersion=nothing, kind=nothing, metadata=nothing, rules=nothing)
        o = new()
        validate_property(IoK8sApiRbacV1alpha1ClusterRole, Symbol("aggregationRule"), aggregationRule)
        setfield!(o, Symbol("aggregationRule"), aggregationRule)
        validate_property(IoK8sApiRbacV1alpha1ClusterRole, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiRbacV1alpha1ClusterRole, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiRbacV1alpha1ClusterRole, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiRbacV1alpha1ClusterRole, Symbol("rules"), rules)
        setfield!(o, Symbol("rules"), rules)
        o
    end
end # type IoK8sApiRbacV1alpha1ClusterRole

const _property_map_IoK8sApiRbacV1alpha1ClusterRole = Dict{Symbol,Symbol}(Symbol("aggregationRule")=>Symbol("aggregationRule"), Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("rules")=>Symbol("rules"))
const _property_types_IoK8sApiRbacV1alpha1ClusterRole = Dict{Symbol,String}(Symbol("aggregationRule")=>"IoK8sApiRbacV1alpha1AggregationRule", Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("rules")=>"Vector{IoK8sApiRbacV1alpha1PolicyRule}")
Base.propertynames(::Type{ IoK8sApiRbacV1alpha1ClusterRole }) = collect(keys(_property_map_IoK8sApiRbacV1alpha1ClusterRole))
Swagger.property_type(::Type{ IoK8sApiRbacV1alpha1ClusterRole }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiRbacV1alpha1ClusterRole[name]))}
Swagger.field_name(::Type{ IoK8sApiRbacV1alpha1ClusterRole }, property_name::Symbol) =  _property_map_IoK8sApiRbacV1alpha1ClusterRole[property_name]

function check_required(o::IoK8sApiRbacV1alpha1ClusterRole)
    true
end

function validate_property(::Type{ IoK8sApiRbacV1alpha1ClusterRole }, name::Symbol, val)
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.6.3
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# * https://twitter.com/asunokibou/status/1457694534016507922
# * https://twitter.com/mathsorcerer/status/1457937816235634691
# * https://twitter.com/sophia8wendel/status/1457980892966576128

# %%
ENV["LINES"] = 100
using Combinatorics
sol = [p for p in permutations(1:9, 7) if sum(i//p[i] for i in 1:7) == 7]

# %%

max_stars_repo_path: 0024/asunokibou puzzle.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0024/asunokibou puzzle.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0024/asunokibou puzzle.jl

max_forks_repo_name: genkuroki/public


text:
@test SolverLogging.getwidth("%0.24e") == -1
@test SolverLogging.getwidth("%01.24e") == 1
@test SolverLogging.getwidth("%0112.24e") == 112
@test SolverLogging.getwidth("%-1.24e") == 1
@test SolverLogging.getwidth("012%-13.24e") == 13
@test SolverLogging.getwidth("%.24e") == -1
@test SolverLogging.getwidth("%e") == -1
@test SolverLogging.getwidth("%5d") == 5
@test SolverLogging.getwidth("%512d") == 512
@test SolverLogging.getwidth("%0512d") == 512


text:
function latexOutput(outfile::String = "")
    latexRegressandTransform(s::String,colmin::Int64,colmax::Int64) = "\\multicolumn{$(colmax-colmin+1)}{c}{$s}"
    latexTableHeader(numberOfResults::Int64, align::String) = "\\begin{tabular}{$align}"
    latexTableFooter(numberOfResults::Int64, align::String) = "\\end{tabular}"
    function latexHeaderRule(headerCellStartEnd::Vector{Vector{Int64}})
        if length(headerCellStartEnd)<2
            error("Invalid headerCellStartEnd: need to have at least two columns.")
        end
        s = ""
        for i in headerCellStartEnd[2:end]
            s = s * "\\cmidrule(lr){$(i[1])-$(i[2])}" * " "
        end
        return s
    end
    toprule = "\\toprule"
    midrule = "\\midrule"
    bottomrule = "\\bottomrule"
    headerrule = latexHeaderRule
    headercolsep = " & "
    colsep = " & "
    linestart = ""
    linebreak = " \\\\ "

    label_fe_yes = "Yes"
    label_fe_no = ""

    label_statistic_n = "\$N\$"
    label_statistic_r2 = "\$R^2\$"
    label_statistic_f = "\$F\$"
    label_statistic_adjr2 = "Adjusted \$R^2\$"
    label_statistic_r2_within = "Within-\$R^2\$"
    label_statistic_p = "\$F\$-test \$p\$ value"
    label_statistic_f_kp = "First-stage \$F\$ statistic"
    label_statistic_p_kp = "First-stage \$p\$ value"
    label_statistic_dof = "Degrees of Freedom"

    label_estimator = "Estimator"
    label_estimator_ols = "OLS"
    label_estimator_iv = "IV"
    label_estimator_nl = "NL"

    foutfile = outfile
    encapsulateRegressand = latexRegressandTransform
    header = latexTableHeader
    footer = latexTableFooter
    return RenderSettings(toprule, midrule, bottomrule, headerrule, headercolsep, colsep, linestart,
        linebreak, label_fe_yes, label_fe_no,
        label_statistic_n, label_statistic_r2, label_statistic_adjr2, label_statistic_r2_within,
        label_statistic_f, label_statistic_p, label_statistic_f_kp, label_statistic_p_kp, label_statistic_dof,
        label_estimator, label_estimator_ols, label_estimator_iv, label_estimator_nl,
        foutfile, encapsulateRegressand, header, footer)
end


text: fun(x,y) = x + y

max_stars_repo_name: waseemssaeed/TetsRepo

max_issues_repo_name: waseemssaeed/TetsRepo

max_forks_repo_name: waseemssaeed/TetsRepo


text:
struct SFloat{T} <: AbstractFloat
    value :: T
end
value(x::SFloat) = x.value
SFloat{T}(x::SFloat) where T = SFloat{T}(T(x.value))

"""
Stochastic Float64

This type represents a double-precision (64-bit) floating-point number, on which
every operation gets randomly rounded upwards or downwards.
"""
const SFloat64 = SFloat{Float64}

"""
Stochastic Float32

This type represents a single-precision (32-bit) floating-point number, on which
every operation gets randomly rounded upwards or downwards.
"""
const SFloat32 = SFloat{Float32}

det_type(::Type{T}) where {T<:Real} = T
det_type(::Type{SFloat{T}}) where {T} = T

sto_type(::Type{T}) where {T} = SFloat{T}

Base.promote_rule(::Type{T1}, ::Type{T2}) where {T1<:SFloat, T2} = sto_type(promote_type(det_type(T1), det_type(T2)))
result_type(x::T1, y::T2) where {T1, T2} = promote_type(T1, T2)

#binary operators
for op = (:+, :-, :*, :/)
    eval(quote
        Base.$op(a::SFloat,b::SFloat) = result_type(a,b)($op(RND, a.value,  b.value))
    end)
end

for op = (:<, :<=)
    eval(quote
        Base.$op(a::SFloat,b::SFloat) = $op(a.value,  b.value)
    end)
end

#unary operators
for op = (:-, :zero, :one, :abs)
    eval(quote
        Base.$op(a::SFloat) = SFloat($op(a.value))
    end)
end

value(x::Number) = x
deref(x) = x
deref(x::Array{T,0}) where T = x[]
macro reliable_digits(expr)
    return quote
        x = [value.($(esc(expr))) for i = 1:10]
        mu = mean(x)
        sigma = std(x)
        s = -log10.((/).(sigma, abs.(mu)))
        zip(mu, s) |> collect |> deref
    end
end


text:
"""
    DynamicsUtils

Utilities for dynamics simulations.
Includes:
* Basic dynamics variables functions
* Density matrix dynamics functions
* Standard callbacks to use during dynamics
* Plotting recipes for outputs
"""
module DynamicsUtils

using NQCDynamics:
    AbstractSimulation,
    Simulation,
    RingPolymerSimulation,
    Calculators,
    masses

"""
    divide_by_mass!(dv, masses)

Divide the contents of `dv` by the `masses`.
Assumes `dv` is an array of size (dofs, atoms) or (dofs, atoms, beads).
`masses` is the vector of masses for each atom that matches length with the second dimension.
"""
divide_by_mass!(dv, masses) = dv ./= masses'

"""
    velocity!(dr, v, r, sim, t)

Write the velocity `v` into `dr`.
Has extra arguments to work with `Dynamical(O/S)DEProblem`s.
"""
velocity!(dr, v, r, sim, t) = dr .= v

"""
    apply_interbead_coupling!(du::DynamicalVariables, u::DynamicalVariables,
                              sim::RingPolymerSimulation)

Applies the force that arises from the harmonic springs between adjacent beads.

Only applies the force for atoms labelled as quantum within the `RingPolymerParameters`.
"""
function apply_interbead_coupling!(dr::AbstractArray{T,3}, r::AbstractArray{T,3}, sim::RingPolymerSimulation) where {T}
    for i in axes(dr, 3)
        iplus = mod1(i+1, sim.beads.n_beads)
        iminus = mod1(i-1, sim.beads.n_beads)
        for j in sim.beads.quantum_atoms
            for k in axes(dr, 1)
                dr[k,j,i] -= sim.beads.ω_n² * (2r[k,j,i] - r[k,j,iplus] - r[k,j,iminus])
            end
        end
    end
    return nothing
end

function classical_hamiltonian(sim, u) end

function classical_kinetic_energy(sim::AbstractSimulation, u)
    classical_kinetic_energy(sim,  DynamicsUtils.get_velocities(u))
end

function classical_kinetic_energy(sim::Simulation, v::AbstractMatrix)
    kinetic = zero(eltype(v))
    for i in axes(v, 2)
        for j in axes(v, 1)
            kinetic += masses(sim, i) * v[j,i]^2
        end
    end
    return kinetic / 2
end

function classical_kinetic_energy(sim::RingPolymerSimulation, v::AbstractArray{T,3}) where {T}
    kinetic = zero(eltype(v))
    for k in axes(v, 3)
        for i in axes(v, 2)
            for j in axes(v, 1)
                kinetic += masses(sim, i) * v[j,i,k]^2
            end
        end
    end
    return kinetic / 2
end

function classical_potential_energy(sim::AbstractSimulation, u)
    classical_potential_energy(sim,  DynamicsUtils.get_positions(u))
end

function classical_potential_energy(sim::Simulation, r::AbstractMatrix)
    Calculators.get_potential(sim.calculator, r)
end

function classical_potential_energy(sim::RingPolymerSimulation, r::AbstractArray{T,3}) where {T}
    V = Calculators.get_potential(sim.calculator, r)
    sum(V) + RingPolymers.get_spring_energy(sim.beads, masses(sim), r)
end

include("dynamics_variables.jl")
include("callbacks.jl")
export CellBoundaryCallback
export TerminatingCallback

include("density_matrix_dynamics.jl")
include("plot.jl")

end # module


text:
using ReinforcementLearning
using StableRNGs
using Flux
using Flux.Losses

function build_dueling_network(network::Chain)
    lm = length(network)
    if !(network[lm] isa Dense) || !(network[lm-1] isa Dense)
        error("The Qnetwork provided is incompatible with dueling.")
    end
    base = Chain([deepcopy(network[i]) for i=1:lm-2]...)
    last_layer_dims = size(network[lm].weight, 2)
    val = Chain(deepcopy(network[lm-1]), Dense(last_layer_dims, 1))
    adv = Chain([deepcopy(network[i]) for i=lm-1:lm]...)
    return DuelingNetwork(base, val, adv)
end

function RL.Experiment(
    ::Val{:JuliaRL},
    ::Val{:DQN},
    ::Val{:CartPole},
    ::Nothing;
    seed = 123,
)
    rng = StableRNG(seed)
    env = CartPoleEnv(; T = Float32, rng = rng)
    ns, na = length(state(env)), length(action_space(env))

    agent = Agent(
        policy = QBasedPolicy(
            learner = DQNLearner(
                approximator = NeuralNetworkApproximator(
                    model = DuelingNetwork(
                        base = Chain(
                            Dense(ns, 128, relu; init = glorot_uniform(rng)),
                            Dense(128, 128, relu; init = glorot_uniform(rng)),
                        ),
                        val = Dense(128, 1; init = glorot_uniform(rng)),
                        adv = Dense(128, na; init = glorot_uniform(rng)),
                    ),
                    optimizer = ADAM(),
                ) |> gpu,
                target_approximator = NeuralNetworkApproximator(
                    model = DuelingNetwork(
                        base = Chain(
                            Dense(ns, 128, relu; init = glorot_uniform(rng)),
                            Dense(128, 128, relu; init = glorot_uniform(rng)),
                        ),
                        val = Dense(128, 1; init = glorot_uniform(rng)),
                        adv = Dense(128, na; init = glorot_uniform(rng)),
                    ),
                ) |> gpu,
                loss_func = huber_loss,
                stack_size = nothing,
                batch_size = 32,
                update_horizon = 1,
                min_replay_history = 100,
                update_freq = 1,
                target_update_freq = 100,
                rng = rng,
            ),
            explorer = EpsilonGreedyExplorer(
                kind = :exp,
                ϵ_stable = 0.01,
                decay_steps = 500,
                rng = rng,
            ),
        ),
        trajectory = CircularArraySARTTrajectory(
            capacity = 1000,
            state = Vector{Float32} => (ns,),
        ),
    )

    stop_condition = StopAfterStep(10_000, is_show_progress=!haskey(ENV, "CI"))
    hook = TotalRewardPerEpisode()
    Experiment(agent, env, stop_condition, hook, "")
end

using Plots
pyplot() #hide
ex = E`JuliaRL_DQN_CartPole`
run(ex)
plot(ex.hook.rewards)
savefig("assets/JuliaRL_DQN_CartPole.png") #hide

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
using PyPlot
using WignerFamilies
using HalfIntegers
# compute all j₂=100, j₃=60, m₂=70, m₃=-55, m₁=-m₂-m₃
# bottom of page 14 of Raynal et al.
# On the zeros of 3j coefficients: polynomial degree versus recurrence order

X₀, Y₀ = 43, 25
# X₀, Y₀ = 109, 63
X(n) = iszero(n) ? X₀ : 7X(n-1) + 12Y(n-1)
Y(n) = iszero(n) ? Y₀ : 4(n-1) + 7Y(n-1)
n_seq = 2
(2X(n_seq)+1)/6, (X(n_seq)+2)/6
##
# j₂, j₃, m₂, m₃ = 100, 60, 70, -55
j₂, j₃, m₂, m₃ = (2X(n_seq)+1)/6, (X(n_seq)+2)/6, HalfInt(3/2), HalfInt(3/2)
# j₂, j₃, m₂, m₃ = 4, 4, 2, -2

# j₂, j₃, m₂, m₃ = HalfInt(23/2), HalfInt(21/2), HalfInt(3/2), HalfInt(3/2)

m₁ = -m₂ - m₃

w = WignerF(BigFloat, j₂, j₃, m₂, m₃)

nmid = Int(ceil((w.nₘᵢₙ + w.nₘₐₓ)/2))
nmid = iseven(nmid) ? nmid : nmid + 1

w3j = get_wigner_array(w)

n₊ = WignerFamilies.rψ!(w, nmid, w3j)
n₋ = WignerFamilies.sψ!(w, nmid, w3j)

for k in (n₊+1):w.nₘₐₓ  # product the ratios upwards
    w3j[k] *= w3j[k-1]
end
for k in (n₋-1):-1:w.nₘᵢₙ  # product the ratios downwards
    w3j[k] *= w3j[k+1]
end

ψn₋ = w3j[n₋]
ψn₊ = w3j[n₊]

WignerFamilies.ψauxplus!(w, n₋, n₊, w3j)

if n₋ > w.nₘᵢₙ
    scale_middle = ψn₋ / w3j[n₋]
    w3j[w.nₘᵢₙ:(n₋+1)] .*= scale_middle
end
if n₊ < w.nₘₐₓ
    scale_end = w3j[n₊] / ψn₊
    w3j[(n₊+1):w.nₘₐₓ] .*= scale_end
end

norm = WignerFamilies.normalization(w, w3j)
w3j ./= norm
if sign(w3j[w.nₘₐₓ]) != WignerFamilies.f_jmax_sgn(w)
    w3j .*= -1
end

##

plt.clf()
using WignerSymbols
exact = [wigner3j(BigFloat, n, w.j₂, w.j₃, m₁, m₂) for n in w.nₘᵢₙ:w.nₘₐₓ]
plt.plot(w.nₘᵢₙ:w.nₘₐₓ, Float64.(w3j))

plt.plot(w.nₘᵢₙ:w.nₘₐₓ, Float64.(exact), alpha=0.5)
plt.xlabel(raw"$j$")
plt.ylabel(raw"$f(j) \, \times \, 10^3$")
plt.tight_layout()
# plt.xlim(100,120)
gcf()

##
plt.clf()
plt.plot( Float64.(exact .- w3j.parent))
plt.gcf()

##
w = WignerF(BigFloat, 100, 100, 0, 0)
w3j = get_wigner_array(w)
@time WignerFamilies.classical_wigner3j_m0!(w, w3j)
##


text:
using Test
using ForwardDiff
using JLD2
using Symbolics
using StaticArrays
using LinearAlgebra
using Random
using SparseArrays
using BenchmarkTools
using ContactImplicitMPC

# Solver
include("solver/qr.jl")
include("solver/lu.jl")
include("solver/ldl.jl")
include("solver/schur.jl")

# Dynamics
include("dynamics/model.jl")
include("dynamics/particle.jl")
include("dynamics/quadruped.jl")

# Simulator
include("simulator/rotations.jl")
include("simulator/environment.jl")
include("simulator/particle.jl")
include("simulator/hopper_3D.jl")
include("simulator/quadruped.jl")
include("simulator/open_loop.jl")

# Controller
include("controller/objective.jl")
include("controller/linearized_step.jl")
include("controller/implicit_dynamics.jl")
include("controller/linearized_solver.jl")
include("controller/newton.jl")
include("controller/newton_structure_solver.jl")
include("controller/trajectory.jl")

# MPC examples
include("controller/mpc_quadruped.jl")
include("controller/mpc_flamingo.jl") #TODO fix newton structure solver test


text:
function load_humangem_mat_model()
    human_file = rawdir(Human1, [ "Human1_Publication_Data_Scripts/tINIT_GEMs/data" ], "HumanGEM.mat")
    mat_model = MAT.matread(human_file)["ihuman"];
    mat_model = MetNets.to_symbol_dict(mat_model)
    MetNets.reshape_mat_dict(mat_model)
end

function load_humangem_raw_model(model_dict = load_humangem_mat_model())
    return MetNets.MetNet(model_dict; reshape = true)
end

save_humangem_base_model(model; verbose = true) =
    sdat(HumanGEM, MetNets.compressed_model(model),
        "HumanGEM_base_model", ".jls";
        verbose
    )

function load_humangem_base_model(;uncompress = false)
    base_model = ldat(HumanGEM,
        "HumanGEM_base_model", ".jls";
        verbose = false
    )
    return uncompress ? MetNets.uncompressed_model(base_model) : base_model
end


text:
fname = :TRKOp
c = Combi(hessian_operator, PDataK, solve_modelKTR, preprocessKTR, decreaseKTR, TparamsKTR())
include("Template.jl")


text:
struct Params{T <: Real, A, W}
                      dt :: T
                       g :: T
             density_ice :: T
           density_ocean :: T
               gas_const :: T
           sec_per_year  :: T
       default_thickness :: T
       default_viscosity :: T
     default_temperature :: T
          default_damage :: T
       accumulation_rate :: A
glen_a_activation_energy :: T
              glen_a_ref :: T
    glen_temperature_ref :: T
                  glen_n :: T
    glen_reg_strain_rate :: T
              weertman_c :: W
              weertman_m :: T
      weertman_reg_speed :: T
     sea_level_wrt_geoid :: T
       minimum_thickness :: T
           evolveShelves :: Bool
                smallHAF :: T
end

"""
Params(; <kwargs>)

Construct a WAVI.jl parameters object for holding physical parameters.

Keyword arguments
=================
- 'dt': model timestep (NB: simulation timestep set in timestepping params, this value is updated when model embedded to the value specified in timestepping_params when passed to simulation)
- 'g': gravitational acceleration (m^2 / s)
- 'density_ice': ice density (kg / m^3)
- 'density_ocean' ocean water density (kg / m^3)
- 'gas_const': gas constant in glen b calculation
- 'sec_per_year': seconds per year (s)
- 'default_thickness': thickness value reverted to if no initial thickness passed (m)
- 'default_viscosity': viscosity value reverted to if no initial thickness passed (Pa s)
- 'default_temperature': temperature value reverted to if no initial thickness passed (K)
- 'default_damage': damage value reverted to if no initial thickness passed (dimensionless)
- 'accumulation_rate': uniform accumulation_rate (m/yr)
- 'glen_a_activation_energy': activation energy in glen b calculation
- 'glen_a_ref': glen a reference value used in glen b calculation
- 'glen_temperature_ref': reference temperature using in glen b calculation
- 'glen_n': exponent in glen b calculation
- 'glen_reg_strain_rate': strain rate regularization value
- 'weertman_c': basal sliding field of coefficients
- 'weertman_m': sliding law exponent
- 'weertman_reg_speed': regularization speed, used to prevent bed speed going to zero
- 'sea_level_wrt_geoid': reference sea level
- 'minimum_thickness': minimum ice thickness on model domain
- 'evolveShelves': flag for turning on and off the evolution of the shelves in the forward run_simulation
- 'smallHAF': small value of HAF used within update_thickness when not evolving shelves
"""
function Params(; g = 9.81,
                  density_ice = 918.0,
                  density_ocean = 1028.0,
                  gas_const = 8.314,
                  sec_per_year =3600*24*365.25,
                  default_thickness= 100.,
                  default_viscosity= 1.0e7,
                  default_temperature=265.700709,
                  default_damage= 0.0,
                  accumulation_rate= 0.0,
                  glen_a_activation_energy = 5.8631e+04,
                  glen_a_ref= 4.9e-16 *sec_per_year * 1.0e-9,
                  glen_temperature_ref= 263.15,
                  glen_n = 3.0,
                  glen_reg_strain_rate = 1.0e-5,
                  weertman_c = 1.0e4,
                  weertman_m  = 3.0,
                  weertman_reg_speed = 1.0e-5,
                  sea_level_wrt_geoid  = 0.0,
                  minimum_thickness = 50.0,
                  evolveShelves = true,
                  smallHAF = 1.0)

  #defualt the timestep to 1.0 (will be updated when the model is embedded in a simulation)
  dt = 1.0

  return Params(
                  dt,
                  g,
                  density_ice,
                  density_ocean,
                  gas_const,
                  sec_per_year,
                  default_thickness,
                  default_viscosity,
                  default_temperature,
                  default_damage,
                  accumulation_rate,
                  glen_a_activation_energy,
                  glen_a_ref,
                  glen_temperature_ref,
                  glen_n,
                  glen_reg_strain_rate,
                  weertman_c,
                  weertman_m,
                  weertman_reg_speed,
                  sea_level_wrt_geoid,
                  minimum_thickness,
                  evolveShelves,
                  smallHAF
                  )
end


text:
# Support for using a matrix as a AbstractVectorOfArray, especially iteration on A[i], etc.
# Note that a proper implementation should probably work with an arbitrary Array, so that a 3-tensor
# appears as a vector of matrices, etc.
struct MatrixVectorOfArray{T,N,A} <: RecursiveArrayTools.AbstractVectorOfArray{T,N,A}
    u::A # A <: AbstractMatrix{T} where size(u, 1) == N
end

Base.Array(VA::MatrixVectorOfArray{T,N,A}) where {T,N,A} = VA.u

MatrixVectorOfArray(mat::AbstractMatrix) = MatrixVectorOfArray{Float64,2,typeof(mat)}(mat)

# Interface for the linear indexing. This is just a view of the underlying nested structure
@inline Base.firstindex(VA::MatrixVectorOfArray) = firstindex(VA.u, 2)
@inline Base.lastindex(VA::MatrixVectorOfArray) = lastindex(VA.u, 2)

@inline Base.length(VA::MatrixVectorOfArray) = size(VA.u, 2)
@inline Base.eachindex(VA::MatrixVectorOfArray) = Base.OneTo(size(VA.u, 2))
@inline Base.IteratorSize(VA::MatrixVectorOfArray) = Base.HasLength()
# Linear indexing will be over the container elements, not the individual elements
# unlike an true AbstractArray
Base.@propagate_inbounds function Base.getindex(VA::MatrixVectorOfArray{T,N}, I::Int) where {T,N}
    return view(VA.u, :, I)
end
Base.@propagate_inbounds function Base.getindex(VA::MatrixVectorOfArray{T,N}, I::Colon) where {T,N}
    return view(VA.u, :, I)
end
Base.@propagate_inbounds function Base.getindex(VA::MatrixVectorOfArray{T,N}, i::Int,
                                                ::Colon) where {T,N}
    return [VA.u[i, j] for j in 1:length(VA)]
end
Base.@propagate_inbounds function Base.getindex(VA::MatrixVectorOfArray{T,N},
                                                ii::CartesianIndex) where {T,N}
    ti = Tuple(ii)
    i = last(ti)
    jj = CartesianIndex(Base.front(ti))
    return VA.u[jj, i]
end
Base.@propagate_inbounds function Base.setindex!(VA::MatrixVectorOfArray{T,N}, v,
                                                 I::Int) where {T,N}
    return VA.u[:, I] = v
end
Base.@propagate_inbounds function Base.setindex!(VA::MatrixVectorOfArray{T,N}, v,
                                                 I::Colon) where {T,N}
    return VA.u[:, I] = v
end
Base.@propagate_inbounds function Base.setindex!(VA::MatrixVectorOfArray{T,N}, v,
                                                 I::AbstractArray{Int}) where {T,N}
    return VA.u[:, I] = v
end
Base.@propagate_inbounds function Base.setindex!(VA::MatrixVectorOfArray{T,N}, v, i::Int,
                                                 ::Colon) where {T,N}
    for j in 1:length(VA)
        VA.u[i, j] = v[j]
    end
    return v
end
Base.@propagate_inbounds function Base.setindex!(VA::MatrixVectorOfArray{T,N}, x,
                                                 ii::CartesianIndex) where {T,N}
    ti = Tuple(ii)
    i = last(ti)
    jj = CartesianIndex(Base.front(ti))
    return VA.u[jj, i] = x
end

# Interface for the two-dimensional indexing, a more standard AbstractArray interface
@inline Base.size(VA::MatrixVectorOfArray) = size(VA.u)
Base.@propagate_inbounds function Base.getindex(VA::MatrixVectorOfArray{T,N}, I::Int...) where {T,N}
    return VA.u[I...]
end
Base.@propagate_inbounds function Base.setindex!(VA::MatrixVectorOfArray{T,N}, v,
                                                 I::Int...) where {T,N}
    return VA.u[I...] = v
end

# The iterator will be over the subarrays of the container, not the individual elements
# unlike an true AbstractArray
# using view just as eachcol(VA.u) would on the undrelying matrix
function Base.iterate(VA::MatrixVectorOfArray, state = 1)
    return state >= size(VA.u, 2) + 1 ? nothing : (view(VA, :, state), state + 1)
end

# function Base.append!(VA::AbstractVectorOfArray{T,N}, new_item::AbstractVectorOfArray{T,N}) where {T,N}
#     for item in copy(new_item)
#         push!(VA, item)
#     end
#     return VA
# end

# Tools for creating similar objects
@inline function Base.similar(VA::MatrixVectorOfArray, ::Type{T} = eltype(VA)) where {T}
    return MatrixVectorOfArray(similar(VA.u))
end

# fill!
# For DiffEqArray it ignores ts and fills only u
Base.fill!(VA::MatrixVectorOfArray, x) = fill!(VA.u, x)

# make it show just like its data
function Base.show(io::IO, m::MIME"text/plain", x::MatrixVectorOfArray)
    return (println(io, summary(x), ':'); show(io, m, x.u))
end
Base.summary(A::MatrixVectorOfArray) = string("MatrixVectorOfArray{", eltype(A), ",", ndims(A), "}")

Base.map(f, A::MatrixVectorOfArray) = map(f, eachcol(A.u))
function Base.mapreduce(f, op, A::MatrixVectorOfArray)
    return mapreduce(f, op, (mapreduce(f, op, x) for x in eachcol(A.u)))
end

# Promotion utilities
promote_vv(vec::AbstractVectorOfArray) = vec
promote_vv(mat::AbstractMatrix) = MatrixVectorOfArray(mat)
promote_vv(vec::AbstractVector) = VectorOfArray(vec)
promote_vv(::Nothing) = nothing


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ParameterDefinitionsValue <: SwaggerModel
    type::Any # spec type: Union{ Nothing, String } # spec name: type
    allowedValues::Any # spec type: Union{ Nothing, Vector{Any} } # spec name: allowedValues
    defaultValue::Any # spec type: Union{ Nothing, Any } # spec name: defaultValue
    metadata::Any # spec type: Union{ Nothing, Dict{String, Any} } # spec name: metadata

    function ParameterDefinitionsValue(;type=nothing, allowedValues=nothing, defaultValue=nothing, metadata=nothing)
        o = new()
        validate_property(ParameterDefinitionsValue, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        validate_property(ParameterDefinitionsValue, Symbol("allowedValues"), allowedValues)
        setfield!(o, Symbol("allowedValues"), allowedValues)
        validate_property(ParameterDefinitionsValue, Symbol("defaultValue"), defaultValue)
        setfield!(o, Symbol("defaultValue"), defaultValue)
        validate_property(ParameterDefinitionsValue, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        o
    end
end # type ParameterDefinitionsValue

const _property_map_ParameterDefinitionsValue = Dict{Symbol,Symbol}(Symbol("type")=>Symbol("type"), Symbol("allowedValues")=>Symbol("allowedValues"), Symbol("defaultValue")=>Symbol("defaultValue"), Symbol("metadata")=>Symbol("metadata"))
const _property_types_ParameterDefinitionsValue = Dict{Symbol,String}(Symbol("type")=>"String", Symbol("allowedValues")=>"Vector{Any}", Symbol("defaultValue")=>"Any", Symbol("metadata")=>"Dict{String, Any}")
Base.propertynames(::Type{ ParameterDefinitionsValue }) = collect(keys(_property_map_ParameterDefinitionsValue))
Swagger.property_type(::Type{ ParameterDefinitionsValue }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ParameterDefinitionsValue[name]))}
Swagger.field_name(::Type{ ParameterDefinitionsValue }, property_name::Symbol) =  _property_map_ParameterDefinitionsValue[property_name]

const _allowed_ParameterDefinitionsValue_type = ["String", "Array", "Object", "Boolean", "Integer", "Float", "DateTime"]

function check_required(o::ParameterDefinitionsValue)
    true
end

function validate_property(::Type{ ParameterDefinitionsValue }, name::Symbol, val)
    if name === Symbol("type")
        Swagger.validate_param(name, "ParameterDefinitionsValue", :enum, val, _allowed_ParameterDefinitionsValue_type)
    end
end


text:
module AffyCelFiles

#http://media.affymetrix.com/support/developer/powertools/changelog/gcos-agcc/cel.html

struct DataHeader
    isDefined::Bool
    dataTypeIdentifier::String
    uniqueID::String
    creationDateTime::String
    locale::String
    nParameters::Int
    names::Array{String,1}
    rawvalues::Array{String,1}
    values::Array{String,1}
    types::Array{String,1}
    nParentFileHeaders::Int
    parentFileHeaders::Array{DataHeader,1}
    function DataHeader(isDefined=false)
        new(
            isDefined
        )
    end
    function DataHeader(args...)
        new(
            args...
        )
    end
end

struct DataSet
    isDefined::Bool
    dataSetName::String
    nParameters::Int
    parameterNames::Array{String,1}
    parameterRawvalues::Array{String,1}
    parameterValues::Array{String,1}
    parameterTypes::Array{String,1}
    dataSetColNames::Array{String,1}
    dataSetColTypeCodes::Array{Int8,1}
    dataSetColTypeSizes::Array{Int32,1}
    dataSetColTypes::Array{Type,1}
    dataSetColumns::Dict{Int,AbstractArray}
    function DataSet(isDefined=false)
        new(
            isDefined
        )
    end
    function DataSet(args...)
        new(
            args...
        )
    end
end

struct DataGroup
    isDefined::Bool
    nDataSets::Int
    dataSets::Array{DataSet,1}
    function DataGroup(isDefined=false)
        new(
            isDefined
        )
    end
    function DataGroup(args...)
        new(
            args...
        )
    end
end

struct Cel
    isDefined::Bool
    nGroups::Int
    dataHeader::DataHeader
    dataGroups::Array{DataGroup,1}
    function Cel(isDefined=false)
		new(
			isDefined
		)
	end
    function Cel(args...)
		new(
			args...
		)
	end
end

function cel_read(file::AbstractString)
	io=open(file,"r")
	cel=cel_read(io)
	close(io)
	cel
end

function cel_read(io::IO)::Cel
	seekstart(io)
    magic=read(io,Int8)
    if magic == 59
        #Command Console generic data file format
        #  http://media.affymetrix.com/support/developer/powertools/changelog/gcos-agcc/generic.html
        return cel_read_generic(io)
    end

	seekstart(io)
    magic=read(io,Int32)
    if magic == 64
        #Version 4 Format
        #  http://media.affymetrix.com/support/developer/powertools/changelog/gcos-agcc/cel.html
        return cel_read_v4(io)
    end

	@warn "Unknown CEL file version"
	return Cel()
end

function cel_read_v4(io::IO)
    #not implemented, need example file
    @warn "CEL file version 4 not implemented, need example file"
	return Cel()
end

function cel_read_generic(io::IO)
    version=read(io,Int8)
    if version != 1
        @warn "Version number of Command Console generic data file is "*string(version)*". Expected 1."
    end
    nGroups=ntoh(read(io,Int32))
    firstGroupPos=ntoh(read(io,UInt32))

    dataHeader=cel_read_generic_header(io)
    #println(firstGroupPos)
    #println(position(io))

    if firstGroupPos != position(io)
        @warn "current io position is "*string(position(io))*", should be first data group at "*string(firstGroupPos)*", performing seek to "*string(firstGroupPos)
        seek(io,firstGroupPos)
    end

    dataGroups=Array{DataGroup}(undef,nGroups)
    dataGroupsIndex=1
    nextGroupPos=firstGroupPos
    while nextGroupPos > 0 && dataGroupsIndex <= nGroups
        nextGroupPos=ntoh(read(io,UInt32))
        nextDataSet=ntoh(read(io,UInt32))
        nDataSets=ntoh(read(io,Int32))
        groupName=read_wstring(io)
        dataSets=Array{DataSet}(undef,nDataSets)
        dataSetsIndex=1
        while dataSetsIndex <= nDataSets
            if nextDataSet != position(io)
                #@warn "current io position is "*string(position(io))*", should be next dataset at "*string(nextDataSet)*", performing seek to "*string(nextDataSet)
                seek(io,nextDataSet)
            end
            dataElementPos=ntoh(read(io,UInt32))
            nextDataSet=ntoh(read(io,UInt32))
            dataSetName=read_wstring(io)
            nParameters=ntoh(read(io,Int32))
            names=Array{String}(undef,nParameters)
            rawvalues=Array{String}(undef,nParameters)
            values=Array{String}(undef,nParameters)
            types=Array{String}(undef,nParameters)
            for i in 1:nParameters
                names[i]=read_wstring(io)
                rawvalues[i]=read_value(io)
                types[i]=read_wstring(io)
                values[i]=calvin_convert(rawvalues[i], Val(Symbol(types[i])) )
            end
            nColumns=ntoh(read(io,UInt32))
            colNames=Array{String}(undef,nColumns)
            colTypeCodes=Array{Int8}(undef,nColumns)
            colTypeSizes=Array{Int32}(undef,nColumns)
            colTypes=Array{Type}(undef,nColumns)
            rowSize=0
            for i in 1:nColumns
                colNames[i]=read_wstring(io)
                colTypeCodes[i]=read(io,Int8)
                colSize=ntoh(read(io,Int32))
                colTypeSizes[i]=colSize
                rowSize+=colSize
                colTypes[i]=decode_type(colTypeCodes[i],colSize)
            end
            nRows=ntoh(read(io,UInt32))
            if dataElementPos != position(io)
                @warn "current io position is "*string(position(io))*", should be next data element at "*string(dataElementPos)*", performing seek to "*string(dataElementPos)
                seek(io,dataElementPos)
            end
            dataSetColumns=Dict{Int,AbstractArray}()
            buffer=zeros(UInt8,nRows*rowSize)
            readbytes!(io,buffer,nRows*rowSize)
            #permBuffer=PermutedDimsArray(reshape(buffer,(rowSize,nRows)),(2,1))
            #startColIndex=0
            #endColIndex=0
            #for i in 1:nColumns
            #    startColIndex=endColIndex+1
            #    endColIndex=startColIndex+colTypeSizes[i]-1
            #    if colTypeCodes[i] <= 6
            #        column=map( row -> reinterpret(colTypes[i],reverse(row[startColIndex:endColIndex]))[1], eachrow(permBuffer))
            #    elseif colTypeCodes[i] == 7 || colTypeCodes[i] == 8
            #        @warn "type string not implemented"
            #        column=Array{String}(undef,0)
            #    else
            #        @error "unknown type code, code="*string(colTypeCodes[i])
            #        column=Array{String}(undef,0)
            #    end
            #    dataSetColumns[i]=column
            #end
            for i in 1:nColumns
                if colTypeCodes[i] <= 6
                    dataSetColumns[i]=Array{colTypes[i]}(undef,nRows)
                elseif colTypeCodes[i] == 7 || colTypeCodes[i] == 8
                    @warn "type string not implemented"
                    column=Array{String}(undef,0)
                else
                    @error "unknown type code, code="*string(colTypeCodes[i])
                    column=Array{String}(undef,0)
                end
            end
            for j in 1:nRows
                startColIndex=(j-1)*rowSize
                endColIndex=startColIndex
                for i in 1:nColumns
                    startColIndex=endColIndex+1
                    endColIndex=startColIndex+colTypeSizes[i]-1
                    if colTypeCodes[i] <= 6
                        dataSetColumns[i][j]=reinterpret(colTypes[i],reverse(buffer[startColIndex:endColIndex]))[1]
                    end
                end
            end
            dataSets[dataSetsIndex]=DataSet(true,dataSetName,nParameters,names,rawvalues,values,types,colNames,colTypeCodes,colTypeSizes,colTypes,dataSetColumns)
            dataSetsIndex+=1
        end
        dataGroups[dataGroupsIndex]=DataGroup(true,nDataSets,dataSets)
        dataGroupsIndex+=1
    end
    return Cel(true,nGroups,dataHeader,dataGroups)
end

function decode_type(c::Int8,s::Int32)
    if c==0 && s==1
        return Int8
    elseif  c==1 && s==1
        return UInt8
    elseif  c==2 && s==2
        return Int16
    elseif  c==3 && s==2
        return UInt16
    elseif  c==4 && s==4
        return Int32
    elseif  c==5 && s==4
        return UInt32
    elseif  c==6 && s==2
        return Float16
    elseif  c==6 && s==4
        return Float32
    elseif  c==6 && s==8
        return Float64
    elseif  c==7 || c==8
        return String
    else
        @error "unknown type code and size combination, code="*string(c)*" and size="*string(s)*", falling back to String"
    end
    return String
end

function cel_read_generic_header(io::IO)
    #read generic data header
    data_type_identifier=read_string(io)
    guid=read_string(io)
    creation=read_string(io)
    locale=read_wstring(io)
    nParameters=ntoh(read(io,Int32))
    names=Array{String}(undef,nParameters)
    rawvalues=Array{String}(undef,nParameters)
    values=Array{String}(undef,nParameters)
    types=Array{String}(undef,nParameters)
    for i in 1:nParameters
        #println(position(io))
        names[i]=read_wstring(io)
        rawvalues[i]=read_value(io)
        types[i]=read_wstring(io)
        values[i]=calvin_convert(rawvalues[i], Val(Symbol(types[i])) )
    end
    nParentFileHeaders=ntoh(read(io,Int32))
    parentFileHeaders=Array{DataHeader}(undef,nParentFileHeaders)
    if nParentFileHeaders > 0
        for i in 1:nParentFileHeaders
            parentFileHeaders[i]=cel_read_generic_header(io)
        end
    end

    return DataHeader(
        true,
        data_type_identifier,
        guid,
        creation,
        locale,
        nParameters,
        names,
        rawvalues,
        values,
        types,
        nParentFileHeaders,
        parentFileHeaders
    )
end

function calvin_convert(v::String,::Val{Symbol("text/x-calvin-integer-8")})
    return string(reinterpret(Int8,UInt16.(collect(v)[2:-1:1]))[1])
end

function calvin_convert(v::String,::Val{Symbol("text/x-calvin-unsigned-integer-8")})
    return string(reinterpret(UInt8,UInt16.(collect(v)[2:-1:1]))[1])
end

function calvin_convert(v::String,::Val{Symbol("text/x-calvin-integer-16")})
    return string(reinterpret(Int16,UInt16.(collect(v)[2:-1:1]))[1])
end

function calvin_convert(v::String,::Val{Symbol("text/x-calvin-unsigned-integer-16")})
    return string(reinterpret(UInt16,UInt16.(collect(v)[2:-1:1]))[1])
end

function calvin_convert(v::String,::Val{Symbol("text/x-calvin-integer-32")})
    return string(reinterpret(Int32,UInt16.(collect(v)[2:-1:1]))[1])
end

function calvin_convert(v::String,::Val{Symbol("text/x-calvin-unsigned-integer-32")})
    return string(reinterpret(UInt32,UInt16.(collect(v)[2:-1:1]))[1])
end

function calvin_convert(v::String,::Val{Symbol("text/x-calvin-float")})
    return string(reinterpret(Float32,UInt16.(collect(v)[2:-1:1]))[1])
end

function calvin_convert(v::String,::Val{Symbol("text/plain")})
    chars=collect(v)
    eos=findfirst(c->c=='\0',chars)
    return isnothing(eos) ? v : v[1:(eos-1)]
end

function calvin_convert(v::String,::Val{Symbol("text/ascii")})
    chars=Char.(reinterpret(UInt8,UInt16.(collect(v))))
    eos=findfirst(c->c=='\0',chars)
    return isnothing(eos) ? join(chars) : join(chars[1:(eos-1)])
end

function read_string(io::IO)
    length=ntoh(read(io,Int32))
    c=read(io,length)
    return join(Array{Char}(c))
end

function read_wstring(io::IO)
    length=2*ntoh(read(io,Int32))
    b=zeros(UInt8,length)
    readbytes!(io,b,length)
    b=ntoh.(reinterpret(UInt16,b))
    return transcode(String,b)
end

function read_value(io::IO)
    length=ntoh(read(io,Int32))
    bufferLength=length
    if !iseven(length)
        bufferLength+=1
    end
    b=zeros(UInt8,bufferLength)
    readbytes!(io,b,length)
    b=ntoh.(reinterpret(UInt16,b))
    return transcode(String,b)
end

end # module

max_stars_repo_name: oheil/AffyCelFiles

max_issues_repo_name: oheil/AffyCelFiles

max_forks_repo_name: oheil/AffyCelFiles


text:
#=
Implement locking in a binary tree. A binary tree node can be locked or unlocked only if all of its descendants or ancestors are not locked.

Design a binary tree node class with the following methods:

    is_locked, which returns whether the node is locked
    lock, which attempts to lock the node. If it cannot be locked, then it should return false. Otherwise, it should lock it and return true.
    unlock, which unlocks the node. If it cannot be unlocked, then it should return false. Otherwise, it should unlock it and return true.

You may augment the node to add parent pointers or any other property you would like. You may assume the class is used in a single-threaded program, so there is no need for actual locks or mutexes. Each method should run in O(h), where h is the height of the tree.
=#

mutable struct Node
    key::Int
    left::Union{Node, Nothing}
    right::Union{Node, Nothing}
    locked::Bool
end

function is_locked(node::Node)
    return node.locked
end

function is_lockable(node::Node)
    if isnothing(node.right) && isnothing(node.left)
        return true
    elseif isnothing(node.right) && is_locked(node.left)
        return is_lockable(node.left)
    elseif is_locked(node.right) && isnothing(node.left)
        return is_lockable(node.right)
    elseif is_locked(node.right) && is_locked(node.left)
        return is_lockable(node.right) && is_lockable(node.left)
    else
        return false
    end
end

function lock(node::Node)
    if is_lockable(node)
        setfield!(node, :locked, true)
        return true
    else
        return is_lockable(node)
    end
end

# for the unlock function, I need a way to lock for the ancestors....

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
mutable struct NoObjModel <: AbstractModel
    settings::Settings
    constraints::Vector{Constraint}
    obj::NoObjective
    output::Output

    NoObjModel(settings, constraints, obj, output) = new(settings, constraints, obj, output)
    NoObjModel() = new(Settings(), Constraint[], NoObjective(), Output())
end


text:
using StructuralCausalModels
import StructuralCausalModels.SymbolList

import Base.-

-(x::SymbolList, y::SymbolList) = (type=:line, x=x, y=y)

↔(x, y) = (type=:darrow, x=x, y=y)

→(x, y) = (type=:rarrow, x=x, y=y)

println(:b - :A)
println(:b ↔ [:A, :D])
println([:X, :T] → :A)

m = (:b - :A, [:X, :T] → :A, :b ↔ [:A, :D])

m[3] |> display


text:
using Test
import ClimaCore: Domains, Meshes, Topologies
import ClimaCore.Geometry: Geometry
using ClimaCore.Meshes: equispaced_rectangular_mesh
using StaticArrays
using IntervalSets

function rectangular_grid(
    n1,
    n2,
    x1periodic,
    x2periodic;
    x1min = 0.0,
    x1max = 1.0,
    x2min = 0.0,
    x2max = 1.0,
)
    domain = Domains.RectangleDomain(
        Geometry.XPoint(x1min)..Geometry.XPoint(x1max),
        Geometry.YPoint(x2min)..Geometry.YPoint(x2max),
        x1periodic = x1periodic,
        x2periodic = x2periodic,
        x1boundary = x1periodic ? nothing : (:west, :east),
        x2boundary = x2periodic ? nothing : (:south, :north),
    )
    #=
    domain = Domains.RectangleDomain(
        x1min..x1max,
        x2min..x2max,
        x1periodic = x1periodic,
        x2periodic = x2periodic,
        x1boundary = x1periodic ? nothing : (:west, :east),
        x2boundary = x2periodic ? nothing : (:south, :north),
    )
    =#
    mesh = equispaced_rectangular_mesh(domain, n1, n2)
    grid_topology = Topologies.Grid2DTopology(mesh)
    return (domain, mesh, grid_topology)
end

@testset "simple rectangular grid opposing face" begin

    @testset "assert correct element numbering" begin
        _, _, grid_topology = rectangular_grid(1, 1, true, true)
        @test_throws AssertionError Topologies.opposing_face(
            grid_topology,
            0,
            1,
        )
        @test_throws AssertionError Topologies.opposing_face(
            grid_topology,
            2,
            1,
        )
    end

    @testset "1×1 element quad mesh with all periodic boundries" begin
        _, _, grid_topology = rectangular_grid(1, 1, true, true)
        @test Topologies.opposing_face(grid_topology, 1, 1) == (1, 2, false)
        @test Topologies.opposing_face(grid_topology, 1, 2) == (1, 1, false)
        @test Topologies.opposing_face(grid_topology, 1, 3) == (1, 4, false)
        @test Topologies.opposing_face(grid_topology, 1, 4) == (1, 3, false)
    end

    @testset "1×1 element quad mesh with 1 periodic boundary" begin
        _, _, grid_topology = rectangular_grid(1, 1, true, false)
        @test Topologies.opposing_face(grid_topology, 1, 1) == (1, 2, false)
        @test Topologies.opposing_face(grid_topology, 1, 2) == (1, 1, false)
        @test Topologies.opposing_face(grid_topology, 1, 3) == (0, 3, false)
        @test Topologies.opposing_face(grid_topology, 1, 4) == (0, 4, false)
    end

    @testset "1×1 element quad mesh with non-periodic boundaries" begin
        _, _, grid_topology = rectangular_grid(1, 1, false, false)
        @test Topologies.opposing_face(grid_topology, 1, 1) == (0, 1, false)
        @test Topologies.opposing_face(grid_topology, 1, 2) == (0, 2, false)
        @test Topologies.opposing_face(grid_topology, 1, 3) == (0, 3, false)
        @test Topologies.opposing_face(grid_topology, 1, 4) == (0, 4, false)
    end

    @testset "2×2 element quad mesh with non-periodic boundaries" begin
        _, _, grid_topology = rectangular_grid(2, 2, false, false)
        @test Topologies.opposing_face(grid_topology, 1, 1) == (0, 1, false)
        @test Topologies.opposing_face(grid_topology, 1, 2) == (2, 1, false)
        @test Topologies.opposing_face(grid_topology, 1, 3) == (0, 3, false)
        @test Topologies.opposing_face(grid_topology, 1, 4) == (3, 3, false)
        @test Topologies.opposing_face(grid_topology, 2, 1) == (1, 2, false)
        @test Topologies.opposing_face(grid_topology, 2, 2) == (0, 2, false)
        @test Topologies.opposing_face(grid_topology, 2, 3) == (0, 3, false)
        @test Topologies.opposing_face(grid_topology, 2, 4) == (4, 3, false)
    end
end

@testset "simple rectangular grid interior faces iterator" begin
    @testset "1×1 element quad mesh with all periodic boundries" begin
        _, _, grid_topology = rectangular_grid(1, 1, true, true)
        @test length(Topologies.interior_faces(grid_topology)) == 2
        faces = collect(Topologies.interior_faces(grid_topology))
        @test faces[1] == (1, 1, 1, 2, false)
        @test faces[2] == (1, 3, 1, 4, false)
    end
    @testset "1×1 element quad mesh with 1 periodic boundary" begin
        _, _, grid_topology = rectangular_grid(1, 1, true, false)
        @test length(Topologies.interior_faces(grid_topology)) == 1
        faces = collect(Topologies.interior_faces(grid_topology))
        @test faces[1] == (1, 1, 1, 2, false)
    end
    @testset "1×1 element quad mesh with non-periodic boundaries" begin
        _, _, grid_topology = rectangular_grid(1, 1, false, false)
        @test length(Topologies.interior_faces(grid_topology)) == 0
        faces = collect(Topologies.interior_faces(grid_topology))
        @test isempty(faces)
    end
    @testset "2×2 element quad mesh with non-periodic boundaries" begin
        _, _, grid_topology = rectangular_grid(2, 2, false, false)
        @test length(Topologies.interior_faces(grid_topology)) == 4
        faces = collect(Topologies.interior_faces(grid_topology))
        @test faces[1] == (2, 1, 1, 2, false) || faces[1] == (1, 2, 2, 1, false)
        @test faces[2] == (3, 3, 1, 4, false) || faces[2] == (3, 2, 4, 1, false)
        @test faces[3] == (4, 1, 3, 2, false) || faces[3] == (1, 4, 3, 3, false)
        @test faces[4] == (4, 3, 2, 4, false) || faces[4] == (2, 4, 4, 3, false)
    end
end

@testset "simple rectangular grid boundary faces iterator" begin
    @testset "1×1 element quad mesh with all periodic boundries" begin
        _, _, grid_topology = rectangular_grid(1, 1, true, true)
        @test length(Topologies.boundary_faces(grid_topology, 1)) == 0
        @test length(Topologies.boundary_faces(grid_topology, 2)) == 0
        @test length(Topologies.boundary_faces(grid_topology, 3)) == 0
        @test length(Topologies.boundary_faces(grid_topology, 4)) == 0
        @test isempty(collect(Topologies.boundary_faces(grid_topology, 1)))
        @test isempty(collect(Topologies.boundary_faces(grid_topology, 2)))
        @test isempty(collect(Topologies.boundary_faces(grid_topology, 3)))
        @test isempty(collect(Topologies.boundary_faces(grid_topology, 4)))
    end
    @testset "1×1 element quad mesh with 1 periodic boundary" begin
        _, _, grid_topology = rectangular_grid(1, 1, true, false)
        @test length(Topologies.boundary_faces(grid_topology, 1)) == 0
        @test length(Topologies.boundary_faces(grid_topology, 2)) == 0
        @test length(Topologies.boundary_faces(grid_topology, 3)) == 1
        @test length(Topologies.boundary_faces(grid_topology, 4)) == 1
        @test isempty(collect(Topologies.boundary_faces(grid_topology, 1)))
        @test isempty(collect(Topologies.boundary_faces(grid_topology, 2)))
        @test collect(Topologies.boundary_faces(grid_topology, 3)) == [(1, 3)]
        @test collect(Topologies.boundary_faces(grid_topology, 4)) == [(1, 4)]
    end
    @testset "1×1 element quad mesh with non-periodic boundaries" begin
        _, _, grid_topology = rectangular_grid(1, 1, false, false)
        @test length(Topologies.boundary_faces(grid_topology, 1)) == 1
        @test length(Topologies.boundary_faces(grid_topology, 2)) == 1
        @test length(Topologies.boundary_faces(grid_topology, 3)) == 1
        @test length(Topologies.boundary_faces(grid_topology, 4)) == 1
        @test collect(Topologies.boundary_faces(grid_topology, 1)) == [(1, 1)]
        @test collect(Topologies.boundary_faces(grid_topology, 2)) == [(1, 2)]
        @test collect(Topologies.boundary_faces(grid_topology, 3)) == [(1, 3)]
        @test collect(Topologies.boundary_faces(grid_topology, 4)) == [(1, 4)]
        @test Topologies.boundary_tag(grid_topology, :west) == 1
        @test Topologies.boundary_tag(grid_topology, :east) == 2
        @test Topologies.boundary_tag(grid_topology, :south) == 3
        @test Topologies.boundary_tag(grid_topology, :north) == 4
    end
    @testset "2×3 element quad mesh with non-periodic boundaries" begin
        _, _, grid_topology = rectangular_grid(2, 3, false, false)
        @test length(Topologies.boundary_faces(grid_topology, 1)) == 3
        @test length(Topologies.boundary_faces(grid_topology, 2)) == 3
        @test length(Topologies.boundary_faces(grid_topology, 3)) == 2
        @test length(Topologies.boundary_faces(grid_topology, 4)) == 2

        @test collect(Topologies.boundary_faces(grid_topology, 1)) ==
              [(1, 1), (3, 1), (5, 1)]
        @test collect(Topologies.boundary_faces(grid_topology, 2)) ==
              [(2, 2), (4, 2), (6, 2)]
        @test collect(Topologies.boundary_faces(grid_topology, 3)) ==
              [(1, 3), (2, 3)]
        @test collect(Topologies.boundary_faces(grid_topology, 4)) ==
              [(5, 4), (6, 4)]
    end
end

@testset "simple rectangular grid vertex iterator" begin
    @testset "1×1 element quad mesh with all periodic boundries" begin
        _, _, grid_topology = rectangular_grid(1, 1, true, true)
        @test length(Topologies.vertices(grid_topology)) == 1
        V = collect(Topologies.vertices(grid_topology))
        @test V[1] isa Topologies.Vertex
        @test length(V[1]) == 4
        @test collect(V[1]) == [(1, 1), (1, 2), (1, 3), (1, 4)]
    end

    @testset "1×1 element quad mesh with 1 periodic boundary" begin
        _, _, grid_topology = rectangular_grid(1, 1, true, false)
        @test length(Topologies.vertices(grid_topology)) == 2
        V = collect(Topologies.vertices(grid_topology))
        @test V[1] isa Topologies.Vertex
        @test length(V[1]) == 2
        @test collect(V[1]) == [(1, 1), (1, 2)]
        @test length(V[2]) == 2
        @test collect(V[2]) == [(1, 3), (1, 4)]
    end

    @testset "1×1 element quad mesh with non-periodic boundaries" begin
        _, _, grid_topology = rectangular_grid(1, 1, false, false)
        @test length(Topologies.vertices(grid_topology)) == 4
        V = collect(Topologies.vertices(grid_topology))
        @test V[1] isa Topologies.Vertex
        @test length(V[1]) == 1
        @test collect(V[1]) == [(1, 1)]
        @test length(V[2]) == 1
        @test collect(V[2]) == [(1, 2)]
        @test length(V[3]) == 1
        @test collect(V[3]) == [(1, 3)]
        @test length(V[4]) == 1
        @test collect(V[4]) == [(1, 4)]
    end

    @testset "2×3 element quad mesh with non-periodic boundaries" begin
        _, _, grid_topology = rectangular_grid(2, 3, false, false)
        @test length(Topologies.vertices(grid_topology)) == 3 * 4
        V = collect(Topologies.vertices(grid_topology))
        @test length(V[1]) == 1
        @test collect(V[1]) == [(1, 1)]
    end

    @testset "2×3 element quad mesh with periodic boundaries" begin
        _, _, grid_topology = rectangular_grid(2, 3, true, true)
        @test length(Topologies.vertices(grid_topology)) == 2 * 3
        V = collect(Topologies.vertices(grid_topology))
        @test length(V) == 6
        @test collect(V[1]) == [(1, 1), (2, 2), (5, 3), (6, 4)]
        @test collect(V[2]) == [(2, 1), (1, 2), (6, 3), (5, 4)] ||
              collect(V[2]) == [(1, 2), (2, 1), (5, 4), (6, 3)]
        @test collect(V[3]) == [(3, 1), (4, 2), (1, 3), (2, 4)] ||
              collect(V[3]) == [(1, 3), (2, 4), (3, 1), (4, 2)]
        @test collect(V[4]) == [(4, 1), (3, 2), (2, 3), (1, 4)] ||
              collect(V[4]) == [(1, 4), (2, 3), (3, 2), (4, 1)]
        @test collect(V[5]) == [(5, 1), (6, 2), (3, 3), (4, 4)] ||
              collect(V[5]) == [(3, 3), (4, 4), (5, 1), (6, 2)]
        @test collect(V[6]) == [(6, 1), (5, 2), (4, 3), (3, 4)] ||
              collect(V[6]) == [(3, 4), (4, 3), (5, 2), (6, 1)]
    end
end

@testset "simple rectangular grid coordinates" begin
    @testset "1×1 element quad mesh with all periodic boundries" begin
        _, _, grid_topology = rectangular_grid(1, 1, true, true)
        c1, c2, c3, c4 = Topologies.vertex_coordinates(grid_topology, 1)
        @test c1 == Geometry.XYPoint(0.0, 0.0)
        @test c2 == Geometry.XYPoint(1.0, 0.0)
        @test c3 == Geometry.XYPoint(0.0, 1.0)
        @test c4 == Geometry.XYPoint(1.0, 1.0)

        _, _, grid_topology = rectangular_grid(
            1,
            1,
            false,
            false;
            x1min = -1.0,
            x1max = 1.0,
            x2min = -1.0,
            x2max = 1.0,
        )
        c1, c2, c3, c4 = Topologies.vertex_coordinates(grid_topology, 1)
        @test c1 == Geometry.XYPoint(-1.0, -1.0)
        @test c2 == Geometry.XYPoint(1.0, -1.0)
        @test c3 == Geometry.XYPoint(-1.0, 1.0)
        @test c4 == Geometry.XYPoint(1.0, 1.0)
    end

    @testset "2×4 element quad mesh with non-periodic boundaries" begin
        _, _, grid_topology = rectangular_grid(2, 4, false, false)
        c1, c2, c3, c4 = Topologies.vertex_coordinates(grid_topology, 1)
        @test c1 == Geometry.XYPoint(0.0, 0.0)
        @test c2 == Geometry.XYPoint(0.5, 0.0)
        @test c3 == Geometry.XYPoint(0.0, 0.25)
        @test c4 == Geometry.XYPoint(0.5, 0.25)

        c1, c2, c3, c4 = Topologies.vertex_coordinates(grid_topology, 8)
        @test c1 == Geometry.XYPoint(0.5, 0.75)
        @test c2 == Geometry.XYPoint(1.0, 0.75)
        @test c3 == Geometry.XYPoint(0.5, 1.0)
        @test c4 == Geometry.XYPoint(1.0, 1.0)
    end

    @testset "check coordinate type accuracy" begin
        _, _, grid_topology = rectangular_grid(
            3,
            1,
            false,
            false,
            x1min = big(0.0),
            x1max = big(1.0),
            x2min = big(0.0),
            x2max = big(1.0),
        )
        c1, c2, c3, c4 = Topologies.vertex_coordinates(grid_topology, 1)
        @test eltype(c2) == BigFloat
        @test getfield(c2, 1) ≈ big(1.0) / big(3.0) rtol = eps(BigFloat)
        @test getfield(c2, 2) == 0.0
    end
end