text:
using Random

@testset "Max det completion" begin
    # -- Test 1: matrix with known completion --
    # From VA Figure 4.2
    n = 17
    ijs = [(3,1), (3,2), (4,1), (4,2), (4,3), (5,1), (5,3), (5,4), (8,7),
           (9,5), (9,6), (9,7), (9,8), (11,10), (13,10), (13,11), (13,12),
           (14,10), (14,11), (14,12), (14,13),
           (15, 1), (15,3), (15,4), (15,5), (15,7), (15,8), (15,9),
           (16,5), (16,6), (16,9), (16,12), (16,13), (16,14), (16,15),
           (17,10), (17,11), (17,12), (17,13), (17,14), (17,15), (17,16)]
    append!(ijs, [(i,i) for i in 1:n])
    II, JJ = Chordal.unzip(ijs)
    A = sparse(II, JJ, ones(length(II)))
    A = (A + A') / 2
    sp = sparsity_pattern(A)

    # True solution
    #=
    X = Semidefinite(n)
    p = maximize(logdet(X), [X[II[x],JJ[x]] == A[II[x],JJ[x]] for x in 1:length(II)])
    Convex.solve!(p, () -> Mosek.Optimizer())
    Xstar = X.value
    =#
    Xstar = [1.0 0.3333370598035662 0.5 0.5 0.5 0.22222287873041388 0.27778132246744036 0.2777811910239065 0.33333530246961796 0.1712975423472041 0.17129709225708564 0.20370126720855525 0.2037019856666862 0.20370171663098346 0.5 0.33333088023046814 0.27777638852035724; 0.3333370598035662 1.0 0.5 0.5 0.33333305125898244 0.14814707689070564 0.1851867428728771 0.1851869926336214 0.22222250153752335 0.11419886412281514 0.11419870887817024 0.135800056328276 0.13580093153877904 0.13580062045920008 0.33333343444912716 0.22221930265929504 0.18518415224568624; 0.5 0.5 1.0 0.5 0.5 0.22222196163744942 0.27778076043784233 0.2777809794496988 0.33333475836820464 0.1712991354296111 0.1712988996359925 0.2037013235821753 0.2037028777637407 0.2037025691566877 0.5 0.33333017624316663 0.2777768088531805; 0.5 0.5 0.5 1.0 0.5 0.2222219730839787 0.27778081718375297 0.27778099169886683 0.3333347721005158 0.171299355067919 0.1712991188490877 0.20370132886164274 0.20370290806545718 0.20370259032235122 0.5 0.3333303500560772 0.27777693617839855; 0.5 0.33333305125898244 0.5 0.5 1.0 0.33333571187398225 0.33333659387305287 0.33333587588587565 0.5 0.2222299016007586 0.222228389459611 0.27777554232420304 0.2777775941596431 0.27777738847461536 0.5 0.5 0.3333356606726846; 0.22222287873041388 0.14814707689070564 0.22222196163744942 0.2222219730839787 0.33333571187398225 1.0 0.2777807210172411 0.27777952955889934 0.5 0.19907461348157693 0.1990723384116564 0.25925704288495866 0.25925593341309217 0.2592559179558187 0.3333352130700529 0.5 0.2777786944527559; 0.27778132246744036 0.1851867428728771 0.27778076043784233 0.27778081718375297 0.33333659387305287 0.2777807210172411 1.0 0.5 0.5 0.1712997337016747 0.17129908722574652 0.2037044400165765 0.2037058562144702 0.20370577454105926 0.5 0.3333317798520487 0.2777766860348415; 0.2777811910239065 0.1851869926336214 0.2777809794496988 0.27778099169886683 0.33333587588587565 0.27777952955889934 0.5 1.0 0.5 0.17129896845287737 0.1712983380154316 0.20370409144796173 0.2037052452025856 0.2037051577978674 0.5 0.3333308462064361 0.27777608804056464; 0.33333530246961796 0.22222250153752335 0.33333475836820464 0.3333347721005158 0.5 0.5 0.5 0.5 1.0 0.22222720796785625 0.22222501887881926 0.2777758212711955 0.2777770394826086 0.27777696102252825 0.5 0.5 0.3333344233472463; 0.1712975423472041 0.11419886412281514 0.1712991354296111 0.171299355067919 0.2222299016007586 0.19907461348157693 0.1712997337016747 0.17129896845287737 0.22222720796785625 1.0 0.5 0.3750025210974436 0.5 0.5 0.2916717850653076 0.37500249564668897 0.5; 0.17129709225708564 0.11419870887817024 0.1712988996359925 0.1712991188490877 0.222228389459611 0.1990723384116564 0.17129908722574652 0.1712983380154316 0.22222501887881926 0.5 1.0 0.3750032717610717 0.5 0.5 0.29167218644030213 0.3750017383599886 0.5; 0.20370126720855525 0.135800056328276 0.2037013235821753 0.20370132886164274 0.27777554232420304 0.25925704288495866 0.2037044400165765 0.20370409144796173 0.2777758212711955 0.3750025210974436 0.3750032717610717 1.0 0.5 0.5 0.33333681129144077 0.5 0.5; 0.2037019856666862 0.13580093153877904 0.2037028777637407 0.20370290806545718 0.2777775941596431 0.25925593341309217 0.2037058562144702 0.2037052452025856 0.2777770394826086 0.5 0.5 0.5 1.0 0.5 0.333338503723327 0.5 0.5; 0.20370171663098346 0.13580062045920008 0.2037025691566877 0.20370259032235122 0.27777738847461536 0.2592559179558187 0.20370577454105926 0.2037051577978674 0.27777696102252825 0.5 0.5 0.5 0.5 1.0 0.3333385764072932 0.5 0.5; 0.5 0.33333343444912716 0.5 0.5 0.5 0.3333352130700529 0.5 0.5 0.5 0.2916717850653076 0.29167218644030213 0.33333681129144077 0.333338503723327 0.3333385764072932 1.0 0.5 0.5; 0.33333088023046814 0.22221930265929504 0.33333017624316663 0.3333303500560772 0.5 0.5 0.3333317798520487 0.3333308462064361 0.5 0.37500249564668897 0.3750017383599886 0.5 0.5 0.5 0.5 1.0 0.5; 0.27777638852035724 0.18518415224568624 0.2777768088531805 0.27777693617839855 0.3333356606726846 0.2777786944527559 0.2777766860348415 0.27777608804056464 0.3333344233472463 0.5 0.5 0.5 0.5 0.5 0.5 0.5 1.0]
    pstar = det(Xstar)

    W = CD.maxdet_completion(A)
    @test ≈(A[CartesianIndex.(ijs)], W[CartesianIndex.(ijs)])
    @test det(W) ≈ pstar

    L, D = Chordal.maxdet_completion_etree(A)
    W = inv(D) * inv(LowerTriangular(Matrix(L)))
    W = L' \ W
    @test ≈(A[CartesianIndex.(ijs)], W[CartesianIndex.(ijs)])
    @test det(W) ≈ pstar

    L, D = Chordal.maxdet_completion_factors(A)
    W = inv(Matrix(D)) * inv(LowerTriangular(Matrix(L)))
    W = L' \ W
    @test ≈(A[CartesianIndex.(ijs)], W[CartesianIndex.(ijs)])
    @test det(W) ≈ pstar

    # -- Test 2: ensure pattern matches on random matrices --
    n = 50
    r = 5
    Random.seed!(0)
    M = rand(n,r)
    M = M*M'

    sp = spzeros(n, n)
    for i in 1:n
        block_size = randn() < 1.2 ? 2 : 10
        sp[i:min(i+block_size,n), i:min(i+block_size,n)] .= 1
    end

    M_ = sparse(sp .* M)
    M_comp = CD.maxdet_completion(M_)
    @test ≈(M_, sp.*M_comp, atol=1e-6)

    L, D = CD.maxdet_completion_etree(M_)
    W = inv(D) * inv(LowerTriangular(Matrix(L)))
    W = L' \ W
    @test ≈(M_, sp.*W, atol=1e-6)

    # L, D = Chordal.maxdet_completion_factors(M_)
    # W = inv(Matrix(D)) * inv(LowerTriangular(Matrix(L)))
    # W = L' \ W
    # @test ≈(M_, sp.*W, atol=1e-6)
end

@testset "Min rank completion" begin
   ## Build a test matrix
    # From Xin Jiang's masters thesis
    n = 9
    nonzero_ind_val = [
        (1,1,923.0),
        (1,4,464.0),
        (1,8,1215.0),
        (2,2,1721.0),
        (2,3,1494.0),
        (2,5,1182.0),
        (2,8,1638.0),
        (3,3,1300.0),
        (3,4,726.0),
        (3,5,1016.0),
        (3,8,1432.0),
        (4,4,418.0),
        (4,5,580.0),
        (4,8,770.0),
        (5,5,921.0),
        (5,8,1209.0),
        (5,9,959.0),
        (6,6,677.0),
        (6,7,421.0),
        (6,9,746.0),
        (7,7,270.0),
        (7,8,661.0),
        (7,9,499.0),
        (8,8,1773.0),
        (8,9,1195.0),
        (9,9,1021.0)
    ]
    sp = sparse(CD.unzip(nonzero_ind_val)...)
    sp = sp + triu(sp, 1)'
    s_pat = sparsity_pattern([sp])

    Y = CD.minrank_completion(sp)
    A_comp = Y*Y'
    matched = s_pat .* A_comp
    @test matched ≈ sp

    function bigger_mrc_test(n=50, r=4)
        function sparsity_structure(n)
            sp = spzeros(n, n)
            for i in 1:n
                block_size = randn() < 1 ? 2 : 10
                sp[i:min(i+block_size,n), i:min(i+block_size,n)] .= 1
            end
            return Int.(sp)
        end

        M = rand(n,r)
        M = M*M'

        sp = sparsity_structure(n)

        M_ = sparse(sp .* M)
        U = CD.minrank_completion(M_)
        M_comp = U*U'
        err = abs.(sp.*M_comp - M_)
        rmse = norm(err) / sqrt(length(M_))
        # @show sp
        # println(sp[end-4:end, end-4:end] == ones(Int, 5,5), rmse > 1e-10)
        return rmse
    end

    for N in [10, 20, 50, 200], r in [2, 5, 10, 20]
        r >= N && continue
        @test maximum([bigger_mrc_test(N, r) for i in 1:10]) < sqrt(eps())
    end

end


text:
export iauLtpequ
"""
**  Long-term precession of the equator.
**
**  This function is part of the International Astronomical Union's
**  SOFA (Standards of Fundamental Astronomy) software collection.
**
**  Status:  support function.
**
**  Given:
**     epj     double         Julian epoch (TT)
**
**  Returned:
**     veq     double[3]      equator pole unit vector
**
**  Notes:
**
**  1) The returned vector is with respect to the J2000.0 mean equator
**     and equinox.
**
**  2) The Vondrak et al. (2011, 2012) 400 millennia precession model
**     agrees with the IAU 2006 precession at J2000.0 and stays within
**     100 microarcseconds during the 20th and 21st centuries.  It is
**     accurate to a few arcseconds throughout the historical period,
**     worsening to a few tenths of a degree at the end of the
**     +/- 200,000 year time span.
**
**  References:
**
**    Vondrak, J., Capitaine, N. and Wallace, P., 2011, New precession
**    expressions, valid for long time intervals, Astron.Astrophys. 534,
**    A22
**
**    Vondrak, J., Capitaine, N. and Wallace, P., 2012, New precession
**    expressions, valid for long time intervals (Corrigendum),
**    Astron.Astrophys. 541, C1
**
**  This revision:  2016 February 9
**
**  SOFA release 2018-01-30
**
**  Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauLtpequ(double epj, double veq[3])
function iauLtpequ(epj::Real)
   veq = zeros(Float64, 3)
   ccall((:iauLtpequ, libsofa_c), Cvoid,
         (Cdouble, Ref{Cdouble}),
          convert(Float64, epj),
          veq)

   return SVector{3}(veq)
end

ext: jl


text:
using StudentHousing, JuMP, SDDP, SDDiP, StatsBase, Distributions

srand(32)

# ==============================================================================
# One stage, demo model.
# ==============================================================================

# Since we don't have real data, we'll generate some synthetic houses.

# Express the universe of all houses we can choose from, to make synthetic data:

# Number of bedrooms
nbedrooms_range = collect(1:2)
nbedrooms_frequency = Weights([0.5, 0.5])
# Number of bathrooms
nbathrooms_range = collect(1:2)
nbathrooms_frequency = Weights([0.5, 0.5])
# The price is at most...
prices_range_pp = [800.0, 1000.0]
# The area is at least ...
area_ranges = [0.0]
# Cache the data
market_data = StudentHousing.MarketData(nbedrooms_range, nbedrooms_frequency,
    nbathrooms_range, nbathrooms_frequency, prices_range_pp, area_ranges)

# The budget we have this period:
budget = 1e6

# Generate a structure holding all of our problem data:
problem_data = StudentHousingData(market_data, nhouses = 2, budget = budget, demand_distribution = Uniform())
# We chose to create 2 houses to choose from.

# Have a look at some of the things we stored
println(problem_data.houses)

# julia> problem_data.houses
# 2-element Array{StudentHousing.House,1}:
# StudentHousing.House(2, 1, 634.796, 614.355)
# StudentHousing.House(1, 1, 557.96, 311.02)

# We have also cached all the possible combinations of characteristics
println(problem_data.all_characteristics)
# julia> problem_data.all_characteristics
# 6-element Array{StudentHousing.Characteristic,1}:
#  StudentHousing.Characteristic(1, 1, 1, 1)
#  StudentHousing.Characteristic(1, 1, 2, 1)
#  StudentHousing.Characteristic(2, 1, 1, 1)
#  StudentHousing.Characteristic(2, 1, 2, 1)
#  StudentHousing.Characteristic(2, 2, 1, 1)
#  StudentHousing.Characteristic(2, 2, 2, 1)

# Our houses could be described as
# house 1: 2 bedrooms, 1 bathroom, <800$,  >= 0 sqrt ft (fits characteristic 3)
# house 1: 2 bedrooms, 1 bathroom, <1000$, >= 0 sqrt ft (fits characteristic 4)
# house 2: 1 bedrooms, 1 bathroom, <800$,  >= 0 sqrt ft (fits characteristic 1)
# house 2: 1 bedrooms, 1 bathroom, <1000$, >= 0 sqrt ft (fits characteristic 2)

# Note that the second house fits the descriptions of the first and second set
# of characteristics. Some preference patterns will allow the first only, some
# will allow the first and second. We don't do anything special to prioritize
# one assignment over another at the moment (e.g. assign lowest rent when
# possible). (TODO)

# If we permitted any combination of approving/disapproving a characteristic,
# we would have 2^6 = 64 preference patterns. However we have cached information
# only about preference patterns that make logical sense.

# julia> println(length(problem_data.patterns_allow))
# 13
# So 13 patterns don't allow a person to be OK with characteristic x but not OK
# with characteristic y that dominates x (better in every single rankable
# feature).

# Build the model and solve
m_one_stage = onestagemodel(problem_data)
@assert solve(m_one_stage) == :Optimal

# Have a look at our solutions:
println(getvalue(m_one_stage[:investment]))
# [1.0, 1.0]
println(getvalue(m_one_stage[:assignment]))
# [1, 6] = 1.0 <------------------
# [1, 7] = 0.0
# [1,11] = 1.0 <------------------
# [1,12] = 0.0
# [1,13] = 0.0
# [2, 3] = 1.0 <------------------
# [2, 4] = 0.0
# [2, 5] = 0.0
# [2, 6] = -0.0
# [2, 7] = 0.0
# [2, 8] = 0.0
# [2, 9] = 0.0
# [2,10] = 0.0
# [2,11] = -0.0
# [2,12] = 0.0
# [2,13] = 0.0

# We have chosen to invest in both houses and assigned the first to one
# entity with preference pattern 6, and one with preference pattern 1.
# We assigned the second house to someone with preference pattern 3.

# Sanity check:
p = problem_data.patterns_allow[6]
# 3-element Array{Int64,1}:
#  1
#  3
#  5
p = problem_data.patterns_allow[11]
# 4-element Array{Int64,1}:
#  1
#  2
#  3
#  5
# ... shows us that preference pattern 6 allows characteristics set 3, and
# preference pattern 11 allows characteristic 3 (which house 1 fits).
p = problem_data.patterns_allow[3]
# 1-element Array{Int64,1}:
#  1
# ... shows us that preference pattern 3 allows characteristics set 1.

# ==============================================================================
# One stage, larger model.
# ==============================================================================
# Let's see how much we can scale up without any large-scale techniques,
# and investigate how much we lose by solving a linear relaxation of our
# problem.

# To scale up, we'll just increase the number of bedrooms allowed for now.

function bedrooms_scale(i::Int)
    # Number of possible bedrooms = 1 -> i
    nbedrooms_range = collect(1:i)
    # Some synthetic sampling probabilities
    weights = 2.0.^(-collect(1:i))
    weights[1] = 1 - sum(weights[2:end])
    nbedrooms_frequency = Weights(weights)
    # Return the bedroom data
    nbedrooms_range, nbedrooms_frequency
end

nbathrooms_range = collect(1:2)
nbathrooms_frequency = Weights([0.75; 0.25])
prices_range_pp = [800.0, 1000.0]
area_ranges = [0.0, 700.0, 850.0]
budget = 1000.0

for i = 1:2
    srand(1)
    nbedrooms_range, nbedrooms_frequency = bedrooms_scale(i)
    @time market_data = StudentHousing.MarketData(nbedrooms_range,
        nbedrooms_frequency, nbathrooms_range, nbathrooms_frequency,
        prices_range_pp, area_ranges)
    @time problem_data = StudentHousingData(market_data, nhouses = 50, budget = budget,
        demand_distribution = Uniform(0.0, 100.0))
    P = StudentHousing.get_npatterns(problem_data)
    println("Using $P patterns, solving MIP:")
    @time m_one_stage = onestagemodel(problem_data)
    # println("Assigned students = ", getvalue(sum(sum(m_one_stage[:assignment]))))
    tic()
    @assert solve(m_one_stage) == :Optimal
    toc()
    println("MIP objective = ", getobjectivevalue(m_one_stage))
    @assert solve(m_one_stage, relaxation = true) == :Optimal
    println("LP objective = ", getobjectivevalue(m_one_stage))
end

# 0.000004 seconds (1 allocation: 64 bytes)
# 0.000155 seconds (327 allocations: 319.641 KiB)
# Using 9 patterns, solving MIP:
# 0.000444 seconds (3.66 k allocations: 209.766 KiB)
# elapsed time: 0.016997735 seconds
# MIP objective = 474.0
# LP objective = 473.5585096078103
# 0.000005 seconds (1 allocation: 64 bytes)
# 1.222919 seconds (1.05 M allocations: 1.084 GiB, 25.72% gc time)
# Using 329 patterns, solving MIP:
# 0.011981 seconds (78.52 k allocations: 4.183 MiB)
# elapsed time: 0.045735724 seconds
# MIP objective = 16616.0
# LP objective = 16615.10611513201

# 0.000017 seconds (1 allocation: 64 bytes)
# too long outside of MIP on PC, 2^48 potential patterns to process

# Even when 2^ncharacteristics is too large to be stored, the
# number of legal patterns we get is only in the hundrends, and the MIP takes
# milliseconds to solve. The bottleneck is processing the patterns before we
# even get to the MIP.

# Let's focus on the multistage problem for a bit.

# ==============================================================================
# Multistage, medium model.
# ==============================================================================
srand(1)
nbedrooms_range = collect(1:2)
nbedrooms_frequency = Weights([0.75, 0.25])

market_data = StudentHousing.MarketData(nbedrooms_range,
    nbedrooms_frequency, nbathrooms_range, nbathrooms_frequency,
    prices_range_pp, area_ranges)

problem_data = StudentHousingData(market_data, nhouses = 50, budget = budget,
    demand_distribution = Uniform(0.0, 100.0), nstages = 10, nnoises = 10)

m_multi_stage = multistagemodel(problem_data)
solve(m_multi_stage, max_iterations = 5)

# converges extremely quickly... and Benders vs Lagrangian bounds v close
# IP: 164.660K LP: 164.592K

# r = simulate(m_multi_stage, 50, [:investment, :assignment, :shortage])
# julia> mean(r[i][:objective] for i = 1:50)
# 164633.34


text:
using RHEOS
using Test

const tol = (eps(RHEOS.RheoFloat))^(0.125)

println("===============================================")
println("Testing RHEOS")
println("===============================================")
println("|")

include("symbols.jl")
include("definitions.jl")
include("base.jl")
include("datagen.jl")
include("IO.jl")
include("processing.jl")
include("interface.jl")

include(joinpath(@__DIR__, "models", "elements.jl"))
include(joinpath(@__DIR__, "models", "maxwell.jl"))
include(joinpath(@__DIR__, "models", "kelvinvoigt.jl"))
include(joinpath(@__DIR__, "models", "zener.jl"))
include(joinpath(@__DIR__, "models", "poynting-thomson.jl"))


text:
#Ost - "E6" - восток, Nord - "N8" - север, West - "W4" - запад, Sud - "S2" - юг
function main_function8!(r)
    side = Ost
    while isborder(r,Nord)==true
        putmarker!(r)
        go_mark(r,side)
        side=inverse(side)
    end
end

function go_mark(r, side)
    while ismarker(r)==true
        move!(r,side)
    end
end
function inverse(side)
    return HorizonSide(mod(Int(side)+2, 4))
end


text:
const GLE_DEFAULTS = LittleDict{Symbol,Any}(
    :palette => [ # imitated from tableau 10 - 2
        RGB(0.33, 0.47, 0.64),
        RGB(0.90, 0.57, 0.26),
        RGB(0.82, 0.37, 0.36),
        RGB(0.51, 0.70, 0.69),
        RGB(0.42, 0.62, 0.35),
        RGB(0.91, 0.79, 0.37),
        RGB(0.66, 0.49, 0.62),
        RGB(0.95, 0.63, 0.66),
        RGB(0.59, 0.46, 0.38),
        RGB(0.72, 0.69, 0.67)
    ],
    #
    :textstyle_font  => "texcmss",
    :textstyle_hei   => 0.35,
    :textstyle_color => nothing,
    #
    :linestyle_lstyle => 0,
    :linestyle_lwidth => nothing,
    :linestyle_smooth => nothing,
    :linestyle_color  => nothing,
    #
    :markerstyle_marker => nothing,
    :markerstyle_msize  => nothing,
    :markerstyle_color  => nothing,
    #
    :barstyle_color => nothing,
    :barstyle_fill  => c"white",
    #
    :fillstyle_fill   => c"cornflowerblue",
    #
    :boxplotstyle_bwidth         => 0.6,
    :boxplotstyle_wwidth         => 0.3,
    :boxplotstyle_wrlength       => 1.5,
    :boxplotstyle_blstyle_lstyle => 1,
    :boxplotstyle_blstyle_lwidth => 0.0,
    :boxplotstyle_blstyle_color  => c"seagreen",
    :boxplotstyle_mlstyle_lstyle => 1,
    :boxplotstyle_mlstyle_lwidth => 0.0,
    :boxplotstyle_mlstyle_color  => c"seagreen",
    :boxplotstyle_mmstyle_marker => "fdiamond",
    :boxplotstyle_mmstyle_msize  => 0.4,
    :boxplotstyle_mmstyle_color  => c"dodgerblue",
    :boxplotstyle_omstyle_marker => "+",
    :boxplotstyle_omstyle_msize  => 0.5,
    :boxplotstyle_omstyle_color  => c"tomato",
    :boxplotstyle_mshow          => true,
    :boxplotstyle_oshow          => true,
    #
    :text2d_position  => "cc",
    :text2d_textstyle => :default_TextStyle,
    #
    :straightline2d_linestyle => :default_LineStyle,
    #
    :box2d_position  => "bl",
    :box2d_nobox     => true,
    :box2d_linestyle => :default_LineStyle,
    :box2d_fillstyle => nothing,
    #
    :colorbar_size     => nothing,
    :colorbar_pixels   => 100,
    :colorbar_nobox    => true,
    :colorbar_position => "right",
    :colorbar_offset   => (0.3, 0.0),
    #
    :title_textstyle => :default_TextStyle,
    :title_dist      => nothing,
    #
    :tickslabels_textstyle => :default_TextStyle,
    :tickslabels_angle     => nothing,
    :tickslabels_format    => nothing,
    :tickslabels_shift     => nothing,
    :tickslabels_dist      => nothing,
    :tickslabels_off       => false,
    #
    :ticks_linestyle => :default_LineStyle,
    :ticks_length    => nothing,
    :ticks_symticks  => false,
    :ticks_off       => false,
    :ticks_grid      => false,
    #
    :legend_position  => nothing,
    :legend_textstyle => :default_TextStyle,
    :legend_offset    => (0.0, 0.0),
    :legend_bgcolor   => nothing,
    :legend_margins   => nothing,
    :legend_nobox     => false,
    :legend_off       => false,
    #
    :axis_ticks     => :default_Ticks,
    :axis_textstyle => :default_TextStyle,
    :axis_title     => nothing,
    :axis_base      => nothing,
    :axis_lwidth    => nothing,
    :axis_min       => nothing,
    :axis_max       => nothing,
    :axis_off       => false,
    :axis_log       => false,
    #
    :axes2d_title  => nothing,
    :axes2d_size   => nothing,
    :axes2d_legend => nothing,
    :axes2d_origin => nothing,
    :axes2d_math   => false,
    :axes2d_off    => false,
    :axes2d_scale  => "auto",
    #
    :scatter2d_linestyles_lstyle => 0,
    :scatter2d_linestyles_lwidth => 0.05,
    :scatter2d_linestyles_smooth => true,
    #
    :hist2d_horiz    => false,
    :hist2d_bins     => nothing,
    :hist2d_scaling  => "none",
    :hist2d_label    => "",
    #
    :bar2d_stacked => false,
    :bar2d_horiz   => false,
    :bar2d_bwidth  => nothing,
    :bar2d_labels  => String[],
    #
    :boxplot_horiz => false,
    #
    :heatmap_cmap      => colormap("RdBu", 10),
    :heatmap_cmiss     => c"white",
    :heatmap_transpose => false,
    #
    :scatter3d_linestyle => :default_LineStyle,
    :scatter3d_markerstyle => :default_MarkerStyle,
    :scatter3d_label => "",
)

"""
    default(T)

Default constructor for an object.
"""
function default(::Type{T}, a...) where T
    tname = T.name.name
    return T(a...,
        (
            begin
                s = Symbol(lowercase("$(tname)_$(fn)"))
                v = get(GLE_DEFAULTS, s, nothing)
                if v isa Symbol
                    sv = string(v)
                    if startswith(sv, "default")
                        T2 = Symbol(split(sv, '_')[2])
                        v  = eval(quote
                            $T2()
                        end)
                    end
                end
                v
            end
            for fn in fieldnames(T)[length(a)+1:end]
        )...
    )
end


text:
module ThieleFourPointInterpolation

export thiele

include("twopoints.jl")
include("threepoints.jl")
include("fourpoints.jl")

end # module


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

# This file contains average properties of the human eye and optical human eye models

# luminance (cd/m2)	Multiple	Value	Item
# 10−6	µcd/m2	1 µcd/m2	Absolute threshold of vision[1]
# 10−5
# 10−4		400 µcd/m2	Darkest sky[2]
# 10−3	mcd/m2	1 mcd/m2	Night sky[3]
# 1.4   mcd/m2	Typical photographic scene lit by full moon[4]
# 5     mcd/m2	Approximate scotopic/mesopic threshold[5]
# 10−2		40 mcd/m2	Phosphorescent markings on a watch dial after 1 h in the dark[6][7]
# 10−1
# 100	cd/m2	2 cd/m2	Floodlit buildings, monuments, and fountains[8]
# 5     cd/m2	Approximate mesopic/photopic threshold[5]
# 101		25 cd/m2	Typical photographic scene at sunrise or sunset[4]
# 30    cd/m2	Green electroluminescent source[2]
# 102		250 cd/m2	Peak luminance of a typical LCD monitor[10][11]
# 700   cd/m2	Typical photographic scene on overcast day[4][8][11]
# 103	kcd/m2	2 kcd/m2	Average cloudy sky[2]
# 5     kcd/m2	Typical photographic scene in full sunlight[4][8]

"""
# Pupil diameter as a function of scene luminance
https://jov.arvojournals.org/article.aspx?articleid=2279420
https://en.wikipedia.org/wiki/Orders_of_magnitude_(luminance)

Pupil diameter is approximately 2.8mm at 100cd/m^2. A typical overcast day is 700cd/m^2
"""

"""computes pupil diameter as a function of scene luminance `L`, in cd/m², and the angular area, `a`, over which this luminance is presented to the eye."""
𝐃sd(L,a) = 7.75 - 5.75 * ((L * a / 846)^.41) / ((L * a / 846)^.41 + 2) # the first letter of this function name is \bfD not D.
export 𝐃sd

eyeradius() = 12mm
export eyeradius

"""Posterior focal length, i.e., optical distance from entrance pupil to the retina. Focal length will change depending on accomodation. This value is for focus at ∞. When the eye is focused at 25cm focal length will be ≈ 22mm. Because the index of refraction of the vitreous humor is approximately 1.33 the physical distance from the entrance pupil to the retina will be 24mm/1.33 = 18mm."""
eyefocallength() = 24mm
export eyefocallength

vc_epupil() = 3mm #distance from vertex of cornea to entrance pupil

""" distance from vertex of cornea to center of rotation"""
cornea_to_eyecenter() = 13.5mm
export cornea_to_eyecenter

entrancepupil_to_retina() = 22.9mm

"""distance from entrance pupil to center of rotation."""
entrancepupil_to_eyecenter() = entrancepupil_to_retina() - eyeradius()
export entrancepupil_to_eyecenter

"""average angle, in degrees, the eye will rotate before users will turn their head"""
comfortable_eye_rotation_angle() = 20°

"""Average one sided translation of the entrance pupil associated with comfortable eye rotation. If you are using this to define an eyebox multiply this value by 2"""
comfortable_entrance_pupil_translation() = sin(comfortable_eye_rotation_angle())*entrancepupil_to_eyecenter()
export comfortable_entrance_pupil_translation


text:
using Test

using MarkovChains.Graphs
@testset "digraph" begin
    @testset "search_dfs" begin
        dg = Digraph()
        n1 = add_node!(dg)
        n2 = add_node!(dg)
        n3 = add_node!(dg)
        n4 = add_node!(dg)
        add_edge!(dg, n1, n2)
        add_edge!(dg, n2, n3)
        add_edge!(dg, n1, n4)
        order = []
        dfs(dg, n1, (p, c) -> begin
            push!(order, c)
        end)
        @test order == [1, 2, 3, 4]
    end
    @testset "scc_one_comp" begin
        dg = Digraph()
        n1 = add_node!(dg)
        n2 = add_node!(dg)
        n3 = add_node!(dg)
        add_edge!(dg, n1, n2)
        add_edge!(dg, n2, n3)
        add_edge!(dg, n3, n1)
        comps = strongly_connected_components(dg)
        @test length(comps) == 1
        @test comps[1].members == [3,2,1]
        @test comps[1].is_bottom == true
    end
    @testset "scc_acyclic" begin
        dg = Digraph()
        n1 = add_node!(dg)
        n2 = add_node!(dg)
        n3 = add_node!(dg)
        add_edge!(dg, n1, n2)
        add_edge!(dg, n2, n3)
        comps = strongly_connected_components(dg)
        @test length(comps) == 3
        @test comps[1].members == [3]
        @test comps[1].is_bottom == true
        @test comps[2].members == [2]
        @test comps[2].is_bottom == false
        @test comps[3].members == [1]
        @test comps[3].is_bottom == false
    end
    @testset "scc_two_comp" begin
        dg = Digraph()
        n1 = add_node!(dg)
        n2 = add_node!(dg)
        n3 = add_node!(dg)
        n4 = add_node!(dg)
        add_edge!(dg, n1, n2)
        add_edge!(dg, n2, n1)
        add_edge!(dg, n2, n3)
        add_edge!(dg, n3, n4)
        add_edge!(dg, n4, n3)
        comps = strongly_connected_components(dg)
        @test length(comps) == 2
        @test comps[1].members == [4, 3]
        @test comps[1].is_bottom == true
        @test comps[2].members == [2, 1]
        @test comps[2].is_bottom == false
    end
end


text:
export Sector, SectoralPrism

#= Sector (for a base shape) =#

mutable struct Sector{D,T} <: Shape{2,4,D,T}  # M = 2K
    c::SVector{2,T}  # center of circle
    r::T  # radius of circle
    ϕ₀::T  # center angle bisecting sector: -π ≤ ϕ₀ < π  (π excluded)
    ∆ϕ2::T  # "radius" in angle dimension: 0 ≤ ∆ϕ2 ≤ π (sector spans from ϕ₀ - ∆ϕ2 to ϕ₀ + ∆ϕ2)
    data::D  # auxiliary data
    Sector{D,T}(c,r,ϕ₀,∆ϕ,data) where {D,T<:Real} = new(c,r,ϕ₀,∆ϕ,data)  # suppress default outer constructor
end

function Sector(c::SVector{2,<:Real}, r::Real, ϕ::Real, ∆ϕ::Real, data::D=nothing) where {D}
    r≥0 || throw(ArgumentError("r = $r must be nonnegative."))
    -2π≤∆ϕ≤2π  || throw(ArgumentError("∆ϕ = $∆ϕ must be between -2π and 2π, inclusive."))

    ϕ₀ = rem(ϕ + ∆ϕ/2, 2π, RoundNearest)  # put ϕ₀ in [-π, π)
    ∆ϕ2 = abs(∆ϕ/2)

    return Sector{D}(c, r, ϕ₀, ∆ϕ2, data)
end

Sector(c::AbstractVector{<:Real},  # center of circle
       r::Real,  # radius of circle
       ϕ::Real,  # start angle
       ∆ϕ::Real,  # width in angle dimension; can be negative (-2π ≤ ∆ϕ ≤ 2π)
       data=nothing) =
    Sector(SVector{2}(c), r, ϕ, ∆ϕ, data)

Base.:(==)(s1::Sector, s2::Sector) = s1.c==s2.c && s1.r==s2.r && s1.ϕ₀==s2.ϕ₀ && s1.∆ϕ2==s2.∆ϕ2 && s1.data==s2.data
Base.isapprox(s1::Sector, s2::Sector) = s1.c≈s2.c && s1.r≈s2.r && s1.ϕ₀≈s2.ϕ₀ && s1.∆ϕ2≈s2.∆ϕ2 && s1.data==s2.data
Base.hash(s::Sector, h::UInt) = hash(s.c, hash(s.r, hash(s.ϕ₀, hash(s.∆ϕ2, hash(s.data, hash(:Sector, h))))))

distangle(ϕ::Real, ϕ₀:: Real) = rem(ϕ-ϕ₀, 2π, RoundNearest)  # ϕ measured from ϕ₀; result within [-π, π)

function Base.in(x::SVector{2,<:Real}, s::Sector)
    d = x - s.c
    ld = norm(d)
    ϕ = ld==0 ? s.ϕ₀ : atan(d[2], d[1])  # angle to x with respect to c

    return ld ≤ s.r && abs(distangle(ϕ, s.ϕ₀)) ≤ s.∆ϕ2
end

function surfpt_nearby(x::SVector{2,<:Real}, s::Sector{D,T}) where {D,T<:Real}
    # Basically mimic the same function for Prism, but proceeds in the (ρ,ϕ) domain.
    d = x - s.c
    ld = norm(d)

    # Calculate the closest point in the ρ dimension and outward normal direciton there.
    r2 = s.r / 2
    ρ = ld - r2  # positive if closer to arc; negative if closer to center
    d̂ = ld ≤ Base.rtoldefault(T) * r2  ? SVector(cos(s.ϕ₀),sin(s.ϕ₀)) : normalize(d)

    surfρ = ρ<0 ? 0.0 : s.r  # scalar: closest point to x between center and perimeter point
    noutρ = copysign(1.0,ρ) * d̂  # SVector{2}: outward direction normal at surfρ

    absρ = abs(ρ)
    abs∆ρ = abs(r2 - absρ)  # radial distance between x and either center or perimeter, whichever closer to x

    onbndρ = abs∆ρ ≤ Base.rtoldefault(T) * r2  # basically r2 ≈ ρ but faster
    isoutρ = (r2 < absρ) || onbndρ

    # Calculate the closest point in the ϕ dimension and outward normal direciton there.
    ϕ = ld==0 ? 0.0 : distangle(atan(d[2], d[1]), s.ϕ₀)  # positive if closer to end side; negative if closer to start side

    ϕtemp = ϕ<0 ? s.ϕ₀ - s.∆ϕ2 : s.ϕ₀ + s.∆ϕ2
    surfϕ = SVector(cos(ϕtemp), sin(ϕtemp))  # SVector{2}: closest point in ϕ dimension
    cosθ = d̂⋅surfϕ  # cosine of angle between d and surfϕ
    ldcosθ = ld*cosθ  # ldconθ .* surfϕ is surface point
    ldsinθ = ld*sqrt(1-cosθ^2)  # always positive

    ϕtemp += copysign(π/2,ϕ)
    noutϕ = SVector(cos(ϕtemp), sin(ϕtemp))  # SVector{2}: outward direction normal at surfϕ

    absϕ = abs(ϕ)
    abs∆ϕ = abs(s.∆ϕ2 - absϕ)  # angular distance between x and closer side of sector

    onbndϕ = abs∆ϕ ≤ Base.rtoldefault(T) * s.∆ϕ2  # basically ∆ϕ2 ≈ ϕ but faster
    isoutϕ = (s.∆ϕ2 < absϕ) || onbndϕ

    # Pick the surface point and outward direction normal depending on the location of x.
    if isoutρ && isoutϕ  # x outside in both ρ and ϕ dimensions
        surf = surfρ .* surfϕ  # one of two end points of arc
        nout = (onbndρ && onbndϕ) ? (noutρ + noutϕ) : (d - surf)
        nout = normalize(nout)
    elseif !isoutρ && isoutϕ  # x inside in ρ dimension, but outside in ϕ dimension
        (surf, nout) = (absϕ < s.∆ϕ2+π/2) ? (ldcosθ .* surfϕ, noutϕ) : (SVector(0.0,0.0), d̂)
    elseif isoutρ && !isoutϕ  # x outside in ρ dimension, but inside in ϕ dimension
        (surf, nout) = (surfρ .* d̂, noutρ)
    else  # !isoutρ && !isoutϕ: x strictly inside sector
        (surf, nout) = (abs∆ρ ≤ ldsinθ) ? (surfρ .* d̂, noutρ) : (ldcosθ .* surfϕ, noutϕ)
    end

    return surf+s.c, nout
end

translate(s::Sector{D}, ∆::SVector{2,T}) where {D,T<:Real} = Sector{D,T}(s.c+∆, s.r, s.ϕ₀, s.∆ϕ2, s.data)

function bounds(s::Sector)
    # Find the minimum and maximum coordinates among the center and two ends of the arc.
    ϕ = SVector(s.ϕ₀ - s.∆ϕ2, s.ϕ₀ + s.∆ϕ2)  # [start angle, end angle]
    v = s.r .* [SVector(0.0,0.0) [cos.(ϕ) sin.(ϕ)]']  # [center, start point of arc, end point of arc]

    # Consider using the code below once https://github.com/JuliaArrays/StaticArrays.jl/issues/498
    # is resolved:
    # l = minimum(v, dims=Val(2))[:,1]
    # u = maximum(v, dims=Val(2))[:,1]

    xs = v[1,:]
    ys = v[2,:]

    xmin = minimum(xs)
    xmax = maximum(xs)

    ymin = minimum(ys)
    ymax = maximum(ys)

    # Consider the extreme points on the arc in the Cartesian directions.
    abs(distangle(0, s.ϕ₀)) ≤ s.∆ϕ2 && (xmax = max(xmax,s.r))  # sector contains +x-direction from center
    abs(distangle(π/2, s.ϕ₀)) ≤ s.∆ϕ2 && (ymax = max(ymax, s.r))  # sector contains +y-direction from center
    abs(distangle(π, s.ϕ₀)) ≤ s.∆ϕ2  && (xmin = min(xmin, -s.r))  # sector contains -x-direction from center
    abs(distangle(3π/2, s.ϕ₀)) ≤ s.∆ϕ2  && (ymin = min(ymin, -s.r))  # sector contains -y-direction from center

    return (SVector(xmin,ymin)+s.c, SVector(xmax,ymax)+s.c)
end

#= Sectoral prism =#
const SectoralPrism = Prism{Sector{Nothing}}

# Below, if we called SectoralPrism(c, ...) in the function body, it would call the inner
# constructor Prism{Sector{Nothing}}(c, ...) because SectoralPrism = Prism{Sector{Nothing}},
# which is not what we want.
# To call the outer constructor of Prism, we should call Prism(c, ...) instead of SectoralPrism(c, ...).
SectoralPrism(c::SVector{3,<:Real},
              r::Real,
              ϕ::Real,
              ∆ϕ::Real,
              h::Real=Inf,
              a::SVector{3,<:Real}=SVector(0.0,0.0,1.0),
              data=nothing) where {K} =
    (â = normalize(a); Prism(c, Sector(SVector(0.0,0.0),r,ϕ,∆ϕ), h, [orthoaxes(â)... â], data))

SectoralPrism(c::AbstractVector{<:Real},  # center of prism
              r::Real,  # radius of sectoral base
              ϕ::Real,  # start angle of sectoral base: 0 ≤ ϕₛ < 2π  (2π excluded)
              ∆ϕ::Real,  # end angle sectoral base: 0 ≤ ϕₑ-ϕₛ ≤ 2π
              h::Real=Inf,  # height of prism
              a::AbstractVector{<:Real}=[0.0,0.0,1.0],  # axis direction of prism
              data=nothing) =
    SectoralPrism(SVector{3}(c), r, ϕ, ∆ϕ, h, SVector{3}(a), data)

function bounds_ctrcut(s::SectoralPrism)
    ax = inv(s.p)  # prism axes: columns are not only unit vectors, but also orthogonal
    if ax ≈ I  # prism axes are aligned with Cartesian directions (this covers most usage)
        l, u = bounds(s.b)  # (SVector{2}, SVector{2})
        a₁₂ = ax[:,SVector(1,2)]  # SMatrix{3,2}

        return a₁₂*l, a₁₂*u
    else  # prism axes are not aligned with Cartesian directions
        b = s.b
        r = b.r

        el = Ellipsoid(SVector(0.0,0.0,0.0), SVector(r,r,0.0), ax)  # center is set at origin to return bounds with respect to prism center
        bp = boundpts(el)  # SMatrix{3,3}: boundary points; all three columns of b are free of NaN because no column of ax is aligned with Cartesian directions
        bp′ = s.p * bp  # SMatrix{3,3}: bp in prism coordinates
        bp2′ = bp′[SVector(1,2),:]  # SMatrix{2,3}: in each column of bp′, third entry is in axis dimension, so must be zero mathematically

        ϕ = atan.(bp2′[2,:], bp2′[1,:])  # SVector{3}: angles of boundary points in base plane
        ϕsym = [ϕ; ϕ.+π]  # SVector{6}: include symmetric points with respect to base center
        ϕall = [ϕsym; SVector(b.ϕ₀-b.∆ϕ2, b.ϕ₀+b.∆ϕ2)]  # SVector{8}: angles of all boundary point candidates, except base center

        bpall′ = r .* [cos.(ϕall) sin.(ϕall) @SVector(zeros(8))]'  # SMatrix{3,8}: each column is boundary point candidate in prism coordinates
        bpall = ax * bpall′  # SMatrix{3,8}: each column is boundary point candidate in external coordinates
        bpallc = [bpall SVector(0.0,0.0,0.0)]  # SMatrix{3,9}: include base center
        ind = abs.(distangle.(ϕsym, b.ϕ₀)) .≤ b.∆ϕ2  # SVector{6}: indices of angles contained in base arc
        indallc = [ind; SVector(true,true,true)]  # SVector{9}: include arc ends and base center

        xs = bpallc[1,:]
        ys = bpallc[2,:]
        zs = bpallc[3,:]

        # It might be possible to implement the following more efficiently using reduce().
        xmin = ymin = zmin = Inf
        xmax = ymax = zmax = -Inf
        for i = 1:9
            if indallc[i]
                xmin = min(xmin, xs[i])
                ymin = min(ymin, ys[i])
                zmin = min(zmin, zs[i])

                xmax = max(xmax, xs[i])
                ymax = max(ymax, ys[i])
                zmax = max(zmax, zs[i])
            end
        end

        return (SVector(xmin,ymin,zmin), SVector(xmax,ymax,zmax))
    end
end


text:
# FunctionGraph

@testset "ex-fg" begin
    b = board("brd", shownavigation=true, showscreenshot=true)
    b ++ slider("a", [[1,-1],[5,-1],[0,1.5,3]])
    b ++ @jsf foo(x) = val(a)*x^2
    b ++ functiongraph(foo, dash=2)
    s = str(b, preamble=false)
    @test isapproxstr(s, """(function(){
        function foo(x){return (val(a)*Math.pow(x,2))};
        var brd = JXG.JSXGraph.initBoard("brd",
                {"boundingbox":[-10,10,10,-10],
                 "axis":false,
                 "showcopyright":false,
                 "shownavigation":true,
                 "showscreenshot":true});
        var a = brd.create(
                'slider',
                [[1.0,-1.0],[5.0,-1.0],[0.0,1.5,3.0]],
                {});
        brd.create(
                'functiongraph',
                [function(t){return foo(t);}, -10, 10],
                {"dash":2});})();
        """)
end

# ParametricCurve
@testset "ex-pc" begin
    b = board("brd", xlim=[-1, 15], ylim=[-0.5, 2.5])
    @jsf f1(t) = t - sin(t)
    @jsf f2(t) = 1 - cos(t)
    sT = slider("T", [[0,2.1],[3,2.1],[0,π,5π]])
    @jsf fb() = val(T)
    @jsf pa() = f1(val(T))
    @jsf pb() = f2(val(T))
    b ++ (f1, f2, sT, fb, pa, pb)
    b ++ plot(f1, f2; a=0, b=fb, dash=2)
    b ++ point(pa, pb, withlabel=false)
    s = str(b, preamble=false)
    @test isapproxstr(s, """(function(){
        function f1(t){return (t-sin(t))};
        function f2(t){return (1-cos(t))};
        function fb(){return val(T)};
        function pa(){return f1(val(T))};
        function pb(){return f2(val(T))};
        var brd=JXG.JSXGraph.initBoard("brd",
                {"boundingbox":[-1.0,2.5,15.0,-0.5],
                "axis":false,
                "showcopyright":false,
                "shownavigation":false});
        var T=brd.create(
                'slider',
                [[0.0,2.1],[3.0,2.1],[0.0,3.141592653589793,15.707963267948966]], {});
        brd.create(
                'curve',
                [function(t){return f1(t);}, function(t){return f2(t);},
                  0, function(t){return fb(t);}],
                {"dash":2});
        brd.create(
                'point',
                [function(t){return pa(t);}, function(t){return pb(t);}],
                {"withlabel":false});})();
        """)
end

# DataPlot
@testset "ex-dp" begin
    b = board("b", xlim=[0, 1], ylim=[0, 1])
    x = [1,2,3,4]
    y = [1,2,3,4]
    b ++ plot(x, y)
    s = str(b, preamble=false)
    @test isapproxstr(s, """(function(){
        var b=JXG.JSXGraph.initBoard("b"
                ,
                {"boundingbox":[0,1,1,0],
                "axis":false,
                "showcopyright":false,
                "shownavigation":false});
        b.create('curve', [[1,2,3,4], [1,2,3,4]], {});})();
        """)
end

# Lissajou
@testset "ex-li" begin
    b = board("brd", xlim=[-12, 12], ylim=[-10,10])
    b ++ (
        slider("a", [[2,8],[6,8],[0,3,6]]),
        slider("b", [[2,7],[6,7],[0,2,6]]),
        slider("A", [[2,6],[6,6],[0,3,6]]),
        slider("B", [[2,5],[6,5],[0,3,6]]),
        slider("delta", [[2,4],[6,4],[0,0,π]], name="δ")
        )
    @jsf f1(t) = val(A)*sin(val(a)*t+val(delta))
    @jsf f2(t) = val(B)*sin(val(b)*t)
    b ++ plot(f1, f2, a=0, b=2π, strokecolor="#aa2233", strokewidth=3)
    s = str(b, preamble=false)
    @test isapproxstr(s, """(function(){
        var brd=JXG.JSXGraph.initBoard("brd"
                ,
                {"boundingbox":[-12,10,12,-10],
                "axis":false,
                "showcopyright":false,
                "shownavigation":false});
        var a=brd.create('slider', [[2,8],[6,8],[0,3,6]], {});
        var b=brd.create('slider', [[2,7],[6,7],[0,2,6]], {});
        var A=brd.create('slider', [[2,6],[6,6],[0,3,6]], {});
        var B=brd.create('slider', [[2,5],[6,5],[0,3,6]], {});
        var delta=brd.create('slider', [[2.0,4.0],[6.0,4.0],[0.0,0.0,3.141592653589793]], {"name":"δ"});
        function f1(t){return (val(A)*sin(((val(a)*t)+val(delta))))};
        function f2(t){return (val(B)*sin((val(b)*t)))};
        brd.create(
                'curve',
                [function(t){return f1(t);}, function(t){return f2(t);},
                 0, 6.283185307179586],
                {"strokecolor":"#aa2233",
                "strokewidth":3});})();
        """)
end


text:
for (func, elty) in [(:CLBlastShad, Float32), (:CLBlastDhad, Float64),
                     (:CLBlastChad, ComplexF32), (:CLBlastZhad, ComplexF64)]
    #TODO: (:CLBlastHhad, Float16)

    @eval function $func(n::Integer, alpha::$elty,
                         x_buffer::cl.CL_mem, x_offset::Integer, x_inc::Integer,
                         y_buffer::cl.CL_mem, y_offset::Integer, y_inc::Integer,
                         beta::$elty,
                         z_buffer::cl.CL_mem, z_offset::Integer, z_inc::Integer,
                         queue::cl.CmdQueue, event::cl.Event)
        err = ccall(
            ($(string(func)), libCLBlast),
            cl.CL_int,
            (Csize_t, $elty, Ptr{Cvoid}, Csize_t, Csize_t, Ptr{Cvoid}, Csize_t, Csize_t, $elty, Ptr{Cvoid}, Csize_t, Csize_t, Ptr{Cvoid}, Ptr{Cvoid}),
            n, alpha, x_buffer, x_offset, x_inc, y_buffer, y_offset, y_inc, beta, z_buffer, z_offset, z_inc, Ref(queue), Ref(event)
        )
        if err != cl.CL_SUCCESS
            println(stderr, "Calling function $(string($func)) failed!")
            throw(cl.CLError(err))
        end
        return err
    end

    @eval function had!(n::Integer, α::Number,
                        x::cl.CLArray{$elty}, x_inc::Integer,
                        y::cl.CLArray{$elty}, y_inc::Integer,
                        β::Number,
                        z::cl.CLArray{$elty}, z_inc::Integer;
                        queue::cl.CmdQueue=cl.queue(z))
        # output event
        event::cl.Event = cl.Event(C_NULL)
        alpha = convert($elty, α)
        beta = convert($elty, β)

        $func(Csize_t(n), alpha,
              pointer(x), Csize_t(0), Csize_t(x_inc),
              pointer(y), Csize_t(0), Csize_t(y_inc),
              beta,
              pointer(z), Csize_t(0), Csize_t(z_inc),
              queue, event)

        # wait for kernel
        cl.wait(event)

        z
    end

end


text:
using RiskAdjustedLinearizations, Test, JLD2
using SparseArrays, SparseDiffTools, FiniteDiff
include(joinpath(dirname(@__FILE__), "..", "..", "examples", "wachter_disaster_risk", "wachter.jl"))

# Load in guesses and true solutions
detout = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "reference", "det_ss_output.jld2"), "r")

# Set up RiskAdjustedLinearization
m = WachterDisasterRisk()
ral = inplace_wachter_disaster_risk(m) # has guesses for z and y already
z = copy(ral.z)
y = copy(ral.y)
Ψ = copy(ral.Ψ)

# Solve!
@test_throws AssertionError solve!(ral, ral.z, ral.y, ral.Ψ; algorithm = :deterministic, verbose = :high, autodiff = :central) # first w/finite diff nlsolve Jacobian
@info "The following series of print statements are expected."
ral.z .= z .* 1.001
ral.y .= y .* 1.001
solve!(ral, ral.z, ral.y; algorithm = :deterministic, verbose = :high, autodiff = :central) # first w/finite diff nlsolve Jacobian
@test maximum(abs.(ral.z - detout["z"])) < 1e-6
@test maximum(abs.(ral.y - detout["y"])) < 1e-6
ral.z .= z .* 1.001
ral.y .= y .* 1.001
RiskAdjustedLinearizations.deterministic_steadystate!(ral, vcat(ral.z, ral.y);
                                                      verbose = :none, autodiff = :central) # first w/finite diff nlsolve Jacobian
@test maximum(abs.(ral.z - detout["z"])) < 1e-6
@test maximum(abs.(ral.y - detout["y"])) < 1e-6

update!(ral, z, y, Ψ) # now autodiff Jacobian
ral.z .= vec(detout["z"]) * 1.001
ral.y .= vec(detout["y"]) * 1.001
solve!(ral, ral.z, ral.y; algorithm = :deterministic, verbose = :high, autodiff = :forward) # now autodiff nlsolve Jacobian
@test maximum(abs.(ral.z - detout["z"])) < 1e-6
@test maximum(abs.(ral.y - detout["y"])) < 1e-6
ral.z .= vec(detout["z"]) * 1.001
ral.y .= vec(detout["y"]) * 1.001
RiskAdjustedLinearizations.deterministic_steadystate!(ral, vcat(ral.z, ral.y);
                                                      verbose = :none, autodiff = :forward)
@test maximum(abs.(ral.z - detout["z"])) < 1e-6
@test maximum(abs.(ral.y - detout["y"])) < 1e-6

# Check sparse Jacobian
sparsity, colorvec = compute_sparsity_pattern(ral, :deterministic; sparsity_detection = false)
jac_cache = preallocate_jac_cache(ral, :deterministic; sparsity_detection = false)
RiskAdjustedLinearizations.deterministic_steadystate!(ral, 1.001 .* vcat(ral.z, ral.y);
                                                      verbose = :none,
                                                      sparse_jacobian = true,
                                                      sparsity = sparsity,
                                                      colorvec = colorvec)
@test maximum(abs.(ral.z - detout["z"])) < 1e-6
@test maximum(abs.(ral.y - detout["y"])) < 1e-6
RiskAdjustedLinearizations.deterministic_steadystate!(ral, 1.001 .* vcat(ral.z, ral.y);
                                                      verbose = :none,
                                                      autodiff = :forward,
                                                      sparse_jacobian = true,
                                                      sparsity = sparsity,
                                                      colorvec = colorvec)
@test maximum(abs.(ral.z - detout["z"])) < 1e-6
@test maximum(abs.(ral.y - detout["y"])) < 1e-6
RiskAdjustedLinearizations.deterministic_steadystate!(ral, 1.001 .* vcat(ral.z, ral.y);
                                                      verbose = :none,
                                                      sparse_jacobian = true,
                                                      jac_cache = jac_cache)
@test maximum(abs.(ral.z - detout["z"])) < 1e-6
@test maximum(abs.(ral.y - detout["y"])) < 1e-6
RiskAdjustedLinearizations.deterministic_steadystate!(ral, 1.001 .* vcat(ral.z, ral.y);
                                                      verbose = :none,
                                                      sparse_jacobian = true,
                                                      sparsity_detection = false)
@test maximum(abs.(ral.z - detout["z"])) < 1e-6
@test maximum(abs.(ral.y - detout["y"])) < 1e-6

# Using SparsityDetection.jl fails
@test_broken RiskAdjustedLinearizations.deterministic_steadystate!(ral, 1.001 .* vcat(ral.z, ral.y);
                                                                   verbose = :none,
                                                                   sparse_jacobian = true,
                                                                   sparsity_detection = true)


text:
# test IOftuart type
using FTD2XX
using Base.Test

# assume at least ont device is attached
deviceinfo = ft_getdeviceinfolist()[1]

# open by description
io = open(FT_Description(deviceinfo.description))
close(io)

# open by serial number
io = open(FT_SerialNumber(deviceinfo.serialnumber))
close(io)

# open by location
io = open(FT_Location(deviceinfo.locid))
close(io)

# open by index (zero based)
uartconfig1 = UARTConfiguration(921600,8,1,"n",readtimeout = 1000,writetimeout = 1000)
uartconfig2 = UARTConfiguration(9600,8,1,"n",readtimeout = 100,writetimeout = 100)
io = open(FT_DeviceIndex(0),uartconfig1)

wait(Timer(0.1))  # why do I need this?

# assume external loopback is connected
for i in 1:50
  write(io, 0x55)
  @test read(io, UInt8) == 0x55

  write(io, "hello world")
  @test readuntil(io, "d") == "hello world"  # carefull, eof == true always

  testvector = [0x44, 0x67, 0x88, 0x01, 0xab]
  byteswritten = write(io, testvector)
  @test byteswritten == length(testvector)
  readbuffer = zeros(UInt8,length(testvector))
  read!(io, readbuffer)
  @test testvector == readbuffer
  print(i%10)
end
println()

close(io)


text:
# Functions for creating batch files to run on NIH Biowulf super computer

using DelimitedFiles

transientgenelist = ["MYC", "DUSP5", "TRIB1", "PMAIP1", "SERPINE1", "SOX9", "ERRFI1", "NR4A2", "JUN", "GEM"]

function makeswarm(G::Int,infolder,swarmfile::String,inlabel,label,nsets,datafolder::String;thresholdlow::Float64=0.,thresholdhigh::Float64=100000000.,conds::String="DMSO-AUXIN",result::String= "2021-03-11",batchsize=1000,maxtime = 3600. * 2,nchains::Int = 2,runcycle::Bool=true,transient::Bool=false,fittedparam::String="",fixedeffects::String="",juliafile::String="/home/carsonc/StochasticGene/runfiles/fitscript.jl",root="/home/carsonc/scrna/")
    if occursin.("-",conds)
        cond = string.(split(conds,"-"))[1]
    else
        cond = conds
    end
    genes = checkgenes(root,cond,datafolder,thresholdlow,thresholdhigh)
    makeswarm(genes,G,infolder,swarmfile,inlabel,label,nsets,datafolder,conds=conds,result=result,batchsize=batchsize,maxtime=maxtime,nchains=nchains,runcycle=runcycle,transient=transient,fittedparam=fittedparam,fixedeffects=fixedeffects,juliafile=juliafile,root=root)
end

function makeswarm(genes::Vector,G::Int,infolder,swarmfile::String,inlabel,label,nsets,datafolder::String;conds::String="DMSO-AUXIN",result::String="testout",batchsize=1000,maxtime=60.,nchains::Int=2,runcycle::Bool=true,transient::Bool=false,fittedparam::String="",fixedeffects="",juliafile::String="/home/carsonc/StochasticGene/runfiles/fitscript.jl",root="/home/carsonc/scrna/")

    resultfolder = joinpath("Results",result)
    infolder = joinpath("Results",infolder)
    ngenes = length(genes)
    println(ngenes)
    println(runcycle)

    if ngenes > batchsize
        nbatches = div(ngenes,batchsize)
        batches = Vector{Vector{String}}(undef,nbatches+1)
        println(batchsize," ",nbatches)
        for i in 1:nbatches
            batches[i] = genes[batchsize*(i-1)+1:batchsize*(i)]
        end
        batches[end] = genes[batchsize*nbatches+1:end]
        for batch in eachindex(batches)
            sfile = swarmfile * "_" * conds * "$batch" * ".swarm"
            writegenes(sfile,batches[batch],nchains,juliafile,conds,G,maxtime,infolder,resultfolder,datafolder,inlabel,label,nsets,runcycle,transient,nchains,fittedparam,fixedeffects)
        end
    else
        sfile = swarmfile * "_" * conds * ".swarm"
        f = open(sfile,"w")
        writegenes(sfile,genes,nchains,juliafile,conds,G,maxtime,infolder,resultfolder,datafolder,inlabel,label,nsets,runcycle,transient,nchains,fittedparam,fixedeffects)
    end
end

function getbatches(genes,ngenes,batchsize)
    nbatches = div(ngenes,batchsize)
    batches = Vector{Vector{String}}(undef,nbatches+1)
    println(batchsize," ",nbatches)
    for i in 1:nbatches
        batches[i] = genes[batchsize*(i-1)+1:batchsize*(i)]
    end
    batches[end] = genes[batchsize*nbatches+1:end]
    return batches
end

function writegenes(sfile,genes,nprocs,juliafile,cond,G,maxtime,infolder,resultfolder,datafolder,inlabel,label,nsets,runcycle,transient,nchains,fittedparam,fixedeffects)
    f = open(sfile,"w")
    for gene in genes
        writedlm(f,["julia -p" nprocs juliafile gene cond G maxtime infolder resultfolder datafolder inlabel label nsets runcycle transient nchains fittedparam fixedeffects])
    end
    close(f)
end

# function checkgenes(cond,datafolder,thresholdlow,thresholdhigh,root = "/Users/carsonc/Box/scrna",halflives="data/HCT116_all_cells_histograms_and_half_lives_March_2021/Starved_Rep7_half_lives.csv")
#     datapath = joinpath(root,datafolder)
#     # genes = readdlm(joinpath(root,"data/HCT_HEK_Half_life_top_below_1h.csv"),',')[2:end,end]
#     hlfile = joinpath(root,halflives)
#     contents = readdlm(joinpath(root,"data/HCT116_all_cells_histograms_and_half_lives_March_2021/Starved_Rep7_half_lives.csv"),',')
#     genedecay = contents[2:end,2]
#     genenames = contents[2:end,1]
#     genes1 = Vector{String}(undef,0)
#     for i in eachindex(genedecay)
#         if typeof(genedecay[i]) <: Number
#             if thresholdlow <= genedecay[i] < thresholdhigh
#                 push!(genes1,genenames[i])
#             end
#         end
#     end
#     println(length(genes1))
#     genes2 = readdlm(joinpath(root,"data/HCT116_alleles_number.txt"))[2:end,1]
#     println(length(genes2))
#     genes = intersect(genes1,genes2)
#     complement = setdiff(genes1,genes2)
#     println(length(genes))
#     println(length(complement))
#     list = Vector{String}(undef,0)
#     for gene in genes
#         if cond == ""
#             datafile = joinpath(datapath,gene * ".txt")
#         else
#             datafile = joinpath(datapath,gene * "_" * cond * ".txt")
#         end
#         if isfile(datafile)
#             push!(list,gene)
#         end
#     end
#     println(length(list))
#     return list
# end
#
#
# function checkgenes(cond,datafolder,root = "/Users/carsonc/Box/scrna")
#     datapath = joinpath(root,datafolder)
#     # genes = readdlm(joinpath(root,"data/HCT_HEK_Half_life_top_below_1h.csv"),',')[2:end,end]
#     genedecay = readdlm(joinpath(root,"data/HCT116_all_cells_histograms_and_half_lives_March_2021/Starved_Rep7_half_lives.csv"),',')
#     # genes1 = readdlm(joinpath(root,"data/HCT116_all_cells_histograms_and_half_lives_March_2021/Starved_Rep7_half_lives.csv"),',')[2:end,1]
#     genes1 = Vector{String}(undef,0)
#     for i in eachindex(genedecay[2:end,2])
#         if typeof(genedecay[i,2]) <: Number
#             push!(genes1,genedecay[i,1])
#         end
#     end
#     genes2 = readdlm(joinpath(root,"data/HCT116_alleles_number.txt"))[2:end,1]
#     println(length(genes1))
#     println(length(genes2))
#     genes = intersect(genes1,genes2)
#     list = Vector{String}(undef,0)
#     for gene in genes
#         datafile = joinpath(datapath,gene * "_" * cond * ".txt")
#         if isfile(datafile)
#             push!(list,gene)
#         end
#     end
#     println(length(list))
#     return list
# end
#
# function checkgenes(datafolder,root)
#     datapath = joinpath(root,datafolder)
#     genedecay = readdlm(joinpath(root,"data/HCT116_all_cells_histograms_and_half_lives_March_2021/Starved_Rep7_half_lives.csv"),',')
#     genes1 = Vector{String}(undef,0)
#     for i in eachindex(genedecay[2:end,2])
#         if typeof(genedecay[i,2]) <: Number
#             push!(genes1,genedecay[i,1])
#         end
#     end
#     genes2 = readdlm(joinpath(root,"data/HCT116_alleles_number.txt"))[2:end,1]
#     println(length(genes1))
#     println(length(genes2))
#     genes = intersect(genes1,genes2)
#     list = Vector{String}(undef,0)
#     for gene in genes
#         datafile = joinpath(datapath,gene * ".txt")
#         if isfile(datafile)
#             push!(list,gene)
#         end
#     end
#     println(length(list))
#     return list
# end

function checkgenes(root,cond,datafolder,thresholdlow::Float64=0.,thresholdhigh::Float64=1e8,hlfile="data/HCT116_all_cells_histograms_and_half_lives_March_2021/Starved_Rep7_half_lives.csv",allelefile="data/HCT116_alleles_number.txt")
    genes = intersect(get_halflives(root,hlfile,thresholdlow,thresholdhigh), get_genes(root,cond,datafolder))
    if allelefile != ""
        return intersect(genes,get_alleles(root,allelefile))
    else
        return genes
    end
end

function get_genes(root,cond,datafolder)
    genes = Vector{String}(undef,0)
    files = readdir(joinpath(root,datafolder))
    for file in files
        if occursin(cond,file)
            push!(genes,split(file,"_")[1])
        end
    end
    return genes
end

function get_halflives(root,file,thresholdlow=0,thresholdhigh=1e8)
    genes = Vector{String}(undef,0)
    file = joinpath(root,file)
    halflives = readdlm(file,',')
    for row in eachrow(halflives)
        if typeof(row[2]) <: Number
            if thresholdlow <= row[2] < thresholdhigh
                push!(genes,string(row[1]))
            end
        end
    end
    return genes
end

get_alleles(root,allelefile = "data/HCT116_alleles_number.txt") = readdlm(joinpath(root,allelefile))[2:end,1]

fix(folder) = writeruns(fixruns(findjobs(folder)))

function findjobs(folder)
    files = readdir(folder)
    files = files[occursin.("swarm_",files)]
    for (i,file) in enumerate(files)
        files[i] = split(file,"_")[2]
    end
    unique(files)
end

function fixruns(jobs,message="FAILED")
    runlist = Vector{String}(undef,0)
    for job in jobs
        if occursin(message,read(`jobhist $job`,String))
            swarmfile = findswarm(job)
            list = readdlm(swarmfile,',')
            runs =  chomp(read(pipeline(`jobhist $job`, `grep $message`),String))
            runs = split(runs,'\n')
            println(job)
            for run in runs
                linenumber = parse(Int,split(split(run," ")[1],"_")[2])
                a = String(list[linenumber+1])
                println(a)
                push!(runlist,a)
            end
        end
    end
    return runlist
end

function writeruns(runs,outfile="fitfix.swarm")

    f = open(outfile,"w")
    for run in runs
        writedlm(f,[run],quotes=false)
    end
    close(f)

end

function findswarm(job)
    sc = "Swarm Command"
    line = read(pipeline(`jobhist $job`, `grep $sc`),String)
    list = split(line," ")
    list[occursin.(".swarm",list)][1]
end

function collate(folders,type="rates_scRNA_T120")
    genes = Vector{Vector}(undef,0)
    for folder in folders
        files = readdir(folder)
        for file in files
            if occursin(type,file)
                contents=readdlm(file,',')
                for row in eachrow(contents)
                    row[1] .== genes
                end
            end
        end
    end
end

large_deviance(measurefile,threshold) = filter_gene(measurefile,"Deviance",threshold)

function filter_gene(measurefile,measure,threshold)
    genes = Vector{String}(undef,0)
    measures,header = readdlm(measurefile,',',header=true)
    println(length(measures[:,1]))
    col = findfirst(header[1,:] .== measure)
    for d in eachrow(measures)
        if d[col] > threshold || isnan(d[col])
            push!(genes,d[1])
        end
    end
    println(length(genes))
    return genes
end

function filter_gene_nan(measurefile,measure)
    genes = Vector{String}(undef,0)
    measures,header = readdlm(measurefile,',',header=true)
    println(length(measures[:,1]))
    col = findfirst(header[1,:] .== measure)
    for d in eachrow(measures)
        if isnan(d[col])
            push!(genes,d[1])
        end
    end
    println(length(genes))
    return genes
end

function get_missing_genes(folder,type,label,cond,model)
    genes = checkgenes("DMSO",folder,0.,100000000.,"../")
    genes1=StochasticGene.getgenes(folder,type,label,cond,model)
    union(setdiff(genes1,genes),setdiff(genes,genes1))
end

function get_missing_genes(genes,folder,type,label,cond,model)
    genes1=StochasticGene.getgenes(folder,type,label,cond,model)
    union(setdiff(genes1,genes),setdiff(genes,genes1))
end

function scan_swarmfiles(jobid,folder=".")
    if ~(typeof(jobid) <: String)
        jobid = string(jobid)
    end
    genes = Array{String,1}(undef,0)
    files = readdir(folder)
    files = files[occursin.(jobid,files)]
    for file in files
        genes = vcat(genes,scan_swarmfile(file))
    end
    return genes
end

function scan_swarmfile(file)
    genes = Array{String,1}(undef,0)
    contents = readdlm(file,'\t')
    lines = contents[occursin.("[\"",string.(contents))]
    for line in lines
        push!(genes,split.(string(line)," ")[1])
    end
    return genes
end

function scan_fitfile(file,folder=".")
    genes = Array{String,1}(undef,0)
    joinpath(folder,file)
    file = readdlm(file,'\t')
    for line in eachrow(file)
        push!(genes,line[4])
    end
    return genes
end

function replace_yield(G,folder1,folder2,cond1,cond2,outfolder)
    if typeof(G) <: Number
        G = string(G)
    end
    if ~isdir(outfolder)
        mkpath(outfolder)
    end
    files1 = getratefile(folder1,G,cond1)
    files2 = getratefile(folder2,G,cond2)
    for file1 in files1
        gene = StochasticGene.getgene(file1)
        file2 = getratefile(files2,gene)
        outfile = joinpath(outfolder,file2)
        r1 = StochasticGene.readrates(joinpath(folder1,file1))
        r2 = StochasticGene.readrates(joinpath(folder2,file2))
        r2[end] = r1[end]
        f = open(outfile,"w")
        writedlm(f,[r2],',')
        close(f)
    end

end

"""
assemble_r(G,folder1,folder2,cond1,cond2,outfolder)

Combine rates from two separate fits into a single rate vector

"""

function assemble_r(G,folder1,folder2,cond1,cond2,outfolder)
    if typeof(G) <: Number
        G = string(G)
    end
    if ~isdir(outfolder)
        mkpath(outfolder)
    end
    files1 = getratefile(folder1,G,cond1)
    files2 = getratefile(folder2,G,cond2)
    for file1 in files1
        gene = StochasticGene.getgene(file1)
        file2 = getratefile(files2,gene)
        if file2 != 0
            file2=joinpath(folder2,file2)
        else
            file2=joinpath(folder1,file1)
        end
        name = replace(file1, cond1 => cond1 * "-" * cond2)
        outfile = joinpath(outfolder,name)
        assemble_r(joinpath(folder1,file1),file2,outfile)
    end
end

function  assemble_r(ratefile1,ratefile2,outfile)
    r1 = StochasticGene.readrates(ratefile1,2)
    r2 = StochasticGene.readrates(ratefile2,2)
    r1[end] = clamp(r1[end],eps(Float64),1-eps(Float64))
    r = vcat(r1[1:end-1],r2[1:end-1],r1[end])
    f = open(outfile,"w")
    writedlm(f,[r],',')
    writedlm(f,[r],',')
    writedlm(f,[r],',')
    writedlm(f,[r],',')
    close(f)
end

function assemble_r(gene,G,folder1,folder2,cond1,cond2,outfolder)
    file1 = getratefile(gene,G,folder1,cond1)[1]
    file2 = getratefile(gene,G,folder2,cond2)[1]
    name = replace(file1, cond1 => cond1 * "-" * cond2)
    println(name)
    outfile = joinpath(outfolder,name)
    println(outfile)
    assemble_r(joinpath(folder1,file1),joinpath(folder2,file2),outfile)
end

function getratefile(files,gene)
    files = files[occursin.("_"*gene*"_",files)]
    if length(files) > 0
        return files[1]
    else
        # println(gene)
        return 0
    end
end

function getratefile(folder,G,cond)
    files = readdir(folder)
    files = files[occursin.("rates_",files)]
    files = files[occursin.("_"*cond*"_",files)]
    files[occursin.("_"*G*"_",files)]
end

getratefile(gene,G,folder,cond) = getfile("rate",gene,G,folder,cond)

function getfile(type,gene::String,G::String,folder,cond)
    files = readdir(folder)
    files = files[occursin.(type,files)]
    files = files[occursin.("_"*gene*"_",files)]
    files = files[occursin.("_"*G*"_",files)]
    files[occursin.("_"*cond*"_",files)]
end

function change_name(folder,oldname,newname)
    files = readdir(folder)
    files = files[occursin.(oldname,files)]
    for file in files
        newfile = replace(file, oldname => newname)
        mv(joinpath(folder,file),joinpath(folder,newfile),force=true)
    end
end

function make_halflife(infile,outfile,col=4)
    f = open(outfile,"w")
    writedlm(f,["Gene" "Halflife"],',')
    contents,rows = readdlm(infile,',',header=true)
    for row = eachrow(contents)
        gene = string(row[1])
        gene = strip(gene,'*')
        h1 = float(row[col])
        h2 = float(row[col+1])
        if h1 > 0 || h2 > 0
            h = (h1 + h2)/(float(h1>0) + float(h2>0))
            writedlm(f,[gene h],',')
        end
    end
    nothing
end

function make_datafiles(infolder,outfolder,label)
    histograms = readdir(infolder)
    if ~isdir(outfolder)
        mkpath(outfolder)
    end
    for file in histograms
        newfile = replace(file,label => "")
        cp(joinpath(infolder,file),joinpath(outfolder,newfile))
    end
end

max_stars_repo_name: nih-niddk-mbs/StochasticGene

max_issues_repo_name: nih-niddk-mbs/StochasticGene

max_forks_repo_name: nih-niddk-mbs/StochasticGene


text:
using Pkg;
using CSV, Plots;
using Statistics, SparseArrays, LinearAlgebra;
cd(@__DIR__) #src

include("func_slp.jl")

mutable struct inputcom
  y        :: Array
  x        :: Array
  w        :: Array
  H_min    :: Int
  H_max    :: Int
  type     :: String
  r        :: Int8
  λ        :: Int8
end

df = CSV.read("data.csv");
df = convert(Array{Float16}, df[:, 2:4]);
T, k = size(df);

## Parameterization
P = 4; # Number of lags used in LP - controlled variable

# start LP at H_min=0 or 1 (H_min=1 if impose no contemporanous impact)
H_max = 20;
ind_response = 1; # Endogenous variable - response
ind_shock    = 3; # Endogenous variable related to the shock

ind_lp = "reg"; # Type of linear regression

ind_contempo = filter(x -> x ≠ ind_shock, ind_all);
#,(P+1):end
y  = df[:, ind_response]; # endogenous variable
x  = df[:, ind_shock]; # endoegnous variable related to the shock

# control variables (contemporaneous vars, lagged vars)
w  = [ df[:, ind_contempo]  matlag( df , P ) ];
w[.!isfinite.(w)] .= 0;
# Packaing everything as input
r = 2; #(r-1)=order of the limit polynomial (so r=2 implies the IR is shrunk towards a line )
λ = 100; # value of penalty

testing₀ = inputcom(y, x, w, H_min, H_max, "reg", 3,4);
testing₁ = inputcom(y, x, w, H_min, H_max, "smooth", r, λ);

lp₀ = smoothlocalprojection(testing₀);

lp₁ = smoothlocalprojection(testing₁);
lp₁.δ
plot(lp₀.IR)
plot!(lp₁.IR)


text:
@hsa_kernel function vcopy(a,b)
    idx = get_global_id(Int32(0)) + 1

    b[idx] = a[idx]
    return nothing
end

@hsa_kernel function vadd(a,b,c)
    i = get_global_id(Int32(0)) + 1
    c[i] = a[i] + b[i]
    return nothing
end

function plus(a,b)
    return a + b
end

@hsa_kernel function vadd_fun(a,b,c)
    i = get_global_id(Int32(0)) + 1
    c[i] = plus(a[i], b[i])
    return nothing
end

@hsa_kernel function mmul(a,b,c,acols)
    # one kernel invocation per column of the result matrix
    arows = get_global_size(Int32(0))
    # i = col
    i = get_global_id(Int32(0))
    # j = row
    for j = 1:arows
        c_ij = 0.0
        for k = 1:acols
            a_idx = (k-1) * arows + j
            b_idx = i * acols + k
            c_ij += a[a_idx] * b[b_idx]
        end

        c_idx = i * arows + j
        c[c_idx] = c_ij
    end
    return nothing
end

@hsa_kernel function mmul2d(a,b,c,acols)
    # one kernel invocation per cell of the result matrix
    arows = get_global_size(Int32(0))
    # i = col
    i = get_global_id(Int32(1))
    # j = row
    j = get_global_id(Int32(0)) + 1

    c_ij = 0.0
    for k = 1:acols
        a_idx = (k-1) * arows + j
        b_idx = i * acols + k
        c_ij += a[a_idx] * b[b_idx]
    end

    c_idx = i * arows + j
    c[c_idx] = c_ij

    return nothing
end

@hsa_kernel function mmul_array(a,b,c)
    # one kernel invocation per cell of the result matrix
    # i = col
    i = get_global_id(Int32(0)) + 1
    # j = row
    j = get_global_id(Int32(1)) + 1

    c_ij = 0.0
    for k = 1:size(a)[1]
        c_ij += a[k,j] * b[k,i]
    end

    c[j,i] = c_ij

    return nothing
end


text:
"""
Delegates loading functions & traits for the NEID spectrograph

Author: Eric Ford and collaborators
Created: August 2020
"""

"""
Module providing types and traits and customized functions for the NEID Spectrograph.
- [NEID blog](https://neid.psu.edu/)
"""
module NEID
using RvSpectMLBase
import RvSpectMLBase: AbstractInstrument, AbstractInstrument1D, AbstractInstrument2D
import ..EchelleInstruments
import ..EchelleInstruments: read_manifest, read_metadata_from_fits, check_metadata_fields_expected_present
#import ..EchelleInstruments: read_header
import ..EchelleInstruments: default_Δv_to_avoid_tellurics

using CSV, DataFrames, Query
using FITSIO
using Interpolations
using NaNMath, Missings

""" Trait for 1D Extracted spectra from NEID """
struct NEID1D <: AbstractInstrument1D end

""" Trait for 2D Extracted spectra from NEID """
struct NEID2D <: AbstractInstrument2D end

""" Trait to specify any 1D or 2D Extracted spectra from NEID """
const AnyNEID = Union{NEID1D,NEID2D}

export NEID, NEID1D, NEID2D, AnyNEID

# traits.jl imports from RvSpectMBase on its own
include("traits.jl")
export min_order, max_order, min_pixel_in_order, max_pixel_in_order
export orders_to_use_default, min_col_default, max_col_default
export metadata_symbols_default, metadata_strings_default
export default_ccf_mask_v_width
export get_inst_module
export bad_col_ranges

include("io.jl")
# not exported, since don't have parameter that specializes to instrument

#import RvSpectMLBase: continuum_normalize_spectrum, continuum_normalize_spectra  # Not yet implemented
#import RvSpectMLBase: filter_line_list
#import RvSpectMLBase: find_worst_telluric_in_each_chunk # Not yet implemented
import RvSpectMLBase: make_clean_line_list_from_tellurics, choose_obs_idx_for_init_guess
include("util.jl") # TODO
#export continuum_normalize_spectrum!, continuum_normalize_spectra! # Not yet implemented
export filter_line_list
export make_clean_line_list_from_tellurics
export make_λ_list_for_bad_columns
#export choose_obs_idx_for_init_guess

end


text:
using SimpleSudoku
using Test

@testset "SimpleSudoku.jl" begin

    # New York Times 2022 02 14 Hard
    grid_str = """
            300007400
            000010083
            500020700
            980000006
            000005001
            010000020
            620070000
            009600000
            050400000"""

    grid = mapreduce(vcat, split(grid_str, "\n")) do str
        reshape(parse.(Int8, split(str, "")), 1, 9)
    end

    solved_grid, _ = solve_sudoku(grid)

    solved_grid
    #  3  6  1  9  8  7  4  5  2
    #  2  9  7  5  1  4  6  8  3
    #  5  4  8  3  2  6  7  1  9
    #  9  8  5  2  4  1  3  7  6
    #  7  3  2  8  6  5  9  4  1
    #  4  1  6  7  3  9  5  2  8
    #  6  2  4  1  7  3  8  9  5
    #  1  7  9  6  5  8  2  3  4
    #  8  5  3  4  9  2  1  6  7

end


text:
### Parsing utilities

include("lib/result.jl")

#=
A relic from the past

macro chk1(expr,label=:error)
    quote
        x = $(esc(expr))
        if isnull(x[1])
            @goto $label
        else
           x[1].value, x[2]
        end
    end
end
=#

macro chk2(expr,label=:error)
    @assert expr.head == :(=)
    lhs, rhs = expr.args

    @assert lhs.head == :tuple
    res, state = lhs.args
    quote
        x = $(esc(rhs))
        $(esc(state)) = x[2] # bubble error location
        if isnull(x[1])
            $(esc(:(@goto $label)))
        else
            $(esc(res)) = x[1].value
        end
    end
end

@inline _isdigit(c::Char) = isdigit(c)

@inline function parse_uint_and_stop(str, i, len, n::T) where {T <: Integer}
    ten = T(10)
    # specialize handling of the first digit so we can return an error
    max_without_overflow = div(typemax(T)-9,10) # the larg
    y1 = iterate(str, i)
    y1===nothing && return n, false, i
    c = y1[1]
    if _isdigit(c) && n <= max_without_overflow
        n *= ten
        n += T(c-'0')
    else
        return n, false, i
    end
    i = y1[2]

    y2 = iterate(str, i)
    while y2!==nothing && n <= max_without_overflow
        c = y2[1]
        if _isdigit(c)
            n *= ten
            n += T(c-'0')
        else
            return n, true, i
        end
        i = y2[2]

        y2 = iterate(str, i)
    end
    return n, true, i
end

# slurp up extra digits
@inline function read_digits(str, i, len)
    y = iterate(str, i)
    while y!==nothing
        c = y[1]
        if !_isdigit(c) # do nothing
            return i
        end
        i = y[2]
        y = iterate(str, i)
    end
    return i
end

@inline function tryparsenext_base10_digit(T,str,i, len)
    y = iterate(str,i)
    y===nothing && @goto error
    c = y[1]; ii = y[2]
    '0' <= c <= '9' || @goto error
    return convert(T, c-'0'), ii

    @label error
    return nothing
end

Base.@pure maxdigits(::Type{T}) where {T} = ndigits(typemax(T))
Base.@pure min_with_max_digits(::Type{T}) where {T} = convert(T, T(10)^(maxdigits(T)-1))

@inline function tryparsenext_base10(T, str,i,len)
    i0 = i
    R = Nullable{T}
    y = tryparsenext_base10_digit(T,str,i, len)
    y===nothing && return R(), i
    r = y[1]; i = y[2]

    # Eat zeros
    while r==0
        y2 = tryparsenext_base10_digit(T,str,i, len)
        y2 === nothing && return R(convert(T, 0)), i
        r = y2[1]; i = y2[2]
    end

    digits = 1
    ten = T(10)
    while true
        y2 = tryparsenext_base10_digit(T,str,i,len)
        y2===nothing && break
        digits += 1
        d = y2[1]; i = y2[2]
        r = r*ten + d
    end

    max_digits = maxdigits(T)

    # Checking for overflow
    if digits > max_digits
        # More digits than the max value we can hold, this is certainly
        # an overflow
        return R(), i0
    elseif digits == max_digits && r < min_with_max_digits(T)
        # Same digits as the max digits we can hold. If the number we computed
        # is now smaller than the smallest number with the same number of
        # digits as the typemax number, we must have overflown, so we
        # again return a parsing failure
        return R(), i0
    end

    return R(convert(T, r)), i
end

@inline function tryparsenext_sign(str, i, len)
    R = Nullable{Int}

    y = iterate(str, i)
    if y===nothing
        return return R(), i
    else
        c = y[1]; ii = y[2]
        if c == '-'
            return R(-1), ii
        elseif c == '+'
            return R(1), ii
        else
            return R(1), i
        end
    end
end

@inline function isspace(c::Char)
    c == ' '
end

@inline function isnewline(c::Char)
    c == '\n' || c == '\r'
end

@inline function eatwhitespaces(str, i=1, l=lastindex(str))
    y = iterate(str, i)
    while y!==nothing
        c = y[1]; ii = y[2]
        if isspace(c)
            i=ii
        else
            break
        end
        y = iterate(str, i)
    end
    return i
end

function eatnewlines(str, i=1, l=lastindex(str))
    count = 0
    y = iterate(str, i)
    while y!==nothing
        c = y[1]; ii = y[2]
        if c == '\r'
            i=ii
            y2 = iterate(str, i)
            if y2!==nothing
                c = y2[1]
                ii = y2[2]
                if c == '\n'
                    i=ii
                end
            end
            count += 1
        elseif c == '\n'
            i=ii
            y3 = iterate(str, i)
            if y3!==nothing
                c = y3[1]
                ii = y3[2]
                if c == '\r'
                    i=ii
                end
            end
            count += 1
        else
            break
        end
        y = iterate(str, i)
    end

    return i, count
end

# Move past consecutive lines that start with commentchar.
# Return a tuple of the new pos in str and the amount of comment lines moved past.
function eatcommentlines(str, i=1, l=lastindex(str), commentchar::Union{Char, Nothing}=nothing)
    commentchar === nothing && return i, 0

    count = 0
    while i <= l && str[i] == commentchar
        i = getlineend(str, i)
        y = iterate(str, i)
        y === nothing && return i, count
        i = y[2]
        i, lines = eatnewlines(str, i)
        count += lines
    end
    return i, count
end

function stripquotes(x)
    x[1] in ('\'', '"') && x[1] == x[end] ?
        strip(x, x[1]) : x
end

function getlineend(str, i=1, l=lastindex(str))
    y = iterate(str, i)
    while y!==nothing
        c = y[1]; ii = y[2]
        isnewline(c) && break
        i = ii
        y = iterate(str, i)
    end

    return prevind(str, i)
end

# This is similar to getlineend, but ignores line ends inside
# quotes
function getrowend(str, i, len, opts, delim)
    i0 = i
    i = eatwhitespaces(str, i, len)
    y = iterate(str, i)
    while y!==nothing
        c = y[1]; i = y[2]
        if c==Char(opts.quotechar)
            # We are now inside a quoted field
            y2 = iterate(str, i)
            while y2!==nothing
                c = y2[1]; i = y2[2]
                if c==Char(opts.escapechar)
                    y3 = iterate(str, i)
                    if y3===nothing
                        if c==Char(opts.quotechar)
                            return prevind(str, i)
                        else
                            error("Parsing error, quoted string never terminated.")
                        end
                    else
                        c2 = y3[1]; ii = y3[2]
                        if c2==Char(opts.quotechar)
                            i = ii
                        elseif c==Char(opts.quotechar)
                            break
                        end
                    end
                elseif c==Char(opts.quotechar)
                    break;
                end
                y2 = iterate(str, i)
                if y2===nothing
                    error("Parsing error, quoted string never terminated.")
                end
            end
            i = eatwhitespaces(str, i, len)
            y4 = iterate(str, i)
            if y4!==nothing
                c = y4[1]; i4 = y4[2]
                if isnewline(c)
                    return prevind(str, i)
                elseif c!=Char(delim)
                    error("Invalid line")
                end
            else
                return prevind(str, i)
            end
        else
            # We are now inside a non quoted field
            while y!==nothing
                c = y[1]; i = y[2]
                if c==Char(delim)
                    i = eatwhitespaces(str, i)
                    break
                elseif isnewline(c)
                    return prevind(str, i, 2)
                end
                y = iterate(str, i)
            end
        end
        y = iterate(str, i)
    end
    return prevind(str, i)
end

### Testing helpers

unwrap(xs) = (get(xs[1]), xs[2:end]...)
failedat(xs) = (@assert isnull(xs[1]); xs[2])

# String speedup hacks

# StrRange
# This type is the beginning of a hack to avoid allocating 3 objects
# instead of just 1 when using the `tryparsenext` framework.
# The expression (Some{String}("xyz"), 4) asks the GC to track
# the string, the nullable and the tuple. Instead we return
# (Some{StrRange}(StrRange(0,3)), 4) which makes 0 allocations.
# later when assigning the column inside `tryparsesetindex` we
# create the string. See `setcell!`
struct StrRange
    offset::Int
    length::Int
    escapecount::Int
end

function getlineat(str, i)
    l = lastindex(str)
    if i <= l
        ii = prevind(str, i)
    else
        ii = l
    end
    line_start = i
    while ii > 0 && !isnewline(str[ii])
        line_start = ii
        ii = prevind(str, line_start)
    end

    # TODO Handle nothing case
    c, ii = iterate(str, line_start)
    line_end = line_start
    while !isnewline(c) && ii <= l
        line_end = ii
        # TODO Handle nothing case
        c, ii = iterate(str, ii)
    end

    line_start:line_end
end

_widen(::Type{UInt8}) = UInt16 # fix for bad Base behavior
_widen(::Type{Int8}) = Int16
_widen(T) = widen(T)


text:
module EnsembleInference

using Distributions: Distributions
using LinearAlgebra
using Manifolds: Manifolds, SpecialEuclidean, SpecialOrthogonal, TranslationGroup
using PDMats: PDMats
using Random
using Statistics: Statistics
using StatsBase: StatsBase

const SO3 = SpecialOrthogonal(3)
const SE3 = SpecialEuclidean(3)

# this type is defined to avoid Manifolds's complicated type vector space type hierarchy and
# to avoid choosing a point representation
struct LieAlgebra{G}
    group::G
end

const so3 = LieAlgebra(SO3)
const se3 = LieAlgebra(SE3)

struct BasisVector{B,I}
    basis::B
end
BasisVector(basis, i::Integer) = BasisVector{typeof(basis),i}(basis)

# basis vectors needed for SO(3) (E1-E3) and SE(3) (E1-E6)
const E1 = BasisVector(Manifolds.DefaultOrthogonalBasis(), 1)
const E2 = BasisVector(Manifolds.DefaultOrthogonalBasis(), 2)
const E3 = BasisVector(Manifolds.DefaultOrthogonalBasis(), 3)
const E4 = BasisVector(Manifolds.DefaultOrthogonalBasis(), 4)
const E5 = BasisVector(Manifolds.DefaultOrthogonalBasis(), 5)
const E6 = BasisVector(Manifolds.DefaultOrthogonalBasis(), 6)

const SHAPESPEC_SE3 = Manifolds.ShapeSpecification(
    Manifolds.StaticReshaper(), Manifolds.base_manifold(SpecialEuclidean(3)).manifolds...
)

"""
    inversion(d::Distribution)

Given a distribution ``d`` of elements ``g ∈ G``, where ``G`` is a group, the inversion of
``d`` is the corresponding distribution of ``g^{-1}``.
"""
inversion

se3_array(x) = Manifolds.ProductArray(SHAPESPEC_SE3, x)

include("utility.jl")
include("dirac.jl")
include("haar.jl")
include("diffusion_normal.jl")
include("representation.jl")

export SO3, so3, SE3, se3
export DiffusionNormal, DiracDelta, Haar
export inversion, se3_array

end # module


text:
# SPDX-License-Identifier: X11
# 2020-11-09 (Timed)
# Akari (500pt)

function main()
  H,W,N,M = parse.(Int, (split ∘ readline)())
  G = zeros(Int, (H, W))
  #==
   = 0: empty
   = 1: bulb
   = 2: block
   =#
  for i ∈ 1:N
    a,b = parse.(Int, (split ∘ readline)())
    G[a,b] = 1
  end
  for i ∈ 1:M
    c,d = parse.(Int, (split ∘ readline)())
    G[c,d] = 2
  end
  I = zeros(Int, (H, W))
  #==
   = 0: dark
   = 1: light
   =#

  for h ∈ 1:H, w ∈ 1:W
    G[h,w] == 2 && continue

    lit = false

    for i ∈ h:-1:1
      if G[i,w] == 1
        I[h,w] = 1
        lit = true
        break
      end
      G[i,w] == 2 && break
    end
    lit && continue

    for i ∈ h + 1:H
      if G[i,w] == 1
        I[h,w] = 1
        lit = true
        break
      end
      G[i,w] == 2 && break
    end
    lit && continue

    for i ∈ w:-1:1
      if G[h,i] == 1
        I[h,w] = 1
        lit = true
        break
      end
      G[h,i] == 2 && break
    end
    lit && continue

    for i ∈ w + 1:W
      if G[h,i] == 1
        I[h,w] = 1
        lit = true
        break
      end
      G[h,i] == 2 && break
    end
  end

  (println ∘ sum)(I)
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
using Pkg

ssh = true
if ssh
    git = "git@github.com:SebastianAment/"
else
    git = "https://github.com/SebastianAment/"
end

add(x) = Pkg.add(Pkg.PackageSpec(url = git * x * ".git"))

add("LazyInverse.jl")
add("LazyLinearAlgebra.jl")
add("LinearAlgebraExtensions.jl")
add("KroneckerProducts.jl")
add("WoodburyIdentity.jl")
add("CovarianceFunctions.jl")
add("OptimizationAlgorithms.jl")
add("BackgroundSubtraction.jl")
println("Trying to add GaussianDistributions.jl via ssh. If this fails, please add the package locally by typing \"]add path/to/GaussianDistributions.jl\" in the REPL.")
add("GaussianDistributions.jl") # needs ssh OR add path to local package


text:
using Revise
using Knet
using IterTools
include("data.reader.jl")
using Base.Iterators: flatten
import .Iterators: cycle, Cycle, take
using Statistics
include("utils.jl")
using Plots; default(fmt=:png, ls=:auto)
include("models.jl")
include("training.jl")
include("plots.jl")


text:
using Test

@time include("ffunctors.jl")


text:
@doc doc"""

```math
dS = μSdt + \sqrt{v}SdW_1
dv = κ(Θ-v)dt + σ\sqrt{v}dW_2
dW_1 dW_2 = ρ dt
```

"""
function HestonProblem(μ, κ, Θ, σ, ρ, u0, tspan; seed = UInt64(0), kwargs...)
    f = function (du, u, p, t)
        du[1] = μ * u[1]
        du[2] = κ * (Θ - u[2])
    end
    g = function (du, u, p, t)
        du[1] = √u[2] * u[1]
        du[2] = Θ * √u[2]
    end
    Γ = [1 ρ; ρ 1] # Covariance Matrix
    noise_rate_prototype = nothing

    if seed == 0
        seed = rand(UInt64)
    end
    noise = CorrelatedWienerProcess!(Γ, tspan[1], zeros(2), zeros(2),
                                     rng = Xorshifts.Xoroshiro128Plus(seed))

    sde_f = SDEFunction{true}(f, g)
    SDEProblem(sde_f, g, u0, tspan, noise = noise, seed = seed, kwargs...)
end

@doc doc"""

``d \ln S(t) = (r(t) - q(t) - \frac{Θ(t,S)^2}{2})dt + σ dW_t``

Solves for ``log S(t)``.
"""
function GeneralizedBlackScholesProblem(r, q, Θ, σ, u0, tspan; kwargs...)
    f = function (u, p, t)
        r(t) - q(t) - Θ(t, exp(u))^2 / 2
    end
    g = function (u, p, t)
        σ
    end
    SDEProblem{false}(f, g, u0, tspan; kwargs...)
end

@doc doc"""

``d \ln S(t) = (r(t) - \frac{Θ(t,S)^2}{2})dt + σ dW_t``

Solves for ``log S(t)``.
"""
mutable struct BlackScholesProblem{uType, tType, tupType, isinplace, NP, F, F2, C, ND, MM
                                   } <:
               DiffEqBase.AbstractSDEProblem{uType, tType, isinplace, ND}
    r::tType
    Θ::tType
    σ::tType
    u0::uType
    tspan::tupType
    p::Nothing
    f::F
    g::F2
    noise::NP
    callback::C
    noise_rate_prototype::ND
    mass_matrix::MM
    seed::UInt64
end

function BlackScholesProblem(r, Θ, σ, u0, tspan; callback = CallbackSet(),
                             noise_rate_prototype = nothing, seed = UInt64(0))
    GeneralizedBlackScholesProblem(r, (t) -> 0, Θ, σ, u0, tspan, callback = callback,
                                   seed = seed)
end

@doc doc"""

``dx = a(b(t)-x)dt + σ dW_t``

"""
function ExtendedOrnsteinUhlenbeckProblem(a, b, σ, u0, tspan; kwargs...)
    f = function (u, p, t)
        a * (b(t) - u)
    end
    g = function (u, p, t)
        σ
    end
    SDEProblem{false}(f, g, u0, tspan; kwargs...)
end

@doc doc"""

``dx = a(r-x)dt + σ dW_t``

"""
function OrnsteinUhlenbeckProblem(a, r, σ, u0, tspan; kwargs...)
    f = function (u, p, t)
        a * (r - u)
    end
    g = function (u, p, t)
        σ
    end
    SDEProblem{false}(f, g, u0, tspan; kwargs...)
end

@doc doc"""

``dx = μ dt + σ dW_t``

"""
function GeometricBrownianMotionProblem(μ, σ, u0, tspan; kwargs...)
    f = function (u, p, t)
        μ
    end
    g = function (u, p, t)
        σ
    end
    SDEProblem{false}(f, g, u0, tspan; kwargs...)
end

@doc doc"""

``dx = σ(t)e^{at} dW_t``

"""
function MfStateProblem(a, σ, u0, tspan; kwargs...)
    f = function (u, p, t)
        0
    end
    g = function (u, p, t)
        σ(t) * exp(a * t)
    end
    SDEProblem{false}(f, g, u0, tspan; kwargs...)
end


text:
# Solvers for mixed complementarity problems

macro reformulate(df)
    esc(quote
        if reformulation == :smooth
            rf = mcp_smooth($df, lower, upper)
        elseif reformulation == :minmax
            rf = mcp_minmax($df, lower, upper)
        else
            throw(ArgumentError("Unknown reformulation $reformulation"))
        end
    end)
end

function mcpsolve(df::OnceDifferentiable,
                  lower::AbstractArray{T},
                  upper::AbstractArray{T},
                  initial_x::AbstractArray{T};
                  method::Symbol = :trust_region,
                  reformulation::Symbol = :smooth,
                  xtol::Real = zero(T),
                  ftol::Real = convert(T,1e-8),
                  iterations::Integer = 1_000,
                  store_trace::Bool = false,
                  show_trace::Bool = false,
                  extended_trace::Bool = false,
                  linesearch = LineSearches.BackTracking(),
                  factor::Real = one(T),
                  autoscale::Bool = true) where T

    @reformulate df
    nlsolve(rf,
            initial_x, method = method, xtol = xtol, ftol = ftol,
            iterations = iterations, store_trace = store_trace,
            show_trace = show_trace, extended_trace = extended_trace,
            linesearch = linesearch, factor = factor, autoscale = autoscale)
end

function mcpsolve(f,
                  j,
                  lower::AbstractArray{T},
                  upper::AbstractArray{T},
                  initial_x::AbstractArray{T};
                  method::Symbol = :trust_region,
                  reformulation::Symbol = :smooth,
                  xtol::Real = zero(T),
                  ftol::Real = convert(T,1e-8),
                  iterations::Integer = 1_000,
                  store_trace::Bool = false,
                  show_trace::Bool = false,
                  extended_trace::Bool = false,
                  linesearch = LineSearches.BackTracking(),
                  factor::Real = one(T),
                  autoscale = true,
                  inplace = true) where T
    if inplace
        df = OnceDifferentiable(f, initial_x, initial_x)
    else
        df = OnceDifferentiable(not_in_place(f, j)..., initial_x, initial_x)
    end
    @reformulate df
    nlsolve(rf,
            initial_x, method = method, xtol = xtol, ftol = ftol,
            iterations = iterations, store_trace = store_trace,
            show_trace = show_trace, extended_trace = extended_trace,
            linesearch = linesearch, factor = factor, autoscale = autoscale)
end

function mcpsolve(f,
                  lower::AbstractArray{T},
                  upper::AbstractArray{T},
                  initial_x::AbstractArray{T};
                  method::Symbol = :trust_region,
                  reformulation::Symbol = :smooth,
                  xtol::Real = zero(T),
                  ftol::Real = convert(T,1e-8),
                  iterations::Integer = 1_000,
                  store_trace::Bool = false,
                  show_trace::Bool = false,
                  extended_trace::Bool = false,
                  linesearch = LineSearches.BackTracking(),
                  factor::Real = one(T),
                  autoscale::Bool = true,
                  autodiff = :central,
                  inplace = true) where T
    if inplace
        df = OnceDifferentiable(f, initial_x, initial_x, autodiff)
    else
        df = OnceDifferentiable(not_in_place(f), initial_x, initial_x, autodiff)
    end

    @reformulate df
    nlsolve(rf,
            initial_x, method = method, xtol = xtol, ftol = ftol,
            iterations = iterations, store_trace = store_trace,
            show_trace = show_trace, extended_trace = extended_trace,
            linesearch = linesearch, factor = factor, autoscale = autoscale)
end


text:
push!(LOAD_PATH, joinpath(homedir(),"GitHub/KenyaCoV/src"))
using Plots,Parameters,Distributions,DifferentialEquations,JLD2,DataFrames
using Revise
import KenyaCoV

"""
Load a completely susceptible population differentiated by county and urban/rural.
Also, load the optimised mixing matrix for spatial transmission
"""
# #Load data and completely susceptible Population
# u0,P,transport_matrix = model_ingredients_from_data("./src/2009_National_Estimates_of_Rural_and_Urban_Populations_by_County.csv",
#Load data completely susceptible Population
u0,P,transport_matrix = KenyaCoV.model_ingredients_from_data("data/combined_population_estimates.csv",
                                                             "data/optimal_transition_matrix.jld2",
                                                            "data/optimal_movement_matrix.jld2",
                                                            "data/flight_numbers.csv",
                                                            "data/projected_global_prevelance.csv")

"""
Example of methods that modify underlying parameters
"""
#This method modifies the parameter set for changing the mixing structure
# KenyaCoV.transportstructure_params!(P,[0.001 for i = 1:KenyaCoV.n],transport_matrix)
#You can modify other parameters directly
P.τ = 0. #e.g. no treatment rate
for (i,p) in enumerate(P.global_prev)
    P.global_prev[i] = 0.
end
# P.global_prev
P.μ₁ = 0.
P.ϵ = 1.
P.β = 2.2*P.γ
#Define initial conditions by modifying the completely susceptible population
P.dt = 0.25

u0[30,4,1] += 10#One diseased in Nairobi

prob = KenyaCoV.create_KenyaCoV_non_neg_prob(u0,(0.,365.),P)
# prob_tl = KenyaCoV.create_KenyaCoV_prob(u0,(0.,365.),P)
# prob_ode = KenyaCoV.create_KenyaCoV_ode_prob(u0,(0.,365.),P)
# sol_ode = solve(prob_ode,Tsit5())
# @time sol_tl = solve(prob_tl,SimpleTauLeaping(),dt = P.dt)

@time sol = solve(prob,FunctionMap(),dt = P.dt)


text:
@testset "length - tests" begin
    @testset "length(r)" begin
        for (r,b) in ((OneToMRange(10), OneTo(10)),
                      (OneToSRange(UInt(10)), OneTo(UInt(10))),
                      (UnitMRange(1, 10), UnitRange(1, 10)),
                      (UnitSRange(1.,10.), UnitRange(1.0, 10.0)),
                      (StepMRange(1, 2, 10), StepRange(1, 2, 10)),
                      (StepSRange(UInt32(1), UInt32(2), UInt32(10)), StepRange(UInt32(1), UInt32(2), UInt32(10))),
                      (mrange(1.0, step=2.0, stop=10.0), 1.0:2.0:10.0),
                      (srange(UInt32(1), step=UInt32(2), stop=UInt32(10)), UInt32(1):UInt32(2):UInt32(10)),
                     )
            @test @inferred(length(r)) == length(b)
            @test @inferred(length(r)) == length(b)
            if b isa StepRangeLen
                @test @inferred(stephi(r)) == stephi(b)
                @test @inferred(steplo(r)) == steplo(b)
            end
        end
    end

    @testset "can_set_length" begin
        @test @inferred(can_set_length(LinRange)) == false
        @test @inferred(can_set_length(1:10)) == false
        @test @inferred(can_set_length(LinMRange)) == true
        @test @inferred(can_set_length(StepMRangeLen)) == true
        @test @inferred(can_set_length(StepMRange)) == true
        @test @inferred(can_set_length(UnitMRange)) == true
        @test @inferred(can_set_length(OneToMRange)) == true
    end

    @testset "set_length!" begin
        @test @inferred(set_length!(LinMRange(1, 10, 5), 10)) == LinMRange(1, 10, 10)
        @test @inferred(set_length!(LinMRange(1, 10, 5), UInt32(10))) == LinMRange(1, 10, 10)
        @test @inferred(set_length!(LinMRange(1,1,0), 1)) == LinMRange(1,1,1)
        @test @inferred(set_length!(StepMRangeLen(1, 1, 10), 11)) == StepMRangeLen(1, 1, 11)
        @test @inferred(set_length!(StepMRangeLen(1, 1, 10), UInt32(11))) == StepMRangeLen(1, 1, 11)
        @test @inferred(set_length!(StepMRange(1, 1, 10), UInt32(11))) == StepMRange(1, 1, 11)
        @test @inferred(set_length!(OneToMRange(10), UInt32(11))) == OneToMRange(11)
    end

    @testset "set_lediv!" begin
        @test @inferred(set_lendiv!(LinMRange(1, 10, 5), 9)) == LinMRange(1, 10, 10)
        @test @inferred(set_lendiv!(LinMRange(1, 10, 5), UInt32(9))) == LinMRange(1, 10, 10)
    end

    for frange in (mrange, srange)
        @testset "length($(frange))" begin
            @test length(frange(.1, step=.1, stop=.3)) == 3
            @test length(frange(1.1, step=1.1, stop=3.3)) == 3
            @test length(frange(1.1, step=1.3, stop=3)) == 2
            @test length(frange(1, step=1, stop=1.8)) == 1
            @test length(frange(1, step=.2, stop=2)) == 6
            @test length(frange(1., step=.2, stop=2.)) == 6
            @test length(frange(2, step=-.2, stop=1)) == 6
            @test length(frange(2., step=-.2, stop=1.)) == 6
            @test length(frange(2, step=.2, stop=1)) == 0
            @test length(frange(2., step=.2, stop=1.)) == 0

            @test length(frange(1, 0)) == 0
            @test length(frange(0.0, -0.5)) == 0
            @test length(frange(1, step=2, stop=0)) == 0
            @test length(frange(Char(0), Char(0x001fffff))) == 2097152
            @test length(frange(typemax(UInt64)//one(UInt64), step=1, stop=typemax(UInt64)//one(UInt64))) == 1
        end

        @testset "length($(frange)) with typemin/typemax" begin
            let r = frange(typemin(Int64), step=2, stop=typemax(Int64)), s = frange(typemax(Int64), step=-2, stop=typemin(Int64))
                @test first(r) == typemin(Int64)
                @test last(r) == (typemax(Int64)-1)
                #@test_throws OverflowError length(r)

                @test first(s) == typemax(Int64)
                @test last(s) == (typemin(Int64)+1)
                #@test_throws OverflowError length(s)
            end

            @test length(frange(typemin(Int64), step=3, stop=typemax(Int64))) == 6148914691236517206
            @test length(frange(typemax(Int64), step=-3, stop=typemin(Int64))) == 6148914691236517206

            #= No static type for big
            for s in 3:100
                @test length(frange(typemin(Int), step=s, stop=typemax(Int))) == length(big(typemin(Int)):big(s):big(typemax(Int)))
                @test length(frange(typemax(Int), step=-s, stop=typemin(Int))) == length(big(typemax(Int)):big(-s):big(typemin(Int)))
            end
            =#

            #= TODO
            @test length(frange(UInt(1), step=UInt(1), stop=UInt(0))) == 0
            @test length(frange(typemax(UInt), step=UInt(1), stop=(typemax(UInt)-1))) == 0
            @test length(frange(typemax(UInt), step=UInt(2), stop=(typemax(UInt)-1))) == 0
            @test length((frange(typemin(Int)+3, step=5, stop=(typemin(Int)+1)))) == 0
            =#
        end

        # avoiding intermediate overflow (#5065)
        @test length(frange(1, step=4, stop=typemax(Int))) == div(typemax(Int),4) + 1

        # issue #6364
        @test length((frange(1, 64))*(pi/5)) == 64

        @testset "overflow in length($frange)" begin
            Tset = Int === Int64 ? (Int,UInt,Int128,UInt128) :
                                   (Int,UInt,Int64,UInt64,Int128, UInt128)
            for T in Tset
                @test_throws OverflowError length(zero(T):typemax(T))
                @test_throws OverflowError length(typemin(T):typemax(T))
                @test_throws OverflowError length(zero(T):one(T):typemax(T))
                @test_throws OverflowError length(typemin(T):one(T):typemax(T))
                if T <: Signed
                    @test_throws OverflowError length(-one(T):typemax(T)-one(T))
                    @test_throws OverflowError length(-one(T):one(T):typemax(T)-one(T))
                end
            end
        end
    end
end