text:
using Revise
using LarSurf
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm, SparseArrays
using Pandas
using Seaborn
using Dates
using Logging

using TimerOutputs

const to = TimerOutput()
# Compilation
# block_size=[4,4,4]
#
# segmentation = LarSurf.generate_segmentation567(1)
# Flin1, larmodel1 = LarSurf.grid_get_surf_Fvec_larmodel(segmentation)
# Flin2, larmodel2 = LarSurf.grid_get_surf_Fvec_larmodel_old(segmentation)
# Flin_block2 = LarSurf.__grid_get_surface_Fchar_per_block_old_implementation(segmentation, block_size)
# V1, top1 = larmodel1
# V2, top2 = larmodel2
# (VV1, EV1, FV1, CV1) = top1
# (VV2, EV2, FV2, CV2) = top2

segmentation = LarSurf.generate_almost_cube(100)
println("segmentation size: ", size(segmentation))
block_size=[16,16,16]

function one_test_iteration(to)
    print("=")
    @timeit to "surf" LarSurf.get_surface_grid(segmentation)
    @timeit to "surf old" LarSurf.get_surface_grid_old(segmentation)
    print(".")
    @timeit to "surf block" LarSurf.get_surface_grid_per_block(segmentation, block_size)
    @timeit to "surf block full" LarSurf.get_surface_grid_per_block_full(segmentation, block_size)
    print(".")
    @timeit to "surf block FVr" LarSurf.get_surface_grid_per_block_FVreduced(segmentation, block_size)
    @timeit to "surf block Vr FVr" LarSurf.get_surface_grid_per_block_Vreduced_FVreduced(segmentation, block_size)
    print(".")
    @timeit to "surf block Vr FVr fixed" LarSurf.get_surface_grid_per_block_Vreduced_FVreduced_fixed_block_size(segmentation, block_size)
    @timeit to "surf block Vr FVr old" LarSurf.get_surface_grid_per_block_Vreduced_FVreduced_old(segmentation, block_size)
    print(".")
    @timeit to "surf block VR FVr parallel" LarSurf.get_surface_grid_per_block_Vreduced_FVreduced_parallel(segmentation, block_size)
    print(".")
    return to
end

println("Warming...")
one_test_iteration(to)
println("Evaluation")

reset_timer!(to)

@timeit to "size20" begin
    for i=1:2
        one_test_iteration(to)
    end
end

display(to)
#


text:
include("computstate.jl")
include("graphstate.jl")
include("chpstate.jl")

ext: jl

max_stars_repo_name: Roger-luo/openQC

max_issues_repo_name: Roger-luo/openQC

max_forks_repo_name: Roger-luo/openQC


text:
using Test
using Bigsimr
using Distributions
import LinearAlgebra: PosDefException
import Bigsimr: ValidCorrelationError

@testset "Random Vector Utilities" begin

    @testset "Fast Standard Random Normal" begin
        # Must work for each bits floating point type
        for T in (Float64, Float32, Float16)
            @test_nowarn Bigsimr._randn(T, 4, 4)
            @test eltype(Bigsimr._randn(T, 4, 4)) === T
        end

        @test_nowarn Bigsimr._randn(Float64, 4.0, 4.0)
        @test_nowarn Bigsimr._randn(Float64, 4.0, 4)
        @test_nowarn Bigsimr._randn(Float64, 4,   4.0)

        @test_nowarn Bigsimr._randn(4.0, 4.0)
        @test_nowarn Bigsimr._randn(4.0, 4)
        @test_nowarn Bigsimr._randn(4  , 4.0)

        @test_throws InexactError Bigsimr._randn(4.5, 5.5)
        @test_throws InexactError Bigsimr._randn(4.5, 5)
        @test_throws InexactError Bigsimr._randn(4  , 5.5)
    end

    @testset "Random Multivariate Normal Generation" begin
        r_negdef = [
            1.00 0.82 0.56 0.44
            0.82 1.00 0.28 0.85
            0.56 0.28 1.00 0.22
            0.44 0.85 0.22 1.00
        ]

        # Must fail for negative semidefinite matrices
        @test_throws PosDefException Bigsimr._rmvn(10, r_negdef)

        # Must work for each bits floating point type
        for T in (Float64, Float32, Float16)
            r = cor_randPD(T, 4)
            @test_nowarn Bigsimr._rmvn(10, r)
            @test eltype(Bigsimr._rmvn(10, r)) === T
        end
    end

    @testset "User Random Multivariate Normal Generation" begin
        # Must work for each bits floating point type
        for T in (Float64, Float32, Float16)
            r = cor_randPD(T, 4)
            @test_nowarn rmvn(10, r)
            @test eltype(rmvn(10, r)) === T
        end

        r = cor_randPD(4)
        @test_nowarn rmvn(10.0, r)
        @test_throws InexactError rmvn(10.5, r)
    end

    @testset "Normal to Marginal" begin
        # Must work for scalars, vectors, and matrices of type Float64
        D = Binomial(10, 0.5)

        x = randn(Float64)
        y = randn(Float64, 2)
        z = randn(Float64, 2, 2)

        ω = randn(Float64, 2, 2, 2)
        γ = randn(Float32)

        @test_nowarn Bigsimr.normal_to_margin(D, x)
        @test_nowarn Bigsimr.normal_to_margin(D, y)
        @test_nowarn Bigsimr.normal_to_margin(D, z)

        @test_nowarn Bigsimr.normal_to_margin(D, ω)
        @test_nowarn Bigsimr.normal_to_margin(D, γ)

        # Standard normal to standard normal should be invariant
        z = rand(Normal(0, 1), 100000)
        @test z ≈ Bigsimr.normal_to_margin(Normal(0, 1), z)

        # Estimated parameters must be close to true parameters
        d1 = Binomial(20, 0.2)
        d2 = Poisson(3)
        d3 = Normal(12, π)

        x1 = Bigsimr.normal_to_margin(d1, z)
        x2 = Bigsimr.normal_to_margin(d2, z)
        x3 = Bigsimr.normal_to_margin(d3, z)

        f1 = fit_mle(Binomial, 20, x1)
        f2 = fit_mle(Poisson, x2)
        f3 = fit_mle(Normal, x3)

        @test all(isapprox.(params(d1), params(f1), rtol=0.01))
        @test all(isapprox.(params(d2), params(f2), rtol=0.01))
        @test all(isapprox.(params(d3), params(f3), rtol=0.01))
    end

end

@testset "Random Vector Generation" begin

    @testset "rvec" begin
        # Must throw an error if a margin is not a univariate distribution
        r = cor_randPD(2)
        m = [Binomial(10, 0.2), MvNormal(zeros(2), r)]
        @test_throws MethodError rvec(2, r, m)

        # Must throw an error if r is not a valid correlation matrix
        m = [Binomial(10, 0.3), Gamma(10, 3)]
        r = Float64[1.0 2.33333; 0.333333 1.0] # Not positive definite
        c = Float64[2 4; 4 100]                # Is covariance, not correlation
        @test_throws ValidCorrelationError rvec(3, r, m)
        @test_throws ValidCorrelationError rvec(3, c, m)

        # Must throw an arror if the dimensions of r do not match the number of margins
        m = [Binomial(10, 0.3), Gamma(10, 3)]
        r = cor_randPD(3)
        @test_throws DimensionMismatch rvec(4, r, m)

        # Distributions.jl only returns random Float64 types, but `rvec` should
        # still accept correlation matrices of different floating point types
        for T in (Float64, Float32, Float16)
            r = cor_randPD(T, 2)
            @test_nowarn rvec(10, r, m)
            @test_nowarn rvec(10.0, r, m)
            @test_throws InexactError rvec(10.5, r, m)
        end
    end

end


text:
module t

using BioSequences

struct Unsafe end

const N_AA = length(alphabet(AminoAcid))

struct CodonSet <: AbstractSet{RNACodon}
    x::UInt64

    CodonSet(x::UInt64, ::Unsafe) = new(x)
end
CodonSet() = CodonSet(UInt64(0), Unsafe())
CodonSet(itr) = foldl(push, itr, init=CodonSet())

function Base.iterate(x::CodonSet, s::UInt64=x.x)
    iszero(s) ? nothing : (reinterpret(RNACodon, trailing_zeros(s)), s & (s-1))
end

function push(s::CodonSet, x::RNACodon)
    CodonSet(s.x | (UInt64(1) << (reinterpret(UInt64, x) & 63)), Unsafe())
end

Base.length(x::CodonSet) = count_ones(x.x)
Base.in(c::RNACodon, s::CodonSet) = isodd(s.x >>> (reinterpret(UInt64, c) & 63))
delete(s::CodonSet, x::RNACodon) = setdiff(s, CodonSet((x,)))
Base.issubset(a::CodonSet, b::CodonSet) = isempty(setdiff(a, b))
Base.filter(f, s::CodonSet) = CodonSet(Iterators.filter(f, s))
Base.setdiff(a::CodonSet, b::Vararg{CodonSet}) = CodonSet(a.x & ~(union(b...).x), Unsafe())

for (name, f) in [(:union, |), (:intersect, &), (:symdiff, ⊻)]
    @eval function Base.$(name)(a::CodonSet, b::Vararg{CodonSet})
        CodonSet(mapreduce(i -> i.x, $f, b, init=a.x), Unsafe())
    end
end

struct ReverseGeneticCode <: AbstractDict{AminoAcid, CodonSet}
    name::String
    sets::NTuple{N_AA-1, CodonSet}
end

function ReverseGeneticCode(x::BioSequences.GeneticCode)
    ind(aa::AminoAcid) = reinterpret(UInt8, aa) + 1

    sets = fill(CodonSet(), N_AA-1)
    for i in Int64(0):Int64(63)
        aa = x.tbl[i + 1]
        sets[ind(aa)] = push(sets[ind(aa)], reinterpret(RNACodon, i))
    end

    # Ambiguous amino acids
    for (n, (a, b)) in [(AA_B, (AA_D, AA_N)), (AA_J, (AA_I, AA_L)), (AA_Z, (AA_E, AA_Q))]
        sets[ind(n)] = sets[ind(a)] ∪ sets[ind(b)]
    end

    # AA_X codes for all amino acids
    sets[ind(AA_X)] = CodonSet(typemax(UInt64), Unsafe())

    # Pyrrolysine and selenocysteine are never part of the "forward" genetic
    # code, but can be unambiguously resolved in the reverse genetic code.
    sets[ind(AA_U)] = CodonSet((mer"UGA"r,))
    sets[ind(AA_O)] = CodonSet((mer"UAG"r,))

    ReverseGeneticCode(x.name, Tuple(sets))
end

for symbol in names(BioSequences, all=true)
    Base.isdeprecated(BioSequences, symbol) && continue
    isdefined(BioSequences, symbol) || continue
    thing = getproperty(BioSequences, symbol)
    thing isa BioSequences.GeneticCode || continue
    @eval const $(Symbol(:rev_, symbol)) = ReverseGeneticCode($thing)
end

function Base.getindex(s::ReverseGeneticCode, a::AminoAcid)
    if reinterpret(UInt8, a) > (N_AA - 2) # cannot translate gap
        error("Cannot reverse translate element: ", a)
    end
    @inbounds s.sets[reinterpret(UInt8, a) + 1]
end

function reverse_translate!(
    v::Vector{CodonSet},
    seq::AminoAcidSeq,
    code=rev_standard_genetic_code
)
    resize!(v, length(seq))
    @inbounds for i in eachindex(v)
        v[i] = code[seq[i]]
    end
    v
end

function reverse_translate(seq::AminoAcidSeq, code=rev_standard_genetic_code)
    reverse_translate!(Vector{CodonSet}(undef, length(seq)), seq, code)
end

end # module

max_stars_repo_name: jakobnissen/play

max_issues_repo_name: jakobnissen/play

max_forks_repo_name: jakobnissen/play


text:
using Base.Test
import Ogre
import WaterData

module test_integrator_resources
import Ogre
import WaterData

# Simple ODEs
# dy/dt = 6  -->  y = 6t (y0=0)
ode_1(t, y::Real) = 6
ode_1(t, y::Vector) = [6]
ode_1_sol(t) = 6.*t
# dy/dt = 2t -->  t = t^2 + y0 (y0=0)
ode_2(t, y) = 2.*t
ode_2_sol(t) = t.^2
# dy/dt = y  -->  y = y0 e^t (y0=1)
ode_3(t, y) = y
ode_3_sol(t) = exp.(t)
# dy1/dt = -y2, dy2/dt = y1  --> oscillating solution in y and w
ode_4(t, y) = [-y[2], y[1]]
ode_4_sol(t) = hcat(cos.(t)-2*sin.(t), 2*cos.(t) + sin.(t))

simple_odes = [ode_1, ode_2, ode_3, ode_4]
simple_ode_solutions = [ode_1_sol, ode_2_sol, ode_3_sol, ode_4_sol]

# Simple no-temperature EOSes
fe = let
    vinet = Ogre.Vinet(8.30e3, 156.2e9, 6.08, [1e3, 1e6]...)
    tfd = Ogre.TFD(26, 55.845)
    Ogre.PressurePiecewiseEOS([vinet, tfd], [0, 2.09e13, Inf])
end
mgsio3 = let
    bme = Ogre.BME3(4.10e3, 247e9, 3.97, [1e3, 1e6]...)
    tfd = Ogre.TFD([12, 14, 8], [24.305, 28.0855, 15.9994], [1., 1., 3.])
    Ogre.PressurePiecewiseEOS([bme, tfd], [0, 1.35e13, Inf])
end
h2o = let
    bme = Ogre.BME3(1460., 23.7e9, 4.15, [1e3, 1e6]...)
    tfd = Ogre.TFD([1, 8], [1.00794, 15.9994], [2., 1.])
    Ogre.PressurePiecewiseEOS([bme, tfd], [0, 7686e9, Inf])
end

# Simple with-temperature EOS
type JustAPressureEOS <: Ogre.EOS end
Ogre.@addEOSCall JustAPressureEOS
(::JustAPressureEOS)(P) = 1000 + P
Ogre.istempdependent(::JustAPressureEOS) = false

type NotReallyATemperatureEOS <: Ogre.EOS end
Ogre.@addEOSCall NotReallyATemperatureEOS
(::NotReallyATemperatureEOS)(P, T) = 1000 + P
Ogre.istempdependent(::NotReallyATemperatureEOS) = true

type ReallyATemperatureEOS <: Ogre.EOS end
Ogre.@addEOSCall ReallyATemperatureEOS
(::ReallyATemperatureEOS)(P, T) = 1000 + 1/T
Ogre.istempdependent(::ReallyATemperatureEOS) = true

# Planet structures
M_earth = Ogre.M_earth
R_earth = Ogre.R_earth
atmospheric_pressure = 1e5
surface_temperature = 300
solution_grid = linspace(M_earth, 0, 100)
r_bracket = [0, 2]*M_earth
dual_layer_eos = let
    eoses = [mgsio3, fe]
    transitions = [0, 2/3, 1] * M_earth
    Ogre.MassPiecewiseEOS(eoses, transitions)
end
tri_layer_eos = let
    eoses = [h2o, mgsio3, fe]
    transitions = [0, 1/6, 2/3, 1] * M_earth
    Ogre.MassPiecewiseEOS(eoses, transitions)
end
tri_layer_system = let
    trial_surface = Ogre.BoundaryValues(M_earth, R_earth, atmospheric_pressure)
    system = Ogre.PlanetSystem(M_earth, tri_layer_eos, trial_surface,
                               solution_grid)
    Ogre.converge!(system)
    system
end
tri_layer_radius = Ogre.currentradiusguess(tri_layer_system)

# Full temperature structure
withtemp = let
    heatcap = WaterData.ConstantHeatCapacity(4200)
    eos = ReallyATemperatureEOS()
    surface = Ogre.BoundaryValues(M_earth, R_earth, atmospheric_pressure,
                                  surface_temperature)
    system = Ogre.PlanetSystem(M_earth, eos, heatcap, surface, solution_grid,
                               r_bracket)
    Ogre.converge!(system)
    system
end
end

@testset "Integrator tests" begin
    mass = Ogre.mass
    radius = Ogre.radius
    pressure = Ogre.pressure
    temperature = Ogre.temperature
    centre = Ogre.centre
    surface = Ogre.surface
    res = test_integrator_resources

    @testset "Standard RK4 integrator test cases" begin
        t = collect(Float64, 0:.1:10)

        # Four different integrator tests
        odes = res.simple_odes
        soln_funcs = res.simple_ode_solutions
        solns = map(f -> f(t), soln_funcs)

        @testset "The integrator initialises and steps properly" begin
            integrator = Ogre.RK4
            s1 = integrator(odes[1], 0., t)
            s2 = integrator(odes[1], [0.], t)
            s3 = integrator(odes[4], [1., 2.], t)
            @test Ogre.y0(s1) == 0.0
            @test start(s1) == 1
            @test next(s1, 1)[2] == 2
            @test Ogre.y0(s2) == [0.0]
            @test start(s2) == 1
            @test next(s2, 1)[2] == 2
            @test Ogre.y0(s3) == [1., 2.]
            @test start(s3) == 1
            @test next(s3, 1)[2] == 2
        end

        @testset "The dense solutions are as expected" begin
            f = Ogre.ode4_dense
            @test isapprox(f(odes[1], 0., t), solns[1])
            @test isapprox(f(odes[2], 0., t), solns[2])
            @test isapprox(f(odes[3], 1., t), solns[3], rtol=5e-3)
            @test isapprox(f(odes[4], [1., 2.], t), solns[4], rtol=5e-3)
        end
    end

    @testset "Planetary integrator test cases" begin
        @testset "Solving for a radius (high-level funcs)" begin
            @testset "We can match Madhu and Sara's results to within 1%" begin
                test_masses = [0.5, 1.0, 5.0]
                fe = res.fe
                mgsio3 = res.mgsio3
                target_fe_radii = [0.628896, 0.76727, 1.17721]
                target_mgsio3_radii = [0.842621, 1.04239, 1.65182]

                fe_radii = Ogre.R(test_masses, fe, in_earth_units=true)
                mgsio3_radii = Ogre.R(test_masses, mgsio3, in_earth_units=true)

                @test isapprox(fe_radii, target_fe_radii, rtol=0.01)
                @test isapprox(mgsio3_radii, target_mgsio3_radii, rtol=0.01)
            end

            @testset "Dual- and tri-layer solutions work" begin
                # accept a pretty large error since we're really just looking
                # to make sure that it produces the right type of solution
                @test isapprox(Ogre.R(res.M_earth, res.dual_layer_eos),
                               res.R_earth, rtol=0.05)
                @test isapprox(Ogre.R(res.M_earth, res.tri_layer_eos),
                               res.R_earth, rtol=0.05)
            end
        end

        @testset "Solving for internal structure (tri-layer)" begin
            system = res.tri_layer_system
            soln = Ogre.solve(system)

            @testset "Solution outputs are the correct type" begin
                @test isa(soln, Ogre.PlanetStructure{Ogre.NoTemp})
            end

            @testset "Solution boundaries match the boundary conditions" begin
                @test mass(centre(soln)) == 0
                @test mass(surface(soln)) == res.M_earth
                @test pressure(centre(soln)) > pressure(surface(soln))
                @test pressure(surface(soln)) == res.atmospheric_pressure
                @test radius(centre(soln)) < 1000
                @test radius(centre(soln)) > 0
                @test radius(surface(soln)) == res.tri_layer_radius
            end
        end

        @testset "Solving a single-layer solution two different ways" begin
            eos_notemp = res.JustAPressureEOS()
            eos_withtemp = res.NotReallyATemperatureEOS()
            heatcap = WaterData.ConstantHeatCapacity(1000.)
            local radius = Ogre.R(res.M_earth, eos_notemp)
            tempdep_radius = Ogre.R(res.M_earth, eos_withtemp, heatcap)
            @test tempdep_radius == radius
        end

        @testset "Temperature dependence" begin
            system = res.withtemp
            soln = Ogre.solve(system)

            @test mass(centre(soln)) == 0
            @test mass(surface(soln)) == res.M_earth
            @test pressure(centre(soln)) > pressure(surface(soln))
            @test pressure(surface(soln)) == res.atmospheric_pressure
            @test temperature(centre(soln)) > temperature(surface(soln))
            @test temperature(surface(soln)) == 300
            @test radius(centre(soln)) < 1000
            @test isa(soln, Ogre.PlanetStructure{Ogre.WithTemp})
        end
    end
end


text:
################################################################################
# NodeMap data structure.
################################################################################
# Keeps track of where nodes in the taskgraph are mapped to the architecture
# as well as any applicable metadata.
const PPLC = Union{Path{Link},Path{Port},Path{Component}}

"""
Record of how [`Taskgraph`](@ref)s and [`TaskgraphEdge`](@ref)s in a
[`Taskgraph`](@ref) map to a [`TopLevel`](@ref).
"""
mutable struct Mapping
    """
    `Dict{String, Path{Component}}` - Takes a node name and returns the path
    to the [`Component`](@ref) where that node is mapped.
    """
    nodes::Dict{String,Path{Component}}

    """
    `Vector{SparseDiGraph{Union{Path{Link},Path{Port},Path{Component}}}}` -
    Takes a integer index for an edge in the parent [`Taskgraph`](@ref) and
    returns a graph whose node types `Path`s to architectural compoennts.o

    Edge connectivity in the graph describes how the [`TaskgraphEdge`](@ref) is
    routed through the [`TopLevel`](@ref)
    """
    edges::Vector{SparseDiGraph{PPLC}}
end

getpath(m::Mapping, nodename::String) = m.nodes[nodename]
getpath(m::Mapping, i::Integer) = m.edges[i]
Base.setindex!(m::Mapping, v, i::Integer) = setindex!(m.edges, v, i)

"""
Top level data structure. Summary of parameters:

* `T` - The [`RuleSet`](@ref) used to control placement and routing.
* `D` - The number of dimensions in the architecture (will usually be 2 or 3).
"""
mutable struct Map{D,T<:RuleSet}
    "[`RuleSet`](@ref) for assigning `taskgraph` to `toplevel`."
    ruleset::T

    "[`TopLevel{A,D}`](@ref) - The `TopLevel` to be used for the mapping."
    toplevel::TopLevel{D}

    "The [`Taskgraph`](@ref) to map to the `toplevel`."
    taskgraph::Taskgraph
    options::Dict{Symbol,Any}

    "How `taskgraph` is mapped to `toplevel`."
    mapping::Mapping
    metadata::Dict{String,Any}
end

function Map(
    ruleset::T,
    toplevel::TopLevel{D},
    taskgraph::Taskgraph;
    options = Dict{Symbol,Any}(),
    metadata = Dict{String,Any}(),
) where {T<:RuleSet,D}

    # Create a new Mapping data type for the new map
    # Get the node names
    nodes = Dict(n => Path{Component}() for n in nodenames(taskgraph))
    edges = [SparseDiGraph{PPLC}() for i in 1:num_edges(taskgraph)]
    mapping = Mapping(nodes, edges)
    return Map(ruleset, toplevel, taskgraph, options, mapping, metadata)
end

rules(map::Map) = map.ruleset

function Base.show(io::IO, map::Map{T,D}) where {T,D}
    print(io, "Map{$T,$D} with $(num_nodes(map.taskgraph)) nodes")
    return nothing
end

################################################################################
# Methods for interacting with the Map.
################################################################################
getpath(m::Map, nodename::String) = getpath(m.mapping, nodename)
getpath(m::Map, i::Integer) = getpath(m.mapping, i)
getaddress(map::Map, nodename::String) = getaddress(map.toplevel, getpath(map, nodename))

function isused(m::Map{D}, addr::CartesianIndex{D}) where {D}
    for path in values(m.mapping.nodes)
        getaddress(m.toplevel, path) == addr && return true
    end
    return false
end

function gettask(m::Map{D}, addr::CartesianIndex{D}) where {D}
    for (taskname, path) in m.mapping.nodes
        if getaddress(m.toplevel, path) == addr
            return getnode(m.taskgraph, taskname)
        end
    end
    return nothing
end

################################################################################
# Saving/loading
################################################################################

function makejls(filepath)
    dir, file = splitdir(filepath)

    filename = "$(first(split(file, "."))).jls"
    finalpath = joinpath(dir, filename)
    return finalpath
end

function save(m::Map, filepath)
    f = open(makejls(filepath), "w")
    serialize(f, m.mapping)
    return close(f)
end

function load(m::Map, filepath)
    f = open(makejls(filepath), "r")
    mapping = deserialize(f)
    close(f)
    # bind deserialized result
    return m.mapping = mapping
end


text:
@testset "test of substitution model plot, 2 binary traits" begin

m3 = TwoBinaryTraitSubstitutionModel([2.0,1.2,1.1,2.2,1.0,3.1,2.0,1.1],
     ["carnivory", "noncarnivory", "wet", "dry"]);
@test_logs plot(m3)

end


text:
# Functions focusing on finding relative orderings between datasets and reordering those datasets accordingly
function get_orders(order1::Array{String,1}, order2::Array{String,1})
    common_ids = intersect(order1, order2)
    new_order1 = get_id_positions(common_ids, order1)
    new_order2 = get_id_positions(common_ids, order2)
    return new_order1, new_order2
end
#
function get_id_positions(ids::Array{String,1}, old_order::Array{String,1})
    new_order = zeros(Int64, size(ids,1))
    for i in 1:size(ids,1)
        for j in 1:size(old_order,1)
            if ids[i] == old_order[j]
                new_order[i] = j
                break
            end
        end
    end
    return new_order
end
#
function get_phenotype_intersect!(phenotypes::Phenotypes, new_order::Array{Int64,1})
    phenotypes.trait = phenotypes.trait[new_order,:]
    phenotypes.IID = phenotypes.IID[new_order]
end
# Orders snps and also converts to a BinMarker object from BinaryMarkers
function get_genotype_intersect_bin(markers::BinaryMarkers, new_order::Array{Int64,1})
    return BinMarkers(markers.markerInfo, markers.IID[new_order], markers.binarygenotypes[new_order, :])
end

function get_weight_intersect!(weights::Weights, new_order::Array{Int64,1})
    weights.rinverse = weights.rinverse[new_order]
    weights.IID = weights.IID[new_order]
end

max_stars_repo_name: egmreynolds/MANA

max_issues_repo_name: egmreynolds/MANA

max_forks_repo_name: egmreynolds/MANA


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Julia 1.3.1
#     language: julia
#     name: julia-1.3
# ---

# # Argo Profile Data : _a North Pacific Example_
#
# Here we plot Argo temperature anomalies over the years in a North Pacific box, where we see possible hints of `NP blob 2.0` & its origin.
#
# _Algorithm:_
#
# - 1. get Argo data set from `GDAC` (ftp)
# - 2. apply QC & interpolate to standard levels (`MITprof`)
# - 3. collocate with climatology & uncertainty field (`MITprof`)
# - 4. subtract climatology & bin average over region (this notebook)

# + {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# _Data Source & Format:_
#
# - Roemmich, et al, 2019: On the Future of Argo: A Global, Full-Depth, Multi-Disciplinary Array. Frontiers in Marine Science, 6. https://doi.org/10.3389/fmars.2019.00439
# - Forget, G., J.-M. Campin, P. Heimbach, C. N. Hill, R. M. Ponte, and C. Wunsch, 2015: ECCO version 4: an integrated framework for non-linear inverse modeling and global ocean state estimation. Geoscientific Model Development, 8, 3071-3104, <https://doi.org/10.5194/gmd-8-3071-2015>
#

# + {"slideshow": {"slide_type": "skip"}, "cell_type": "markdown"}
# ### Set up tools
#

# + {"cell_style": "center", "slideshow": {"slide_type": "skip"}}
using ArgoData, Plots, MAT, Dates, Statistics
argo_T = matread("argo_T.mat")
argo_S = matread("argo_S.mat")

# + {"slideshow": {"slide_type": "skip"}, "cell_type": "markdown"}
# ### Depth-Time Anomaly Function

# + {"slideshow": {"slide_type": "skip"}}
years=collect(2006:2019); ny=length(years); nd=Int(argo_T["nr"])
tim=[y+m/12 for m in 1:12, y in years][:]
dep=-argo_T["prof_depth"]

LL=50.
xrng=(2005.,2020.)
δt=Millisecond(86400. *1000. *30.) # 30 <-> +/-30 days in miliseconds
δl=2.5 # 5 <-> +/-5 degree ranges
argo_time = julian2datetime.( datetime2julian(DateTime(0)) .+ argo_T["prof_date"] )

function TimeDepthMedian(x,λ=40.)
    anom=Array{Float64,2}(undef,ny*12,nd)
    for y=1:ny, m=1:12, d=1:nd
        tt=findall( (abs.(argo_time.-DateTime(years[y],m,15)).<δt) .&
                    (abs.(argo_T["prof_lat"].-λ).<δl) )
        tt=[tt[i][1] for i in 1:length(tt)]
        xx=(x)[tt,d]
        xx=xx[findall((!isnan).(xx))]
        ~isempty(xx) ? xm=median(xx) : xm=NaN
        anom[m+(y-1)*12,d]=xm
    end

    return anom,λ
end

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### Depth-Time Sections

# + {"cell_style": "center", "slideshow": {"slide_type": "subslide"}}
anom,λ=TimeDepthMedian(argo_T["monclim"],LL); v="Tclim"
contourf(tim,vec(dep)[end:-1:1],transpose(anom[:,end:-1:1]),title="median( $v ) in $λ +/- $δl °N",xlim=xrng,ylim=(-300.,0.))

# + {"cell_style": "center", "slideshow": {"slide_type": "subslide"}}
anom,λ=TimeDepthMedian(argo_T["prof"],LL); v="Tobs"
contourf(tim,vec(dep)[end:-1:1],transpose(anom[:,end:-1:1]),title="median( $v ) in $λ +/- $δl °N",xlim=xrng,ylim=(-300.,0.))

# + {"cell_style": "center", "slideshow": {"slide_type": "subslide"}}
anom,λ=TimeDepthMedian(argo_S["monclim"],LL); v="Sclim"
contourf(tim,vec(dep)[end:-1:1],transpose(anom[:,end:-1:1]),title="median( $v ) in $λ +/- $δl °N",xlim=xrng,ylim=(-300.,0.))

# + {"cell_style": "center", "slideshow": {"slide_type": "subslide"}}
anom,λ=TimeDepthMedian(argo_S["prof"],LL); v="Sobs"
contourf(tim,vec(dep)[end:-1:1],transpose(anom[:,end:-1:1]),title="median( $v ) in $λ +/- $δl °N",xlim=xrng,ylim=(-300.,0.))

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### Depth-Time Anomaly Plots
#

# + {"cell_style": "center", "slideshow": {"slide_type": "subslide"}}
anom,λ=TimeDepthMedian(argo_T["prof"]-argo_T["monclim"],LL); v="Tobs-Tclim"
contourf(tim,vec(dep)[end:-1:1],transpose(anom[:,end:-1:1]),title="median( $v ) in $λ +/- $δl °N",xlim=xrng,ylim=(-300.,0.),clim=(-1.5,1.5))

# + {"cell_style": "center", "slideshow": {"slide_type": "subslide"}}
anom,λ=TimeDepthMedian(argo_S["prof"]-argo_S["monclim"],LL); v="Sobs-Sclim"
contourf(tim,vec(dep)[end:-1:1],transpose(anom[:,end:-1:1]),title="median( $v ) in $λ +/- $δl °N",xlim=xrng,ylim=(-300.,0.),clim=(-0.3,0.3))

max_stars_repo_name: gaelforget/MarineEcosystemNotebooks

max_forks_repo_name: gaelforget/MarineEcosystemNotebooks


text:
#=
This examples can only be run
when a Nextion display is connected
to USB TTL converter

Correct test.tft file should be previously
send to display using uSD card
=#

using Nextion
using Test

#portname = "/dev/ttyUSB0"

@testset "Examples" begin
    #include("../examples/hello_serialports.jl")  # this is broken using LibSerialPort.jl instead
    include("../examples/hello_libserial.jl")
    include("../examples/hello_nextion.jl")
    include("../examples/raw_cmd_valid.jl")
    include("../examples/raw_cmd_valid.jl")
    include("../examples/page.jl")
    include("../examples/text.jl")
    include("../examples/draw.jl")
    include("../examples/scroll.jl")
    include("../examples/number.jl")
    include("../examples/pbar.jl")
    include("../examples/picture.jl")
    include("../examples/gauge.jl")
    include("../examples/slider.jl")
#    include("../examples/waveform.jl")
    #waveform(10)  # to avoid looping infinitely
    include("../examples/checkbox.jl")
    include("../examples/radio.jl")
    #include("../examples/button.jl")  # to avoid looping infinitely
end


text:
import Base.append!
append!(gp, x::VecF64, y::Float64) = append!(gp, reshape(x, :, 1), [y])
function append!(gp::GPE{X,Y,M,K,P,D}, x::MatF64, y::VecF64) where {X,Y,M,K,P <: ElasticPDMat, D}
    size(x, 2) == length(y) || error("$(size(x, 2)) observations, but $(length(y)) targets.")
    newcov = [cov(gp.kernel, gp.x, x); cov(gp.kernel, x, x) + (exp(2*gp.logNoise) + 1e-5)*I]
    append!(gp.data, gp.kernel, gp.x, x)
    append!(gp.x, x)
    append!(gp.cK, newcov)
    gp.nobs += length(y)
    append!(gp.y, y)
    update_target!(gp, kern = false, noise = false)
end

wrap_cK(cK::ElasticPDMat, Σbuffer, chol) = cK
mat(cK::ElasticPDMat) = view(cK.mat)
cholfactors(cK::ElasticPDMat) = view(cK.chol).factors

function ElasticGPE(dim; mean::Mean = MeanZero(), kernel = SE(0.0, 0.0),
                    logNoise::Float64 = -2.0, kwargs...)
    x = ElasticArray(Array{Float64}(undef, dim, 0))
    y = ElasticArray(Array{Float64}(undef, 0))
    ElasticGPE(x, y, mean, kernel, logNoise; kwargs...)
end
function ElasticGPE(x::MatF64, y::VecF64, mean::Mean, kernel::Kernel,
                    logNoise::Float64 = -2.0;
                    capacity = 10^3, stepsize = 10^3)
    data = ElasticKernelData(kernel, x, capacity = capacity, stepsize = stepsize)
    N = length(y)
    gp = GPE(ElasticArray(x), ElasticArray(y), mean, kernel, data,
             ElasticPDMat(Σ_default(x, kernel, data, logNoise),
                          capacity = capacity, stepsize = stepsize),
             logNoise)
    initialise_target!(gp)
end
export ElasticGPE

function prepareappend!(kd, Xnew)
    dim, nobs_new = size(Xnew)
    nobs = kd.dims[1]
    if nobs + nobs_new > kd.capacity[1]
        kd.capacity = kd.capacity .+ kd.stepsize
        ElasticPDMats.resize!(kd)
    end
    kd, dim, nobs, nobs_new
end

function ElasticKernelData(k::Isotropic, X::MatF64; capacity = 10^3, stepsize = 10^3)
    kerneldata = IsotropicData(AllElasticArray(2; capacity = (capacity, capacity), stepsize = (stepsize, stepsize)))
    nobs = size(X, 2)
    distance!(view(kerneldata.R, 1:nobs, 1:nobs), k, X)
    setdimension!(kerneldata.R, nobs, 1:2)
    kerneldata
end
function append!(kerneldata::IsotropicData{<:AllElasticArray}, k::Isotropic, X::MatF64, Xnew::MatF64)
    kd, dim, nobs, nobs_new = prepareappend!(kerneldata.R, Xnew)
    distance!(view(kd, 1:nobs, nobs + 1:nobs + nobs_new), k, X, Xnew)
    copyto!(view(kd, nobs + 1:nobs + nobs_new, 1:nobs), transpose(view(kd, 1:nobs, nobs + 1:nobs + nobs_new)))
    distance!(view(kd, nobs + 1:nobs + nobs_new, nobs + 1:nobs + nobs_new), k, Xnew)
    setdimension!(kd, nobs + nobs_new, 1:2)
    kerneldata
end
append!(kerneldata, k, X, Xnew::VecF64) = append!(kerneldata, k, X, reshape(Xnew, :, 1))

function ElasticKernelData(k::StationaryARD, X::MatF64; capacity = 10^3, stepsize = 10^3)
    dim, nobs = size(X)
    dist_stack = AllElasticArray(3; capacity = (capacity, capacity, size(X, 1)),
                                    stepsize = (stepsize, stepsize, 0))
    for d in 1:dim
        grad_ls = view(dist_stack, 1:nobs, 1:nobs, d)
        distance!(grad_ls, SqEuclidean(), view(X, d:d, :))
    end
    setdimension!(dist_stack, nobs, 1:2)
    setdimension!(dist_stack, dim, 3)
    StationaryARDData(dist_stack)
end
function append!(kerneldata::StationaryARDData, kernel::StationaryARD, X::MatF64, Xnew::MatF64)
    kd, dim, nobs, nobs_new = prepareappend!(kerneldata.dist_stack, Xnew)
    for d in 1:dim
        grad_ls = view(kd, 1:nobs, nobs + 1:nobs + nobs_new, d)
        distance!(grad_ls, SqEuclidean(), view(X, d:d, :), view(Xnew, d:d, :))
        copyto!(view(kd, nobs + 1:nobs + nobs_new, 1:nobs, d),
                transpose(grad_ls))
        distance!(view(kd, nobs + 1:nobs + nobs_new,
                       nobs + 1:nobs + nobs_new, d), SqEuclidean(), view(Xnew, d:d, :))
    end
    setdimension!(kd, nobs + nobs_new, 1:2)
    kerneldata
end

function ElasticKernelData(k::LinArd, X::MatF64; capacity = 10^3, stepsize = 10^3)
    dim, nobs = size(X)
    XtX_d = AllElasticArray(3; capacity = (capacity, capacity, size(X, 1)),
                                    stepsize = (stepsize, stepsize, 0))
    @inbounds @simd for d in 1:dim
        for i in 1:nobs
            for j in 1:i
                XtX_d[i, j, d] = XtX_d[j, i, d] = X[d, i] * X[d, j]
            end
        end
    end
    setdimension!(XtX_d, nobs, 1:2)
    setdimension!(XtX_d, dim, 3)
    LinArdData(XtX_d)
end
function append!(kerneldata::LinArdData, kernel::LinArd, X::MatF64, Xnew::MatF64)
    kd, dim, nobs, nobs_new = prepareappend!(kerneldata.XtX_d, Xnew)
    @inbounds @simd for d in 1:dim
        for i in 1:nobs
            for j in 1:nobs_new
                kd[i, nobs + j, d] = kd[nobs + j, i, d] = X[d, i] * Xnew[d, j]
            end
        end
        for i in 1:nobs_new
            for j in 1:i
                kd[nobs + i, nobs + j, d] = kd[nobs + j, nobs + i, d] = Xnew[d, i] * Xnew[d, j]
            end
        end
    end
    setdimension!(kd, nobs + nobs_new, 1:2)
    kerneldata
end

function ElasticKernelData(k::LinIso, X::MatF64; capacity = 10^3, stepsize = 10^3)
    dim, nobs = size(X)
    XtX = AllElasticArray(2, capacity = (capacity, capacity), stepsize = (stepsize, stepsize))
    @inbounds @simd for d in 1:dim
        for i in 1:nobs
            for j in 1:i
                XtX[i, j] = XtX[j, i] += X[d, i] * X[d, j]
            end
        end
    end
    setdimension!(XtX, nobs, 1:2)
    LinIsoData(XtX)
end
function append!(kerneldata::LinIsoData, kernel::LinIso, X::MatF64, Xnew::MatF64)
    kd, dim, nobs, nobs_new = prepareappend!(kerneldata.XtX, Xnew)
    @inbounds @simd for d in 1:dim
        for i in 1:nobs
            for j in 1:nobs_new
                kd[i, j + nobs] = kd[j + nobs, i] += X[d, i] * Xnew[d, j]
            end
        end
        for i in 1:nobs_new
            for j in 1:i
                kd[i + nobs, j + nobs] = kd[j + nobs, i + nobs] += Xnew[d, i] * Xnew[d, j]
            end
        end
    end
    setdimension!(kd, nobs + nobs_new, 1:2)
    kerneldata
end


text:
using LARS

function get_t_stat(y, w, x_i)
    design_matrix = hcat(ones(size(w, 1)), w[:, 2:end], x_i[size(w, 2):end])
    coefficients = inv(design_matrix'*design_matrix)*(design_matrix'*y)  # OLS estimate
    residuals = y - design_matrix * coefficients
    # calcualte variance-covariance matrix according to White(1980)
    coefficient_covariance = inv(design_matrix'design_matrix)*(design_matrix'*diagm(residuals.^2)*design_matrix)*inv(design_matrix'design_matrix)
    t_stat = coefficients./sqrt(diag(coefficient_covariance))
    last(t_stat)
end

function targeted_predictors(y_index::Int64, x::Matrix{Float64}, number_of_lags::Int64, thresholding::String="hard"; number_of_steps_in_lars=30, significance_level=0.05)
    # Bai and Ng (2008) -> regress y on w and each x_i and keep only the predictors x which are significant
    #predictors = DataFrame(hcat(y, w, x))  # unfortunately using formulae and the GLM package is a bit clumsy here
    w = lagged_matrix(x[:, y_index], [0:number_of_lags])
    y = w[:, y_index]
    if thresholding == "hard"  # TODO: hard thresholding does not work properly, only admits one variable plus the white covariance has negative diagonal elements
        critical_t_stat = quantile(TDist(size(w, 1) - (size(w, 2) + 1)), 1-significance_level/2)  # number of parameters is constant + lag number + x_i
        t_stats = Float64[get_t_stat(y, w, x[:, i]) for i in [1:size(x, 2)]]  # t_stat for y regressed on itself will be huge!
        return abs(t_stats) .> critical_t_stat
    end
    if thresholding == "soft"  # lasso coefficients which are non 0 are used
        design_matrix = hcat(x[size(w, 2):end, [1:size(x,2)] .!= y_index], w[:, 2:end])
        #res=lars(design_matrix, y, intercept=true, standardize=true, use_gram=true, maxiter=typemax(Int), lambda_min=0.0, verbose=true)
        res = lars(design_matrix, y)  # adds an intercept by default
        in_the_model = falses(size(x, 2))
        betas_added = Int[step.added for step in res.steps]
        betas_added = betas_added[betas_added.<size(x, 2)]  # only set the values for beta coefficients to true (not for coefficients of w)
        betas_added = betas_added[betas_added.>0][1:number_of_steps_in_lars]  # in some steps no variable is added
        betas_added[betas_added.>y_index] = betas_added[betas_added.>y_index] + 1  # all indices after y_index need to be shifted by 1 to account for y
        in_the_model[betas_added] = true
        in_the_model[y_index] = true
        return in_the_model
    end
end


text:
push!(LOAD_PATH,"../../src/")
import UNSflow

alphadef = UNSflow.ConstDef(5. *pi/180)

hdef = UNSflow.ConstDef(0.)

udef = UNSflow.ConstDef(1.)

full_kinem = UNSflow.KinemDef(alphadef, hdef, udef)

pvt = 0.25

geometry = "FlatPlate"

lespc = [10.15;]

#The value of leading-edge radius (rho) must be specified according to geometry
surf = UNSflow.TwoDSurf(geometry, pvt, full_kinem, lespc, rho=0.016)

curfield = UNSflow.TwoDFlowField()

dtstar = UNSflow.find_tstep(alphadef)

t_tot = 3.

nsteps =Int(round(t_tot/dtstar))+1

startflag = 0

writeflag = 1

writeInterval = t_tot/10.

#delvort = delSpalart(500, 12, 1e-5)
delvort = UNSflow.delNone()

mat, surf, curfield = UNSflow.ldvmLin(surf, curfield, nsteps, dtstar,startflag, writeflag, writeInterval, delvort)

UNSflow.makeVortPlots2D()

UNSflow.makeForcePlots2D()

UNSflow.cleanWrite()

#Edge velocity at end of simulation
q_u, q_l = UNSflow.calc_edgeVel(surf, [curfield.u[1], curfield.w[1]])

using PyPlot
plot(surf.x, q_u)
plot(surf.x, q_l)


text:
using Plots
using Random
pyplot()
include("N_plots.jl")
model, obj, circles = Dynamics.dubinscar_escape
Random.seed!(2);
N=101
opts = SolverOptions()
opts.verbose = false
opts.cost_tolerance = 1e-6
opts.cost_tolerance_intermediate = 0.01
opts.constraint_tolerance = 1e-4
opts.resolve_feasible = true
opts.outer_loop_update_type = :default
opts.penalty_scaling = 50
opts.penalty_initial = 10
opts.R_infeasible = 1
opts.square_root = true
ipopt_options = Dict("tol"=>opts.cost_tolerance,"constr_viol_tol"=>opts.constraint_tolerance)

# Constrained
solver = Solver(model,obj,N=N,opts=opts)
n,m,N = get_sizes(solver)

U0 = ones(m,N-1) + rand(m,N-1)/3
X_guess = [2.5 2.5 0.;4. 5. .785;5. 6.25 0.;7.5 6.25 -.261;9 5. -1.57;7.5 2.5 0.]
X0 = TrajectoryOptimization.interp_rows(N,obj.tf,Array(X_guess'))
# plot(X0[1,:],X0[2,:])

res, stats = solve(solver,U0)
res_inf, stats_inf = solve(solver,X0,U0)
@btime solve($solver,$X0,$U0)
stats_inf["iterations"]
stats_inf["runtime"]
evals(solver,:f) / stats_inf["iterations"]

res_i, stats_i = solve_dircol(solver,X0,U0, options=ipopt_options)
@btime solve_dircol($solver,$X0,$U0, options=$ipopt_options)
stats_i["iterations"]
evals(solver,:f)/stats_i["iterations"]
stats_i["runtime"]

plt = plot(title="",aspect_ratio=:equal,size=(500,300),xlim=[-0.5,10.5],ylim=[-0.5,6.6],grid=:off,bg=:white)
plot_obstacles(circles,:grey)
scatter!(X_guess[:,1],X_guess[:,2],marker=:circle,markersize=7,markerstrokecolor=:gold,color=:gold,label="")
scatter!((X_guess[end,1],X_guess[end,2]),marker=:circle,markersize=7,markerstrokecolor=:gold,color=:gold,label="Initial Guess")
plot_trajectory!(X0,style=:dot,color=:black,width=2,label="Interpolation")
plot_trajectory!(res_i.X,width=3,color=:blue,label="DIRCOL")
plot_trajectory!(to_array(res.X),width=2,color=:darkorange2,label="ALTRO",aspect_ratio=:equal,legend=:topleft)
plot_trajectory!(to_array(res_inf.X),width=2,color=:darkorange2,style=:dash,label="ALTRO (inf)",aspect_ratio=:equal)
scatter!([obj.x0[1]],[obj.x0[2]],label=:start,color=:red,markerstrokecolor=:red,markersize=12,markershape=:rtriangle)
scatter!([obj.xf[1]],[obj.xf[2]],label=:goal,color=:green,markerstrokecolor=:green,markersize=12,markershape=:rtriangle,axis=:off)

savefig(joinpath(IMAGE_DIR,"escape_traj.svg"))

# Newton solve comparison
Random.seed!(2);
opts = SolverOptions()
opts.verbose = false
opts.cost_tolerance = 1e-6
opts.cost_tolerance_intermediate = 0.01
opts.constraint_tolerance = 1e-8
opts.resolve_feasible = true
opts.outer_loop_update_type = :default
opts.penalty_scaling = 50
opts.penalty_initial = 10
opts.R_infeasible = 1
opts.square_root = true
ipopt_options = Dict("tol"=>opts.cost_tolerance,"constr_viol_tol"=>opts.constraint_tolerance)

solver = Solver(model,obj,N=N,opts=opts)
n,m,N = get_sizes(solver)

U0 = ones(m,N-1) + rand(m,N-1)/3
X_guess = [2.5 2.5 0.;4. 5. .785;5. 6.25 0.;7.5 6.25 -.261;9 5. -1.57;7.5 2.5 0.]
X0 = TrajectoryOptimization.interp_rows(N,obj.tf,Array(X_guess'))

res_inf, stats_inf = solve(solver,X0,U0)
stats_inf["iterations"]
stats_inf["runtime"]
evals(solver,:f) / stats_inf["iterations"]

res_i, stats_i = solve_dircol(solver,X0,U0, options=ipopt_options)
stats_i["iterations"]
evals(solver,:f)/stats_i["iterations"]
stats_i["runtime"]

# Newton results
solver.opts.resolve_feasible = false
solver.opts.feasible_projection = false
idx_newton = 13
solver.opts.iterations = idx_newton # cut out at this constraint violation
@btime solve(solver,$X0,$U0)
res_n, stats_n = solve(solver,X0,U0)
max_violation(res_n)
stats_n["c_max"][end]

length(stats_n["c_max"])

using TrajectoryOptimization: gen_usrfun_newton, NewtonVars, gen_newton_functions, newton_projection
t_start = time_ns()
V_ = newton_projection(solver,res_n,eps=1e-12,verbose=false)
res_ = ConstrainedVectorResults(solver,V_.Z.X,V_.Z.U)
backwardpass!(res_,solver)
rollout!(res_,solver,0.0)
# max_violation(res_)
t_newton = float(time_ns()-t_start)/1e9

t_ilqr = 0.617336
t_dircol = 9.08
t_inf = stats_inf["runtime"]
time_newton = collect(range(0,t_ilqr,length=stats_n["iterations"]))
push!(time_newton,t_ilqr + t_newton)
time_dircol = range(0,t_dircol,length=stats_i["iterations"])
time_ilqr = range(0,length=stats_inf["iterations"],step=time_newton[2]-time_newton[1])
c_max = [stats_n["c_max"]; max_violation(res_)]
stats_inf["c_max"][31]

markersize = 3
width = 1.5
fontsize=8
gr()
pyplot()
p = plot(time_dircol,stats_i["c_max"][2:end],yscale=:log10,label="DIRCOL",color=:blue,width=width,
    markershape=:circle,markerstrokecolor=:blue,markersize=markersize)
plot!(time_ilqr,stats_inf["c_max"],yscale=:log10,label="ALTRO",width=width,color=:darkorange2,
    markershape=:circle,markerstrokecolor=:darkorange2,markersize=markersize)
plot!(time_newton,c_max,yscale=:log10,ylim=[1e-9,1e-1],xlim=[0,10],label="ALTRO*",legend=:topright,xlabel="Time (s)",ylabel="Max Constraint Violation",color=:green,width=width,
    markershape=:circle,markerstrokecolor=:green,markersize=markersize,dpi=400,size=(500,250))
plot_vertical_lines!(p,[t_ilqr],label="projected Newton")
plot_vertical_lines!(p,[time_ilqr[stats_inf["max_mu_iteration (infeasible)"]]],label="max penalty",color=:red,style=:dot,legend=:bottom,
    tickfontsize=fontsize,legendfontsize=fontsize,guidefontsize=fontsize)

savefig(p,joinpath(IMAGE_DIR,"escape_newton.png"))

cost(solver,res_)
stats_i["cost"][end-3]
stats_inf["cost"][end]

# solver_truth = Solver(model, obj, N=601)
#
# Ns = [101,145,201,251,281]
# group = "escape"
# run_step_size_comparison(model, obj, U0, group, Ns, opts=opts, integrations=[:rk3,:ipopt],benchmark=false, infeasible=true, X0=X0, dt_truth=solver_truth.dt)
# plot_stat("runtime",group,legend=:topleft,["rk3","ipopt"],title="Escape",color=[:blue :darkorange2],size=(500,250),ylim=(0,50))
# savefig(joinpath(IMAGE_DIR,"escape_runtime.eps"))
# plot_stat("iterations",group,legend=:topleft,["rk3","ipopt"],title="Escape",size=(500,250))
# savefig(joinpath(IMAGE_DIR,"escape_iters.eps"))
# plot_stat("error",group,yscale=:log10,legend=:right,["rk3","ipopt"],title="Escape")
#
# Ns, data = load_data("runtime",["rk3","ipopt"],"escape")
# Ns, err = load_data("std",["rk3","ipopt"],"escape")
# plot(Ns,data[1],yerr=err[1],label="ALTRO",markershape=:circle,markersize=6,color=:darkorange2,markerstrokecolor=:darkorange2)
# plot!(Ns,data[2],yerr=err[2],label="Ipopt",markershape=:circle,markersize=6,color=:blue,markerstrokecolor=:blue)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Singular_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Singular")
JLLWrappers.@generate_main_file("Singular", UUID("43d676ae-4934-50ba-8046-7a96366d613b"))
end  # module Singular_jll


text:
mutable struct Moon
    pos
    vel
end

function Moon(pos)
    Moon(pos, zeros(Int64, 3))
end

mutable struct System
    moons
end

function apply_gravity_pair!(moon1, moon2)
    velocity_change_1 = sign.(moon2.pos - moon1.pos)
    moon1.vel += velocity_change_1
    moon2.vel -= velocity_change_1
end

function apply_gravity!(moons)
    for i=1:length(moons)
        for j=i+1:length(moons)
            apply_gravity_pair!(moons[i], moons[j])
        end
    end
end

function apply_velocity_single!(moon)
    moon.pos += moon.vel
end

function apply_velocity!(moons)
    for moon in moons
        apply_velocity_single!(moon)
    end
end

function step!(moons)
    apply_gravity!(moons)
    apply_velocity!(moons)
end

function run!(system, steps)
    for i=1:steps
        step!(system.moons)
    end
end

function total_energy(system::System)
    sum(total_energy(moon) for moon in system.moons)
end

function potential_energy(moon::Moon)
    sum(abs.(moon.pos))
end

function kinetic_energy(moon::Moon)
    sum(abs.(moon.vel))
end

function total_energy(moon::Moon)
    potential_energy(moon) * kinetic_energy(moon)
end

test_moons = [Moon([-1,0,2]), Moon([2,-10,-7]), Moon([4,-8,8]), Moon([3,5,1])]

input_system = System([
    Moon([15, -2, -6]),
    Moon([-5, -4, -11]),
    Moon([0, -6, 0]),
    Moon([5, 9, 6]),
])

max_stars_repo_name: st--/advent-of-code

max_issues_repo_name: st--/advent-of-code

max_forks_repo_name: st--/advent-of-code


text:
# Random Maclaurin Feature Maps
# from: Kar, Karnick, "Random Feature Maps for Dot Product Kernels", AISTATS 2012
#
# Danny Perry (dperry@cs.utah.edu)
# Nov 2014
#
# Polynomial Kernel Approximation:  (alpha * x'*y + c)^d

# @param alpha - polynomail kernel parameter multiplier
# @param c - polynomail kernel parameter offset
# @param d - polynomail kernel parameter power
# @param dimension - dimension of input data
# @param num_functions - how many functions to generate
function GenerateFunctionsPolynomial(sigma, input_dimension, num_functions)
	# draw directions from a normal distribution with appropriate sigma corresponding to the kernel:
	w = randn(num_functions, input_dimension) * (1/sigma)
	# draw random offsets from a uniform distribution in [0,2*pi]
	b = rand(num_functions) * 2 * pi

	# construct the random function struct:
	RandomMaclaurinFunctions(w,b)
end


text:
function check_server_updating(configuration)
    datadir = joinpath(ENV["HOME"], "ft2k", "Software", "Watchman", "data")
    result = readlines(`stat $datadir/ctimelog.csv`)
    tm = match(r"Modify: ([0-9/-]+ [0-9:.]+)", result[6])[1] # [0-9:.]
    dm = DateTime(split(tm," ")[1], "yyyy-mm-dd")
    tm = split(split(split(tm, " ")[2],".")[1], ":")
    tm = Meta.parse.(tm)
    timemodified = Dates.DateTime(Dates.Year(dm).value, Dates.Month(dm).value, Dates.Day(dm).value,tm[1], tm[2], tm[3])
    println("ctimelog changed last $timemodified")
    timediff = (now() - timemodified)
    println("timediff is $timediff")
    if timediff > convert(Dates.Millisecond, Dates.Minute(40)) # 20 min
        pull_success = false
        subject = "server might be down"
        msg = "ctimelog unchanged for longer than 20 min, server might be down"
        println(msg)
        email_user(subject, msg,configuration)
        return false
    end
    return true
end

function email_user(subject, msg,configuration)
    domain = Ft2k.conf["Watchman"]["domain"]
    #user_email = configuration["user_email"]
    emailcontents =
"""EHLO $domain
AUTH LOGIN
QUtJQUk0Nlg0TFVBNjNOSjI2VFE=
QW16eHpHc2MxUnRyYm9VbVFPN1EwUUQrOGJNQWxYT3FSWG5WWnRINHFlNmk=
MAIL FROM: ft2k.sysmon@gmail.com
RCPT TO: andrew.s.mckay@gmail.com
DATA
From: ft2k.sysmon@gmail.com
To: andrew.s.mckay@gmail.com
Subject: $subject

$msg
.
QUIT"""
    f = open("mail_out_watchman.txt","w")

        write(f,emailcontents )

    close(f)
    A = `openssl s_client -crlf -quiet -starttls smtp -connect email-smtp.us-west-2.amazonaws.com:587`
    B = pipeline(A, stdin = "mail_out_watchman.txt") # /usr/sbin/ssmtp mail_out_watchman.txt
    run(B)
    #rm("mail_out_watchman.txt")
end

# pull and configure data
function pull_configure_data(configuration)
    # locating the files
    ctimefile = joinpath(ENV["HOME"],Ft2k.conf["Data"]["ctimelog"])
    cfeedfile = joinpath(ENV["HOME"],Ft2k.conf["Data"]["cfeedlog"])
    sqlitedb = joinpath(ENV["HOME"], Ft2k.conf["Data"]["dbsqlite"])

    #datafolder = joinpath(pwd(), "data")
    # loading time_db
    time_db = Array{Any,1}(undef,2)
    time_db[1] = readdlm(ctimefile, header = false,',')
    time_db[1] = time_db[1][:,1:4]
    x = [time_db[1][i,3]+Dates.datetime2unix(DateTime(time_db[1][i,4])) for i in 1:size(time_db[1])[1]]
    time_db[1] = hcat(time_db[1], x)
    time_db[1][:,1] = string.(time_db[1][:,1])
    time_db[2] = ["id" "v" "sec" "date" "unixtime"]
    #time_db = loadtable(joinpath(datafolder, "ctimelog.csv"), header_exists = false, colnames = ["id" "v" "sec" "date" "na"])
    feed_db = Array{Any,1}(undef, 2)
    feed_db[1] =readdlm(cfeedfile, header = false, ',')
    feed_db[1][:,1] = string.(feed_db[1][:,1])
    feed_db[1] = feed_db[1][:,1:6]
    feed_db[1] = [feed_db[1][i,:] for i in 1:size(feed_db[1])[1] if feed_db[1][i,6] != "2085f78900"]
    feed_db[1] = map(p -> permutedims(p[:,:], [2,1]), feed_db[1])

    feed_db[1] = vcat(feed_db[1]...)
    feed_db[1][:,6] = Dates.Date.(feed_db[1][:,6], "yyyy-mm-dd")
    y = [feed_db[1][i,2]+Dates.datetime2unix(DateTime(feed_db[1][i,6])) for i in 1:size(feed_db[1])[1]]
    feed_db[1] = hcat(feed_db[1], y)
    feed_db[2] = ["id" "sec" "out" "back" "cal" "date"]
    #feed_db = loadtable(joinpath(datafolder, "cfeedlog.csv"), header_exists = false, colnames = ["id" "sec" "out" "back" "cal" "date"])
    #db = SQLite.DB(joinpath(ENV["HOME"], "ft2k", "Software", "Server", "db.sqlite"))
    db = SQLite.DB(sqlitedb)
    sqlite_db = Array{Any,1}(undef,2)
    sqlite_db[1] = DataFrame(SQLite.Query(db, "SELECT * FROM clients;"))
    sqlite_db[1] = convert(Array, sqlite_db[1])
    sqlite_db[2] = ["id" "feedings" "cal" "diff"]
    #sqlite_db = table(SQLite.query(db, "SELECT * FROM clients;"))
    return time_db, feed_db, sqlite_db
end

function findmaxarray(array)
    da = Dict(zip(unique(array[:,1]), zeros(size(array[:,1]))))
    for i in 1:size(array)[1]
        if da[array[i,1]] < array[i,5]
            da[array[i,1]] = array[i,5]
        end
    end
    return da
end

# check whether any feeders have mysteriously fallen off the grid
function check_man_down(tdict)
    max_device_sleep_sec = Meta.parse(Ft2k.conf["Watchman"]["max_sleep_time_sec"])

    #time_db2 = pushcol(time_db, :unixtime, map(p -> p.sec + Dates.datetime2unix(DateTime(p.date)), time_db))
    #most_recent = groupreduce(max, time_db2, :id, select= :unixtime)

    most_recent = Dict()
    for (k,v) in tdict
        most_recent[k] = maximum(v)
    end
    nowtime = now(Dates.UTC)-Hour(7) # NOTE: this needs to change for daylight savings
    println("Now is $nowtime")
    nowtime = Dates.datetime2unix(nowtime) # needed to keep this on the right time zone
    feedersDown = []
    for (k,v) in most_recent
         if (nowtime -(2*max_device_sleep_sec)) < v < (nowtime- max_device_sleep_sec)
             push!(feedersDown, (k, v))
        end
    end
    # feedersDown = filter(p -> (nowtime-(2*max_device_sleep_sec)) < p.max < (nowtime-max_device_sleep_sec) , most_recent)
    return feedersDown
end

function get_feedings(sqlite_db, id)
    return sqlite_db[1][sqlite_db[1][:,1] .== id, :]
end

function get_last_putative_feeding(array)
    da = Dict(zip(array[:,1], array[:,2]))
    if now(Dates.UTC) - now() < Dates.Millisecond(1000)
        nowtime = now()-Dates.Hour(7)
    else
        nowtime = now()
    end
    now_sec = Dates.Second(nowtime).value + 60*Dates.Minute(nowtime).value + 60*60*Dates.Hour(nowtime).value - 20*60 # add on 20 min buffer

    last_feeding = Dict()
    for (k,v) in da
        feedings = Meta.parse.(split(v, ','))
        feedings = feedings[feedings .< now_sec]
        try
            last_feeding[k] = maximum(feedings)
        catch
            last_feeding[k] = maximum(Meta.parse.(split(v, ',')))
        end
    end
    return last_feeding
end

# check whether any feedings have been missed outside given window
function check_missed_feedings(fdict, sqlite_db)
    # set to 120 seconds for now, in config.toml
    feeding_threshold_sec = Meta.parse(Ft2k.conf["Watchman"]["feeding_threshold_sec"])

    last_feeding = get_last_putative_feeding(sqlite_db[1])
    #sqlite_db2 = pushcol(sqlite_db, :most_recent_scheduled, map(p -> get_last_putative_feeding(p.orders), sqlite_db))
    nowtime  = now(Dates.UTC) - Hour(7)
    todayfeedings= Dict()
    for (k,v) in fdict
        todayfeedings[k] = filter(p -> p[5] == Date(nowtime), v)
    end
    #feed_db[1][feed_db[1][:,6] .== Date(nowtime),:]
    mostrecent_feeding = Dict()
    for (k,v) in todayfeedings
        if length(v) > 0
            mostrecent_feeding[k] = maximum(v)
        end
    end

    missedfeedings = []
    for (k,v) in last_feeding # the last putative feeding
        if haskey(mostrecent_feeding, k) && ((60*10)> (v - mostrecent_feeding[k][1]) > feeding_threshold_sec)
           # println(k, " ", v, " ", mostrecent_feeding[k][1], " ", ((60*10)+ v - mostrecent_feeding[k][1]))
            push!(missedfeedings, (k, v, Date(nowtime)))
        end
    end

   return missedfeedings
end

function checkAndEmail(logged, feedersdown, missedfeedings)
    # log both
    formail = []
    for m in missedfeedings
        if !in(m, logged)
           # read, write, create, append
            mm = join(m, ",")
            f = open("logs/feedingsmissed", "a+")
            write(f, mm)
            close(f)
            push!(formail, m)
            push!(logged, m)
        end
    end
    # email remaining
    if length(formail) > 0
        mout = join([p[1] for p in formail], ",")
        msg = "Feedings were missed for $mout"
        subject = "Missed feedings"
        email_user(subject, msg,configuration)
    end

    formail = []
    for m in feedersdown
        # read, write, create, append
        f = open("logs/feedersdown", "a+")
        write(f, join(m, ","))
        close(f)
        push!(formail, m)
        push!(logged, m)
    end
    if length(formail) > 0
        mout = join([p[1] for p in formail], ",")
        msg = "Feeders have not checked in: $mout"
        subject = "Dropped feeders"
        email_user(subject, msg,configuration)
    end
    return logged
end

function dailycheck(tdict, fdict, ids, sqlite_db)
    todayis = Date(now()-Hour(7))
    unixtoday = Dates.datetime2unix(DateTime(todayis))

    # filter out non-feeders
    #fdict = filter((k,v) -> length(v) > 5, fdict)

    # filter dictionaries to only have today's values
    for (k,v) in tdict
        tdict[k] = filter(p -> p > unixtoday, v)
    end

    for (k,v) in fdict
        fdict[k] = filter(p -> p[6] > unixtoday, v)
    end

    # checked in over last day
    tdict = filter((k,v) -> length(v) > 0, tdict)
    fdict = filter((k,v) -> length(v) > 0, fdict)

    # number of check ins?

    # number of putative feedings over last day
    putative_feed = Dict(zip(sqlite_db[1][:,1] , length.(split.(sqlite_db[1][:,2],','))))

    outa = Array{Any}[]

    println(length(fdict))
    println(length(tdict))
    println(length(putative_feed))
    for (k,v) in tdict
        differencev = minimum(map(p -> p[3] - p[2], fdict[k]))
        push!(outa, [k length(v) length(fdict[k]) (putative_feed[k]) (length(fdict[k]) - putative_feed[k]) differencev])
    end

    outa = vcat(outa...)

    outa = outa[sortperm(outa[:,5]), :]

    subject = "Daily Check"
    msg= DataFrame(outa)
    #names!(msg, Dict(:x1 => :id, :x2 => :check_ins, :x3 => :feedingstoday, :x4 => :putative_feedings, :x5 => :feeding_diff, :x6 => :measured_diff))
    email_user(subject, msg,Ft2k.conf)
end


text:
#TO DO
# Doc string
function maskedscore(m::Langmodel, temp_lm::DefaultDict, word, context)
   score(m, temp_lm, lookup(m.vocab, [word])[1], lookup(m.vocab, [context])[1])
end

function logscore(m::Langmodel, temp_lm::DefaultDict, word, context)
    log2(maskedscore(m, temp_lm, word, context))
end

function entropy(m::Langmodel, lm::DefaultDict, text_ngram)
    local log_set=Float64[]
    for ngram in text_ngram
        ngram = split(ngram)
        push!(log_set, logscore(m, lm, ngram[end], join(ngram[1:end-1], " ")))
        #println(logscore(m,lm,ngram[end],ngram[1:end-1]))
    end
    return(sum(log_set)/length(log_set))
end

function perplexity(m::Langmodel, lm::DefaultDict, text_ngram)
    return(2^(entropy(m, lm, text_ngram)))
end


text:
function info_for_logger()
    return (; timestamp=Dates.format(now(), dateformat"yyyy-mm-dd HH:MM:SS"),
            worker_id=myid())
end

"""
    jsonable(obj::Any) = obj

Add methods to `BeaconK8sUtilities.jsonable` in order to specify how [`json_logger`](@ref) should
serialize types that themselves do not already have a `StructTypes.StructType` defined.

## Example

For example, if we want to log `TimerOutput`s from TimerOutputs.jl, we can define the
following.

```julia
using TimerOutputs

BeaconK8sUtilities.jsonable(to::TimerOutput) = TimerOutputs.todict(to)
```

This pirates `BeaconK8sUtilities.jsonable` but since it is only used for the logger,
it is better than pirating `StructTypes.StructType` which could be used more widely.
"""
jsonable(obj::Any) = obj

# https://github.com/JuliaLogging/LoggingExtras.jl/pull/59/
function make_maxlog_logger(logger)
    counts = Dict{Any,Int}()
    return ActiveFilteredLogger(logger) do log
        maxlog = get(log.kwargs, :maxlog, nothing)
        maxlog === nothing && return true # no limit
        c = get(counts, log.id, 0)
        if c < maxlog
            # log this message and update the count
            counts[log.id] = c + 1
            return true
        else
            return false
        end
    end
end

"""
    json_logger(level=Logging.Info; info_for_logger=BeaconK8sUtilities.info_for_logger, io=stderr)

Provides a logger which:

* Emits logs as JSON (via `LoggingFormats.JSON(; recursive=true)`)
* Adds extra information to every log message via `info_for_logger()`, which can be any function that evaluates
to a `NamedTuple`. The default provides a timestamp and the current worker id (`Distributed.myid()`).
* the respects `maxlog` convention (see <https://docs.julialang.org/en/v1/stdlib/Logging/#Logging.@logmsg>)

Set the positional argument to the minimum-enabled logging level, and `io` to the IO handle where the logs should be emitted.

See also [`BeaconK8sUtilities.jsonable`](@ref) for how to customize how objects are converted to JSON.
"""
function json_logger(level=Logging.Info; info_for_logger=info_for_logger, io=stderr)
    t = TransformerLogger(FormatLogger(LoggingFormats.JSON(; recursive=true), io)) do log
        transformed_kwargs = map(jsonable, NamedTuple(log.kwargs))
        return merge(log, (; kwargs=merge(info_for_logger(), transformed_kwargs)))
    end
    return MinLevelLogger(make_maxlog_logger(t), level)
end


text:
using Test, QuakeML
using Dates: DateTime

@testset "Construction" begin
    @testset "Convenience" begin
        @test QuakeML.RealQuantity(1) == QuakeML.RealQuantity(value=1)
        @test QuakeML.IntegerQuantity(2) == QuakeML.IntegerQuantity(value=2)
        @test convert(QuakeML.IntegerQuantity, 3) == QuakeML.IntegerQuantity(value=3)
        @test QuakeML.Origin(time=DateTime(2000), longitude=1, latitude=2,
            public_id="smi:QuakeML.jl/origin/a").time.value == DateTime(2000)
    end

    @testset "Schema version" begin
    end

    @testset "ResourceReference" begin
        # Not a proper URI
        @test_throws ArgumentError QuakeML.ResourceReference(value="bad URI")
        # String too long
        @test_throws ArgumentError QuakeML.ResourceReference(value="smi:local/" * "a"^246)
        @test QuakeML.ResourceReference("quakeml:QuakeML.jl/refA") isa QuakeML.ResourceReference
        @test QuakeML.ResourceReference("smi:QuakeML.jl/refB") isa QuakeML.ResourceReference
        @test QuakeML.ResourceReference("smi:local/a") isa QuakeML.ResourceReference
    end

    @testset "WhitespaceOrEmptyString" begin
        @test_throws ArgumentError QuakeML.WhitespaceOrEmptyString("x")
        @test QuakeML.WhitespaceOrEmptyString(" \t  \n ").value == " \t  \n "
        @test isempty(QuakeML.WhitespaceOrEmptyString("").value)
    end
end

max_stars_repo_name: anowacki/QuakeML


text:
"""
# KM3NeT Data Definitions v1.3.1-8-gae7dfb2
https://git.km3net.de/common/km3net-dataformat
"""

module Reconstruction
  const JPP_RECONSTRUCTION_TYPE = 4000
  const JMUONBEGIN = 0
  const JMUONPREFIT = 1
  const JMUONSIMPLEX = 2
  const JMUONGANDALF = 3
  const JMUONENERGY = 4
  const JMUONSTART = 5
  const JLINEFIT = 6
  const JMUONEND = 99
  const JSHOWERBEGIN = 100
  const JSHOWERPREFIT = 101
  const JSHOWERPOSITIONFIT = 102
  const JSHOWERCOMPLETEFIT = 103
  const JSHOWER_BJORKEN_Y = 104
  const JSHOWERENERGYPREFIT = 105
  const JSHOWERPOINTSIMPLEX = 106
  const JSHOWERDIRECTIONPREFIT = 107
  const JSHOWEREND = 199
  const DUSJ_RECONSTRUCTION_TYPE = 200
  const DUSJSHOWERBEGIN = 200
  const DUSJSHOWERPREFIT = 201
  const DUSJSHOWERPOSITIONFIT = 202
  const DUSJSHOWERCOMPLETEFIT = 203
  const DUSJSHOWEREND = 299
  const AANET_RECONSTRUCTION_TYPE = 101
  const AASHOWERBEGIN = 300
  const AASHOWERFITPREFIT = 302
  const AASHOWERFITPOSITIONFIT = 303
  const AASHOWERFITDIRECTIONENERGYFIT = 304
  const AASHOWEREND = 399
  const JUSERBEGIN = 1000
  const JMUONVETO = 1001
  const JMUONPATH = 1003
  const JMCEVT = 1004
  const JUSEREND = 1099
  const RECTYPE_UNKNOWN = -1
  const RECSTAGE_UNKNOWN = -1
end


text:
@testset "1638.count-substrings-that-differ-by-one-character.jl" begin
    @test count_substrings("ab", "bb") == 3
    @test count_substrings("a", "a") == 0
    @test count_substrings("aba", "baba") == 6
    @test count_substrings("abe", "bbc") == 10
end


text:
struct TimesteppingParams{T <: Real, N <: Integer, TO, C, P}
            niter0 :: N      #starting iteration number
                dt :: T      #timestep
          end_time :: T      #end time of this simulation
                t0 :: T      #start time of this simulation
        chkpt_freq :: T      #temporary checkpoint frequency
       pchkpt_freq :: T      #permanent checkpoint frequency
       chkpt_path  :: String #path to location of permanent and tempoary checkpoint output
      n_iter_total :: TO     #total number of timesteps counting from zero
      n_iter_chkpt :: C      #number of iterations per temporary checkpoint
     n_iter_pchkpt :: P      #number of iterations per permanent checkpoint
    step_thickness :: Bool   #toggle whether to step the thickness at each timestep or not (coupling control)
end

"""
TimesteppingParams(;
                    niter0 = 0,
                    dt = 1.0,
                    end_time = 1.0,
                    t0 = nothing,
                    chkpt_freq = Inf,
                    pchkpt_freq = Inf,
                    chkpt_path = './',
                    step_thickness = true)

Construct a WAVI.jl TimesteppingParams object.
TimesteppingParams stores information relating to timestepping.

Keyword arguments
=================
- 'niter0': Iteration number of the first timestep. niter0 = 0 corresponds to a new simulation, while niter0 > 0 (positive integer) corresponds to a pickup.
- 'dt': Model timestep
- 'end_time': Simulation termination time
- 't0': Starting time of the simulation
- 'chkpt_freq': Frequency of outputting temporary checkpoints
- 'pchkpt_freq': Frequecy with which permanent checkpoints are pass
- 'chkpt_path' : Path to location checkpoint output
- 'step_thickness': Toggle whether to update the ice thickness (true) or not (false) at each timestep
"""
function TimesteppingParams(;
                        niter0 = 0,
                        dt = 1.0,
                        end_time = nothing,
                        n_iter_total = nothing,
                        t0 = nothing,
                        chkpt_freq = Inf,
                        pchkpt_freq = Inf,
                        chkpt_path = "./",
                        step_thickness = true)

    #initialize t0 (really you should read start time from pickup file)
    t0 = niter0 > 0 ? niter0 * dt : 0
    t0 = map(typeof(dt), t0)

    #check compatibility of n_iter_total and end_time, and compute them
    end_time, n_iter_total = compute_iterations_and_end_time(end_time, n_iter_total,dt)

    #compute number of timesteps checkpoint number of timesteps
    chkpt_freq == Inf ? n_iter_chkpt = Inf : n_iter_chkpt  = round(Int, chkpt_freq/dt)
    pchkpt_freq == Inf ? n_iter_pchkpt = Inf : n_iter_pchkpt = round(Int, pchkpt_freq/dt)

    #check the output path ends in '/' and exists
    endswith(chkpt_path, "/") || (chkpt_path = string(chkpt_path, "/"))
    if ~isdir(chkpt_path)
        @warn string("Did not find output path ", chkpt_path, ". Any outputs will go to the working directory", pwd())
        chkpt_path = "./"
    end

    return TimesteppingParams(niter0, dt, end_time, t0, chkpt_freq, pchkpt_freq,
                            chkpt_path,n_iter_total, n_iter_chkpt, n_iter_pchkpt, step_thickness)
end

"""
    compute_iterations_and_end_time(end_time, n_iter_total, dt)

Check input end time and total iterations, and compute those not passed.
"""
function compute_iterations_and_end_time(end_time, n_iter_total, dt)
    (isa(dt, Real) && (dt > 0))  || throw(ArgumentError("timestep dt must be a positive number"))

    if (~(n_iter_total === nothing) && ~(end_time === nothing)) #if both passed, throw error if incompatible
        (end_time ≈ (n_iter_total * dt)) ||  throw(ArgumentError("You have specified both end time (end_time) and total iterations (n_iter_total), but their values are incompatible: end time mustequal n_iter_total * dt"))
    elseif ((n_iter_total === nothing) && ~(end_time === nothing)) #if only end time passed, n_iter_total is the nearest integer
        end_time == Inf ? n_iter_total = Inf : n_iter_total  = round(Int, end_time/dt)
    elseif (~(n_iter_total === nothing) && (end_time === nothing)) #if only number of iterations
        n_iter_total == Inf ? end_time = Inf : end_time  = dt*n_iter_total
    elseif (n_iter_total === nothing) && (end_time === nothing) #if neither is passed
        throw(ArgumentError("You must pass at least one of end_time or n_iter_total"))
    end

    #check both are positive numbers, now that both end_time and n_iter_total assigned
    (isa(n_iter_total, Integer) && (n_iter_total > 0)) || throw(ArgumentError("n_iter_total must be a positive integer"))
    (isa(end_time, Real) && (end_time > 0))|| throw(ArgumentError("end_time must be a positive number"))
    return end_time, n_iter_total
end


text:
module TestBasic

using GarishPrint
using Test

struct T1
    x::Int
end

struct T2
    x::Int
    y::Float64
end

struct T3
    x::Int
    y::T1
    z::T2
end

struct T4{T}
    x::T
    y::T
end

struct T5{T}
    x::T
end

struct TTTTTTTTTTTTTTTTTTTT{A,B}
    x::A
    y::B
end

pprint(T1(1))

pprint(T2(1, 2.0))

pprint(T3(1, T1(2), T2(1, 2.0)))

pprint(T4(1+2im, 2im))
pprint(T4(2 + 0im, 2im))

pprint(T5(Any[1, 2, 3]))
pprint(stdout, fill(undef))

T4(T5([1, 2, 3]), T5([1, 2, 3]))|>pprint

pprint(T5(π))

pprint(T5(rand(100)))
pprint(T5(Dict(1=>2)))
pprint(T5(Dict()))
pprint(T5(π))
pprint(T5(ℯ))
pprint(T5(fill(π)))

pprint(T5(0.10f0))
pprint(T5("0.10f0"))
pprint(T5(Dict("a"=>2.0, "b"=>2im)))
pprint(T5(Dict("a"=>T4(T5([1, 2, 3]), T5([1, 2, 3])), "b"=>2im)))

pprint(T5(Dict("a"=>(1, 2, 3), "b"=>Any)))

a = TTTTTTTTTTTTTTTTTTTT(1, 2)
b = TTTTTTTTTTTTTTTTTTTT(a, a)
c = TTTTTTTTTTTTTTTTTTTT(b, b)

pprint(c)

end


text:
module NVRTC

using ..CUDA

include("../lib/$(CUDA.major).$(CUDA.minor)/libnvrtc.jl")
include("../lib/$(CUDA.major).$(CUDA.minor)/libnvrtc_types.jl")

if is_windows()
    const libnvrtc = Libdl.find_library(["nvrtc64_80","nvrtc64_75"])
else
    const libnvrtc = Libdl.find_library(["libnvrtc"])
end
isempty(libnvrtc) && error("NVRTC library cannot be found.")

function check_nvrtcresult(status)
    status == NVRTC_SUCCESS && return nothing
    warn("NVRTC error triggered from:")
    Base.show_backtrace(STDOUT, backtrace())
    throw(unsafe_string(nvrtcGetErrorString(status)))
end

function getlog(prog::Ptr{Void})
    logsize = Csize_t[0]
    nvrtcGetProgramLogSize(prog, logsize)
    log = Array(UInt8, logsize[1])
    nvrtcGetProgramLog(prog, log)
    unsafe_string(log)
end

function compile(code::String; headers=(), include_names=())
    ref = Ptr{Void}[0]
    headers = Ptr{Void}[pointer(h) for h in headers]
    include_names = Ptr{Void}[pointer(n) for n in include_names]
    nvrtcCreateProgram(ref, code, C_NULL, length(headers), headers, include_names)

    prog = ref[1]
    options = ["--gpu-architecture=compute_30"]
    nvrtcCompileProgram(prog, length(options), options)
    ptxsize = Csize_t[0]
    nvrtcGetPTXSize(prog, ptxsize)
    ptx = Array(UInt8, ptxsize[1])
    nvrtcGetPTX(prog, ptx)
    #log = getlog(prog)
    #throw(log)
    nvrtcDestroyProgram(ref)
    String(ptx)
end

end


text:
"""
    @einsum (i,j...) -> expr
    @einsum expr

Conducts tensor computation based on [Einstein summation convention](https://en.wikipedia.org/wiki/Einstein_notation).
The arguments of the anonymous function are regard as **free indices**.
If arguments are not given, they are guessed based on the order that indices appears from left to right.

# Examples
```jldoctest einsum
julia> A = rand(Mat{3,3})
3×3 Tensor{Tuple{3, 3}, Float64, 2, 9}:
 0.590845  0.460085  0.200586
 0.766797  0.794026  0.298614
 0.566237  0.854147  0.246837

julia> B = rand(Mat{3,3})
3×3 Tensor{Tuple{3, 3}, Float64, 2, 9}:
 0.579672   0.066423  0.112486
 0.648882   0.956753  0.276021
 0.0109059  0.646691  0.651664

julia> @einsum (i,j) -> A[i,k] * B[k,j]
3×3 Tensor{Tuple{3, 3}, Float64, 2, 9}:
 0.643225  0.609151  0.32417
 0.962977  1.00373   0.500018
 0.885164  1.01445   0.460311

julia> @einsum A[i,k] * B[k,j] # same as above
3×3 Tensor{Tuple{3, 3}, Float64, 2, 9}:
 0.643225  0.609151  0.32417
 0.962977  1.00373   0.500018
 0.885164  1.01445   0.460311

julia> @einsum A[i,j] * B[i,j]
2.454690093453888
```
"""
macro einsum(expr)
    freeinds, body = anonymous_args_body(expr)
    einex = einsum_instantiate(body)
    freeinds === nothing && return einex.ex
    isempty(freeinds) && return einex.ex
    perm = find_perm(einex.freeinds => freeinds)
    :(permutedims($(einex.ex), $(ValTuple(perm...))))
end

ValTuple(x...) = Val(x)

function anonymous_args_body(func::Expr)
    if Meta.isexpr(func, :->)
        lhs = func.args[1]
        body = func.args[2]
        if Meta.isexpr(lhs, :tuple)
            freeinds = lhs.args
        elseif lhs isa Symbol
            freeinds = [lhs]
        else
            throw(ArgumentError("wrong arguments in anonymous function expression"))
        end
        if Meta.isexpr(body, :block)
            body = only([x for x in body.args if !(x isa LineNumberNode)])
        end
        freeinds, body
    else
        nothing, func
    end
end

function find_perm((inds, freeinds)::Pair)::Vector{Int}
    map(freeinds) do index
        I = findall(==(index), inds)
        @assert I !== nothing
        only(I)
    end
end

function find_freeindices(allinds::Vector)
    freeinds = []
    for index in unique(allinds)
        x = findall(==(index), allinds)
        length(x) > 2 && error("@einsum: index $index appears more than twice")
        length(x) == 1 && push!(freeinds, index)
    end
    freeinds
end

struct EinsumExpr
    ex::Any
    freeinds::Vector
    allinds::Vector
    function EinsumExpr(ex, freeinds, allinds)
        find_freeindices(allinds) # check dummy indices
        new(ex, freeinds, allinds)
    end
end

isscalarexpr(x::EinsumExpr) = isempty(x.freeinds)

function einsum_instantiate(expr)
    if Meta.isexpr(expr, :call)
        if expr.args[1] == :*
            einex = einsum_instantiate(expr.args[2])
            for ex in expr.args[3:end]
                einex = einsum_instantiate_contraction(einex, einsum_instantiate(ex))
            end
            return einex
        elseif expr.args[1] == :/
            lhs = einsum_instantiate(expr.args[2])
            rhs = einsum_instantiate(expr.args[3])
            return einsum_instantiate_division(lhs, rhs)
        elseif expr.args[1] == :+ || expr.args[1] == :-
            einex = einsum_instantiate(expr.args[2])
            for ex in expr.args[3:end]
                einex = einsum_instantiate_addition(expr.args[1], einex, einsum_instantiate(ex))
            end
            return einex
        end
    elseif Meta.isexpr(expr, :ref)
        ex = esc(expr.args[1])
        allinds = expr.args[2:end]
        freeinds = find_freeindices(allinds)
        return einsum_instantiate_tensor(EinsumExpr(ex, freeinds, allinds))
    end
    EinsumExpr(expr, [], [])
end

# ref case
function einsum_instantiate_tensor(einex::EinsumExpr)
    if isscalarexpr(einex) # handle `A[i,i]`
        ex = :($einsum_contraction(Val(()), ($(einex.ex),), ($(ValTuple(einex.allinds...)),)))
        return EinsumExpr(ex, einex.freeinds, einex.allinds)
    else
        return einex
    end
end

# division
function einsum_instantiate_division(lhs::EinsumExpr, rhs::EinsumExpr)
    @assert isscalarexpr(rhs)
    ex = Expr(:call, :/, lhs.ex, rhs.ex)
    EinsumExpr(ex, lhs.freeinds, lhs.allinds) # is it ok to ignore indices of `rhs`?
end

# summation
function einsum_instantiate_addition(op::Symbol, lhs::EinsumExpr, rhs::EinsumExpr)
    @assert Set(lhs.freeinds) == Set(rhs.freeinds)
    perm = find_perm(rhs.freeinds => lhs.freeinds)
    ex = Expr(:call, op, lhs.ex, :(permutedims($(rhs.ex), $(ValTuple(perm...)))))
    EinsumExpr(ex, lhs.freeinds, lhs.freeinds) # reset allinds
end

# contraction
function einsum_instantiate_contraction(lhs::EinsumExpr, rhs::EinsumExpr)
    if isscalarexpr(lhs)
        ex = Expr(:call, :*, lhs.ex, rhs.ex)
        return EinsumExpr(ex, rhs.freeinds, vcat(lhs.allinds, rhs.allinds))
    elseif isscalarexpr(rhs)
        ex = Expr(:call, :*, lhs.ex, rhs.ex)
        return EinsumExpr(ex, lhs.freeinds, vcat(lhs.allinds, rhs.allinds))
    else
        freeinds = find_freeindices(vcat(lhs.freeinds, rhs.freeinds))
        allinds = vcat(lhs.allinds, rhs.allinds)
        ex = :($einsum_contraction($(ValTuple(freeinds...)), ($(lhs.ex), $(rhs.ex)), ($(ValTuple(lhs.freeinds...)), $(ValTuple(rhs.freeinds...)))))
        return EinsumExpr(ex, freeinds, allinds)
    end
end

# for dummy indices
@inline function sumargs(x, ys...)
    ret = prod(x)
    @simd for i in eachindex(ys)
        @inbounds ret += prod(ys[i])
    end
    ret
end

function einsum_contraction_expr(freeinds::Tuple, tensors::Tuple{Vararg{Any, N}}, tensorinds::Tuple{Vararg{AbstractVector, N}}) where {N}
    allinds = vcat([collect(x) for x in tensorinds]...)
    dummyinds = setdiff(allinds, freeinds)
    allinds = [freeinds..., dummyinds...]

    # check dimensions
    dummyaxes = Base.OneTo{Int}[]
    for symbol in dummyinds
        dim = 0
        count = 0
        for (i, inds) in enumerate(tensorinds)
            for I in findall(==(symbol), inds)
                if dim == 0
                    dim = size(tensors[i], I)
                    push!(dummyaxes, axes(tensors[i], I))
                else
                    size(tensors[i], I) == dim || error("@einsum: dimension mismatch")
                end
                count += 1
            end
        end
        count == 2 || error("@einsum: index $symbol appears more than twice")
    end

    # tensor -> global indices (connectivities)
    whichindices = Vector{Int}[]
    for (i, inds) in enumerate(tensorinds)
        length(inds) == ndims(tensors[i]) || error("@einsum: the number of indices does not match the number of dimensions")
        whichinds = map(inds) do index
            I = findall(==(index), allinds)
            @assert I !== nothing
            only(I)
        end
        push!(whichindices, collect(whichinds))
    end

    T = promote_type(map(eltype, tensors)...)
    if isempty(freeinds)
        TT = T
        freeaxes = ()
        tupleinds = Colon()
    else
        perm = map(freeinds) do index
            only(findall(==(index), vcat(tensorinds...)))
        end
        TT = tensortype(_permutedims(otimes(map(Space, tensors)...), Val(perm))){T}
        freeaxes = axes(TT)
        tupleinds = indices(TT)
    end

    sumexps = map(CartesianIndices(freeaxes)) do finds
        xs = map(CartesianIndices(Tuple(dummyaxes))) do dinds
            ainds = [Tuple(finds)..., Tuple(dinds)...]
            exps = map(enumerate(tensors)) do (i, t)
                inds = ainds[whichindices[i]]
                getindex_expr(t, :(tensors[$i]), inds...)
            end
            Expr(:tuple, exps...)
        end
        :($sumargs($(xs...)))
    end

    TT, sumexps[tupleinds]
end

@generated function _einsum_contraction(::Val{freeinds}, tensors::Tuple, tensorinds::Tuple{Vararg{Val}}) where {freeinds}
    TT, exps = einsum_contraction_expr(freeinds, Tuple(tensors.parameters), Tuple(map(p -> collect(p.parameters[1]), tensorinds.parameters)))
    quote
        @_inline_meta
        $TT($(exps...))
    end
end

# call contraction methods if possible
@generated function einsum_contraction(_freeinds::Val{freeinds}, tensors::Tuple, _tensorinds::Tuple{Vararg{Val}}) where {freeinds}
    default = quote
        @_inline_meta
        _einsum_contraction(_freeinds, tensors, _tensorinds)
    end

    length(tensors.parameters) != 2 && return default
    lhsall, rhsall = map(p -> collect(p.parameters[1]), _tensorinds.parameters)

    lhsdummy = Int[]
    rhsdummy = Int[]
    for i in eachindex(lhsall)
        index = lhsall[i]
        I = findall(==(index), rhsall)
        isempty(I) && continue
        push!(lhsdummy, i)
        push!(rhsdummy, only(I))
    end

    issequence(x) = isempty(x) ? false : x == first(x):first(x)+length(x)-1
    if issequence(lhsdummy) && issequence(rhsdummy) &&
       last(lhsdummy) == lastindex(lhsall) && first(rhsdummy) == firstindex(rhsall)
        lhsfree = deleteat!(lhsall, lhsdummy)
        rhsfree = deleteat!(rhsall, rhsdummy)
        if freeinds == tuple(lhsfree..., rhsfree...)
            return quote
                @_inline_meta
                contraction(tensors..., $(Val(length(lhsdummy))))
            end
        end
    end

    default
end


text:
using Colors
using FixedPointNumbers
using Test
using Random
using NeuroCore.ColorChannels
using AxisIndices
using MappedArrays
using TimeAxes

using AxisIndices.PaddedViews: filltype, paddedviews, PaddedView

using NeuroCore.ColorChannels: GenericImage, GenericGrayImage
using NeuroCore.ColorChannels: NumberLike, Pixel

include("channels.jl")
include("views.jl")
include("convert_reinterpret.jl")
include("traits.jl")

using FFTW

@testset "functions" begin
    ag = rand(Gray{Float32}, 4, 5)
    ac = rand(RGB{Float32}, 4, 5)
    for (f, args) in ((fft, (ag,)), (fft, (ag, 1:2)), (plan_fft, (ag,)),
                      (rfft, (ag,)), (rfft, (ag, 1:2)), (plan_rfft, (ag,)),
                      (fft, (ac,)), (fft, (ac, 1:2)), (plan_fft, (ac,)),
                      (rfft, (ac,)), (rfft, (ac, 1:2)), (plan_rfft, (ac,)))
        ret = @test_throws ErrorException f(args...)
        @test occursin("channel_view", ret.value.msg)
        @test occursin(eltype(args[1])<:Gray ? "1:2" : "2:3", ret.value.msg)
    end
    for (a, dims) in ((ag, 1:2), (ac, 2:3))
        @test ifft(fft(channel_view(a), dims), dims) ≈ channel_view(a)
        ret = @test_throws ErrorException rfft(a)
        @test occursin("channel_view", ret.value.msg)
        @test occursin("$dims", ret.value.msg)
        @test irfft(rfft(channel_view(a), dims), 4, dims) ≈ channel_view(a)
    end

    a = [RGB(1,0,0), RGB(0,1,0), RGB(0,0,1)]
    @test a' == [RGB(1,0,0) RGB(0,1,0) RGB(0,0,1)]

    a = [RGB(1,0,0) RGB(0,1,0); RGB(0,0,1) RGB(0, 0, 0)]
    @test a' == [RGB(1,0,0) RGB(0,0,1); RGB(0,1,0) RGB(0, 0, 0)]
end


text:
@testset "Hashing" begin
    Mod.@polyvar x y z
    @test hash(x) != hash(y)
    @test hash(1x) == hash(1.0x)
    @test hash(1x+3y) == hash(1.0x+3.0y)
    @test hash(one(x)) == hash(x^0)
    @test hash(x*y) == hash(polynomial(x*y))
    @test hash(0.0x) == hash(0.0)
    @test hash(1.0x) == hash(x)
    @test hash(x-x) == hash(zero(x))
    @test hash(monovec([z^3, z*x, y])[1:2]) == hash(monovec([z^3, z*x]))
    @test hash(monovec([z^3, x, y])[1:1]) == hash(monovec([z^3, x])[1:1])
    @test hash(emptymonovec(x)) == hash([])
    @test hash(1) == hash(one(x))
    @test hash(1) == hash(constantmonomial(x * y))
    @test hash(2) != hash(constantmonomial(x * y))
    @test hash(0.0) == hash(x-x)
end


text:
using POMDPs
using ContinuousPOMDPTreeSearchExperiments
using POMDPToolbox
using BasicPOMCP
using POMCPOW
using ParticleFilters
using CPUTime
using DataFrames
using ARDESPOT
using Distributions
using DiscreteValueIteration
using VDPTag2
using QMDP
using MCTS

@everywhere using VDPTag2
@everywhere using POMDPs

pomdp = VDPTagPOMDP(mdp=VDPTagMDP(barriers=CardinalBarriers(0.2, 2.8)))

@show max_time = 1.0
@show max_depth = 10
@show RO = RandomSolver

function gen_sims(x::Vector{Float64}, n, k, seed)
    m = round(Int, x[1])
    @assert m >= 1
    c = x[2]
    @assert c >= 0.0
    k_act = x[3]
    @assert k_act >= 1.0
    inv_alpha_act = x[4]
    @assert inv_alpha_act >= 0.1
    k_obs = x[5]
    @assert k_obs >= 1.0
    inv_alpha_obs = x[6]
    @assert inv_alpha_obs >= 0.1

    sims = []

    for i in 1:n
        rng = MersenneTwister(0)

        # node_updater = ObsAdaptiveParticleFilter(deepcopy(pomdp),
        #                                    LowVarianceResampler(m),
        #                                    0.05, rng)

        # ro = ToNextML(mdp(pomdp), rng)
        ro = RandomSolver(rng)::RO
        ev = SampleRollout(solve(ro, pomdp), rng)
        solver = DPWSolver(n_iterations=typemax(Int),
                           exploration_constant=c,
                           depth=max_depth,
                           max_time=max_time,
                           k_action = k_act,
                           alpha_action = 1/inv_alpha_act,
                           k_state = k_obs,
                           alpha_state = 1/inv_alpha_obs,
                           check_repeat_state=false,
                           check_repeat_action=false,
                           estimate_value=ev,
                           next_action=RootToNextMLFirst(rng),
                           rng=rng
                          )

        belief_mdp = MeanRewardBeliefMDP(pomdp,
                                         LowVarianceResampler(m),
                                         0.05
                                        )

        planner = solve(solver, belief_mdp)

        policy = PolicyWrapper(planner) do planner, s
            local a::TagAction
            try
                a = action(planner, s)
            catch ex
                @show typeof(ex)
                warn("error in policy evaluation; using default")
                a = TagAction(false, 0.0)
            end
            return a::TagAction
        end

        filter = ObsAdaptiveParticleFilter(deepcopy(pomdp),
                                           LowVarianceResampler(10_000),
                                           0.05,
                                           MersenneTwister(i+10000*k))

        srand(policy, i+40000*k)
        sim = Sim(deepcopy(pomdp),
                  policy,
                  filter,
                  rng=MersenneTwister(i+50_000*k),
                  max_steps=50,
                  metadata=Dict(:i=>i, :k=>k)
                 )
        push!(sims, sim)
    end

    return sims
end

# start_mean = [100.0, 2.0, 10.0]
# start_cov = diagm([100.0^2, 10.0^2, 20.0^2])
start_mean =        [15.0,   85.0,   20.0,   20.0,   8.0,   60.0]
start_cov = diagm([8.0^2, 60.0^2,  10.0^2, 10.0^2, 5.0^2, 30.0^2])
d = MvNormal(start_mean, start_cov)
K = 160 # 60 # number of parameter samples
n = 20  # 100 # number of evaluation simulations
m = 40  # 15 # number of elite samples
max_iters = 100

for i in 1:max_iters
    sims = []
    params = Vector{Vector{Float64}}(K)
    print("creating $K simulation sets")
    for k in 1:K
        p = rand(d)
        p[1] = max(1.0, p[1])
        p[2] = max(1.0, p[2])
        p[3] = max(1.0, p[3])
        p[4] = max(1.0, p[4])
        p[5] = max(1.0, p[5])
        p[6] = max(1.0, p[6])
        params[k] = p
        k_sims = gen_sims(p, n, k, i)
        print(".")
        append!(sims, k_sims)
    end
    println()
    results = run_parallel(sims)
    # results = run(sims)
    combined = by(results, :k) do df
        DataFrame(mean_reward=mean(df[:reward]))
    end
    @show mean(combined[:mean_reward])
    order = sortperm(combined[:mean_reward])
    elite = params[combined[:k][order[K-m+1:end]]]
    elite_matrix = Matrix{Float64}(length(start_mean), m)
    for k in 1:m
        elite_matrix[:,k] = elite[k]
    end
    try
        d = fit(typeof(d), elite_matrix)
    catch ex
        if ex isa Base.LinAlg.PosDefException
            println("pos def exception")
            d = fit(typeof(d), elite_matrix += 0.01*randn(size(elite_matrix)))
        else
            rethrow(ex)
        end
    end
    println("iteration $i")
    @show mean(d)
    @show det(cov(d))
    @show ev = eigvals(cov(d))
    for j in 1:length(ev)
        @show eigvecs(cov(d))[:,j]
    end
end


text:
doc"""
    tanh(x::Var)

Hyperbolic tangent function.
"""
Base.tanh(x::Var) = Var(tanh(x.data), ∇tanh!, (x,))
Base.tanh(x::Array) = tanh.(x)
Base.tanh(x::CuArray) = CUDNN.tanh(x)
Base.tanh(x::Node) = Node(tanh, (x,))

function ∇tanh!(y::Var, x::Var)
    isnothing(x.grad) && return
    ∇tanh!(y.data, y.grad, x.data, x.grad)
end

function ∇tanh!(y::Array{T}, gy::Array{T}, x::Array{T}, gx::Array{T}) where T
    @inbounds for i = 1:length(gx)
        gx[i] += gy[i] * (T(1) - y[i] * y[i])
    end
end

∇tanh!(y::CuArray, gy::CuArray, x::CuArray, gx::CuArray) = CUDNN.∇tanh!(y, gy, x, gx)


text:
# NOTE: FilterHelper has its own implementation in filter.jl
"""
"""
function process_args{V}(::Type{V}, exs, index, primary)::Tuple{Set{Symbol}, Expr}
    helpers_ex = Expr(:ref, V)
    # Each source gets its own symbol set, which we use to build the
    # mappings/reverse mappings (see build_f_ex!)
    srcs_used = Set{Symbol}()
    for ex in exs
        # println(1)
        process_arg!(srcs_used, helpers_ex, V, ex, index, primary)
    end

    return srcs_used, helpers_ex
end

function process_arg! end


text:
"""
Utilities for constructing interpolant-basis polynomial sum-of-squares models.
"""
module PolyUtils

using DocStringExtensions

using LinearAlgebra
# delete later, affects qr. see https://github.com/JuliaLang/julia/pull/40623
if VERSION < v"1.7.0-DEV.1188"
    const ColumnNorm = Val{true}
end

import Combinatorics
import SpecialFunctions: gamma_inc

include("realdomains.jl")
include("realinterp.jl")
include("complex.jl")

end


text:
@testset "CIS square example with $factory" for factory in blp_factories
    include(dirname(dirname(pathof(HybridSystems))), "examples", "square.jl")
    horizontal_jump_example(Polyhedra.DefaultLibrary{Float64}(), false)
    p = getis(hs, factory)
    @test p[1].c ≈ [0, 0] atol=1e-5
    @test p[1].Q ≈ [1 0; 0 1] rtol=1e-5
end


text:
struct DBLaboratory
    database::SQLite.DB
    pkey::Int
    name::String
    contact::DBPerson
end

Base.show(io::IO, lab::DBLaboratory) = print(io, lab.name)

function Base.write(db::SQLite.DB, ::Type{DBLaboratory}, name, contactkey::Int)::Int
    stmt1 = SQLite.Stmt(db, "INSERT INTO LABORATORY ( NAME, CONTACT ) VALUES ( ?, ? );")
    r = DBInterface.execute(stmt1, ( name, contactkey))
    res = DBInterface.lastrowid(r)
    write(db, DBMember, Int(res), contactkey)
    return res
end

Base.write(db::SQLite.DB, ::Type{DBLaboratory}, name, contact::DBPerson)::Int =
    write(db, DBLaboratory, name, contact.pkey)

function Base.read(db::SQLite.DB, ::Type{DBLaboratory}, pkey::Int)::DBLaboratory
    stmt1 = SQLite.Stmt(db, "SELECT * FROM LABORATORY WHERE PKEY=?;")
    q = DBInterface.execute(stmt1, ( pkey, ))
    if SQLite.done(q)
        error("No known laboratory with key '$(pkey)'.")
    end
    r = SQLite.Row(q)
    contact = read(db, DBPerson, r[:CONTACT])
    return DBLaboratory(db, r[:PKEY], r[:NAME], contact)
end

function Base.findall(db::SQLite.DB,::Type{DBLaboratory})::Vector{DBLaboratory}
    stmt1 = SQLite.Stmt(db, "SELECT * FROM LABORATORY;")
    q = DBInterface.execute(stmt1)
    return [ DBLaboratory(r[:PKEY], r[:NAME], read(db, DBPerson, r[:CONTACT])) for r in q]
end

struct DBMember
    person::DBPerson
    laboratory::DBLaboratory
end

function Base.write(db::SQLite.DB, ::Type{DBMember}, labkey::Int, personkey::Int)::Int
    stmt1 = SQLite.Stmt(db, "INSERT INTO LABMEMBER ( LABKEY, PERSONKEY ) VALUES ( ?, ? );")
    r = DBInterface.execute(stmt1, ( labkey, personkey ))
    return DBInterface.lastrowid(r)
end

Base.write(db::SQLite.DB, ::Type{DBMember}, lab::DBLaboratory, person::DBPerson)::Int =
    write(db, DBMember, lab.pkey, person.pkey)

"""
    Base.findall(db::SQLite, ::Type{DBLaboratory}, person::DBPerson)::Vector{DBLaboratory}

Find all labs that the specified person is associated with.
"""
function Base.findall(db::SQLite.DB, ::Type{DBLaboratory}, person::DBPerson)::Vector{DBLaboratory}
    stmt1 = SQLite.Stmt(db, "SELECT LABKEY, PERSONKEY FROM LABMEMBER WHERE PERSONKEY=?;")
    q = DBInterface.execute(stmt1, ( person.pkey, ))
    return [ read(db, DBLaboratory, r[:LABKEY]) for r in q ]
end

"""
    Base.findall(db::SQLite, ::Type{DBPerson}, lab::DBLaboratory)::Vector{DBPerson}

Find all people associated with the specified lab.
"""
function Base.findall(db::SQLite.DB, ::Type{DBPerson}, lab::DBLaboratory)::Vector{DBPerson}
    stmt1 = SQLite.Stmt(db, "SELECT LABKEY, PERSONKEY FROM LABMEMBER WHERE LABKEY=?;")
    q = DBInterface.execute(stmt1, ( lab.pkey, ))
    return [ read(db, DBPerson, r[:PERSONKEY]) for r in q ]
end


text:
include("makedocs.jl")

deploydocs(
    repo = "github.com/HomodyneCT/MartaCT.jl.git",
    deploy_config = Documenter.GitHubActions()
)


text:
"""
    square(bbox::Vector{T}) where {T <: Real}

Take a bounding box and calculates the minimum square bounding box that would contain the input.

# Examples
```jldoctest
julia> bbox = [-1, 1, 2, 3.5]
4-element Array{Float64,1}:
 -1.0
  1.0
  2.0
  3.5

julia> square(bbox)
4-element Array{Float64,1}:
 -1.0
  0.75
  2.0
  3.75
```
"""
function square(bbox::Vector{T}) where {T <: Real}
    west = bbox[1]
    south = bbox[2]
    east = bbox[3]
    north = bbox[4]

    hDist = distance(float(bbox[1:2]), float([east, south]))
    vDist = distance(float(bbox[1:2]), float([west, north]))

    if hDist >= vDist
        vMidPoint = (south + north) / 2

        return [west, vMidPoint - ((east - west) / 2), east, vMidPoint + ((east - west) / 2)]
    else
        hMidPoint = (west + east) / 2

        return [hMidPoint - ((north - south) / 2), south, hMidPoint + ((north - south) / 2),  north]
    end
end


text:
#	cylinders and toroidal surface by Julia's implementation of PLaSM operators

using ViewerGL; GL = ViewerGL
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

#	1D complex in 2D generated by Lar and PLaSM operators
circle = Lar.circle()()
GL.VIEW([ GL.GLFrame2, GL.GLGrid( circle...,GL.COLORS[1],1 ) ]);

#	1D complex in 1D generated by Lar and PLaSM operators
intervals = Lar.qn(5)([2,-.5])
GL.VIEW([ GL.GLFrame2, GL.GLGrid( intervals...,GL.COLORS[1],1 ) ]);

#	2D complex in 3D generated by Lar and PLaSM operators
cylinders = Lar.larModelProduct([ circle, intervals ])
GL.VIEW([ GL.GLFrame2, GL.GLGrid( cylinders...,GL.COLORS[1],1 ) ]);

#	other visualizations
V,FV = cylinders
FW = [[f] for f in FV]
GL.VIEW( push!( GL.GLExplode( V,FW, 1.2,1.2,1.2, 1,0.2 ), GL.GLFrame2) );
GL.VIEW( push!( GL.GLExplode( V,FW, 1.,1.,1., 99,1 ), GL.GLFrame2) );

# see
function schlegel3D(t::Float64,d::Float64)
    function schlegel3D0(V)
        W = Array{Float64,2}(undef,size(V))
        for k=1:size(V,2)
            W[4,k] = V[4,k] + t # translation
            W[1,k] = V[1,k]*d/W[4,k] # projection
            W[2,k] = V[2,k]*d/W[4,k] # projection
            W[3,k] = V[3,k]*d/W[4,k] # projection
        end
        return W[1:3,:]
    end
    return schlegel3D0
end

# 2D complex in 4D generated as topological product of circles
toroidal = Lar.larModelProduct([ Lar.circle(8)(), Lar.circle(2)() ])
V,FV = toroidal

# removal of 4-th coodinate: ortho-projection 4D -> 3D
GL.VIEW([ GL.GLFrame2, GL.GLGrid( V[1:3,:], FV, GL.COLORS[1],0.2 ) ]);

# schlegel projection 4D -> 3D
W = schlegel3D(-6.,2.)(V)
GL.VIEW([ GL.GLFrame2, GL.GLGrid( W, FV, GL.COLORS[1],0.2 ) ]);

# schlegel projection 4D -> 3D rainbow
meshes = []
W,EW = W, FV
for k=1:length(EW)
	color = GL.COLORS[k%12+1] - (rand(Float64,4)*0.1)
	push!(meshes, GL.GLGrid(W,[EW[k]],color,1) )
end
GL.VIEW(meshes);

#	cylinders and toroidal surface by Julia's implementation of PLaSM operators

using ViewerGL; GL = ViewerGL
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

#	1D complex in 2D generated by Lar and PLaSM operators
circle = Lar.circle()()
GL.VIEW([ GL.GLFrame2, GL.GLGrid( circle...,GL.COLORS[1],1 ) ]);

#	1D complex in 1D generated by Lar and PLaSM operators
intervals = Lar.qn(5)([2,-.5])
GL.VIEW([ GL.GLFrame2, GL.GLGrid( intervals...,GL.COLORS[1],1 ) ]);

#	2D complex in 3D generated by Lar and PLaSM operators
cylinders = Lar.larModelProduct([ circle, intervals ])
GL.VIEW([ GL.GLFrame2, GL.GLGrid( cylinders...,GL.COLORS[1],1 ) ]);

#	other visualizations
V,FV = cylinders
FW = [[f] for f in FV]
GL.VIEW( push!( GL.GLExplode( V,FW, 1.2,1.2,1.2, 1,0.2 ), GL.GLFrame2) );
GL.VIEW( push!( GL.GLExplode( V,FW, 1.,1.,1., 99,1 ), GL.GLFrame2) );

# see
function schlegel3D(t::Float64,d::Float64)
    function schlegel3D0(V)
        W = Array{Float64,2}(undef,size(V))
        for k=1:size(V,2)
            W[4,k] = V[4,k] + t # translation
            W[1,k] = V[1,k]*d/W[4,k] # projection
            W[2,k] = V[2,k]*d/W[4,k] # projection
            W[3,k] = V[3,k]*d/W[4,k] # projection
        end
        return W[1:3,:]
    end
    return schlegel3D0
end

# 2D complex in 4D generated as topological product of circles
toroidal = Lar.larModelProduct([ Lar.circle(8)(), Lar.circle(2)() ])
V,FV = toroidal

# removal of 4-th coodinate: ortho-projection 4D -> 3D
GL.VIEW([ GL.GLFrame2, GL.GLGrid( V[1:3,:], FV, GL.COLORS[1],0.2 ) ]);

# schlegel projection 4D -> 3D
W = schlegel3D(-6.,2.)(V)
GL.VIEW([ GL.GLFrame2, GL.GLGrid( W, FV, GL.COLORS[1],0.2 ) ]);

# schlegel projection 4D -> 3D rainbow
meshes = []
W,EW = W, FV
for k=1:length(EW)
	color = GL.COLORS[k%12+1] - (rand(Float64,4)*0.1)
	push!(meshes, GL.GLGrid(W,[EW[k]],color,1) )
end
GL.VIEW(meshes);

#	cylinders and toroidal surface by Julia's implementation of PLaSM operators

using ViewerGL; GL = ViewerGL
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

#	1D complex in 2D generated by Lar and PLaSM operators
circle = Lar.circle()()
GL.VIEW([ GL.GLFrame2, GL.GLGrid( circle...,GL.COLORS[1],1 ) ]);

#	1D complex in 1D generated by Lar and PLaSM operators
intervals = Lar.qn(5)([2,-.5])
GL.VIEW([ GL.GLFrame2, GL.GLGrid( intervals...,GL.COLORS[1],1 ) ]);

#	2D complex in 3D generated by Lar and PLaSM operators
cylinders = Lar.larModelProduct([ circle, intervals ])
GL.VIEW([ GL.GLFrame2, GL.GLGrid( cylinders...,GL.COLORS[1],1 ) ]);

#	other visualizations
V,FV = cylinders
FW = [[f] for f in FV]
GL.VIEW( push!( GL.GLExplode( V,FW, 1.2,1.2,1.2, 1,0.2 ), GL.GLFrame2) );
GL.VIEW( push!( GL.GLExplode( V,FW, 1.,1.,1., 99,1 ), GL.GLFrame2) );

# see
function schlegel3D(t::Float64,d::Float64)
    function schlegel3D0(V)
        W = Array{Float64,2}(undef,size(V))
        for k=1:size(V,2)
            W[4,k] = V[4,k] + t # translation
            W[1,k] = V[1,k]*d/W[4,k] # projection
            W[2,k] = V[2,k]*d/W[4,k] # projection
            W[3,k] = V[3,k]*d/W[4,k] # projection
        end
        return W[1:3,:]
    end
    return schlegel3D0
end

# 2D complex in 4D generated as topological product of circles
toroidal = Lar.larModelProduct([ Lar.circle(8)(), Lar.circle(2)() ])
V,FV = toroidal

# removal of 4-th coodinate: ortho-projection 4D -> 3D
GL.VIEW([ GL.GLFrame2, GL.GLGrid( V[1:3,:], FV, GL.COLORS[1],0.2 ) ]);

# schlegel projection 4D -> 3D
W = schlegel3D(-6.,2.)(V)
GL.VIEW([ GL.GLFrame2, GL.GLGrid( W, FV, GL.COLORS[1],0.2 ) ]);

# schlegel projection 4D -> 3D rainbow
meshes = []
W,EW = W, FV
for k=1:length(EW)
	color = GL.COLORS[k%12+1] - (rand(Float64,4)*0.1)
	push!(meshes, GL.GLGrid(W,[EW[k]],color,1) )
end
GL.VIEW(meshes);


text:
# This file is part of Kpax3. License is MIT.

include("exceptions.jl")

include("settings.jl")

include("prior_col.jl")
include("prior_row.jl")

include("data.jl")

include("state.jl")
include("state_list.jl")
include("support.jl")