text:
"""
    replace_ref_end!(ex)

Recursively replace occurrences of the symbol :end in a "ref" expression (i.e. A[...]) `ex`
with the appropriate function calls (`lastindex` or `size`). Replacement uses
the closest enclosing ref, so

    A[B[end]]

should transform to

    A[B[lastindex(B)]]

"""
replace_ref_end!(ex) = replace_ref_end_!(ex, nothing)[1]
# replace_ref_end_!(ex,withex) returns (new ex, whether withex was used)
function replace_ref_end_!(ex, withex)
    used_withex = false
    if isa(ex,Symbol) && ex == :end
        withex === nothing && error("Invalid use of end")
        return withex, true
    elseif isa(ex,Expr)
        if ex.head == :ref
            ex.args[1], used_withex = replace_ref_end_!(ex.args[1],withex)
            S = isa(ex.args[1],Symbol) ? ex.args[1]::Symbol : gensym(:S) # temp var to cache ex.args[1] if needed
            used_S = false # whether we actually need S
            # new :ref, so redefine withex
            nargs = length(ex.args)-1
            if nargs == 0
                return ex, used_withex
            elseif nargs == 1
                # replace with lastindex(S)
                ex.args[2], used_S = replace_ref_end_!(ex.args[2],:($lastindex($S)))
            else
                n = 1
                J = lastindex(ex.args)
                for j = 2:J
                    exj, used = replace_ref_end_!(ex.args[j],:($lastindex($S,$n)))
                    used_S |= used
                    ex.args[j] = exj
                    if isa(exj,Expr) && exj.head == :...
                        # splatted object
                        exjs = exj.args[1]
                        n = :($n + length($exjs))
                    elseif isa(n, Expr)
                        # previous expression splatted
                        n = :($n + 1)
                    else
                        # an integer
                        n += 1
                    end
                end
            end
            if used_S && S !== ex.args[1]
                S0 = ex.args[1]
                ex.args[1] = S
                ex = Expr(:let, :($S = $S0), ex)
            end
        else
            # recursive search
            for i = eachindex(ex.args)
                ex.args[i], used = replace_ref_end_!(ex.args[i],withex)
                used_withex |= used
            end
        end
    end
    ex, used_withex
end

"""
    @view A[inds...]

Creates a `SubArray` from an indexing expression. This can only be applied directly to a
reference expression (e.g. `@view A[1,2:end]`), and should *not* be used as the target of
an assignment (e.g. `@view(A[1,2:end]) = ...`).  See also [`@views`](@ref)
to switch an entire block of code to use views for slicing.

```jldoctest
julia> A = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> b = @view A[:, 1]
2-element view(::Array{Int64,2}, :, 1) with eltype Int64:
 1
 3

julia> fill!(b, 0)
2-element view(::Array{Int64,2}, :, 1) with eltype Int64:
 0
 0

julia> A
2×2 Array{Int64,2}:
 0  2
 0  4
```
"""
macro view(ex)
    if Meta.isexpr(ex, :ref)
        ex = replace_ref_end!(ex)
        if Meta.isexpr(ex, :ref)
            ex = Expr(:call, view, ex.args...)
        else # ex replaced by let ...; foo[...]; end
            @assert Meta.isexpr(ex, :let) && Meta.isexpr(ex.args[2], :ref)
            ex.args[2] = Expr(:call, view, ex.args[2].args...)
        end
        Expr(:&&, true, esc(ex))
    else
        throw(ArgumentError("Invalid use of @view macro: argument must be a reference expression A[...]."))
    end
end

############################################################################
# @views macro code:

# maybeview is like getindex, but returns a view for slicing operations
# (while remaining equivalent to getindex for scalar indices and non-array types)
@propagate_inbounds maybeview(A, args...) = getindex(A, args...)
@propagate_inbounds maybeview(A::AbstractArray, args...) = view(A, args...)
@propagate_inbounds maybeview(A::AbstractArray, args::Number...) = getindex(A, args...)
@propagate_inbounds maybeview(A) = getindex(A)
@propagate_inbounds maybeview(A::AbstractArray) = getindex(A)

# _views implements the transformation for the @views macro.
# @views calls esc(_views(...)) to work around #20241,
# so any function calls we insert (to maybeview, or to
# lastindex in replace_ref_end!) must be interpolated
# as values rather than as symbols to ensure that they are called
# from Base rather than from the caller's scope.
_views(x) = x
function _views(ex::Expr)
    if ex.head in (:(=), :(.=))
        # don't use view for ref on the lhs of an assignment,
        # but still use views for the args of the ref:
        lhs = ex.args[1]
        Expr(ex.head, Meta.isexpr(lhs, :ref) ?
                      Expr(:ref, _views.(lhs.args)...) : _views(lhs),
             _views(ex.args[2]))
    elseif ex.head == :ref
        Expr(:call, maybeview, _views.(ex.args)...)
    else
        h = string(ex.head)
        # don't use view on the lhs of an op-assignment a[i...] += ...
        if last(h) == '=' && Meta.isexpr(ex.args[1], :ref)
            lhs = ex.args[1]

            # temp vars to avoid recomputing a and i,
            # which will be assigned in a let block:
            a = gensym(:a)
            i = [gensym(:i) for k = 1:length(lhs.args)-1]

            # for splatted indices like a[i, j...], we need to
            # splat the corresponding temp var.
            I = similar(i, Any)
            for k = 1:length(i)
                if Meta.isexpr(lhs.args[k+1], :...)
                    I[k] = Expr(:..., i[k])
                    lhs.args[k+1] = lhs.args[k+1].args[1] # unsplat
                else
                    I[k] = i[k]
                end
            end

            Expr(:let,
                 Expr(:block,
                      :($a = $(_views(lhs.args[1]))),
                      [:($(i[k]) = $(_views(lhs.args[k+1]))) for k=1:length(i)]...),
                 Expr(first(h) == '.' ? :(.=) : :(=), :($a[$(I...)]),
                      Expr(:call, Symbol(h[1:end-1]),
                           :($maybeview($a, $(I...))),
                           _views.(ex.args[2:end])...)))
        else
            Expr(ex.head, _views.(ex.args)...)
        end
    end
end

"""
    @views expression

Convert every array-slicing operation in the given expression
(which may be a `begin`/`end` block, loop, function, etc.)
to return a view. Scalar indices, non-array types, and
explicit `getindex` calls (as opposed to `array[...]`) are
unaffected.

!!! note
    The `@views` macro only affects `array[...]` expressions
    that appear explicitly in the given `expression`, not array slicing that
    occurs in functions called by that code.

# Examples
```jldoctest
julia> A = zeros(3, 3);

julia> @views for row in 1:3
           b = A[row, :]
           b[:] .= row
       end

julia> A
3×3 Array{Float64,2}:
 1.0  1.0  1.0
 2.0  2.0  2.0
 3.0  3.0  3.0
```
"""
macro views(x)
    esc(_views(replace_ref_end!(x)))
end

max_stars_repo_name: ninjin/julia

max_issues_repo_name: ninjin/julia

max_forks_repo_name: ninjin/julia


text:
# Parameter searching to meet the conditions
using SciMLBase
using SteadyStateDiffEq
using OrdinaryDiffEq
using ModelingToolkit
using Optim
using ThreadsX
using DataFrames

"""
Find parameters that satisfy the boolean conditions in the retrograde (RTG) signalling model using Optimization methods.
"""
function optim_params(
    Model=RtgMTK;
    datafile=joinpath(@__DIR__, "data", "boolean_table_RTG13.csv"),
    knockoutlevel=1e-4,
    proteinlevels=STRESSED,
    steadyStateSolver=DynamicSS(Rodas5()),
    lowerbound=1e-3,
    upperbound=1e3,
    hilllowerbound=2.0,
    hillupperbound=7.0,
    xinit=1.0,
    optimsolver=Optim.SAMIN(),
    optimoptions=Optim.Options(iterations=10^5, show_trace=true, show_every=1000),
    targetratio=10)

    # Boolean conditions
    conds = load_data(datafile)

    @named sys = Model(ONE_SIGNAL; proteinlevels)
    prob = SteadyStateProblem(sys, resting_u0(sys))
    params = copy(prob.p)

    # Indices to parameters in the system
    param2idx = Dict(k => i for (i, k) in enumerate(parameters(sys)))
    iΣRtg1 = param2idx[ΣRtg1]
    iΣRtg2 = param2idx[ΣRtg2]
    iΣRtg3 = param2idx[ΣRtg3]
    iΣMks = param2idx[ΣMks]
    imul_S = param2idx[mul_S]

    # Parameters to be optimized
    params_optim = [k for k in parameters(sys) if !any(isequal(k), (ΣRtg1, ΣRtg2, ΣRtg3, ΣMks, ΣBmh, mul_S))]

    # Mapping indices of the x vector in Optim to params in the system
    xi2params = [param2idx[k] for k in params_optim]
    xidxnS = findfirst(isequal(n_S), params_optim)

    # The cost is bounded by target protein concentration ratio
    scorecap = -log10(targetratio)

    # Cost function
    function cost(x)
        Σcost = 0.0
        for cond in conds
            # params = copy(prob.p)
            # Adjust params according to conditions
            params[iΣRtg1] = cond[:Rtg1] == 0 ? knockoutlevel : proteinlevels[ΣRtg1]
            params[iΣRtg2] = cond[:Rtg2] == 0 ? knockoutlevel : proteinlevels[ΣRtg2]
            params[iΣRtg3] = cond[:Rtg3] == 0 ? knockoutlevel : proteinlevels[ΣRtg3]
            params[iΣMks] = cond[:Mks] == 0 ? knockoutlevel : proteinlevels[ΣMks]
            params[imul_S] = cond[:s]

            # Assign optim vector to ODE parameters
            for i in 1:length(x)
                params[xi2params[i]] = x[i]
            end

            # Calculate cost for this steady state solution
            # Since the objective is nuclear accumulation or not
            # The cost is the logarithm of protein concentration ratios (nuclear vs cytosol)
            sol = solve(remake(prob, p=params), steadyStateSolver)

            if cond[:gfp] == "rtg3"
                score = -log10(rtg3_nucleus(sol) / rtg3_cytosol(sol)) * ifelse(cond[:Trans2Nuc] == 1, 1, -1)
            elseif cond[:gfp] == "rtg1"
                score = -log10(rtg1_nucleus(sol) / rtg1_cytosol(sol)) * ifelse(cond[:Trans2Nuc] == 1, 1, -1)
            else
                score = 0
            end
            score = max(score - scorecap, zero(score))
            Σcost += score
        end

        return Σcost / length(conds)
    end

    # Initial conditions and lower / upper bounds for Optim
    x0 = similar(xi2params, Float64)
    x0 .= xinit
    lb = similar(x0)
    lb .= lowerbound
    ub = similar(x0)
    ub .= upperbound
    ub[xidxnS] = hillupperbound
    lb[xidxnS] = hilllowerbound
    x0 .= clamp.(x0, lb, ub)

    res = Optim.optimize(cost, lb, ub, x0, optimsolver, optimoptions)
    parammap = Dict(params_optim .=> Optim.minimizer(res))

    return (res=res, parammap=parammap)
end

"""Find n parameter sets using multithreading"""
function optim_params_threads(n::Int, Model=RtgMTK; optimoptions=Optim.Options(iterations=10^5), kwargs...)
    parammaps = ThreadsX.map(1:n) do i
        res = optim_params(Model; optimoptions=optimoptions, kwargs...)
        res.parammap
    end
    return DataFrame(parammaps)
end

ext: jl


text:
@testset "utils.jl" begin
    @testset "writing and reading" begin
        S = [rand(ComplexF64, 3) for i = 1:2]
        write_solutions("tmp_sols.txt", S)
        @test read_solutions("tmp_sols.txt") == S

        p = rand(ComplexF64, 14)
        write_parameters("tmp_params.txt", p)
        @test read_parameters("tmp_params.txt") == p
    end

    Q = HC.rand_unitary_matrix(8, ComplexF64)
    @test norm(LinearAlgebra.I - Q' * Q, Inf) < 1e-14

    @test HC.fast_abs(-2.0) ≈ abs(-2.0)
    @test HC.fast_abs(2.0 + 3im) ≈ abs(2.0 + 3im)

    @test HC.nthroot(42.2, 0) == 1.0
    @test HC.nthroot(42.2, 1) == 42.2
    @test HC.nthroot(42.1, 2) == sqrt(42.1)
    @test HC.nthroot(42.1, 3) == cbrt(42.1)
    @test HC.nthroot(42.1, 4) == sqrt(sqrt(42.1))
    @test HC.nthroot(42.1, 5) == 42.1^(1 / 5)

    @testset "SegmentStepper" begin
        seg = HC.SegmentStepper(0, 1)
        @test seg.start == 0
        @test seg.target == 1
        @test seg.abs_Δ == 1.0
        @test seg.forward
        @test !isempty(sprint(show, seg))

        HC.propose_step!(seg, exp2(-60))
        @test seg.Δs == exp2(-60)
        @test seg.s′ == exp2(-60)
        @test seg.t′ == exp2(-60)
        @test seg.t == 0
        @test seg.Δt == exp2(-60)
        HC.step_success!(seg)
        @test seg.s == exp2(-60)
        @test seg.t == exp2(-60)
        @test seg.t′ == exp2(-60)
        HC.propose_step!(seg, exp2(-50))
        @test seg.t == exp2(-60)
        @test seg.s′ == exp2(-60) + exp2(-50)
        @test seg.t′ == exp2(-60) + exp2(-50)

        HC.propose_step!(seg, 4)
        @test seg.s′ == 1.0
        @test seg.t′ == 1.0
        @test !HC.is_done(seg)
        HC.step_success!(seg)
        @test HC.is_done(seg)

        seg = HC.SegmentStepper(im, 0)
        @test seg.start == im
        @test seg.target == 0
        @test seg.abs_Δ == 1.0
        @test !seg.forward

        HC.propose_step!(seg, exp2(-20))
        @test seg.Δs == exp2(-20)
        @test seg.s′ == 1 - exp2(-20)
        @test seg.t′ == im * (1 - exp2(-20))
        @test seg.t == im
        @test seg.Δt == -im * exp2(-20)
        HC.step_success!(seg)
        @test seg.s == 1 - exp2(-20)
        @test seg.t == im * (1 - exp2(-20))
        @test seg.t′ == im * (1 - exp2(-20))
        HC.propose_step!(seg, exp2(-10))
        @test seg.t == im * (1 - exp2(-20))
        @test seg.s′ == 1 - (exp2(-20) + exp2(-10))
        @test seg.t′ == im * (1 - (exp2(-20) + exp2(-10)))

        HC.propose_step!(seg, 2)
        @test seg.s′ == 0.0
        @test seg.t′ == 0.0
        @test !HC.is_done(seg)
        HC.step_success!(seg)
        @test HC.is_done(seg)
    end

    struct Foo_print_fieldnames
        a::Int
        b::Float64
    end

    @test sprint(HC.print_fieldnames, Foo_print_fieldnames(2, 2.2)) ==
          "Foo_print_fieldnames:\n • a → 2\n • b → 2.2\n"
end


text:
#!/usr/bin/julia

run(`clang kaleidoscope.cpp`)


text:
using AoC2020
using Test

@testset "AoC2020.jl" begin
    # Write your tests here.
end


text:
module InterfacesTests

using Test
using TensorPolynomialBases
using StaticArrays

P = SVector{2,Float64}
V = SVector{3,Float64}
G = SMatrix{2,3,Float64,6}
filter(e,order) = sum(e) <= order
order = 1
basis = MonomialBasis{P,V}(filter,order)

x = SVector(2.0,3.0)

cache = ScratchData(basis)
v = zeros(V,length(basis))
g = zeros(G,length(basis))
evaluate!(v,basis,x,cache)
gradient!(g,basis,x,cache)
test_polynomial_basis(basis,x,v,g)

v = evaluate(basis,x)
g = gradient(basis,x)
test_polynomial_basis(basis,x,v,g)

@test P == point_type(basis)
@test V == value_type(basis)
@test G == gradient_type(basis)

end # module InterfacesTests


text:
using Libtask
using LinearAlgebra
using BenchmarkTools

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

function benchmark_driver!(f, x...; f_displayname=string(f))
    x = (x..., nothing)

    println("benchmarking $(f_displayname)...")
    tf = Libtask.TapedFunction(f, x...);

    print("  Run Original Function:")
    @btime $f($(x)...)
    GC.gc()

    print("  Run TapedFunction:")
    @btime $tf($(x)...)
    GC.gc()

    ctf = Libtask.compile(tf)
    print("  Run TapedFunction (compiled):")
    @btime $ctf($(x)...)
    GC.gc()

    print("  Run TapedTask: ")
    x = (x[1:end-1]..., produce);
    # show the number of produce calls inside `f`
    f_task = (f, x; verbose=false) -> begin
        tt = TapedTask(f, x...);
        c = 0
        while consume(tt)!==nothing
            c+=1
        end
        verbose && print("#produce=", c, "; ");
    end
    # Note that we need to pass `f` instead of `tf` to avoid
    #  default continuation in `TapedTask` constructor, see, e.g.
    #  https://github.com/TuringLang/Libtask.jl/pull/135
    f_task(f, x; verbose=true) # print #produce calls
    @btime $f_task($f, $x)
    GC.gc()
end

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

function rosenbrock(x, callback=nothing)
    i = x[2:end]
    j = x[1:end-1]
    ret = sum((1 .- j).^2 + 100*(i - j.^2).^2)
    callback !== nothing && callback(ret)
    return ret
end

x = rand(100000)
benchmark_driver!(rosenbrock, x)

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

function ackley(x::AbstractVector, callback=nothing)
    a, b, c = 20.0, -0.2, 2.0*π
    len_recip = inv(length(x))
    sum_sqrs = zero(eltype(x))
    sum_cos = sum_sqrs
    for i in x
        sum_cos += cos(c*i)
        sum_sqrs += i^2
        callback !== nothing && callback(sum_sqrs)
    end
    return (-a * exp(b * sqrt(len_recip*sum_sqrs)) -
            exp(len_recip*sum_cos) + a + MathConstants.e)
end

x = rand(100000)
benchmark_driver!(ackley, x)

####################################################################
function generate_matrix_test(n)
    return (x, callback=nothing) -> begin
        # @assert length(x) == 2n^2 + n
        a = reshape(x[1:n^2], n, n)
        b = reshape(x[n^2 + 1:2n^2], n, n)
        ret = log.((a * b) + a - b)
        callback !== nothing && callback(ret)
        return ret
    end
end

n = 100
matrix_test = generate_matrix_test(n)
x = collect(1.0:(2n^2 + n))
benchmark_driver!(matrix_test, x; f_displayname="matrix_test")

####################################################################
relu(x) = log.(1.0 .+ exp.(x))
sigmoid(n) = 1. / (1. + exp(-n))

function neural_net(w1, w2, w3, x1, callback=nothing)
    x2 = relu(w1 * x1)
    x3 = relu(w2 * x2)
    ret = sigmoid(LinearAlgebra.dot(w3, x3))
    callback !== nothing && callback(ret)
    return ret
end

xs = (randn(10,10), randn(10,10), randn(10), rand(10))
benchmark_driver!(neural_net, xs...)

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

println("done")

max_stars_repo_name: TuringLang/Libtask

max_issues_repo_name: TuringLang/Libtask

max_forks_repo_name: TuringLang/Libtask


text:
abstract type ASetOperator <: AOperator end
abstract type ASetUnion <: ASetOperator end
struct SetUnion <: ASetUnion
end
∪(x::Node, y::Node) = Node("∪", ASetUnion, (x, y))

abstract type ASetIntersection <: ASetOperator end
struct SetIntersection <: ASetIntersection
end
∩(x::Node, y::Node) = Node("∩", ASetIntersection, (x, y))

abstract type ASetDifference <: ASetOperator end
struct SetDifference <: ASetDifference
end

abstract type AConditionSet <: ASetOperator end
struct ConditionSet <: AConditionSet
end
ConditionSet(x::Node, condition::Node) = Node("|", AConditionSet, (x, condition))  # {x | condition}


text:
using DiffEqJump, Catalyst, DelaySSAToolkit
rn = @reaction_network begin
   C, 0 --> Xₐ
   γ, Xₐ --> 0
   β, Xₐ --> Xᵢ
   γ, Xᵢ --> 0
end C γ β
jumpsys = convert(JumpSystem, rn, combinatoric_ratelaws = false)

u0 = [0, 0]
tf = 30.
saveat = .1
C, γ, β = [2., 0.1, 0.5]
p = [C, γ, β]
tspan = (0.,tf)
dprob = DiscreteProblem(u0, tspan, p)

τ = 15.
delay_trigger_affect! = function (integrator, rng)
   append!(integrator.de_chan[1], τ)
end
delay_trigger = Dict(3=>delay_trigger_affect!)
delay_complete = Dict(1=>[2=>-1])
delay_affect! = function (integrator, rng)
    i = rand(rng, 1:length(integrator.de_chan[1]))
    deleteat!(integrator.de_chan[1],i)
end
delay_interrupt = Dict(4=>delay_affect!)
delaysets = DelayJumpSet(delay_trigger,delay_complete,delay_interrupt)

algo_list =[DelayDirect(),DelayRejection(),DelayMNRM(),DelayDirectCR()]

djprob = DelayJumpProblem(jumpsys, dprob, algo_list[1],  delaysets, de_chan0, save_positions=(false,false))
sol =@benchmark solve(djprob, SSAStepper(),saveat =.1)

djprob = DelayJumpProblem(jumpsys, dprob,algo_list[2],  delaysets, de_chan0, save_positions=(false,false))
sol =@benchmark solve(djprob, SSAStepper(),saveat =.1)

djprob = DelayJumpProblem(jumpsys, dprob,algo_list[3],  delaysets, de_chan0, save_positions=(false,false))
sol =@benchmark solve(djprob, SSAStepper(),saveat =.1)

djprob = DelayJumpProblem(jumpsys, dprob,algo_list[4],  delaysets, de_chan0, save_positions=(false,false))
sol =@benchmark solve(djprob, SSAStepper(),saveat =.1)

meanlist=[]
medianlist=[]
minlist=[]
maxlist=[]
stdlist=[]

for algo in algo_list
    djprob = DelayJumpProblem(jumpsys, dprob,algo,delaysets, de_chan0, save_positions=(false,false))
    a=@benchmark solve(djprob, SSAStepper())
    push!(meanlist,mean(a).time/1e9)
    push!(medianlist,median(a).time/1e9)
    push!(minlist,minimum(a).time/1e9)
    push!(maxlist,maximum(a).time/1e9)
    push!(stdlist,std(a).time/1e9)
    print(djprob.aggregator," OK","\n")
end

df=DataFrame(mean=meanlist,median=medianlist,min=minlist,max=maxlist,std=stdlist)
CSV.write("test/zxytest/results/delay_degradation_save_positions=f.csv",df)

df=CSV.read("C:/Users/86158/Desktop/algotest/dep_gr_save_positions=F.csv",DataFrame)
medianlist=df.median

sa=[string(round(mt,digits=6),"s") for mt in meanlist]
medianvalue=[string(round(mt,digits=6),"s") for mt in medianlist]

algo_name = ["Direct","Rejection","MNRM","DirectCR"]
using Plots
bar(algo_name,meanlist,legend=:false,title="delay_degradation",ylabel="meantime")
scatter!(algo_name, 0.00015 .+ meanlist , markeralpha=0, series_annotations=sa)

bar(algo_name,medianlist,legend=:false,title="delay_degradation",ylabel="mediantime")
scatter!(algo_name, 0.00015 .+ medianlist , markeralpha=0, series_annotations=medianvalue)


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

include(joinpath(dirname(@__DIR__), "src", "EyeTrackingUtils.jl"))

using .EyeTrackingUtils
using Test
using CSV, DataFrames # CSV is currently a test-specific dependency

const data_file = joinpath(dirname(@__DIR__), "data", "word_recognition.csv")

@time @testset "EyeTracking.jl" begin
    word_recognition = DataFrame(CSV.File(data_file))
    # https://github.com/jwdink/eyetrackingR/blob/master/tests/testthat/test_analyze_time_bins.R#L29-L35
    data = make_clean_data(
        word_recognition,
        participant_column = :ParticipantName,
        trial_column = :Trial,
        time_column = :TimeFromTrialOnset,
        trackloss_column = :TrackLoss,
        aoi_columns = [:Animate, :Inanimate],
        treat_non_aoi_looks_as_missing = true,
        item_columns = nothing
        )
    @test true
end


text:
mutable struct MatrixFreeOperator{F,N,S,O} <: AbstractMatrixFreeOperator{F}
  f::F
  args::N
  size::S
  opnorm::O
  ishermitian::Bool
  function MatrixFreeOperator(f::F, args::N;
                              size=nothing, opnorm=true, ishermitian=false) where {F,N}
    @assert (N <: Tuple && length(args) in (1,2)) "Arguments of a "*
    "MatrixFreeOperator must be a tuple with one or two elements"
    return new{F,N,typeof(size),typeof(opnorm)}(f, args, size, opnorm, ishermitian)
  end
end
MatrixFreeOperator(f) = MatrixFreeOperator(f, (nothing,))

# Overloading Base functions
function Base.size(M::MatrixFreeOperator)
  M.size === nothing && error("M.size is nothing, please define size as a tuple of integers")
  return M.size
end
@inline function Base.size(M::MatrixFreeOperator, n)
  M.size === nothing && error("M.size is nothing, please define size as a tuple of integers")
  n <= 0 && error("dimension out of range")
  return n <= length(M.size) ? M.size[n] : 1
end
@inline ==(M1::MatrixFreeOperator, M2::MatrixFreeOperator) = M1.f == M2.f && M1.args == M2.args && M1.size == M2.size && M1.opnorm == M2.opnorm && M1.ishermitian == M2.ishermitian
@inline Base.:*(A::MatrixFreeOperator, X::AbstractVecOrMat) = mul!(similar(X), A, X)

# Overloading LinearAlgebra functions
LinearAlgebra.ishermitian(M::MatrixFreeOperator) = M.ishermitian
function LinearAlgebra.opnorm(M::MatrixFreeOperator, p::Real)
  M.opnorm === nothing && error("""
    M.opnorm is nothing, please define opnorm as a function that takes one
    argument. E.g. `(p::Real) -> p == Inf ? 100 : error("only Inf norm is
    defined")`
  """)
  opn = M.opnorm
  return opn isa Number ? opn : M.opnorm(p)
end

# DiffEqBase
DiffEqBase.numargs(::MatrixFreeOperator) = 4

# Interface
isconstant(M::MatrixFreeOperator) = length(M.args) == 1
function update_coefficients!(M::MatrixFreeOperator, u, p, t)
  !isconstant(M) && (M.args = (p, t))
  return M
end

function (M::MatrixFreeOperator{F,N})(du, u, p, t) where {F,N}
  update_coefficients!(M,u,p,t)
  if isconstant(M)
    M.f(du, u, p)
  else
    M.f(du, u, p, t)
  end
  du
end

function (M::MatrixFreeOperator{F,N})(u, p, t) where {F,N}
  update_coefficients!(M,u,p,t)
  du = similar(u)
  if isconstant(M)
    M.f(du, u, p)
  else
    M.f(du, u, p, t)
  end
  du
end

@inline function mul!(y::AbstractVector, A::MatrixFreeOperator{F,N}, x::AbstractVector) where {F,N}
  if isconstant(A)
    A.f(y, x, A.args[1])
  else
    A.f(y, x, A.args[1], A.args[2])
  end
  return y
end

@inline function mul!(Y::AbstractMatrix, A::MatrixFreeOperator{F,N}, X::AbstractMatrix) where {F,N}
  m,  n = size(Y)
  k, _n = size(X)
  if n != _n
    throw(DimensionMismatch("the second dimension of Y, $N, does not match the second dimension of X, $_n"))
  end
  for i in 1:n
    y = view(Y, :, i)
    x = view(X, :, i)
    mul!(y, A, x)
  end
  Y
end


text:
isvar(n::Number) = false
isvar(n::Symbol) = !endswith(string(n), "_const")
isvar(n::Expr) = n.head == :(::) ? isvar(n.args[1]) : true
isconst(n) = !isvar(n)

const reversenames = Dict(:times => :(*), :plus => :(+), :divide => :(/), :minus => :(-), :power => :(^),
                          :dot_times => :(.*), :dot_plus => :(.+), :dot_divide => :(./), :dot_minus => :(.-), :dot_power => :(.^))

function delta(ops)
    funcname = Symbol("δ$(ops.name)")
    func = :(function $funcname($(ops.inputs...)); end)
    body = func.args[2].args::Vector{Any}
    empty!(body)
    nablas = Dict()
    last_info = [Expr(:line)]
    for line in ops.body
        push_if_changed!(body, last_info, line.info)
        name = replace(string(line.name), r"^\.", "")
        constname = Symbol("δ$(name)_const")
        if line.name == :ref_const2
            nablas[line] = (:ref_const2, line.inputs[2])
            push!(body, :($(toexpr(line.outputs)) = $(line.inputs[1])[$(line.inputs[2])]))
        elseif isdefined(constname) || all(isconst, line.inputs) || all(isconst, line.outputs)
            sname = Symbol(name)
            sname = get(reversenames, sname, sname)
            if sname == :fanout
                temp = gensym(name)
                push!(body, :($temp = $(line.inputs...)))
                [push!(body, :($(line.outputs[k]) = $temp[$k])) for k in 1:length(line.outputs)]
            elseif startswith(string(line.name), ".")
                push!(body, :($(toexpr(line.outputs)) = $sname.($(line.inputs...))))
            else
                push!(body, :($(toexpr(line.outputs)) = $sname($(line.inputs...))))
            end
        else
            nabla = gensym("∇" * name)
            nablas[line] = nabla
            temp = gensym(name)
            funcname = Symbol("δ" * name)
            # push!(body, :(($(line.outputs...), $nabla) = $name($(line.inputs...)))),
            # work around for https://github.com/JuliaLang/julia/issues/15276
            if startswith(string(line.name), ".")
                mcast = Symbol("δ$(name_const("multicast", line.inputs, line.outputs))")
                push!(body, :($temp = $mcast($funcname, $(line.inputs...))))
            else
                push!(body, :($temp = $funcname($(line.inputs...))))
            end
            [push!(body, :($(line.outputs[k]) = $temp[$k])) for k in 1:length(line.outputs)]
            push!(body, :($nabla = $temp[$(length(line.outputs)+1)]))
            # end work around
        end
    end
    push!(body, ∇(ops, nablas))
    push!(body, :($(ops.outputs...), $(Symbol("∇$(ops.name)"))))
    #    @show func
    func
end

function ∇(ops, nablas)
    name = Symbol("∇$(ops.name)")
    inputs = map(topartial, filter(isvar, ops.outputs))
    func = :(function $name($(inputs...)); end)
    if length(inputs) == 1
        func.args[1].args[2] = Expr(:kw, func.args[1].args[2], 1.0)
    end
    body = func.args[2].args::Vector{Any}
    empty!(body)
    dupes = Set(inputs)
    emptys = Set()
    last_info = [Expr(:line)]
    for line in reverse(ops.body)
        push_if_changed!(body, last_info, line.info)
        if haskey(nablas, line)
            nabla = nablas[line]
            if isa(nabla, Tuple) && nabla[1] == :ref_const2
                ins = map(topartial, filter(isvar, line.outputs))
                outs = map(topartial, filter(isvar, line.inputs))
                if outs[1] in dupes
                    push!(body, :($(outs[1])[$(nabla[2])] += $(toexpr(ins))))
                else
                    push!(body, :($(outs[1]) = δzeros($(line.inputs[1]))))
                    push!(body, :($(outs[1])[$(nabla[2])] = $(toexpr(ins))))
                    push!(dupes, outs[1])
                end
            else
                ins = map(topartial, filter(isvar, line.outputs))
                outs = map(topartial, filter(isvar, line.inputs))
                if length(outs) > 0
                    for o in filter(isvar, line.outputs)
                        op = topartial(o)
                        if op in emptys && !(op in dupes)
                            push!(body, :($op = δzeros($o)))
                        end
                    end
                    dedup = [k in dupes ? gensym(k) : (push!(dupes, k); k) for k in outs]
                    push!(body, :($(toexpr(dedup)) = $nabla($(ins...))))
                    [push!(body, :($(outs[k]) += $(dedup[k]))) for k in find(outs .!= dedup)]
                end
            end
        else
            [push!(emptys, topartial(o)) for o in filter(isvar, line.inputs)]
        end
    end
    push!(body, toexpr(map(topartial, filter(isvar, ops.inputs))))
    func
end

topartial(expr::Symbol) = Symbol("∂$expr")
topartial(expr::Expr) = Symbol("∂$(expr.args[1])")
toexpr(symbols) = length(symbols) == 1 ? symbols[1] : Expr(:tuple, symbols...)

function push_if_changed!(body, last_info, info)
    if last_info[1] != info
        last_info[1] = info
        push!(body, info)
    end
end


text:
# Julia wrapper for header: /usr/local/ffmpeg/2.4/include/libavutil/channel_layout.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_get_channel_layout,
    av_get_channel_layout_string,
    av_bprint_channel_layout,
    av_get_channel_layout_nb_channels,
    av_get_default_channel_layout,
    av_get_channel_layout_channel_index,
    av_channel_layout_extract_channel,
    av_get_channel_name,
    av_get_channel_description,
    av_get_standard_channel_layout

function av_get_channel_layout(name)
    ccall((:av_get_channel_layout,libavutil),Uint64,(Ptr{Uint8},),name)
end

function av_get_channel_layout_string(buf,buf_size::Integer,nb_channels::Integer,channel_layout::Uint64)
    ccall((:av_get_channel_layout_string,libavutil),Void,(Ptr{Uint8},Cint,Cint,Uint64),buf,buf_size,nb_channels,channel_layout)
end

function av_bprint_channel_layout(bp,nb_channels::Integer,channel_layout::Uint64)
    ccall((:av_bprint_channel_layout,libavutil),Void,(Ptr{AVBPrint},Cint,Uint64),bp,nb_channels,channel_layout)
end

function av_get_channel_layout_nb_channels(channel_layout::Uint64)
    ccall((:av_get_channel_layout_nb_channels,libavutil),Cint,(Uint64,),channel_layout)
end

function av_get_default_channel_layout(nb_channels::Integer)
    ccall((:av_get_default_channel_layout,libavutil),Int64,(Cint,),nb_channels)
end

function av_get_channel_layout_channel_index(channel_layout::Uint64,channel::Uint64)
    ccall((:av_get_channel_layout_channel_index,libavutil),Cint,(Uint64,Uint64),channel_layout,channel)
end

function av_channel_layout_extract_channel(channel_layout::Uint64,index::Integer)
    ccall((:av_channel_layout_extract_channel,libavutil),Uint64,(Uint64,Cint),channel_layout,index)
end

function av_get_channel_name(channel::Uint64)
    ccall((:av_get_channel_name,libavutil),Ptr{Uint8},(Uint64,),channel)
end

function av_get_channel_description(channel::Uint64)
    ccall((:av_get_channel_description,libavutil),Ptr{Uint8},(Uint64,),channel)
end

function av_get_standard_channel_layout(index::Integer,layout,name)
    ccall((:av_get_standard_channel_layout,libavutil),Cint,(Uint32,Ptr{Uint64},Ptr{Ptr{Uint8}}),index,layout,name)
end


text:
module ModMesh
# Creating the mesh
#
# Hongyang Zhou, hyzhou@umich.edu 05/12/2019

export Mesh, createMesh, createTriMesh

using Triangle

const ϵ = 1e-5
const East_,West_,North_,South_ = 1,2,3,4
const INF = 10000.0 # For determining point inside polygon or not

struct Mesh
   Connectivity::Array{Int,2}
   BClines::Array{Float64,2}
   BCDir::Vector{Int}
   Coords::Array{Float64,2}
   BCnode::Vector{Vector{Int}}
end

"""
Given three colinear points p, q, r, the function checks if
point q lies on line segment 'pr'
"""
function onSegment(p::Vector{Float64}, q::Vector{Float64}, r::Vector{Float64})

   if (q[1] ≤ max(p[1], r[1]) && q[1] ≥ min(p[1], r[1]) &&
      q[2] ≤ max(p[2], r[2]) && q[2] ≥ min(p[2], r[2]))
      return true
   else
      return false
   end

end

"""
To find orientation of ordered triplet (p, q, r).
The function returns following values
0 --> p, q and r are colinear
1 --> clockwise
2 --> counterclockwise
"""
function orientation(p::Vector{Float64}, q::Vector{Float64}, r::Vector{Float64})

   val = (q[2] - p[2])*(r[1] - q[1]) - (q[1] - p[1])*(r[2] - q[2])

   if val == 0
      return 0 # colinear
   else
      return val > 0 ? 1 : 2 # clockwise or counter-clockwise
   end
end

function doIntersect(p1::Vector{Float64},q1::Vector{Float64},
   p2::Vector{Float64}, q2::Vector{Float64})
   # Find the four orientations needed for general and special cases
   o1 = orientation(p1, q1, p2)
   o2 = orientation(p1, q1, q2)
   o3 = orientation(p2, q2, p1)
   o4 = orientation(p2, q2, q1)

   # General case
   if o1 != o2 && o3 != o4 return true end

   # Special Cases
   # p1, q1 and p2 are colinear and p2 lies on segment p1q1
   if o1 == 0 && onSegment(p1, p2, q1) return true end

   # p1, q1 and p2 are colinear and q2 lies on segment p1q1
   if o2 == 0 && onSegment(p1, q2, q1) return true end

   # p2, q2 and p1 are colinear and p1 lies on segment p2q2
   if o3 == 0 && onSegment(p2, p1, q2) return true end

   # p2, q2 and q1 are colinear and q1 lies on segment p2q2
   if o4 == 0 && onSegment(p2, q1, q2) return true end

   return false # Doesn't fall in any of the above cases
end

"Returns true if the point p lies inside the polygon[] with n vertices."
function isInside(point::Vector{Float64},BClines::Array{Float64,2})

   count = 0
   nlines = size(BClines)[2]
   # loop over all vertical boundary lines
   for i=1:nlines
      line = BClines[:,i]
      if doIntersect(line[1:2],line[3:4],point,[INF, point[2]]) count += 1 end
   end

   return isodd(count)
end

# Rules:
# 1. boundaries are ordered counter-clockwisely.
function createMesh(BClines::Array{Float64,2},nI::Int, nJ::Int,
   LatticeType::Int)

   nlines = size(BClines,2)
   BCDir = zeros(Int,nlines)
   Coords = Array{Float64}(undef, 0, 0)

   # Check the orientation of linestyle
   for i = 1:nlines
      if abs(BClines[1,i] - BClines[3,i]) < ϵ
         BClines[2,i] > BClines[4,i] ? BCDir[i] = West_ : BCDir[i] = East_
      else
         BClines[1,i] > BClines[3,i] ? BCDir[i] = North_ : BCDir[i] = South_
      end
   end

   # Find the outermost locations
   xmin,xmax = extrema(BClines[:,1:2:3])
   ymin,ymax = extrema(BClines[:,2:2:4])

   # Shift a little bit to avoid issue at the boundary
   xmesh = range(xmin+ϵ,stop=xmax-ϵ,length=nI)
   ymesh = range(ymin+ϵ,stop=ymax-ϵ,length=nJ)

   dx = xmesh[2] - xmesh[1]
   dy = ymesh[2] - ymesh[1]

   # Check if the point lies inside the region
   for x in xmesh
      for y in ymesh
         if isInside([x,y],BClines)
            # Add to list
            try
               Coords = [Coords [x,y]]
            catch
               Coords = [x,y]
            end
         end
      end
   end

   # Create Connectivity list
   nNode = size(Coords,2)
   if LatticeType == 1 || LatticeType == 2
      Connectivity = fill(-1, (4,nNode))
   else
      Connectivity = fill(-1, (8,nNode))
   end
   for i = 1:nNode
      # East
      iNeiEast =
      findfirst([ abs(Coords[1,i]+dx-Coords[1,j]) < ϵ &&
      abs(Coords[2,i]-Coords[2,j]) < ϵ for j=1:nNode])

      # West
      iNeiWest =
      findfirst([ abs(Coords[1,i]-dx-Coords[1,j]) < ϵ &&
      abs(Coords[2,i]-Coords[2,j]) < ϵ for j=1:nNode])

      # South
      iNeiSouth =
      findfirst([ abs(Coords[1,i]-Coords[1,j]) < ϵ &&
      abs(Coords[2,i]-Coords[2,j]-dy) < ϵ for j=1:nNode])

      # North
      iNeiNorth =
      findfirst([ abs(Coords[1,i]-Coords[1,j]) < ϵ &&
      abs(Coords[2,i]-Coords[2,j]+dy) < ϵ for j=1:nNode])

      # SouthEast
      iNeiSouthEast =
      findfirst([ abs(Coords[1,i]+dx-Coords[1,j]) < ϵ &&
      abs(Coords[2,i]-Coords[2,j]-dy) < ϵ for j=1:nNode])

      # SouthWest
      iNeiSouthWest =
      findfirst([ abs(Coords[1,i]-dx-Coords[1,j]) < ϵ &&
      abs(Coords[2,i]-Coords[2,j]-dy) < ϵ for j=1:nNode])

      # NorthEast
      iNeiNorthEast =
      findfirst([ abs(Coords[1,i]+dx-Coords[1,j]) < ϵ &&
      abs(Coords[2,i]-Coords[2,j]+dy) < ϵ for j=1:nNode])

      # NorthWest
      iNeiNorthWest =
      findfirst([ abs(Coords[1,i]-dx-Coords[1,j]) < ϵ &&
      abs(Coords[2,i]-Coords[2,j]+dy) < ϵ for j=1:nNode])

      if iNeiWest  != nothing Connectivity[1,i] = iNeiWest end
      if iNeiEast  != nothing Connectivity[2,i] = iNeiEast end
      if iNeiSouth != nothing Connectivity[3,i] = iNeiSouth end
      if iNeiNorth != nothing Connectivity[4,i] = iNeiNorth end

      #
      if iNeiSouthEast == nothing &&
         iNeiSouth != nothing && iNeiEast != nothing
         Connectivity[[2,3],i] .= -1
      end
      if iNeiSouthWest == nothing &&
         iNeiSouth != nothing && iNeiWest != nothing
         Connectivity[[1,3],i] .= -1
      end
      if iNeiNorthEast == nothing &&
         iNeiNorth != nothing && iNeiEast != nothing
         Connectivity[[2,4],i] .= -1
      end
      if iNeiNorthWest == nothing &&
         iNeiNorth != nothing && iNeiWest != nothing
         Connectivity[[1,4],i] .= -1
      end
      #
   end

   BCnode = Vector{Vector{Int}}(undef,nlines)
   for i=1:nlines BCnode[i] = [] end

   for i = 1:nNode
      if Connectivity[1,i] == -1
         Coord = [Coords[:,i],[-INF,Coords[2,i]]]
         for iBC = 1:nlines
            if BCDir[iBC]<=2 && doIntersect(
               BClines[1:2,iBC].+[0.0,dy],BClines[3:4,iBC].+[0.0,-dy],
               Coord[1],Coord[2])
               push!(BCnode[iBC],i)
               break
            end
         end
      end
      if Connectivity[2,i] == -1
         Coord = [Coords[:,i],[INF,Coords[2,i]]]
         for iBC = 1:nlines
            if BCDir[iBC]<=2 && doIntersect(
               BClines[1:2,iBC].+[0.0,-dy],BClines[3:4,iBC].+[0.0,dy],
               Coord[1],Coord[2])
               push!(BCnode[iBC],i)
               break
            end
         end
      end
      if Connectivity[3,i] == -1
         Coord = [Coords[:,i],[Coords[1,i],-INF]]
         for iBC = 1:nlines
            if BCDir[iBC] > 2 && doIntersect(
               BClines[1:2,iBC].+[-dx,0.0],BClines[3:4,iBC].+[dx,0.0],
               Coord[1],Coord[2])
               push!(BCnode[iBC],i)
               break
            end
         end
      end
      if Connectivity[4,i] == -1
         Coord = [Coords[:,i],[Coords[1,i],INF]]
         for iBC = nlines:-1:1 # Ugly trick
            if BCDir[iBC] > 2 && doIntersect(
               BClines[1:2,iBC].+[dx,0.0],BClines[3:4,iBC].+[-dx,0.0],
               Coord[1],Coord[2])
               push!(BCnode[iBC],i)
               break
            end
         end
      end
   end

   mesh = Mesh(Connectivity,BClines,BCDir,Coords,BCnode)

   return mesh
end

"""
Create triangle mesh for visualization.
"""
function createTriMesh(mesh::Mesh)
   points = transpose(mesh.Coords)
   # The functions in Triangle does not accept abstract array.
   points = copy(points)
   points_map = Array{Int64,1}(collect(1:1:size(mesh.Coords)[2]))
   # Find edge list: actually the size is known beforehand
   edges_list_1 = Array{Int64}(undef,0)
   edges_list_2 = Array{Int64}(undef,0)

   BCnode = copy(mesh.BCnode)
   BCDir = mesh.BCDir

   nBClines = length(mesh.BCnode)
   push!(BCnode,BCnode[1]) # For easy implementation of the final line

   for i = 1:nBClines
      if BCDir[i] == East_ || BCDir[i] == South_
         for j = 1:size(BCnode[i])[1]
            if j == size(BCnode[i])[1]
               push!(edges_list_1, BCnode[i][j])
               push!(edges_list_2, BCnode[i+1][end])
               continue
            end
            push!(edges_list_1, BCnode[i][j])
            push!(edges_list_2, BCnode[i][j+1])
         end
      else
         for j = size(BCnode[i])[1]:-1:1
            if j == 1
               push!(edges_list_1, BCnode[i][j])
               push!(edges_list_2, BCnode[i+1][end])
               continue
            end
            push!(edges_list_1, BCnode[i][j])
            push!(edges_list_2, BCnode[i][j-1])
         end
      end
   end

   edges_list = [edges_list_1 edges_list_2]

   edge_boundary = [true for i=1:size(edges_list)[1]]

   triMesh = Triangle.constrained_triangulation(points,points_map,edges_list,
   edge_boundary)

   return triMesh
end

end


text:
@testset "Control Bound Constraint" begin

    # Test ControlBoundConstraint
    T = Float64
    dt = 0.1
    n = 4
    m = 6
    X = SVector{n,T}([13.0, 1.0, -12.0, 1.0])
    U = SVector{m,T}([13.0, 1.0, -12.0, 1.0, 2.0, 30.0])
    Z = KnotPoint(X,U,dt)
    u_max = [Inf,   Inf,  -11.0, 15.0, 2.0,  30.0]
    u_min = [-Inf, -10.0, -Inf,  1.0, -2.0, -30.0]

    con = ControlBoundConstraint(m,u_max=u_max, u_min=u_min)
    @test TrajectoryOptimization.evaluate(con,Z) == [-1.0, -14.0, 0.0, 0.0, -11.0, 0.0, -4.0, -60.0]
    con.inds == [3,4,5,6,8,10,11,12]
    @test (@ballocated TrajectoryOptimization.evaluate($con,$Z)) == 0

    function easy_jacobian(con,Z)
        U = control(Z)
        function local_evaluate(U)
            Z_ = KnotPoint(state(Z), U, Z.dt)
            return TrajectoryOptimization.evaluate(con,Z_)
        end
        return ForwardDiff.jacobian(local_evaluate, control(Z))
    end
    ∇c_easy = easy_jacobian(con, Z)
    ∇c = zeros(MMatrix{8,m,T,8*m})
    jacobian!(∇c,con,Z)
    @test norm(∇c - ∇c_easy, 1) < 1e-10
    @test (@ballocated jacobian!($∇c,$con,$Z)) == 0

end


text:
using LinearAlgebra: ldiv!
using SimpleSolvers
using Test

struct LinearSolverTest{T} <: LinearSolver{T} end

test_solver = LinearSolverTest{Float64}()

@test_throws ErrorException factorize!(test_solver)
@test_throws ErrorException ldiv!(test_solver)

A = [[+4.  +5.  -2.]
     [+7.  -1.  +2.]
     [+3.  +1.  +4.]]
x = [+4., -4., +5.]
b = [-14., +42., +28.]

function test_lu_solver(solver, A, b, x)
    for T in (Float64, ComplexF64, Float32, ComplexF32)
        AT = convert(Matrix{T}, A)
        bT = convert(Vector{T}, b)
        xT = convert(Vector{T}, x)

        lu1 = solver(AT)
        x1  = ldiv!(zero(xT), lu1, bT)
        @test x1 ≈ xT atol=8*eps(real(T))

        lu2 = solver(rand(T, size(A)...))
        x2  = zero(xT)
        factorize!(lu2, AT)
        ldiv!(x2, lu2, bT)
        @test x2 ≈ xT atol=8*eps(real(T))
    end
end

test_lu_solver(LUSolverLAPACK, A, b, x)
test_lu_solver(LUSolver, A, b, x)


text:
export
AABBNode,
AABBNodeData,
entity,
query,
aabb

import Base.isless
using Base.Iterators

struct AABBNodeData{T,S}
  aabb::AABB{T}
  entity::S
end

AABBNodeData{T,S}() where {T,S} = AABBNodeData{T,S}(AABB{T}(), nothing)

const AABBNode{T,S} = BinaryTree{AABBNodeData{T,S}}

# Traits #
abstract type QueryTrait  end
struct IntervalTree <: QueryTrait end
struct SearchTree <: QueryTrait end

key(n::AABBNodeData{T,S}) where {T,S} = n.aabb
aabb(n::AABBNodeData{T,S}) where {T,S} = n.aabb

#fullness(::Type{BinaryTree{AABBNodeData{T,S}}}) where {T,S} = FullTree
fullness(::Type{BT}) where {T, S, A<:AABBNodeData{T,S}, BT<:BinaryTree{A}} = FullTree

query_trait(::Type) = SearchTree()
query_trait(::Type{<:Node{AABBNodeData{T,S}}}) where {T,S} = IntervalTree()
query_trait(::Type{<:Node{<:BalanceData{AABBNodeData{T,S}}}}) where {T,S} = IntervalTree()

entity(n::AABBNode{T,S}) where {T,S} = data(n).entity
entity(n::AABBNodeData{T,S}) where {T,S} = n.entity
entity(_::Type{S}, n::AABBNodeData{T,S}) where{T,S} = n.entity
entity(s::Type{S}, n::Nothing) where {S} = zero(S)
data(n::AABBNodeData{T,S}) where {T,S} = n.entity
value(n::AABBNodeData{T,S}) where {T,S} = n.entity

function branch_compare( aabb::AABB{T} ) where {T}
  n -> begin
    if left(n) == nothing || right(n) == nothing ||
      volume(union(aabb,key(left(n)))) < volume(union(aabb,key(right(n))))
      :left
    else
      :right
    end
  end
end

partial(f,x) = y->f(x,y)

function rebuild_parent_data(direction,
                             lcd::AABBNodeData{T,S},
                             rcd::AABBNodeData{T,S},
                             parent::Union{Nothing, AABBNodeData{T,S}}) where {T,S}
  larger_aabb = union(lcd.aabb,rcd.aabb)
  AABBNodeData{T,S}(larger_aabb, entity(S,parent))
end

function rebuild_parent_data(direction,
                             lcd::Nothing,
                             rcd::AABBNodeData{T,S},
                             parent::Union{Nothing, AABBNodeData{T,S}}) where {T,S}
  AABBNodeData{T,S}(rcd.aabb, entity(S,parent))
end

function rebuild_parent_data(direction,
                             lcd::AABBNodeData{T,S},
                             rcd::Nothing,
                             parent::Union{Nothing, AABBNodeData{T,S}}) where {T,S}
  AABBNodeData{T,S}(lcd.aabb, entity(S,parent))
end

function rebuild_parent_data(direction,
                             lcd::Nothing,
                             rcd::Nothing,
                             parent::AABBNodeData{T,S}) where {T,S}
  parent
end

insert( t::AABBNode{T,S}, aabb::AABB{T}, entity::S ) where {T,S} = insert( t, AABBNodeData(aabb, entity) )

query( t::AABBNode{T,S}, qaabb::AABB{T} ) where {T,S} = entity.(Leaves(t,n->overlaps(data(n).aabb, qaabb )))


text:
"""
$(TYPEDEF)

Abstract type which indicates which regions of a [`TriMatrix`](@ref) are stored in
memory. Type parameter is a `Bool` indicating whether the matrix diagonal is
stored.
"""
abstract type TriLayout{D} end

# Check the type parameter of a TriLayout and throw an error if it isn't a Bool
_check_layout_param(D) = D isa Bool || throw(ArgumentError("Type parameter of TriLayout must be a Bool"))

for T in [:TriUpper, :TriLower, :TriSymmetric]
	@eval struct $T{D} <: TriLayout{D}
		$T{D}() where D = (_check_layout_param(D); new{D}())
		$T() = $T{true}()
	end
end

@doc """
$(TYPEDEF)

The upper triangle of the matrix is stored, values beneath the diagonal are
zero.
"""
TriUpper

@doc """
$(TYPEDEF)

The lower triangle of the matrix is stored, values above the diagonal are
zero.
"""
TriLower

@doc """
$(TYPEDEF)

Matrix is symmetric across the diagonal, one value is stored for each pair of
non-diagonal entries.
"""
TriSymmetric

const TriUL{D} = Union{TriUpper{D}, TriLower{D}}

"""
	$(FUNCTIONNAME)(::Type{<:TriLayout})::Bool
	$(FUNCTIONNAME)(layout::TriLayout)::Bool

Check if the given layout/layout type stores values along the diagonal.
"""
hasdiag(::Type{<:TriLayout{D}}) where D = D
hasdiag(::TriLayout{D}) where D = D

"""
	$(FUNCTIONNAME)(::Type{<:TriLayout}, n)::Integer
	$(FUNCTIONNAME)(layout::TriLayout, n)::Integer

Get the number of elements needed to store the data of a `TriMatrix`
with the given layout.
"""
nelems(::Type{L}, n::Integer) where {L<:TriLayout} = trinum(max(hasdiag(L) ? n : n - 1, 0))
nelems(layout::TriLayout, n::Integer) = nelems(typeof(layout), n)

"""
	$(FUNCTIONNAME)(::Type{<:TriLayout})
	$(FUNCTIONNAME)(layout::TriLayout)

Get the `TriLayout` of the transpose of a `TriMatrix` with the given layout.

This maps `TriUpper` and `TriLower` to each other, and `TriSymmetric` to itself.
"""
transpose_layout(::Type{TriUpper{D}}) where D = TriLower{D}
transpose_layout(::Type{TriLower{D}}) where D = TriUpper{D}
transpose_layout(::Type{TriSymmetric{D}}) where D = TriSymmetric{D}
transpose_layout(::Type{TriUpper}) = TriLower
transpose_layout(::Type{TriLower}) = TriUpper
transpose_layout(::Type{TriSymmetric}) = TriSymmetric
transpose_layout(::L) where {L <: TriLayout} = transpose_layout(L)()


text:
using Distributions, SpecialFunctions, Roots, LinearAlgebra
# import SpecialFunctions: lbeta
"""
GNas10: Modified algorithm of Nassar et al., J Neuro, 2010 for Gaussian prior
        for the Gaussian task (section 2.3.1)
"""
struct GNas10 <: GLearner
    χ_0::Float64
    ν_0::Float64
    mu_0::Float64
    sigma_0::Float64 # std of observations
    sigma::Float64 # std of observations
    changeprobability::Float64
    m::Float64
    ρ::Float64
    χ_n::Array{Float64, 1}
    ν_n::Array{Float64, 1}
    r_n::Array{Float64, 1}
    var_n::Array{Float64, 1}
    mu_n::Array{Float64, 1} #parameter we are interested in
    Sgm_n::Array{Float64, 1} # Bayes Factor Surprise
    γ_n::Array{Float64, 1}
    Ssh_n::Array{Float64, 1} # Shannon Surprise
end
function GNas10(; mu_0 = 0., sigma_0 = 1., sigma = 0.1, changeprobability = .001)
    ρ = sigma^2 / sigma_0^2
    m = changeprobability/ (1. - changeprobability)
    ν_0 = calcν(sigma, sigma_0)
    χ_0 = calcχ(mu_0, ν_0, sigma)
    var_n_0 = sigma^2 / ν_0
    mu_n_0 = sigma * χ_0 / ν_0
    GNas10(χ_0, ν_0, mu_0, sigma_0, sigma, changeprobability, m, ρ, [χ_0], [ν_0],
                [1.], [var_n_0], [mu_n_0], [1.], [m / (1. + m)], [1.])
end
export GNas10
function update!(learnerG::Union{GNas10, GNas10original}, y)
    calcSgm!(learnerG, y)
    calcγ!(learnerG) # Ω = learnerG.m * learnerG.Sgm_n[1]/(1. + learnerG.m * learnerG.Sgm_n[1])
    calcShannon!(learnerG, y)
    Δrule_mu_n = learnerG.mu_n[1] + (y - learnerG.mu_n[1])/(learnerG.ρ + learnerG.r_n[1] + 1.)
    Δrule_mu_0 = learnerG.mu_0 + (y - learnerG.mu_0)/(learnerG.ρ + 1.)
    computemu_n!(learnerG, Δrule_mu_n, Δrule_mu_0)
    computer_n!(learnerG)
    computevar_n!(learnerG)
    transformto_χ_ν!(learnerG)
end
export update!
function computemu_n!(learnerG::Union{GNas10, GNas12, GNas10original, GNas12original},
                        Δrule_mu_n, Δrule_mu_0)
    learnerG.mu_n[1] = (1. - learnerG.γ_n[1]) * Δrule_mu_n + learnerG.γ_n[1] * Δrule_mu_0
end
function computer_n!(learnerG::Union{GNas10, GNas10original})
    learnerG.r_n[1] = (1. - learnerG.γ_n[1]) * (learnerG.r_n[1] + 1.) + learnerG.γ_n[1]
end
function computevar_n!(learnerG::Union{GNas10, GNas10original})
    learnerG.var_n[1] = 1. / ((1. / learnerG.sigma_0^2) + (learnerG.r_n[1] / learnerG.sigma^2))
end
function transformto_χ_ν!(learnerG::Union{GNas10, GNas12, GNas10original, GNas12original})
    learnerG.ν_n[1] = calcν(learnerG.sigma, sqrt(learnerG.var_n[1]))
    learnerG.χ_n[1] = calcχ(learnerG.mu_n[1], learnerG.ν_n[1], learnerG.sigma)
end


text:
# Functions for SICR scripts to run set of models in tuple modeltypes on region

using Distributed
using DelimitedFiles
using Dates
using LinearAlgebra
using StatsBase
using Statistics
using Distributions
using DifferentialEquations
using DataFrames
using CSV

const txtstr = ".txt"

include("SICRmodels.jl")
include("SICURmodels.jl")
include("LinearModels.jl")

"""
	function mcmc(data,cols,p,N,total::Int,sicrfn!)

Metropolis-Hastings MCMC of the SIR ODE model assuming a Poisson likelihood
for the cumulative number of cases and deaths, and a noninformative prior
data = DataFrame of
cols = [active cases, deaths, recovered cases]
p[1:end-1] = array of parameters used in the model
p[end] = number infected on day 0 of data (Jan 22, 2020)
pfit = array of parameters to be fit
N = population of data region
total = number of MCMC steps
output = (loglikelihoods, parameters, maximum loglikelihood, max loglikelihood parameters)
sicrfn! = ODE file to use

"""
function mcmc(data::DataFrame,cols::Array,p::Array,pfit::Array,N::Float64,total::Int,sicrfn!,temp::Float64 = 1.,sigmaguess::Float64 = .002)

	chi = ll(data,cols,p,N,sicrfn!)
	accept = 0
	pt = copy(p)
	chiml = chi
	pml = copy(p)

	chiout = Array{Float64,1}(undef,total)
	pout = Array{Float64,2}(undef,total,length(p))

	for step in 1:total
		sampleparams!(pt,p,pfit,sigmaguess)
		chit = ll(data,cols,pt,N,sicrfn!)
		if rand() < exp((chi - chit)/temp)
			chi = chit
			p = copy(pt)
			accept += 1
			if chi < chiml
				chiml = chi
				pml = copy(p)
			end
		end
		chiout[step] = chi
		pout[step,:] = p
	end
	println("maxll: ",chiml)
	return chiout, pout, accept/total, chiml, pml
end

"""
    function simulate(p)

"""

function simulate(sicrfn!,p,firstday,lastday)

    tspan = (firstday,lastday)
    sim = Array{Int,2}(undef,lastday-firstday+1,3)
    uin = uinitial(data,cols,firstday,N)
    prediction,days = model(p[1:end-1],cols,uin,float.(tspan),N,sicrfn!)
    phi = p[end]
    for set in eachindex(cols)
        for i in eachindex(days)
            mu = max(prediction[i,set],eps(Float64))
            p = mu/(mu+phi)
            d = Distributions.NegativeBinomial(phi,p)
            sim[i,set] = rand(d)
        end
    end
    return sim

end

# generate MCMC step guess from current parameter value r using logNormal distirbution for all
# parameters except start day with scale parameter manually adjusted so approximately a third of guesses are accepted
function sampleparams!(pt,p,pfit,sigmaguess)
	for i in pfit
		d = Distributions.LogNormal(log(p[i]),sigmaguess)
		pt[i] = rand(d)
	end
end

# loglikelihood of data given parameters at r using a Poisson likelihood function
function ll(data,cols,p,N,sicrfn!)
	firstday = findfirst(data[:,cols[1]] .> 0)
	lastday = length(data[:,cols[1]])
	tspan = (firstday,lastday)
	uin = uinitial(data,cols,firstday,N)
	prediction,days = model(p,cols,uin,float.(tspan),N,sicrfn!)
	chi = 0
	for set in eachindex(cols)
		for i in eachindex(days)
			d = Distributions.Poisson(max(prediction[i,set],eps(Float64)))
			chi -= loglikelihood(d,[max(data[days[i],cols[set]],0)])
		end
	end
	return chi
end

function null_ll(data,cols,p,N,sicrfn!)
	firstday = findfirst(data[:,cols[1]] .> 0)
	lastday = length(data[:,cols[1]])
	tspan = (firstday,lastday)
	uin = uinitial(data,cols,firstday,N)
	prediction,days = model(p,cols,uin,float.(tspan),N,sicrfn!)
	chi = 0
	ent = 0
	for set in eachindex(cols)
		for i in eachindex(days)
			d = Distributions.Poisson(max(prediction[i,set],eps(Float64)))
			psample = rand(d)
			chi -= loglikelihood(d,[psample])
			ent += entropy(d)
		end
	end
	return chi,ent
end

function uinitial(data,cols,firstday,N)
	uin = Array{Float64,1}(undef,length(cols))
	for i in eachindex(cols)
		uin[i] = data[firstday,cols[i]]/N
	end
	return uin
end

# produce arrays of cases and deaths
function model(pin,cols,uin,tspan,N,sicrfn!)
	sol = modelsol(pin,uin,tspan,N,sicrfn!)
	prediction = Array{Float64,2}(undef,length(sol[1,:]),length(cols))
	prediction[:,1] = pin[2]*sol[1,:]*N
	for i = 2:length(cols)
		prediction[:,i] = pin[i+1]*sol[2,:]*N
	end
	return prediction, round.(Int,sol.t)
end

function modelsol(pin,uin,tspan,N,sicrfn!)
	p = pin[1:end-1]
	# I0 = (p[1]- p[3])/p[2] * uin[1] + pin[end]/N
	u0 = [pin[end]/N,uin[1]]
	prob = ODEProblem(sicrfn!,u0,tspan,p)
	solve(prob,saveat=1.)
end

function makehistograms(out)
	for i in 1:size(out)[2]
		subplot(2,3,i)
		hist(out[:,i],bins=100,density=true)
		title("param$(i)",fontsize=15)
	end
	tight_layout()
end

# returns predictions for C,R,D, matching data, and array of days
# solution array is in same row col order as data
# data columns have been ordered to match solutions.
function modelprediction(data,cols,p,N,sicrfn!)
	firstday = findfirst(data[:,cols[1]] .>0)
	lastday = length(data[:,cols[1]])
	sol = modelprediction(data,cols,p,N,sicrfn!,lastday)
	prediction = Array{Float64,2}(undef,length(sol[1,:]),length(cols))
	prediction[:,1] = p[2]*sol[1,:]*N
	for i = 2:length(cols)
		prediction[:,i] = p[i+1]*sol[2,:]*N
	end
	return prediction, data[firstday:end,cols], Int.(sol.t)
end

# returns solutions to model from day of first case to requested last day
function modelprediction(data,cols,p,N,sicrfn!,lastday)
	firstday = findfirst(data[:,cols[1]] .>0)
	tspan = (firstday,lastday)
	uin = uinitial(data,cols,firstday,N)
	modelsol(p,float.(uin),float.(tspan),N,sicrfn!)
end

# returns statistics for posteriors and other measures
function measures(psamples,pml)
	nparams = size(psamples,2)
	pout = Array{Tuple,1}(undef,nparams)
	for i in 1:nparams
		p = psamples[:,i]
		pout[i] = pml[i],mean(p),std(p),quantile(p,[.025,.5,.975])
	end

	R0 = r0(pml), mean(r0(psamples)),std(r0(psamples)),quantile(r0(psamples),[.025,.5,.975])
	return pout, R0
end

function r0(p::Array{Float64,1})
	beta,sigmac,sigmar,sigmad,q,f = p[1:6]
	sigma = sigmar + sigmad
	return (beta * (sigma + q * sigmac)) / (sigma * (sigmac + f * sigma))
end

function r0(p::Array{Float64,2})
	beta = p[:,1]
	sigmac = p[:,2]
	sigmar = p[:,3]
	sigmad = p[:,4]
	q = p[:,5]
	f = p[:,6]
	sigma = sigmar .+ sigmad
	return (beta .* (sigma .+ q .* sigmac)) ./ (sigma .* (sigmac .+ f .* sigma))
end

function runmcmc(total,datestrin,datestrout,temp=1.,N=1e7,pinstr= "pl",modeltype= ["rqfmt"])

    println(total,' ',temp)
    println(pinstr)

    pathin, pathout = makepaths(datestrin,datestrout)
    datafiles,regions = getdatafiles("data")
    println(datestrin,"->",datestrout)

    pstatstr = pathout * "pstats" * txtstr
    fpstat = open(pstatstr,"w")
    for j in eachindex(datafiles)
        datafile = datafiles[j]
        df = DataFrame(CSV.File("data/" * datafile))
        println(regions[j])
        for i in eachindex(modeltype)
            p,pfit,cols,data = mcmcargs(modeltype[i],regions[j],df,pathin,pinstr)
            # println(modeltype[i])
            outtuple = @spawn mcmc(data,cols,p,pfit,N,total,linearrd!,temp)

            model = modeltype[i]
            region = regions[j]
            writemeasures(outtuple,model,region,pstatstr)
            writecovariances(outtuple,model,region,pathout)
        end
    end
end

function runmcmcMP(total,datestrin,datestrout,temp=1.,N=1e7,pinstr= "pl",modeltype= ["rqfmt"])

    println(total,' ',temp)
    println(pinstr)

    pathin, pathout = makepaths(datestrin,datestrout)
    datafiles,regions = getdatafiles("data")
    println(datestrin,"->",datestrout)

    out = Array{Future,2}(undef,length(modeltype),length(datafiles))

    for j in eachindex(datafiles)
        datafile = datafiles[j]
        df = DataFrame(CSV.File("data/" * datafile))
        println(regions[j])
        for i in eachindex(modeltype)
            p,pfit,cols,data = mcmcargs(modeltype[i],regions[j],df,pathin,pinstr)
            # println(modeltype[i])
            out[i,j] = @spawn mcmc(data,cols,p,pfit,N,total,linearrd!,temp)
        end
    end

    pstatstr = pathout * "pstats" * txtstr
    fpstat = open(pstatstr,"w")
    for j in eachindex(datafiles), i in eachindex(modeltype)
        outtuple = fetch(out[i,j])
        model = modeltype[i]
        region = regions[j]
        writemeasures(outtuple,model,region,pstatstr)
        writecovariances(outtuple,model,region,pathout)
    end
end

# retries list of data files in folder datafolder
function getdatafiles(datafolder)
    datafiles = readdir("data")
    datafiles = datafiles[occursin.("covid_timeseries",datafiles)]
    regions = Array{String,1}(undef,0)
    for datafile in datafiles
        push!(regions,replace(replace(datafile,".csv"=>""),"covid_timeseries_"=>""))
    end
    return datafiles,regions
end

# returns prior for parameters for original r and u models
# param order: p =  beta,sigmac,sigmar,sigmad,q,f,mbase,mlocation,mrate,cmax,c50
function getpfit(modelname)
    if modelname == "r"
        return  [1;2;3;4;12]
    end
    if modelname == "rq"
        return [1;2;3;4;5;12]
    end
    if modelname == "rf"
        return [1;2;3;4;6;12]
    end
    if modelname == "rqf"
        return [1;2;3;4;5;6;12]
    end
    if modelname == "rm"
        return [1;2;3;4;7;8;9;12]
    end
    if modelname == "rqm"
        return [1;2;3;4;5;7;8;9;12]
    end
    if modelname == "rfm"
        return [1;2;3;4;6;7;8;9;12]
    end
    if modelname == "rqfm"
        return [1;2;3;4;5;6;7;8;9;12]
    end
    if modelname == "rmt"
        return [1;2;3;4;5;6;7;8;9;10;11;12]
    end
    if modelname == "rqmt" || modelname == "rfm"
        return [1;2;3;4;5;6;7;8;9;10;11;12]
    end
    if modelname == "rfmt"
        return [1;2;3;4;5;6;7;8;9;10;11;12]
    end
    if modelname == "rqfmt"
        return [1;2;3;4;5;6;7;8;9;10;11;12]
    end
    if modelname == "u"
        return [.25,.1,.1,.1]
    end
    if modelname == "uq" || modelname == "uf"
        return [.25,.1,.1,.1,.1]
    end
    if modelname == "um"
        return [.25,.1,.1,.1,5.,.1]
    end
    if modelname == "uqm" || modelname == "ufm"
        return [.25,.1,.1,.1,.1,5.,.1]
    end
    if modelname == "uqfm"
        return [.25,.1,.1,.1,.1,.1,5.,.1]
    end
end

# name = Dict(sicr! => "r", sicrq! => "rq", sicrf! => "rf", sicrqf! => "rqf", sicrm! => "rm", sicrqm! => "rqm", sicrfm! => "rfm", sicrqfm! => "rqfm",sicu! => "u", sicuq! => "uq", sicuf! => "uf", sicuq! => "uq", sicuqf! => "uqf", sicum! => "um", sicuqm! => "uqm", sicufm! => "ufm", sicuqfm! => "uqfm")

function makepaths(datestrin,datestrout)
    if ~isdir("Results")
        mkdir("Results")
    end
    pathin = "Results/$datestrin/"
    pathout = "Results/$datestrout/"
    if ~isdir(pathout)
        mkdir(pathout)
    end
    return pathin, pathout
end

# write out parameter and R0 statistics
function writemeasures(outtuple,name,region,pstatstr)
    pstats,R0 = measures(outtuple[2],outtuple[5])
    f = open(pstatstr,"a")
    writedlm(f,[region * "_" * name])
    writedlm(f,[outtuple[4] R0[1] R0[2] R0[3] R0[4]'])
    for ptup in pstats
        writedlm(f, [ptup[1] ptup[2] ptup[3] ptup[4]'])
    end
    close(f)
end

# write out max likelihood params
function writemaxll(outtuple,name,region,pathout)
    mlstr = pathout * "maxll" * "_" * region * "_" * name * txtstr
    f=open(mlstr,"w")
    f=open(mlstr,"a")
    writedlm(f,[outtuple[4];outtuple[5]])
    close(f)
end

# write out all loglikelihood mcmc samples
function writellsamples(outtuple,name,region,pathout)
    llstr = pathout * "llsamples" * "_" * region * "_" * name * txtstr
    f=open(llstr,"w")
    f=open(llstr,"a")
    writedlm(f,outtuple[1])
    close(f)
end

# write out last parameter in run
function writelastparam(outtuple,name,region,pathout)
    pstr = pathout * "pl" * "_" * region * "_" * name * txtstr
    f=open(pstr,"w")
    f=open(pstr,"a")
    writedlm(f,[outtuple[4];outtuple[2][end,:]])
    close(f)
end

# write out all parameters parameter in run
function writeparamsamples(outtuple,name,region,pathout)
    pstr = pathout * "paramsamples" * "_" * region * "_" * name * txtstr
    f=open(pstr,"w")
    f=open(pstr,"a")
    for i in 1:size(outtuple[2],1)
        writedlm(f,[outtuple[2][i,:]])
    end
    close(f)
end

function writecovariances(outtuple,name,region,pathout)
    pstr = pathout * "covcor" * "_" * region * "_" * name * ".csv"
    f=open(pstr,"w")
    f=open(pstr,"a")
    writedlm(f,cov(outtuple[2]),',')
    writedlm(f,cor(outtuple[2]),',')
    close(f)
end

# output arguments used in mcmc function
function mcmcargs(name,region,df,pathin,pinstr)
    colsr = [5,6,7]
    colsu = [5,6]
    if occursin("u",name)
        data = sicuprep(df)
        cols = colsu
    else
        data = df
        cols = colsr
    end

    infile = pathin * pinstr * "_" * region * "_" * name * txtstr
    if isfile(infile) && ~isempty(read(infile))
        p = readdlm(infile)[2:end]
        pfit = getpfit(name)
    else
        pfit = getpfit(name)
        p = [.25,.1,.01,.01,.01,1.,.1,1.,1.,.1,10.,1.]
        println(name,": No prior")
    end
    return p,pfit,cols,data
end

max_stars_repo_name: cannin/covid-sicr-test

max_issues_repo_name: cannin/covid-sicr-test

max_forks_repo_name: cannin/covid-sicr-test


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    alr(c)

Additive log-ratio transformation of composition `c`.
"""
function alr(c::Composition{D}) where {D}
  comps = components(c) .+ eps()
  SVector(ntuple(i -> log(comps[i] / comps[D]), D-1))
end

"""
    alrinv(x)

Inverse alr for coordinates `x`.
"""
alrinv(x::SVector{D,T}) where {D,T<:Real} =
  Composition(𝓒([exp.(x); SVector(one(T))]))

"""
    clr(c)

Centered log-ratio transformation of composition `c`.
"""
function clr(c::Composition{D}) where {D}
  comps = components(c) .+ eps()
  gmean = geomean(comps)
  SVector(ntuple(i -> log(comps[i] / gmean), D))
end

"""
    clrinv(x)

Inverse clr for coordinates `x`.
"""
clrinv(x) = Composition(𝓒(exp.(x)))

"""
    ilr(c)

Isometric log-ratio transformation of composition `c`.
"""
function ilr(c::Composition{D}) where {D}
  comps = components(c) .+ eps()
  logs  = log.(comps)
  T = eltype(comps)
  x = MVector(ntuple(i->zero(T), D-1))
  for i in 1:D-1
    s = zero(T)
    sqrtinv = 1/sqrt(i*(i+1))
    for j in 1:i+1
      if j < i+1
        s += - sqrtinv * logs[j]
      elseif j==i+1
        s += i * sqrtinv * logs[j]
      end
    end
    x[i] = s
  end
  SVector(x)
end

"""
    ilrinv(x)

Inverse ilr for coordinates `x`.
"""
function ilrinv(x::SVector{D,T}) where {D,T<:Real}
  z = MVector(ntuple(i->zero(T), D+1))
  for i in 1:D+1
    s = zero(T)
    for j in 1:D
      sqrtinv = 1/sqrt(j*(j+1))
      if i < j+1
        s += - sqrtinv * x[j]
      elseif i == j+1
        s += j * sqrtinv * x[j]
      end
    end
    z[i] = exp(s)
  end
  Composition(𝓒(z))
end


text:
import Base.push!
import Base: push!, isempty

# struct AssignPair
#     lhs::Symbol
#     rhs::Union{Symbol, Expr, Int}
# end
abstract type Stmt end

struct AssignPair <: Stmt
    lhs::Symbol
    rhs::Union{Symbol, Expr, Int}
end

function Base.show(io::IO, p::AssignPair)
    print(io, (lhs(p), rhs(p)))
end

# const CompoundStmt = Array{AssignPair,1}

struct CompoundStmt
    pairs::Array{AssignPair,1}
end

assignments(cs::CompoundStmt) = cs.pairs

push!(cs::CompoundStmt, p::AssignPair) = push!(cs.pairs, p)
isempty(cs::CompoundStmt) = isempty(cs.pairs)
merge(cs1::CompoundStmt, cs2::CompoundStmt) = CompoundStmt([cs1.pairs; cs2.pairs])

CompoundStmt() = CompoundStmt(Array{AssignPair,1}())
CompoundStmt(cs::CompoundStmt...) = CompoundStmt(vcat(assignments.(cs)))

# start(x::CompoundStmt) = false
# next(x::CompoundStmt, state) = (x, true)
# done(x::CompoundStmt, state) = state
# in(x::CompoundStmt, y::CompoundStmt) = x == y

function Base.show(io::IO, cs::CompoundStmt)
    print(io, "Compound(")
    for p in cs.pairs
        print(io, p)
    end
    print(io, ")")
end

# const AssignBlock = Array{AssignPair,1}

struct IfStmt <: Stmt
    tbranch
    ebranch
end

exists_else(stmt::IfStmt) = stmt.ebranch != nothing

IfStmt(b) = IfStmt(b, nothing)

else_branch(is::IfStmt) = is.ebranch
then_branch(is::IfStmt) = is.tbranch

struct AssignBlock <:Stmt
    stmts::Array{Union{CompoundStmt,IfStmt},1}
end

stmts(b::AssignBlock) = b.stmts

function Base.show(io::IO, b::IfStmt)
    print(io, "{THEN} ")
    print(io, then_branch(b))
    print(io, " {ELSE} ")
    print(io, else_branch(b))
end

AssignBlock() = AssignBlock(Array{Union{CompoundStmt,IfStmt},1}())

function Base.show(io::IO, b::AssignBlock)
    # print(io, "BLOCK")
    for stmt in b.stmts
        print(io, stmt)
    end
end

push!(b::AssignBlock, stmt::Union{CompoundStmt,IfStmt}) = push!(b.stmts, stmt)

lhs(a::AssignPair) = a.lhs
rhs(a::AssignPair) = a.rhs

variables(a::CompoundStmt) = lhs.(a.pairs)

struct SymbolicAssign
    fn::SymPy.SymFunction # symbolic function for recursively changed variable
    lc::SymPy.Sym         # symbol for loop counter
    lhs::SymPy.Sym
    rhs::SymPy.Sym
end

#-------------------------------------------------------------------------------

function transform(stmt::CompoundStmt, before::CompoundStmt, after::CompoundStmt)
    [before; stmt; after]
end

function transform(stmt::IfStmt, before::CompoundStmt, after::CompoundStmt)
    if exists_else(stmt)
        [transform(then_branch(stmt), before, after), transform(else_branch(stmt), before, after)]
    else
        [transform(then_branch(stmt), before, after), transform(CompoundStmt(), before, after)]
    end
end

function transform(block::AssignBlock, before::CompoundStmt, after::CompoundStmt)
    transform(stmts(block)..., before, after)
end

function transform(c1::CompoundStmt, stmt::IfStmt, before::CompoundStmt, after::CompoundStmt)
    transform(stmt, merge(before, c1), after)
end

function transform(c1::CompoundStmt, stmt::IfStmt, c2::CompoundStmt, before::CompoundStmt, after::CompoundStmt)
    transform(stmt, merge(before, c1), merge(c1, after))
end

function transform(stmt::IfStmt, c2::CompoundStmt, before::CompoundStmt, after::CompoundStmt)
    transform(stmt, before, merge(c1, after))
end

#-------------------------------------------------------------------------------

extract_assign(expr::LineNumberNode, level) = nothing

function extract_assign(expr::Expr, level)
    h = string(expr.head)
    if h == "while"
        # ignore guard in expr.args[1]
        return extract_assign(expr.args[2], level+2)
    elseif h == "if"
        # ignore guard in expr.args[1]
        return IfStmt([extract_assign(expr.args[i],level+2) for i in 2:length(expr.args)]...)
        # abc = [extract_assign(expr.args[1],level+2), extract_assign(expr.args[2],level+2)]
        # println("Abc: ", abc)
        return IfStmt(abc...)
    elseif h == "block"
        block = AssignBlock()
        cmps = CompoundStmt()
        for arg in expr.args
            stmt = extract_assign(arg, level+2)
            if isa(stmt, AssignPair)
                push!(cmps, stmt)
            elseif isa(stmt, IfStmt)
                if !isempty(cmps)
                    push!(block, cmps)
                end
                cmps = CompoundStmt()
                push!(block, stmt)
            end
        end
        if !isempty(cmps)
            push!(block, cmps)
        end
        return block
    elseif h == "="
        return AssignPair(expr.args[1], expr.args[2])
    end
end

#-------------------------------------------------------------------------------

function symbolic2(lc::Sym, recvars::Array{Symbol,1}, visited::Array{Symbol,1}, pair::AssignPair, pairs::AssignPair...)
    fn = SymFunction(string(lhs(pair)))
    sa = SymbolicAssign(fn, lc, symbolic(lhs(pair), fn, lc+1, recvars, visited, false), symbolic(rhs(pair), fn, lc, recvars, visited, true))
    rs = symbolic2(lc, recvars, [visited; lhs(pair)], pairs...)
    [sa; rs]
end

function symbolic2(lc::Sym, recvars::Array{Symbol,1}, visited::Array{Symbol,1})
    SymbolicAssign[]
end

#-------------------------------------------------------------------------------

function symbolic(loop::CompoundStmt, lc::Sym, recvars::Array{Symbol,1})
    symbolic2(lc, recvars, Symbol[], loop.pairs...)
end

# function symbolic(expr::AssignPair, lc::Sym, recvars::Array{Sym,1})
#     fn = SymFunction(string(lhs(expr)))
#     return SymbolicAssign(fn, lc, symbolic(lhs(expr), fn, lc+1, recvars, visited), symbolic(rhs(expr), fn, lc, recvars, visited))
# end

function symbolic(s::Symbol, f::SymFunction, lc::Sym, recvars::Array{Symbol,1}, visited::Array{Symbol,1}, rhs::Bool)
    if s in recvars
        f = SymFunction(string(s))
        if rhs && s in visited
            return f(lc + 1)
        else
            return f(lc)
        end
    end
    return SymPy.Sym(string(s))
end

function symbolic(i::Int, f::SymFunction, lc::Sym, recvars::Array{Symbol,1}, visited::Array{Symbol,1}, rhs::Bool)
    return Sym(i)
end

function symbolic(expr::Expr, f::SymFunction, lc::Sym, recvars::Array{Symbol,1}, visited::Array{Symbol,1}, rhs::Bool)
    if expr.head == :call && expr.args[1] in (:+, :-, :*, :/, :^)
        return eval(Expr(:call, expr.args[1], [symbolic(expr.args[i], f, lc, recvars, visited, rhs) for i in 2:length(expr.args)]... ))
    else
        error("Not supported rhs in assignment")
    end
end

#-------------------------------------------------------------------------------

function recurrence(expr::SymbolicAssign)
    rec = expr.lhs - expr.rhs
    eq2rec(rec, expr.fn, expr.lc)
end

function eq2rec(eq::Sym, fn::SymFunction, lc::Sym)

    fns = symfunctions(eq)
    w0 = Wild("w0")
    args = [get(match(fn(lc + w0), f), w0, nothing) for f in fns]
    args = filter(x -> x!=nothing, args)
    minidx = convert(Int, minimum(args))
    if minidx < 0
        eq = eq |> subs(lc, lc - minidx)
    else
        minidx = 0
    end
    ord = convert(Int, maximum(args)) - minidx
    coeffs = Sym[]
    for i in 0:ord
        c, eq = coeff_rem(eq, fn(lc + i))
        coeffs = [coeffs; c]
    end
    return CFiniteRecurrence(coeffs, fn, lc, eq)
end

#-------------------------------------------------------------------------------

function canonical!(loops::Array{CompoundStmt,1})
    vars = union(variables.(loops)...)
    for i in 1:length(loops)
        vmiss = setdiff(vars, variables(loops[i]))
        if length(vmiss) > 0
            amiss = [AssignPair(v, v) for v in vmiss]
            push!(loops[i], amiss...)
        end
    end
    loops
end

#-------------------------------------------------------------------------------

function flattenall(a::AbstractArray)
    while any(x->typeof(x)<:AbstractArray, a)
        a = collect(Iterators.flatten(a))
    end
    return a
end

function extract_loop(str::String)
    loops = extract_assign(Meta.parse(str), 0)
    loops = transform(loops, CompoundStmt(), CompoundStmt())
    loops = flattenall(loops)
    loops = filter(x -> !isempty(x), loops)

    if all(isa.(loops, CompoundStmt))
        loops = canonical!(Array{CompoundStmt,1}(loops))
    else
        error("Something went wrong while flattening the loop.")
    end

    recvars = union(variables.(loops)...)

    ls = SingleLoop[]
    for (i, loop) in enumerate(loops)
        lc = Sym("n_$(i)")
        recs = recurrence.(symbolic(loop, lc, recvars))
        loop = SingleLoop(recs, lc, Sym.(string.(variables(loop))))
        push!(ls, loop)
    end

    if length(ls) == 0
        return EmptyLoop()
    elseif length(ls) == 1
        return ls[1]
    else
        return MultiLoop(ls, ls[1].vars)
    end
end


text:
module SBMopt

using ..SBMs

export simpleopt!, fasteropt!, TMopt!, KLopt!
export multiple_trials

"""

    simpleopt!(sbm::SBM, nepoch=100)

Performs the most simple optimization procedure I could think of
(direct hill climbing)
"""
function simpleopt!(sbm::SBM, nepoch=100)

  lls = []

  for n in 1:nepoch
    #print("\rEpoch #$n")
    ll = 0

    for i in shuffle(SBMs.getindices(sbm))

      next_move, ll = SBMs.bestmove(sbm,i; forcemove=false)
      SBMs.enactmove(sbm,i,next_move)

      push!(lls,ll)
    end
  end
  println("")

  lls
end

"""
A faster version of simpleopt! using incrementally updated parameters.
The simple version works for smaller graphs (e.g. 40 nodes), but this scales
much better
"""
function fasteropt!(sbm::SBM, nepoch=100)

  lls = []

  for n in 1:nepoch
    #print("\rEpoch #$n")
    ll = 0
    ps = SBMs.computeparams(sbm)

    for i in shuffle(SBMs.getindices(sbm))

      next_move, ll = SBMs.bestmove(sbm,i,nothing,ps; forcemove=false)

      old_g = SBMs.getgroups(sbm)
      SBMs.enactmove(sbm,i,next_move)
      if old_g[i] != next_move  SBMs.updateparams!(sbm,ps,old_g)  end

      push!(lls,ll)
    end
  end
  println("")

  lls
end

"""

    TMopt!(sbm::SBM, nepoch=100)

Tommy's optimization algorithm
"""
function TMopt!(sbm::SBM, nepoch=100)

  lls = []

  ps = SBMs.computeparams(sbm)
  for n in 1:nepoch

    groups_start_epoch = SBMs.getgroups(sbm)

    for i in shuffle(SBMs.getindices(sbm))

      next_move, ll = SBMs.bestmove(sbm,i,nothing,ps; forcemove=false)

      old_g = SBMs.getgroups(sbm)
      SBMs.enactmove(sbm,i,next_move)

      if old_g[i] != next_move
        SBMs.updateparams!(sbm,ps,old_g)
        push!(lls,ll)
        break
      end

    end #for i in randperm

    groups_end_epoch = SBMs.getgroups(sbm)
    if minimum( groups_start_epoch .== groups_end_epoch )  break  end

  end #for n in 1:nepoch

  lls
end

function KLlikeopt!(sbm::SBM, nepoch=100)

  lls = [];

  for n in 1:nepoch
    print("\rEpoch #$n")

    ps = SBMs.computeparams(sbm)
    #keeping two copies so we can update parameters later
    g  = SBMs.getgroups(sbm)
    last_g = copy(g)

    epoch_moves = Tuple{Int,Int}[];
    epoch_lls = Float64[];

    start_ll = SBMs.loglikelihood(sbm,nothing,ps)

    to_move = IntSet(SBMs.getindices(sbm))

    while !isempty(to_move)

      vertex, move, ll = find_best_move_overall(sbm,g,ps, to_move)

      #enacting move for this epoch
      g[vertex] = move
      SBMs.updateparams!(sbm,ps,last_g,g)
      last_g[vertex] = move

      #bookkeeping
      delete!(to_move, vertex)
      push!(epoch_moves, (vertex,move))
      push!(epoch_lls, ll)

    end #while !empty(to_move)

    max_ll, max_i = findmax(epoch_lls)

    if start_ll >= max_ll  break  end

    #replaying changes up to max_ll over the epoch
    for i in 1:max_i
      v, cl = epoch_moves[i]
      SBMs.enactmove(sbm,v,cl)
      #recording relevant lls
      push!(lls, epoch_lls[i])
    end

  end

  lls
end

function find_best_move_overall(sbm,g,ps, candidates)

  max_ll = -Inf; max_vertex = 0; max_move = 0;

  for i in candidates

    best_move_i, best_ll_i = SBMs.bestmove(sbm,i,g,ps;forcemove=true)

    if best_ll_i > max_ll
      max_ll = best_ll_i
      max_vertex = i
      max_move = best_move_i
    end

  end

  max_vertex, max_move, max_ll
end

function multiple_trials!(sbm::SBM, opt_fn::Function, num_trials=10, nepochs=100)

  best_lls = []
  best_groups = []
  final_lls = []
  best_final_ll = -Inf
  final_groups = []

  for t in 1:num_trials

    if t == 1
      print("\rTrial $t")
    else
      print("\rTrial $t, Last Trial Time: $(round(toq(),3)) seconds")
    end
    tic();

    SBMs.randomize_g!(sbm)

    trial_lls = opt_fn(sbm,nepochs)

    push!(final_lls, trial_lls[end])
    push!(final_groups, SBMs.getgroups(sbm))

    if trial_lls[end] > best_final_ll
      best_lls = trial_lls
      best_groups = SBMs.getgroups(sbm)
    end
  end

  toq();

  final_lls, final_groups, best_lls, best_groups
end

end #SBMopt


text:
import Pkg.instantiate
instantiate()

import Aggregator.store_inferred
store_inferred()

max_stars_repo_name: Open-Reviews/Mangrove

max_issues_repo_name: Open-Reviews/Mangrove

max_forks_repo_name: Open-Reviews/Mangrove


text:
module LaxFriedrichs
export NumFlux
# Lax Friedrichs flux:
function NumFlux(F,a::Float64,b::Float64,alpha::Float64)
    0.5*(F.flux(a)+F.flux(b) - alpha*(b-a))
end
end

max_stars_repo_name: danielcort/benchmarks-python-julia-c

max_issues_repo_name: danielcort/benchmarks-python-julia-c

max_forks_repo_name: danielcort/benchmarks-python-julia-c


text:
module NLCOP

using Kombinator
using JuMP

export ExactNonlinearSolver, ApproximateNonlinearSolver

import Kombinator: solve

include("exact.jl")
include("approx.jl") # The actual NLCOP algorithm.

end


text:
using Random, Images, JLD2, LinearAlgebra
using JOLI, Statistics, FFTW
using JUDI
using Polynomials
using ArgParse
include("../utils/parse_cmd.jl")
parsed_args = parse_commandline()
L = parsed_args["nv"]

JLD2.@load "../models/Compass_tti_625m.jld2"
JLD2.@load "../data/Conc.jld2" conc

extentx = (n[1]-1)*d[1]
extentz = (n[2]-1)*d[2]

v = 1f3.*sqrt.(1f0./m)

idx_wb = find_water_bottom(rho.-rho[1,1])
idx_ucfmt = find_water_bottom((v.-3500f0).*(v.>3500f0))

Kh = zeros(Float32, n)
for i = 1:n[1]
    Kh[i,1:idx_wb[i]-1] .= 1f-4  # water layer
    Kh[i,idx_wb[i]:idx_ucfmt[i]-1] .= 1.03*1f-3*v[i,idx_wb[i]:idx_ucfmt[i]-1].+15f0
    Kh[i,idx_ucfmt[i]:idx_ucfmt[i]+7] .= 1f-4
    Kh[i,idx_ucfmt[i]+8:end] .= 1.03*1f-3*v[i,idx_ucfmt[i]+8:end].+200f0
end

phi = zeros(Float32,n)
for i = 1:n[1]
    for j = 1:n[2]
        p = Polynomial([-0.0314^2*Kh[i,j],2*0.0314^2*Kh[i,j],-0.0314^2*Kh[i,j],1.527^2])
        phi[i,j] = minimum(real(roots(p)[findall(real(roots(p)).== roots(p))]))
    end
    for j = idx_ucfmt[i]:idx_ucfmt[i]+7
        phi[i,idx_ucfmt[i]:idx_ucfmt[i]+7] = Float32.(range(0.056,stop=0.1,length=8))
    end
end

vp = convert(Array{Float32},v)
vs = vp ./ sqrt(3f0)

idx_wb = maximum(find_water_bottom(rho.-rho[1,1]))

function Patchy(sw, vp, vs, rho, phi; bulk_min = 36.6f9, bulk_fl1 = 2.735f9, bulk_fl2 = 0.125f9, ρw = 7.766f2, ρo = 1.053f3)
    bulk_sat1 = rho .* (vp.^2f0 - 4f0/3f0 .* vs.^2f0) * 1f3
    shear_sat1 = rho .* (vs.^2f0) * 1f3

    bulk_min = zeros(Float32,size(bulk_sat1))

    bulk_min[findall(v.>=3350)] .= 50f9   # mineral bulk moduli
    bulk_min[findall(v.<3350)] .= 1.2 * bulk_sat1[findall(v.<3350)] # mineral bulk moduli

    patch_temp = bulk_sat1 ./(bulk_min .- bulk_sat1) - bulk_fl1 ./ phi ./ (bulk_min .- bulk_fl1) + bulk_fl2 ./ phi ./ (bulk_min .- bulk_fl2)

    for i = 1:n[1]
        for j = idx_ucfmt[i]:idx_ucfmt[i]+8
            patch_temp[i,idx_ucfmt[i]:idx_ucfmt[i]+8] = patch_temp[i,idx_ucfmt[i]+9:idx_ucfmt[i]+17]
        end
    end

    bulk_sat2 = bulk_min./(1f0./patch_temp .+ 1f0)
    bulk_sat2[findall(bulk_sat2-bulk_sat1.>0)] = bulk_sat1[findall(bulk_sat2-bulk_sat1.>0)]

    bulk_new = (bulk_sat1+4f0/3f0*shear_sat1).*(bulk_sat2+4f0/3f0*shear_sat1) ./( (1f0.-sw).*(bulk_sat2+4f0/3f0*shear_sat1)
    + sw.*(bulk_sat1+4f0/3f0*shear_sat1) ) - 4f0/3f0*shear_sat1

	bulk_new[:,1:idx_wb] = bulk_sat1[:,1:idx_wb]
    bulk_new[findall(sw.==0)] = bulk_sat1[findall(sw.==0)]
    rho_new = rho + phi .* sw * (ρw - ρo) / 1f3
    rho_new[findall(sw.==0)] = rho[findall(sw.==0)]
    Vp_new = sqrt.((bulk_new+4f0/3f0*shear_sat1)./rho_new/1f3)
    Vp_new[findall(sw.==0)] = vp[findall(sw.==0)]
    Vs_new = sqrt.((shear_sat1)./rho_new/1f3)

    return Vp_new, Vs_new, rho_new
end

indices = Int.(round.(range(1, stop=size(conc,1), length=L)))
sw_stack = [Float32.(imresize([zeros(Float32,19,325);conc[indices[i],:,:]]',n)) for i = 1:L]

vp_stack = [Patchy(sw_stack[i], vp, vs, rho, phi)[1] for i = 1:L]
rho_stack = [Patchy(sw_stack[i], vp, vs, rho, phi)[3] for i = 1:L]

JLD2.@save "../models/timelapsevrho$(L)vint_finer.jld2" vp_stack rho_stack


text:
using Ramnet.Utils: stack
using Ramnet.Partitioners: RandomPartitioner, LinearPartitioner
using Ramnet.Models.Nodes: RegressionNode, GeneralizedRegressionNode
using Ramnet.Models: train!, predict
using Ramnet.Models: Discriminator, BleachingDiscriminator

@testset "Classification Discriminator" begin
    all_active = ones(Bool, 9)
    one_off    = Bool[0, ones(Bool, 8)...]

    d = Discriminator(length(all_active), 3)

    train!(d, all_active)

    @test predict(d, all_active) == 3
    @test predict(d, one_off) == 2

    inputs = stack(all_active, one_off)
    @test predict(d, inputs) == [3, 2]

    # Discriminator with externally-instantiated partitioner
    partitioner   = RandomPartitioner(length(all_active), 3)
    d_partitioner = Discriminator(partitioner)
    train!(d_partitioner, all_active)
    @test predict(d_partitioner, all_active) == 3

    # One-shot instantiation and training
    duplicated_input = stack(all_active, all_active)
    one_shot_d_1 = Discriminator(all_active, partitioner)
    one_shot_d_2 = Discriminator(all_active, 3)
    one_shot_d_3 = Discriminator(duplicated_input, 3)

    @test predict(one_shot_d_1, all_active) == 3
    @test predict(one_shot_d_2, all_active) == 3
    @test predict(one_shot_d_3, all_active) == 3

    # Discriminators with bleaching
    b_d = BleachingDiscriminator(length(all_active), 3)

    train!(b_d, all_active)

    @test predict(b_d, all_active) == 3
    @test predict(b_d, one_off) == 2

    @test predict(b_d, all_active; b=1) == 0
    @test predict(b_d, one_off; b=1) == 0

    train!(b_d, one_off)

    @test predict(b_d, one_off) == 3
    @test predict(b_d, all_active; b=1) == 2
    @test predict(b_d, one_off; b=1) == 2
end

@testset "Original and Discounted Regression Discriminator" begin
    X_train = Bool[
        1 1 0 0 0 0
        1 1 1 0 0 0
        1 1 1 1 0 0
        1 1 1 1 1 0
        1 1 1 1 1 1
    ]

    y_train = Float64[1, 2, 3, 4, 5]

    X_test = Bool[
        0 0 0 0 0 0
        1 0 0 0 0 0
        1 1 0 0 0 0
        1 1 1 1 1 1
        1 1 1 0 0 0
    ]

    # As in the original regression discriminator
    y_test_original = Float64[7 / 4, 7 / 4, 22 / 9, 32 / 9, 23 / 9]

    partitioner = LinearPartitioner(6, 2)

    regressor = Discriminator{LinearPartitioner,RegressionNode}(partitioner; γ=1.0)

    train!(regressor, X_train, y_train)

    @test predict(regressor, X_test) == y_test_original

    # Discounted regression discriminator
    e(n; γ=0.7) = (1 - γ^n) / (1 - γ) # sum of weights
    function s(vs; γ=0.7) # Discounted sum
        su = 0.0
        for (i, v) in Iterators.enumerate(Iterators.reverse(vs))
            su += γ^(i - 1) * v
        end

        su
    end

    y_test_discounted = Float64[
        (s([1]) + s([1,2,3])) / (e(1) + e(3)),
        (s([1]) + s([1,2,3])) / (e(1) + e(3)),
        (s([1,2,3,4,5]) + s([1]) + s([1,2,3])) / (e(5) + e(1) + e(3)),
        (s([1,2,3,4,5]) + s([3,4,5]) + s([5])) / (e(5) + e(3) + e(1)),
        (s([1,2,3,4,5]) + s([2]) + s([1,2,3])) / (e(5) + e(1) + e(3))
    ]

    regressor = Discriminator{LinearPartitioner,RegressionNode}(partitioner; γ=0.7)

    train!(regressor, X_train, y_train)

    @test predict(regressor, X_test) ≈ y_test_discounted
end

@testset "Generalized Regression Discriminator" begin
    X_train = Bool[
        1 1 0 0 0 0
        1 1 1 0 0 0
        1 1 1 1 0 0
        1 1 1 1 1 0
        1 1 1 1 1 1
    ]

    y_train = Float64[1, 2, 3, 4, 5]

    X_test = Bool[
        0 0 0 0 0 0
        1 0 0 0 0 0
        1 1 0 0 0 0
        1 1 1 1 1 1
        1 1 1 0 0 0
    ]

    # Fixed α
    y_test = Float64[5.25 / 4, 5.25 / 4, 13.3125 / 6, 20.8125 / 6, 14.3125 / 6]

    partitioner = LinearPartitioner(6, 2)

    regressor = Discriminator{LinearPartitioner,GeneralizedRegressionNode}(partitioner; α=0.5)

    train!(regressor, X_train, y_train)

    @test predict(regressor, X_test) == y_test
end

max_stars_repo_name: RafaelFK/ramnet

max_issues_repo_name: RafaelFK/ramnet

max_forks_repo_name: RafaelFK/ramnet


text:
"""
    centralized location for physical parameters.
"""

using Unitful

const k = 1.380649e-23 * u"m"^2 * u"kg" / u"s"^2 / u"K"

const Nₐ = 6.022140857e23 / u"mol"

# Gas constant in J mol^-1 K^-1 = m^2 kg mol^-1 s^-2 K^-1
# J = kg m^2 s^-2
const R = k * Nₐ

const e = 1.60217662e-19 * u"C"

# Relative permittivity of air at approx.
# 296.15 K and 101325 Pa and 40% RH
# See https://www.osti.gov/servlets/purl/1504063
# Previously known as the "dielectric constant"
# Often denoted as epsilon
const ε_air_room = 1.000530569 # unitless
# At STP (273.15 K, 1 atm):
# see: https://en.wikipedia.org/wiki/Relative_permittivity
const ε_air_stp = 1.00058986 # unitless

# select one of the two:
const ε = ε_air_room

# permittivity of free space in F/m
# Also known as the electric constant, permittivity of free space
# Often denoted by epsilon_0
const ε₀ = 8.85418782e-12 * u"F" / u"m"

# These values are used to calculate the dynamic viscosity of air
# Here, REF temperature and viscosity are at STP:
# Standard temperature and pressure (273.15 K and 101325 Pa)
const µ₀_air_stp = 1.716e-5 * u"Pa" * u"s"
const T₀_air_stp = 273.15 * u"K"

const Sₖ = 110.4 * u"K"
const M_air = 28.9644 * u"g" / u"mol" |> upreferred


text:
models = [
    rs.actM1(),
    rs.rtgM1(),
    rs.rtgM2(),
    rs.rtgM3(),
    rs.rtgM4(),
]

m4=rs.rtgM4()
m4rand = rs.rtgM4(m4; u=rs.init_u(m4))
rs.getSteady(m4rand)

@show rs.get_protein_lookup(models[end].model)

rs.getSteady.(models)
[rs.knockout(m, [:Rtg1]) for m in models]

@test rs.isValid(rs.rtgM4()) == true

valid, _ = rs.try_conditions(models[end])
@test valid == true


text:
using GUIFrameworkPlayground
using Base.Test

# write your own tests here
include(joinpath("..","examples","example01.jl"))


text:
module DataDragon # LOLToolsDataDragon

include(normpath(@__DIR__, "..", "gen", "locales.jl"))

end # module LOLToolsDataDragon.DataDragon