text:
module materials

using Interpolations
using DelimitedFiles

export LoadMaterial, spesifics

function LoadMaterial(path) # Loads material where complex refractive index is contained in one file
    material = readdlm(path, ',', header=true)[1][:, 1:3]

    n_data = [material[:, 1] .* 1e-6, material[:, 2]]
    k_data = [material[:, 1] .* 1e-6, material[:, 3]]
    n_real = LinearInterpolation(n_data[1], n_data[2])
    n_complex = LinearInterpolation(k_data[1], k_data[2])

    n_total(λ) = n_real(λ) - abs(n_complex(λ)) * 1im
    n_total
end

function LoadMaterial(path, number::Number) # Loads material where several complex refractive inicies are contained in one file
    material = readdlm(path, ',', skipstart=2)

    n_data = [material[:, 1] .* 1e-9, material[:, (number*2)]]
    k_data = [material[:, 1] .* 1e-9, material[:, (number*2 + 1)]]
    n_real = LinearInterpolation(n_data[1], n_data[2])
    n_complex = LinearInterpolation(k_data[1], k_data[2])

    n_total(λ) = n_real(λ) - abs(n_complex(λ)) * 1im
    n_total
end

function LoadMaterial(path_n, path_k, number) # Loads material where n and k are split and dumped by webplotdigitizer
    material_n = readdlm(path_n, ',', skipstart=2)[:, (number*2 - 1):(number*2)]
    material_k = readdlm(path_k, ',', skipstart=2)[:, (number*2 - 1):(number*2)]

    material_n = material_n[(material_n[:, 1] .!= ""), :]
    material_k = material_k[(material_k[:, 1] .!= ""), :]

    material_n = material_n[sortperm(material_n[:, 1]), :]
    material_k = material_k[sortperm(material_k[:, 1]), :]

    n_estimator = LinearInterpolation(material_n[:, 1] .* 1e-9, material_n[:, 2])
    k_estimator = LinearInterpolation(material_k[:, 1] .* 1e-9, material_k[:, 2])

    n_total(λ) = n_estimator(λ) - abs(k_estimator(λ)) * 1im
    n_total
end

module spesifics

import ..LoadMaterial
using ...analyticalmaterials

export Air, Au_Johnson, Au_Werner, Au_11nm_Rosenblatt, Au_21nm_Rosenblatt, Au_44nm_Rosenblatt, Au_Babar, Au_McPeak, Au_OlmonEvaporated, Au_OlmonSingleChrystaline, Au_OlmonTemplateStripped, Pd_Johnson, Pd_Werner, Pd_Palm_2018, Ag, LiF, SiO2_core_Sellmeier, SiO2_thinfilm_Ciprian, Au_unloaded, Pd014_unloaded, Pd034_unloaded, Pd034_loaded, Pd042_unloaded, Pd042_loaded, Pd052_unloaded, Pd052_loaded, Pd073_unloaded, Pd073_loaded, Pd_unloaded, Pd_loaded, Pd000S5, Pd100S20, Pd070S9to17, Pd041S21, Pd051S22, Pd061S23, Pd069S24, Pd064S7and18, Pd053S6, Pd074S8drude, Pd061S19, Pd061S19drude

Air(λ) = 1

Au_11nm_Rosenblatt_estimator = LoadMaterial("materials/refractive_index/Au_11nm_Rosenblatt.csv")
Au_11nm_Rosenblatt(λ) = Au_11nm_Rosenblatt_estimator(λ)

Au_21nm_Rosenblatt_estimator = LoadMaterial("materials/refractive_index/Au_21nm_Rosenblatt.csv")
Au_21nm_Rosenblatt(λ) = Au_21nm_Rosenblatt_estimator(λ)

Au_44nm_Rosenblatt_estimator = LoadMaterial("materials/refractive_index/Au_44nm_Rosenblatt.csv")
Au_44nm_Rosenblatt(λ) = Au_44nm_Rosenblatt_estimator(λ)

Au_Babar_estimator = LoadMaterial("materials/refractive_index/Au_Babar.csv")
Au_Babar(λ) = Au_Babar_estimator(λ)

Au_McPeak_estimator = LoadMaterial("materials/refractive_index/Au_McPeak.csv")
Au_McPeak(λ) = Au_McPeak_estimator(λ)

Au_OlmonEvaporated_estimator = LoadMaterial("materials/refractive_index/Au_OlmonEvaporated.csv")
Au_OlmonEvaporated(λ) = Au_OlmonEvaporated_estimator(λ)

Au_OlmonSingleChrystaline_estimator = LoadMaterial("materials/refractive_index/Au_OlmonSingleChrystaline.csv")
Au_OlmonSingleChrystaline(λ) = Au_OlmonSingleChrystaline_estimator(λ)

Au_OlmonTemplateStripped_estimator = LoadMaterial("materials/refractive_index/Au_OlmonTemplateStripped.csv")
Au_OlmonTemplateStripped(λ) = Au_OlmonTemplateStripped_estimator(λ)

Au_Johnson_estimator = LoadMaterial("materials/refractive_index/Au_Johnson.csv")
Au_Johnson(λ) = Au_Johnson_estimator(λ)

Pd_Johnson_estimator = LoadMaterial("materials/refractive_index/Pd_Johnson.csv")
Pd_Johnson(λ) = Pd_Johnson_estimator(λ)

Au_Werner_estimator = LoadMaterial("materials/refractive_index/Au_Werner.csv")
Au_Werner(λ) = Au_Werner_estimator(λ)

Pd_Werner_estimator = LoadMaterial("materials/refractive_index/Pd_Werner.csv")
Pd_Werner(λ) = Pd_Werner_estimator(λ)

Ag_estimator = LoadMaterial("materials/refractive_index/Ag.csv")
Ag(λ) = Ag_estimator(λ)

LiF_estimator = LoadMaterial("materials/refractive_index/LiF.csv")
LiF(λ) = LiF_estimator(λ)

Pd_Palm_2018_estimator = LoadMaterial("materials/refractive_index/Pd_Palm.csv")
Pd_Palm_2018(λ) = Pd_Palm_2018_estimator(λ)

# F. Downes eq.1
# Note that small wavelengths give imaginary refractive index, this gives unphysical  results
SiO2_core_Sellmeier(λ) = √( 1
    + 0.6961663*(λ*1e6)^2/((λ*1e6)^2-0.0684043)
    + 0.4079426*(λ*1e6)^2/((λ*1e6)^2-0.1162414)
    + 0.8974794*(λ*1e6)^2/((λ*1e6)^2-9.896161) )

# F. Downes eq. 2
# Note that small wavelengths give imaginary refractive index, this gives unphysical  results
SiO2_thinfilm_Ciprian(λ) = √( 1 + 1.1336*(λ*1e6)^2/((λ*1e6)^2-9.261e-2) )

Au_unloaded_estimator = LoadMaterial("materials/palm_PdAu-alloys/unloaded_n.csv", "materials/palm_PdAu-alloys/unloaded_k.csv", 1)
Au_unloaded(λ) = Au_unloaded_estimator(λ)

Pd014_unloaded_estimator = LoadMaterial("materials/palm_PdAu-alloys/unloaded_n.csv", "materials/palm_PdAu-alloys/unloaded_k.csv", 2)
Pd014_unloaded(λ) = Pd014_unloaded_estimator(λ)

Pd034_unloaded_estimator = LoadMaterial("materials/palm_PdAu-alloys/unloaded_n.csv", "materials/palm_PdAu-alloys/unloaded_k.csv", 3)
Pd034_unloaded(λ) = Pd034_unloaded_estimator(λ)

Pd042_unloaded_estimator = LoadMaterial("materials/palm_PdAu-alloys/unloaded_n.csv", "materials/palm_PdAu-alloys/unloaded_k.csv", 4)
Pd042_unloaded(λ) = Pd042_unloaded_estimator(λ)

Pd052_unloaded_estimator = LoadMaterial("materials/palm_PdAu-alloys/unloaded_n.csv", "materials/palm_PdAu-alloys/unloaded_k.csv", 5)
Pd052_unloaded(λ) = Pd052_unloaded_estimator(λ)

Pd073_unloaded_estimator = LoadMaterial("materials/palm_PdAu-alloys/unloaded_n.csv", "materials/palm_PdAu-alloys/unloaded_k.csv", 6)
Pd073_unloaded(λ) = Pd073_unloaded_estimator(λ)

Pd_unloaded_estimator = LoadMaterial("materials/palm_PdAu-alloys/unloaded_n.csv", "materials/palm_PdAu-alloys/unloaded_k.csv", 7)
Pd_unloaded(λ) = Pd_unloaded_estimator(λ)

Pd034_loaded_estimator = LoadMaterial("materials/palm_PdAu-alloys/loaded_n.csv", "materials/palm_PdAu-alloys/loaded_k.csv", 3)
Pd034_loaded(λ) = Pd034_loaded_estimator(λ)

Pd042_loaded_estimator = LoadMaterial("materials/palm_PdAu-alloys/loaded_n.csv", "materials/palm_PdAu-alloys/loaded_k.csv", 4)
Pd042_loaded(λ) = Pd042_loaded_estimator(λ)

Pd052_loaded_estimator = LoadMaterial("materials/palm_PdAu-alloys/loaded_n.csv", "materials/palm_PdAu-alloys/loaded_k.csv", 5)
Pd052_loaded(λ) = Pd052_loaded_estimator(λ)

Pd073_loaded_estimator = LoadMaterial("materials/palm_PdAu-alloys/loaded_n.csv", "materials/palm_PdAu-alloys/loaded_k.csv", 6)
Pd073_loaded(λ) = Pd073_loaded_estimator(λ)

Pd_loaded_estimator = LoadMaterial("materials/palm_PdAu-alloys/loaded_n.csv", "materials/palm_PdAu-alloys/loaded_k.csv", 7)
Pd_loaded(λ) = Pd_loaded_estimator(λ)

Pd000S5_estimator = LoadMaterial("materials/thesisdata/Au_Pd_Multi_sample_9_17.txt", 1)
Pd000S5(λ) = Pd000S5_estimator(λ)

Pd100S20_estimator = LoadMaterial("materials/thesisdata/Au_Pd_Multi_sample_9_17.txt", 2)
Pd100S20(λ) = Pd100S20_estimator(λ)

Pd070S9to17_estimator = LoadMaterial("materials/thesisdata/Au_Pd_Multi_sample_9_17.txt", 3)
Pd070S9to17(λ) = Pd070S9to17_estimator(λ)

Pd041S21_estimator = LoadMaterial("materials/thesisdata/Samples_21_24_B_spline_N_k.txt", 1)
Pd041S21(λ) = Pd041S21_estimator(λ)

Pd051S22_estimator = LoadMaterial("materials/thesisdata/Samples_21_24_B_spline_N_k.txt", 2)
Pd051S22(λ) = Pd051S22_estimator(λ)

Pd061S23_estimator = LoadMaterial("materials/thesisdata/Samples_21_24_B_spline_N_k.txt", 3)
Pd061S23(λ) = Pd061S23_estimator(λ)

Pd069S24_estimator = LoadMaterial("materials/thesisdata/Samples_21_24_B_spline_N_k.txt", 4)
Pd069S24(λ) = Pd069S24_estimator(λ)

Pd064S7and18_estimator = LoadMaterial("materials/thesisdata/Remaining_samples.txt", 1)
Pd064S7and18(λ) = Pd064S7and18_estimator(λ)

Pd053S6_estimator = LoadMaterial("materials/thesisdata/Remaining_samples.txt", 2)
Pd053S6(λ) = Pd053S6_estimator(λ)

Pd074S8drude_estimator = LoadMaterial("materials/thesisdata/Remaining_samples.txt", 3)
Pd074S8drude(λ) = Pd074S8drude_estimator(λ)

Pd061S19_estimator = LoadMaterial("materials/thesisdata/Remaining_samples.txt", 4)
Pd061S19(λ) = Pd061S19_estimator(λ)

Pd061S19drude_estimator = LoadMaterial("materials/thesisdata/Remaining_samples.txt", 5)
Pd061S19drude(λ) = Pd061S19_estimator(λ)

end # spesifics

end # materials


text:
using MATLAB
using Test
using SparseArrays

# test engine

restart_default_msession()

a = [1. 2. 3.; 4. 5. 6.]
b = [2. 3. 4.; 8. 7. 6.]

@mput a b
mat"""
    r1 = a .* b;
    r2 = a + b;
"""
@mget r1 r2

@test isequal(r1, a .* b)
@test isequal(r2, a + b)

@mget r1::Vector
@test isequal(r1, vec(a .* b))

s = sparse([1. 0. 0.; 2. 3. 0.; 0. 4. 5.])
put_variable(:s, s)
s2 = get_variable(:s)
@test isequal(s, s2)

# mxcall

r = mxcall(:plus, 1, a, b)
@test isequal(r, a + b)

(xx, yy) = mxcall(:meshgrid, 2, [1., 2.], [3., 4.])
@test isequal(xx, [1. 2.; 1. 2.])
@test isequal(yy, [3. 3.; 4. 4.])

close_default_msession()

# test for segfault

s = MSession()
close(s)
@test_throws UndefRefError close(s)

# segfault on deleted references
x = mxarray(3.0)
delete(x)
@test_throws UndefRefError delete(x)
@test_throws UndefRefError nrows(x)
@test_throws UndefRefError is_numeric(x)
@test_throws UndefRefError jscalar(x)
@test_throws UndefRefError jvalue(x)

# make sure restart_default_msession() doesn't error on null references on
# default msession
s = get_default_msession()
close(s)
restart_default_msession()

ext: jl


text:
using Test, NLSolvers, LinearAlgebra
@testset "truncated SPR solver" begin
    m = NLSolvers.TCG()

    H = [0.3 0.0; 0.0 0.9]
    g = [0.2, 0.4]

    m(g, H, 0.7, rand(2), 1)

    # Since -H is negative definite the solution is guaranteed
    # to be at the boundary unless g = 0
    for Δ in range(0, 100; step = 0.1)
        @test norm(m(g, -H, Δ, rand(2), 1).p, 2) ≈ Δ
    end

    # Small gradient entries
    g .= [1e-12, 1e-9]
    for Δ in range(0, 100; step = 0.1)
        @test norm(m(g, -H, Δ, rand(2), 1).p, 2) ≈ Δ
    end

    # Mixed gradient entries
    g .= [1e12, 1e-9]
    for Δ in range(0, 100; step = 0.1)
        @test norm(m(g, -H, Δ, rand(2), 1).p, 2) ≈ Δ
    end

    # Zero case
    g = [0.0, 0.0]
    for Δ in range(0, 100; step = 0.1)
        @test norm(m(g, -H, Δ, rand(2), 1).p, 2) == 0
    end
end

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
"""
The Accounting Trinity or [CAL]!

[]: Important (somehow) - (it's objective alright)
(): an Optional read , acts as a comment
me: personal words, notes , to add an overall flavor & meaning to the comment itself - as needed
A/C : Account - a fancy way accountants use for any account 'book'
1.Assets k let's start - here goes nothing

2.Capital :  an abstract entity me: i.e. bank account value- me: has a Direct Relation with either an 'Asset' or a Liability A/C "Account"
3.Assets: a Tangible object, of Value, has an/a owner(s)

below is the explination, non-academic - based on the Personal View -- for your eyes only,ok
# 3.1 thinking about it - I'm a comment after all oh wait maybe I'm not -- who Am ... I ??? - stop talking this is your brain police stop --or what ? - we'll see
#3.2 Ownership:

the actual Owner: for whom the Account belongs to - usualy 1 or more Seperate entities (persons, companies - or mix of the 2 )

1. One Owner: 1 is the lonliest number, if that's you , kudos to you!
2. A couple, partnerships or best Friends (usually Spoils in, or before the end (whatever that means)  you know,  as nothing stands still on its own for good- hate to spoil that on you!) - but sb gotta spill the Beans, right ;) -- sb: somebody
3. 3 ot more: fits the bill for  family-related  members  - [Time: for either long-term or short-term]

3.Liabilities: Debt,  shares US ()(intl.) - how to say in the UK???(if you're authorized to issue that ), given by other people willingly (hopefully they know what they're doing, right) i.e.  Accounts payable A/P  , me: the 'pay me back later, okay?' accounts  - [either long-term or short-term]

together, they make up the following Golden (brown) Equation:

                Assets = Liabilites + Capital
"""
# the Accounting Equation , :
## Assets = Liabilites + Capital
## ---
#existential question the function is already equal to its absolute value -talk about a Self-reference !
#global capital(totalCapital) =abs(totalCapital);

 asset(totalAssets) = totalLiabilities  # no-abs-allowed #TODO:sanity-check
 liability(totalLiabilities) = totalLiabilities # no-abs-allowed #TODO:sanity-check

#nothing, yet - chop chop, gotta go, are you there yet, where you at!

 #syntax: incomplete: "module" at /Users/ahmadlutfi/Desktop/src/src/accountingTrinity.jl:3 requires end

#TODO: working input ref ... # DONE!
capital = passmissing(input("please insert total capital"))  # insert a value only
liability = passmissing( input("please insert total liability"))#throws(ERROR) : redefinition of Constant Liability #TODO:sanity-check
asset = passmissing(input("please insert total asset"))

print("capital:   $capital") # TODO: Ignore Alpha words -- REGEX  - maybe later -- I guessed it
#@assert capital >= 0 #doesn;t matter  # capital can take both positive & Negative values

#export calcCapital, capital,liability
#=
function calcCapital(totalAssets,totalLiabilities)
asset = abs(totalAssets);  liability = abs(totalLiabilites) #edit: I don't recall writing those , conciously -wow
=#
#---
function calcCapital(asset, liability)
#return max(asset, liability) - min(asset, liability) #that's another thing I've been thinking about - to add or not the add -- that is the Question
"""
capital = asset - liability
"""
capital = asset - liability # for some reason i have to check with liability ... #TODO:
return capital # Capital

#---
function calcAsset(capital, liability)
"""

after subtracting all debt (personal debt included)

asset = -  liab
capital = asset  - liability

"""
# asset = asset < 0 : asset = 0 # if any account is less than 0 , adds it as 0 : Reason:  a cadlag process cannot be less than a  0
asset = capital - liability
# capital = asset - liability
return asset  # capital
end

function calcLiability(capital, asset)
"""
Liability Could be in a form of bonds
Liability (Debt) = (short-term) Debt + (Long-term) Debt
"""
    liability  =   capital - asset
    return liability

    end

end

#--- Asset Calculation
#capital Calculation ?

#TODO:calculate Asset  = capital - liability  #Done!
#"""user inserts  totalCapital"""

#TODO:finish
"""
function calcAssets(totalCapital, totalLiabilities)
capital(totalCapital) = tryparse(:Number, input("please insert total Capital"))  # capital(totalCapital)
liability(totalLiabilities) = tryparse(:Number, input("please insert total Liabilities")) #abs(totalLiabilities)
#TODO:
            return max(capital, liability) - min(capital, liability) #
return capital - liability # = Asset
end
end
"""

"
Capital account v.s liability account

Credit / debit   account algebra :

Cr*Cr = Dr

Dr*Dr = Dr

(-): credit
(+): debit

Dr: Debit
Cr: Credit

liability account decrease could be from
paying off Accounts Payable i.e. (-) Accounts Payable A/C [Credit Account] : (-)credit Accounts Payable A/C,
decreasing(-) Accounts Payable A/C - Credit Account #whenever you Credit a   (-)credit account, that means you are Debiting (decreasing, here) the Credit Account i.e. Cr*Cr = Dr,
adding(+) more Accounts Recievable (+) Accounts Recievable A/R [Debit Account] : (+)debit, Accounts Recieveble A/C - Debit Account (+)is debit #whenever you Debit a Debit account, means you
i.e. some one is Paying you off Dr*Dr = Dr

(-) Accounts Payable A/C

"
function calcAsset(totalCapital, totalLiabilities)

return totalCapital - totalLiabilites

end

#--- prompt Area
#=
using Parsers
 function baseline(s,delimiter=' ')
           pieces = split(s, delimiter, keepempty=false)
           map(pieces) do piece
               parse(Float64, piece)
           end
       end

read  = readline()

 #p = Parsers.DELIMITED(" ", ignorerepeated=true);
=#
#--- ....
#= TODO:READ - Someday, ok - getback -- on it
function _sth parse(::Type{T}, s::AbstractString) where T
    y = tryparse(T, s)
    y === nothing && throw(ArgumentError("Cannot parse as $T: \"$s\""))
    return y
end
=#
#--- Checking region ...
# TODO:  prompt #Done @ 6:15 pm
function prompt(message="How are you feeling?") #Working: #getback: prompt's 2 methods prompt() & prompt(message)

print(message*"\n")
line = readline() #read # ok # string
# splits on spaces and casts to T
#result = produce(map(x -> parse(T, x), split(line)))
# return #result
return parse(Complex,line) #ok
end #Done ! #Correct #checked #prompt(message)
prompt() #gets lost: reason(type input variable)- i.e. is it Float64 ? --unsure: the most existential questionyou might ask
#oh my handling complex vars is complex (in itself!)
#what type of work around & remedyyou'll do
#don't know what type of input user may enter  either string or number
# if number, will it be complex - Problem
#but. let  it be a  Float64 - no problem!

using Test
global digits= 10
global datatype = nothing
function inputHandling(T)
 #--- test Region

    for sign in ('-','+'), Im in ("i","j","im"), s1 in (""," "), s2 in (""," "), s3 in (""," "), s4 in (""," ")
         for r in (1,0,-1), i in (1,0,-1),
             n = Complex(r, sign == '+' ? i : -i)
             s = string(s1, r, s2, sign, s3, i, Im, s4)
             @test n === parse(Complex{Int}, s)
             @test Complex(r) === parse(Complex{Int}, string(s1, r, s2))
             @test Complex(0,i) === parse(Complex{Int}, string(s3, i, Im, s4))
             for T in (Float64, BigFloat)
                 nT = parse(Complex{T}, s)
                 @test nT isa Complex{T}
                 @test nT == n
                 @test n == parse(Complex{T}, string(s1, r, ".0", s2, sign, s3, i, ".0", Im, s4))
                 @test n*parse(T,"1e-3") == parse(Complex{T}, string(s1, r, "e-3", s2, sign, s3, i, "e-3", Im, s4))
             end

         end
         for r in (-1.0,-1e-9,Inf,-Inf,NaN), i in (-1.0,-1e-9,Inf,NaN)
    n = Complex(r, sign == '+' ? i : -i)
    s = lowercase(string(s1, r, s2, sign, s3, i, Im, s4))
    @test n === parse(ComplexF64, s)
    @test Complex(r) === parse(ComplexF64, string(s1, r, s2))
    @test Complex(0,i) === parse(ComplexF64, string(s3, i, Im, s4))
end
end
end
for T in (Int, Float64), bad in ("3 + 4*im", "3 + 4", "1+2ij", "1im-3im", "++4im")
    @test_throws ArgumentError parse(Complex{T}, bad)
end
@test_throws ArgumentError parse(Complex{Int}, "3 + 4.2im")
#test passed

#--- T Input  Handling #Re-Opened!

#TODO: sanity check code below
#TODO: isa Useful ?
if(typeof(T) isa Int)
datatype = Int
@inferred parse(Int, "$digits") #throws(error): @inferred not defined
# end
elseif (typeof(T) isa Float64)
    datatype = Float64
@inferred parse(Float64, "$digits") #throws(error): @inferred not defined
#end
elseif(typeof(T) isa Complex{Int})
    datatype = Complex{Int}
     @inferred parse(Complex{Int}, "$digits") #throws(error): @inferred not defined
end

#you're in my head, Always....

max_stars_repo_name: adamwillisXanax/SolvencyPredictor

max_issues_repo_name: adamwillisXanax/SolvencyPredictor

max_forks_repo_name: adamwillisXanax/SolvencyPredictor


text:
# host array

export oneArray

@enum ArrayState begin
  ARRAY_UNMANAGED
  ARRAY_MANAGED
  ARRAY_FREED
end

mutable struct oneArray{T,N} <: AbstractGPUArray{T,N}
  ptr::ZePtr{Nothing}
  dims::Dims{N}

  state::ArrayState

  ctx::ZeContext
  dev::ZeDevice

  function oneArray{T,N}(::UndefInitializer, dims::Dims{N}) where {T,N}
    Base.isbitsunion(T) && error("oneArray does not yet support union bits types")
    Base.isbitstype(T)  || error("oneArray only supports bits types") # allocatedinline on 1.3+
    ctx = context()
    dev = device()
    ptr = allocate(ctx, dev, prod(dims) * sizeof(T), Base.datatype_alignment(T))
    obj = new{T,N}(ptr, dims, ARRAY_MANAGED, ctx, dev)
    finalizer(unsafe_free!, obj)
    return obj
  end
end

function unsafe_free!(xs::oneArray)
  # this call should only have an effect once, becuase both the user and the GC can call it
  if xs.state == ARRAY_FREED
    return
  elseif xs.state == ARRAY_UNMANAGED
    throw(ArgumentError("Cannot free an unmanaged buffer."))
  end

  release(context(xs), device(xs), xs.ptr)
  xs.state = ARRAY_FREED

  # the object is dead, so we can also wipe the pointer
  xs.ptr = ZE_NULL

  return
end

device(a::oneArray) = a.dev
context(a::oneArray) = a.ctx

## alias detection

Base.dataids(A::oneArray) = (UInt(pointer(A)),)

Base.unaliascopy(A::oneArray) = copy(A)

## convenience constructors

oneVector{T} = oneArray{T,1}
oneMatrix{T} = oneArray{T,2}
oneVecOrMat{T} = Union{oneVector{T},oneMatrix{T}}

# type and dimensionality specified, accepting dims as series of Ints
oneArray{T,N}(::UndefInitializer, dims::Integer...) where {T,N} = oneArray{T,N}(undef, dims)

# type but not dimensionality specified
oneArray{T}(::UndefInitializer, dims::Dims{N}) where {T,N} = oneArray{T,N}(undef, dims)
oneArray{T}(::UndefInitializer, dims::Integer...) where {T} =
    oneArray{T}(undef, convert(Tuple{Vararg{Int}}, dims))

# empty vector constructor
oneArray{T,1}() where {T} = oneArray{T,1}(undef, 0)

# do-block constructors
for (ctor, tvars) in (:oneArray => (), :(oneArray{T}) => (:T,), :(oneArray{T,N}) => (:T, :N))
  @eval begin
    function $ctor(f::Function, args...) where {$(tvars...)}
      xs = $ctor(args...)
      try
        f(xs)
      finally
        unsafe_free!(xs)
      end
    end
  end
end

Base.similar(a::oneArray{T,N}) where {T,N} = oneArray{T,N}(undef, size(a))
Base.similar(a::oneArray{T}, dims::Base.Dims{N}) where {T,N} = oneArray{T,N}(undef, dims)
Base.similar(a::oneArray, ::Type{T}, dims::Base.Dims{N}) where {T,N} = oneArray{T,N}(undef, dims)

function Base.copy(a::oneArray{T,N}) where {T,N}
  b = similar(a)
  @inbounds copyto!(b, a)
end

## array interface

Base.elsize(::Type{<:oneArray{T}}) where {T} = sizeof(T)

Base.size(x::oneArray) = x.dims
Base.sizeof(x::oneArray) = Base.elsize(x) * length(x)

## derived types

export oneDenseArray, oneDenseVector, oneDenseMatrix, oneDenseVecOrMat,
       oneDenseArray, oneDenseVector, oneDenseMatrix, oneDenseVecOrMat,
       oneWrappedArray, oneWrappedVector, oneWrappedMatrix, oneWrappedVecOrMat

oneContiguousSubArray{T,N,A<:oneArray} = Base.FastContiguousSubArray{T,N,A}

# dense arrays: stored contiguously in memory
oneDenseReinterpretArray{T,N,A<:Union{oneArray,oneContiguousSubArray}} = Base.ReinterpretArray{T,N,S,A} where S
oneDenseReshapedArray{T,N,A<:Union{oneArray,oneContiguousSubArray,oneDenseReinterpretArray}} = Base.ReshapedArray{T,N,A}
DenseSuboneArray{T,N,A<:Union{oneArray,oneDenseReshapedArray,oneDenseReinterpretArray}} = Base.FastContiguousSubArray{T,N,A}
oneDenseArray{T,N} = Union{oneArray{T,N}, DenseSuboneArray{T,N}, oneDenseReshapedArray{T,N}, oneDenseReinterpretArray{T,N}}
oneDenseVector{T} = oneDenseArray{T,1}
oneDenseMatrix{T} = oneDenseArray{T,2}
oneDenseVecOrMat{T} = Union{oneDenseVector{T}, oneDenseMatrix{T}}

# strided arrays
oneStridedSubArray{T,N,A<:Union{oneArray,oneDenseReshapedArray,oneDenseReinterpretArray},
                  I<:Tuple{Vararg{Union{Base.RangeIndex, Base.ReshapedUnitRange,
                                        Base.AbstractCartesianIndex}}}} = SubArray{T,N,A,I}
oneStridedArray{T,N} = Union{oneArray{T,N}, oneStridedSubArray{T,N}, oneDenseReshapedArray{T,N}, oneDenseReinterpretArray{T,N}}
oneStridedVector{T} = oneStridedArray{T,1}
oneStridedMatrix{T} = oneStridedArray{T,2}
oneStridedVecOrMat{T} = Union{oneStridedVector{T}, oneStridedMatrix{T}}

Base.pointer(x::oneStridedArray{T}) where {T} = Base.unsafe_convert(ZePtr{T}, x)
@inline function Base.pointer(x::oneStridedArray{T}, i::Integer) where T
    Base.unsafe_convert(ZePtr{T}, x) + Base._memory_offset(x, i)
end

# wrapped arrays: can be used in kernels
oneWrappedArray{T,N} = Union{oneArray{T,N}, WrappedArray{T,N,oneArray,oneArray{T,N}}}
oneWrappedVector{T} = oneWrappedArray{T,1}
oneWrappedMatrix{T} = oneWrappedArray{T,2}
oneWrappedVecOrMat{T} = Union{oneWrappedVector{T}, oneWrappedMatrix{T}}

## interop with other arrays

@inline function oneArray{T,N}(xs::AbstractArray{<:Any,N}) where {T,N}
  A = oneArray{T,N}(undef, size(xs))
  copyto!(A, convert(Array{T}, xs))
  return A
end

# underspecified constructors
oneArray{T}(xs::AbstractArray{S,N}) where {T,N,S} = oneArray{T,N}(xs)
(::Type{oneArray{T,N} where T})(x::AbstractArray{S,N}) where {S,N} = oneArray{S,N}(x)
oneArray(A::AbstractArray{T,N}) where {T,N} = oneArray{T,N}(A)

# idempotency
oneArray{T,N}(xs::oneArray{T,N}) where {T,N} = xs

# Level Zero references
oneL0.ZeRef(x::Any) = oneL0.ZeRefArray(oneArray([x]))
oneL0.ZeRef{T}(x) where {T} = oneL0.ZeRefArray{T}(oneArray(T[x]))
oneL0.ZeRef{T}() where {T} = oneL0.ZeRefArray(oneArray{T}(undef, 1))

## conversions

Base.convert(::Type{T}, x::T) where T <: oneArray = x

## interop with C libraries

Base.unsafe_convert(::Type{Ptr{T}}, x::oneArray{T}) where {T} =
  throw(ArgumentError("cannot take the host address of a $(typeof(x))"))
Base.unsafe_convert(::Type{ZePtr{T}}, x::oneArray{T}) where {T} = convert(ZePtr{T}, x.ptr)

## interop with GPU arrays

function Base.unsafe_convert(::Type{oneDeviceArray{T,N,AS.Global}}, a::oneDenseArray{T,N}) where {T,N}
  oneDeviceArray{T,N,AS.Global}(size(a), reinterpret(LLVMPtr{T,AS.Global}, pointer(a)))
end

Adapt.adapt_storage(::KernelAdaptor, xs::oneArray{T,N}) where {T,N} =
  Base.unsafe_convert(oneDeviceArray{T,N,AS.Global}, xs)

# we materialize ReshapedArray/ReinterpretArray/SubArray/... directly as a device array
Adapt.adapt_structure(::KernelAdaptor, xs::oneDenseArray{T,N}) where {T,N} =
  Base.unsafe_convert(oneDeviceArray{T,N,AS.Global}, xs)

## interop with CPU arrays

# We don't convert isbits types in `adapt`, since they are already
# considered GPU-compatible.

Adapt.adapt_storage(::Type{oneArray}, xs::AbstractArray) =
  isbits(xs) ? xs : convert(oneArray, xs)

# if an element type is specified, convert to it
Adapt.adapt_storage(::Type{<:oneArray{T}}, xs::AbstractArray) where {T} =
  isbits(xs) ? xs : convert(oneArray{T}, xs)

Adapt.adapt_storage(::Type{Array}, xs::oneArray) = convert(Array, xs)

Base.collect(x::oneArray{T,N}) where {T,N} = copyto!(Array{T,N}(undef, size(x)), x)

function Base.copyto!(dest::oneArray{T}, doffs::Integer, src::Array{T}, soffs::Integer,
                      n::Integer) where T
  n==0 && return dest
  @boundscheck checkbounds(dest, doffs)
  @boundscheck checkbounds(dest, doffs+n-1)
  @boundscheck checkbounds(src, soffs)
  @boundscheck checkbounds(src, soffs+n-1)
  unsafe_copyto!(context(dest), device(dest), dest, doffs, src, soffs, n)
  return dest
end

Base.copyto!(dest::oneDenseArray{T}, src::Array{T}) where {T} =
    copyto!(dest, 1, src, 1, length(src))

function Base.copyto!(dest::Array{T}, doffs::Integer, src::oneDenseArray{T}, soffs::Integer,
                      n::Integer) where T
  n==0 && return dest
  @boundscheck checkbounds(dest, doffs)
  @boundscheck checkbounds(dest, doffs+n-1)
  @boundscheck checkbounds(src, soffs)
  @boundscheck checkbounds(src, soffs+n-1)
  unsafe_copyto!(context(src), device(src), dest, doffs, src, soffs, n)
  return dest
end

Base.copyto!(dest::Array{T}, src::oneDenseArray{T}) where {T} =
    copyto!(dest, 1, src, 1, length(src))

function Base.copyto!(dest::oneDenseArray{T}, doffs::Integer, src::oneDenseArray{T}, soffs::Integer,
                      n::Integer) where T
  n==0 && return dest
  @boundscheck checkbounds(dest, doffs)
  @boundscheck checkbounds(dest, doffs+n-1)
  @boundscheck checkbounds(src, soffs)
  @boundscheck checkbounds(src, soffs+n-1)
  @assert device(dest) == device(src) && context(dest) == context(src)
  unsafe_copyto!(context(dest), device(dest), dest, doffs, src, soffs, n)
  return dest
end

Base.copyto!(dest::oneDenseArray{T}, src::oneDenseArray{T}) where {T} =
    copyto!(dest, 1, src, 1, length(src))

function Base.unsafe_copyto!(ctx::ZeContext, dev::ZeDevice,
                             dest::oneDenseArray{T}, doffs, src::Array{T}, soffs, n) where T
  GC.@preserve src dest unsafe_copyto!(ctx, dev, pointer(dest, doffs), pointer(src, soffs), n)
  if Base.isbitsunion(T)
    # copy selector bytes
    error("Not implemented")
  end
  return dest
end

function Base.unsafe_copyto!(ctx::ZeContext, dev::ZeDevice,
                             dest::Array{T}, doffs, src::oneDenseArray{T}, soffs, n) where T
  GC.@preserve src dest unsafe_copyto!(ctx, dev, pointer(dest, doffs), pointer(src, soffs), n)
  if Base.isbitsunion(T)
    # copy selector bytes
    error("Not implemented")
  end

  # copies to the host are synchronizing
  synchronize(global_queue(context(src), device(src)))

  return dest
end

function Base.unsafe_copyto!(ctx::ZeContext, dev::ZeDevice,
                             dest::oneDenseArray{T}, doffs, src::oneDenseArray{T}, soffs, n) where T
  GC.@preserve src dest unsafe_copyto!(ctx, dev, pointer(dest, doffs), pointer(src, soffs), n)
  if Base.isbitsunion(T)
    # copy selector bytes
    error("Not implemented")
  end
  return dest
end

## utilities

zeros(T::Type, dims...) = fill!(oneArray{T}(undef, dims...), 0)
ones(T::Type, dims...) = fill!(oneArray{T}(undef, dims...), 1)
zeros(dims...) = zeros(Float64, dims...)
ones(dims...) = ones(Float64, dims...)
fill(v, dims...) = fill!(oneArray{typeof(v)}(undef, dims...), v)
fill(v, dims::Dims) = fill!(oneArray{typeof(v)}(undef, dims...), v)

function Base.fill!(A::oneDenseArray{T}, val) where T
  B = [convert(T, val)]
  unsafe_fill!(context(A), device(A), pointer(A), pointer(B), length(A))
  A
end

## views

@inline function Base.view(A::oneArray, I::Vararg{Any,N}) where {N}
    J = to_indices(A, I)
    @boundscheck begin
        # Base's boundscheck accesses the indices, so make sure they reside on the CPU.
        # this is expensive, but it's a bounds check after all.
        J_cpu = map(j->adapt(Array, j), J)
        checkbounds(A, J_cpu...)
    end
    J_gpu = map(j->adapt(oneArray, j), J)
    Base.unsafe_view(Base._maybe_reshape_parent(A, Base.index_ndims(J_gpu...)), J_gpu...)
end

device(a::SubArray) = device(parent(a))
context(a::SubArray) = context(parent(a))

# contiguous subarrays
function Base.unsafe_convert(::Type{ZePtr{T}}, V::SubArray{T,N,P,<:Tuple{Vararg{Base.RangeIndex}}}) where {T,N,P}
    return Base.unsafe_convert(ZePtr{T}, parent(V)) +
           Base._memory_offset(V.parent, map(first, V.indices)...)
end

# reshaped subarrays
function Base.unsafe_convert(::Type{ZePtr{T}}, V::SubArray{T,N,P,<:Tuple{Vararg{Union{Base.RangeIndex,Base.ReshapedUnitRange}}}}) where {T,N,P}
   return Base.unsafe_convert(ZePtr{T}, parent(V)) +
          (Base.first_index(V)-1)*sizeof(T)
end

## reshape

device(a::Base.ReshapedArray) = device(parent(a))
context(a::Base.ReshapedArray) = context(parent(a))

Base.unsafe_convert(::Type{ZePtr{T}}, a::Base.ReshapedArray{T}) where {T} =
  Base.unsafe_convert(ZePtr{T}, parent(a))

## reinterpret

device(a::Base.ReinterpretArray) = device(parent(a))
context(a::Base.ReinterpretArray) = context(parent(a))

Base.unsafe_convert(::Type{ZePtr{T}}, a::Base.ReinterpretArray{T,N,S} where N) where {T,S} =
  ZePtr{T}(Base.unsafe_convert(ZePtr{S}, parent(a)))


text:
# Infrastructure for value transformation

"""
    clamp01(x) -> y

Produce a value `y` that lies between 0 and 1, and equal to `x` when
`x` is already in this range. Equivalent to `clamp(x, 0, 1)` for
numeric values. For colors, this function is applied to each color
channel separately.

See also: [`clamp01!`](@ref), [`clamp01nan`](@ref).
"""
clamp01(x::Union{N0f8,N0f16}) = x
clamp01(x::Number) = clamp(x, zero(x), oneunit(x))
clamp01(c::Colorant) = mapc(clamp01, c)

"""
    clamp01!(array::AbstractArray)

Restrict values in array to [0, 1], in-place. See also [`clamp01`](@ref).
"""
function clamp01!(img::AbstractArray)
    # slgihtly faster than map!(clamp01, img, img)
    @inbounds for i in eachindex(img)
        img[i] = clamp01(img[i])
    end
    return img
end

"""
    clamp01nan(x) -> y

Similar to `clamp01`, except that any `NaN` values are changed to 0.

See also: [`clamp01nan!`](@ref), [`clamp01`](@ref).
"""
clamp01nan(x) = clamp01(x)
clamp01nan(x::AbstractFloat) = ifelse(isnan(x), zero(x), clamp01(x))
clamp01nan(c::Colorant) = mapc(clamp01nan, c)

"""
    clamp01nan!(array::AbstractArray)

Similar to `clamp01!`, except that any `NaN` values are changed to 0.

See also: [`clamp01!`](@ref), [`clamp01nan`](@ref)
"""
function clamp01nan!(img::GenericImage)
    # slgihtly faster than map!(clamp01nan, img, img)
    @inbounds for i in eachindex(img)
        img[i] = clamp01nan(img[i])
    end
    return img
end

"""
    scale_minmax(min, max) -> f
    scale_minmax(T, min, max) -> f

Return a function `f` which maps values less than or equal to `min` to
0, values greater than or equal to `max` to 1, and uses a linear scale
in between. `min` and `max` should be real values.

Optionally specify the return type `T`. If `T` is a colorant (e.g.,
RGB), then scaling is applied to each color channel.

# Examples

## Example 1

```julia
julia> f = scale_minmax(-10, 10)
(::#9) (generic function with 1 method)

julia> f(10)
1.0

julia> f(-10)
0.0

julia> f(5)
0.75
```

## Example 2

```julia
julia> c = RGB(255.0,128.0,0.0)
RGB{Float64}(255.0,128.0,0.0)

julia> f = scale_minmax(RGB, 0, 255)
(::#13) (generic function with 1 method)

julia> f(c)
RGB{Float64}(1.0,0.5019607843137255,0.0)
```

See also: [`take_map`](@ref).
"""
function scale_minmax(min::T, max::T) where T
    @inline function(x)
        xp, minp, maxp = promote(x, min, max)  # improves performance to promote explicitly
        y = clamp(xp, minp, maxp)
        (y-minp)/(maxp-minp)
    end
end
function scale_minmax(::Type{Tout}, min::T, max::T) where {Tout,T}
    @inline function(x)
        xp, minp, maxp = promote(x, min, max)
        y = clamp(xp, minp, maxp)
        smmconvert(Tout, (y-minp)/(maxp-minp))
    end
end
@inline smmconvert(::Type{Tout}, x) where {Tout} = convert(Tout, x)
# since we know the result will be between 0 and 1, we can use rem to save a check
@inline smmconvert(::Type{Tout}, x) where {Tout<:Normed} = rem(x, Tout)

scale_minmax(min, max) = scale_minmax(promote(min, max)...)
scale_minmax(::Type{T}, min, max) where {T} = scale_minmax(T, promote(min, max)...)

# TODO: use triangular dispatch when we can count on Julia 0.6+
function scale_minmax(::Type{C}, min::T, max::T) where {C<:Colorant, T<:Real}
    return _scale_minmax(C, eltype(C), min, max)
end
function _scale_minmax(::Type{C}, ::Type{TC}, min::T, max::T) where {C<:Colorant, TC<:Real, T<:Real}
    freal = scale_minmax(TC, min, max)
    @inline function(c)
        C(mapc(freal, c))
    end
end
function _scale_minmax(::Type{C}, ::Type{Any}, min::T, max::T) where {C<:Colorant, T<:Real}
    freal = scale_minmax(min, max)
    @inline function(c)
        C(mapc(freal, c))
    end
end

"""
    scale_signed(maxabs) -> f

Return a function `f` which scales values in the range `[-maxabs,
maxabs]` (clamping values that lie outside this range) to the range
`[-1, 1]`.

See also: [`color_signed`](@ref).
"""
function scale_signed(maxabs::Real)
    maxabs > 0 || throw(ArgumentError("maxabs must be positive, got $maxabs"))
    return x -> clamp(x, -maxabs, maxabs)/maxabs
end

"""
    scale_signed(min, center, max) -> f

Return a function `f` which scales values in the range `[min, center]`
to `[-1,0]` and `[center,max]` to `[0,1]`. Values smaller than
`min`/`max` get clamped to `min`/`max`, respectively.

See also: [`color_signed`](@ref).
"""
function scale_signed(min::T, center::T, max::T) where T<:Real
    min <= center <= max || throw(ArgumentError("values must be ordered, got $min, $center, $max"))
    sneg, spos = 1/(center-min), 1/(max-center)
    function(x)
        Δy = clamp(x, min, max) - center
        ifelse(Δy < 0, sneg*Δy, spos*Δy)
    end
end
function scale_signed(min::Real, center::Real, max::Real)
    return scale_signed(promote(min, center, max)...)
end

"""
    color_signed()
    color_signed(colorneg, colorpos) -> f
    color_signed(colorneg, colorcenter, colorpos) -> f

Define a function that maps negative values (in the range [-1,0]) to
the linear colormap between `colorneg` and `colorcenter`, and positive
values (in the range [0,1]) to the linear colormap between
`colorcenter` and `colorpos`.

The default colors are:

- `colorcenter`: white
- `colorneg`: green1
- `colorpos`: magenta

See also: [`scale_signed`](@ref).
"""
color_signed(neg::C, center::C, pos::C) where {C<:Colorant} = function(x)
    y = clamp(x, -oneunit(x), oneunit(x))
    yabs = abs(y)
    C(ifelse(y>0, weighted_color_mean(yabs, pos, center),
                  weighted_color_mean(yabs, neg, center)))
end

function color_signed(colorneg::C, colorpos::C) where C<:Colorant
    return color_signed(colorneg, C(colorant"white"), colorpos)
end

color_signed() = color_signed(colorant"green1", colorant"magenta")

"""
    take_map(f, A) -> fnew
    take_map(f, T, A) -> fnew

Given a value-mapping function `f` and an array `A`, return a
"concrete" mapping function `fnew`. When applied to elements of `A`,
`fnew` should return valid values for storage or display, for example
in the range from 0 to 1 (for grayscale) or valid colorants. `fnew`
may be adapted to the actual values present in `A`, and may not
produce valid values for any inputs not in `A`.

Optionally one can specify the output type `T` that `fnew` should produce.

# Example:
```julia
julia> A = [0, 1, 1000];

julia> f = take_map(scale_minmax, A)
(::#7) (generic function with 1 method)

julia> f.(A)
3-element Array{Float64,1}:
 0.0
 0.001
 1.0
```
"""
take_map

take_map(f, A) = f
take_map(f, ::Type{T}, A) where {T} = x->T(f(x))

function take_map(::typeof(scale_minmax), A::AbstractArray{T}) where T<:Real
    min, max = extrema(A)
    scale_minmax(min, max)
end
function take_map(::typeof(scale_minmax), A::AbstractArray{C}) where C<:Colorant
    min, max = extrema(channel_view(A))
    scale_minmax(C, min, max)
end
function take_map(::typeof(scale_minmax), ::Type{Tout}, A::AbstractArray{T}) where {Tout,T<:Real}
    min, max = extrema(A)
    scale_minmax(Tout, min, max)
end
function take_map(::typeof(scale_minmax), ::Type{Cout}, A::AbstractArray{C}) where {Cout<:Colorant,C<:Colorant}
    min, max = extrema(channel_view(A))
    scale_minmax(Cout, min, max)
end

function take_map(::typeof(scale_signed), A::AbstractArray{T}) where T<:Real
    mn, mx = extrema(A)
    scale_signed(max(abs(mn), abs(mx)))
end


text:
module PopPK

using DataFrames: index
using Base: Integer, String
using Serialization
using Random
using Pumas

include("pre.jl")
include("post.jl")
include("readers.jl")

green = "#91BE33"
turqoise = "#18C8D0"
blue = "#178CCB"
orange = "#F48024"

end # module


text:
using DirectTrajectoryOptimization
const DTO = DirectTrajectoryOptimization

function point_foot_quadruped_dyn(model, env, h, y, x, u, w)

    # dimensions
    nq = model.nq
    nu = model.nu

    # configurations

    q1⁻ = x[1:nq]
    q2⁻ = x[nq .+ (1:nq)]
    q2⁺ = y[1:nq]
    q3⁺ = y[nq .+ (1:nq)]

    # control
    u_control = u[1:nu]
    γ = u[nu .+ (1:4)]
    β = u[nu + 4 .+ (1:16)]

    E = [1.0 0.0 -1.0 0.0;
         0.0 1.0 0.0 -1.0] # friction mapping
    J = J_func(model, env, q2⁺)
    λ = transpose(J[1:12, :]) * [
                        [E * β[0  .+ (1:4)]; γ[1]];
                        [E * β[4  .+ (1:4)]; γ[2]];
                        [E * β[8  .+ (1:4)]; γ[3]];
                        [E * β[12 .+ (1:4)]; γ[4]]
                       ]
    [
     q2⁺ - q2⁻;
     dynamics(model, h, q1⁻, q2⁺, u_control, zeros(model.nw), λ, q3⁺)
    ]
end

function point_foot_quadruped_dyn1(model, env, h, y, x, u, w)
    nx = 2 * model.nq
    [
     point_foot_quadruped_dyn(model, env, h, y, x, u, w);
     y[nx .+ (1:5)] - [u[model.nu .+ (1:4)]; u[end]];
     y[nx + 5 .+ (1:nx)] - x[1:nx];
     y[nx + 5 + nx .+ (1:model.nu)] - u[1:model.nu];
    ]
end

function point_foot_quadruped_dynt(model, env, h, y, x, u, w)
    nx = 2 * model.nq
    [
     point_foot_quadruped_dyn(model, env, h, y, x, u, w);
     y[nx .+ (1:5)] - [u[model.nu .+ (1:4)]; u[end]];
     y[nx + 5 .+ (1:nx)] - x[nx + 5 .+ (1:nx)];
     y[nx + 5 + nx .+ (1:model.nu)] - u[1:model.nu];
    ]
end

function contact_constraints_inequality_1(model, env, h, x, u, w)
    nq = model.nq
    nu = model.nu
    nx = 2nq

    q2 = x[1:nq]
    q3 = x[nq .+ (1:nq)]

    u_control = u[1:nu]
    γ = u[nu .+ (1:4)]
    β = u[nu + 4 .+ (1:16)]
    ψ = u[nu + 4 + 16 .+ (1:4)]
    η = u[nu + 4 + 16 + 4 .+ (1:16)]
    sα = u[nu + 4 + 16 + 4 + 16 .+ (1:1)]

    ϕ = ϕ_func(model, env, q3)[1:4]

    μ = model.μ_world
    fc = μ .* γ[1:4] - [sum(β[0 .+ (1:4)]); sum(β[4 .+ (1:4)]); sum(β[8 .+ (1:4)]); sum(β[12 .+ (1:4)]);]

    [
     -ϕ;
     -fc;
     β .* η .- sα;
     ψ .* fc  .- sα;
    ]
end

function contact_constraints_inequality_t(model, env, h, x, u, w)
    nq = model.nq
    nu = model.nu
    nx = 2nq

    q2 = x[1:nq]
    q3 = x[nq .+ (1:nq)]

    u_control = u[1:nu]
    γ = u[nu .+ (1:4)]
    β = u[nu + 4 .+ (1:16)]
    ψ = u[nu + 4 + 16 .+ (1:4)]
    η = u[nu + 4 + 16 + 4 .+ (1:16)]
    sα = u[nu + 4 + 16 + 4 + 16 .+ (1:1)]

    ϕ = ϕ_func(model, env, q3)[1:4]
    γ⁻ = x[nx .+ (1:4)]
    sα⁻ = x[nx + 4 .+ (1:1)]

    μ = model.μ_world
    fc = μ .* γ[1:4] - [sum(β[0 .+ (1:4)]); sum(β[4 .+ (1:4)]); sum(β[8 .+ (1:4)]); sum(β[12 .+ (1:4)]);]

    [
     -ϕ;
     -fc;
     γ⁻ .* ϕ .- sα⁻;
     β .* η .- sα;
     ψ .* fc  .- sα;
    ]
end

function contact_constraints_inequality_T(model, env, h, x, u, w)
    nq = model.nq
    nx = 2nq

    q2 = x[1:nq]
    q3 = x[nq .+ (1:nq)]

    ϕ = ϕ_func(model, env, q3)[1:4]
    γ⁻ = x[nx .+ (1:4)]
    sα⁻ = x[nx + 4 .+ (1:1)]

    [
     -ϕ;
     γ⁻ .* ϕ .- sα⁻;
    ]
end

function contact_constraints_equality(model, env, h, x, u, w)
    nq = model.nq
    nu = model.nu

    q2 = x[1:nq]
    q3 = x[nq .+ (1:nq)]

    γ = u[nu .+ (1:4)]
    β = u[nu + 4 .+ (1:16)]
    ψ = u[nu + 4 + 16 .+ (1:4)]
    η = u[nu + 4 + 16 + 4 .+ (1:16)]
    sα = u[nu + 4 + 16 + 4 + 16 .+ (1:1)]

    E = [1.0 0.0 -1.0 0.0;
         0.0 1.0 0.0 -1.0]
    v = (q3 - q2) ./ h[1]
    vT = vcat([E' * v[(i-1) * 3 .+ (1:2)] for i = 1:4]...)

    ψ_stack = vcat([ψi * ones(4) for ψi in ψ]...)

    [
     η - vT - ψ_stack;
    ]
end


text:
module TestExamplesDGMulti

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "dgmulti_2d")

@testset "DGMulti 2D" begin
  @trixi_testset "elixir_euler_weakform.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4,4),
      l2 = [0.0013463253573454783, 0.0014235911638070984, 0.0014235911638076622, 0.004721923810347086],
      linf = [0.0015248269221803668, 0.002070690855385582, 0.002070690855385804, 0.004913338290754687]
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (SBP)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4,4),
      approximation_type = SBP(),
      l2 = [0.007465048853079056, 0.005297478943323888, 0.005297478943334149, 0.014701766352428895],
      linf = [0.02149035033755564, 0.0135293837651278, 0.013529383764659286, 0.03269081948752639]
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (Quadrilateral elements)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4,4),
      element_type = Quad(),
      l2 = [0.0002909691660845978, 0.0002811425883546657, 0.0002811425883549579, 0.0010200600240538172],
      linf = [0.0004970396373780162, 0.0004059109438805386, 0.00040591094388231497, 0.0014247618507141624]
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (EC) " begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4,4),
      volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha),
      surface_integral = SurfaceIntegralWeakForm(flux_ranocha),
      l2 = [0.008107539211003132, 0.007464472445092778, 0.007464472445093055, 0.01597648138530006],
      linf = [0.012298218434060981, 0.013874789519390918, 0.013874789519420005, 0.040393065744379175]
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (SBP, EC)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4,4),
      volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha),
      surface_integral = SurfaceIntegralWeakForm(flux_ranocha),
      approximation_type = SBP(),
      l2 = [0.012858228819248307, 0.010649745431713896, 0.010649745431680024, 0.02628727578633061],
      linf = [0.03733928157930677, 0.053088127555369624, 0.05308812755616854, 0.13379093830601718]
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (Quadrilateral elements, SBP, EC)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4,4),
      element_type = Quad(),
      volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha),
      surface_integral = SurfaceIntegralWeakForm(flux_ranocha),
      approximation_type = SBP(),
      l2 = [0.0029320668218228204, 0.0030626480847183927, 0.0030626480847164876, 0.0068440872356226165],
      linf = [0.013699485223367613, 0.012808498911205168, 0.01280849891121405, 0.022408991793810618]
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (convergence)" begin
    mean_convergence = convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"), 2)
    @test isapprox(mean_convergence[:l2], [4.2498632232077815, 4.133717042428824, 4.133717042041395, 4.086223177072245], rtol=0.05)
  end

  @trixi_testset "elixir_euler_weakform_periodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform_periodic.jl"),
      l2 = [0.0014986508075708323, 0.001528523420746786, 0.0015285234207473158, 0.004846505183839211],
      linf = [0.0015062108658376872, 0.0019373508504645365, 0.0019373508504538783, 0.004742686826709086]
    )
  end

  @trixi_testset "elixir_euler_triangulate_pkg_mesh.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_triangulate_pkg_mesh.jl"),
      l2 = [4.664661209491976e-6, 3.7033509525940745e-6, 4.794877426562555e-6, 1.2682723101532175e-5],
      linf = [2.5099852761334418e-5, 2.2683684021362893e-5, 2.6180448559287584e-5, 5.5752932611508044e-5]
    )
  end
end

EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "dgmulti_3d")

@testset "DGMulti 3D" begin

  # 3d tet/hex tests
  @trixi_testset "elixir_euler_weakform.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      l2 = [0.0010029534292051608, 0.0011682205957721673, 0.001072975385793516, 0.000997247778892257, 0.0039364354651358294],
      linf = [0.003660737033303718, 0.005625620600749226, 0.0030566354814669516, 0.0041580358824311325, 0.019326660236036464]
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (EC)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      surface_integral = SurfaceIntegralWeakForm(flux_ranocha),
      volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha),
      l2 = [0.014932088450136542, 0.017080219613061526, 0.016589517840793006, 0.015905000907070196, 0.03903416208587798],
      linf = [0.06856547797256729, 0.08225664880340489, 0.06925055630951782, 0.06913016119820181, 0.19161418499621874]
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (Hexahedral elements)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      element_type = Hex(),
      l2 = [0.00030580190715769566, 0.00040146357607439464, 0.00040146357607564597, 0.000401463576075708, 0.0015749412434154315],
      linf = [0.00036910287847780054, 0.00042659774184228283, 0.0004265977427213574, 0.00042659774250686233, 0.00143803344597071]
    )
  end

  @trixi_testset "elixir_euler_weakform_periodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform_periodic.jl"),
      l2 = [0.0010317074322517949, 0.0012277090547035293, 0.0011273991123913515, 0.0010418496196130177, 0.004058878478404962],
      linf = [0.003227752881827861, 0.005620317864620361, 0.0030514833972379307, 0.003987027618439498, 0.019282224709831652]
    )
  end

  @trixi_testset "elixir_euler_weakform_periodic.jl (Hexahedral elements)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform_periodic.jl"),
      element_type = Hex(),
      l2 = [0.00034230612468547436, 0.00044397204714598747, 0.0004439720471461567, 0.0004439720471464591, 0.0016639410646990126],
      linf = [0.0003674374460325147, 0.0004253921341716982, 0.0004253921340786615, 0.0004253921340831024, 0.0014333414071048267]
    )
  end

end

end # module


text:
include( "../src/argminProb.jl" )
using Test

  println( "running argminProb unit tests." );

  vec = repeat( [2.0], outer=[3 1] )[:,1];
  append!( vec, rand( 4.0:20.0, 100 ) );
  nrRuns = 100000;
  res    = zeros( Float64, 103 );
  for i in 1:1:nrRuns
    index = argminUniProb.argminProb( vec );
    res[index] += 1;
  end
  minVal = minimum( res[1:3] );
  maxVal = maximum( res[1:3] );

  @test maxVal / minVal > 0.95
  @test maxVal / minVal < 1.05

  len        = 100;
  lmatrix    = repeat( [2.0], outer=[len 3] );
  rmatrix    = rand( 4.0:20.0, len, 3 );
  vmatrix    = hcat( lmatrix, rmatrix );
  minIndices = argminUniProb.argminProb( vmatrix );

  @test all( minIndices .>= 1 )
  @test all( minIndices .<= 3 )

  minIndices = zeros( Int64, size( vmatrix, 2 ) );
  argminUniProb.argminProb!( vmatrix, minIndices );

  @test all( minIndices .>= 1 )
  @test all( minIndices .<= 3 )


text:
# include("./plot_utils.jl")
# GPU = GeneralPlotUtils
# LU = LabelUtils
# include("./plot_utils.jl")
# GPU = GeneralPlotUtils
# LU = LabelUtils

module Plotting
# using Plots; gr()
using Plots; pyplot()
using Reproduce
using FileIO
using Statistics
using ProgressMeter
using JLD2
using LaTeXStrings
using StatsPlots
# Plots.scalefontsizes(2)
# import ..GeneralPlotUtils
# import ..LabelUtils
# include("./plot_utils.jl")
include("./plot_utils.jl")
GPU = GeneralPlotUtils
LU = LabelUtils

 # default(titlefont = (20, "times"), legendfontsize = 18, guidefont = (18, :darkgreen), tickfont = (12, :orange), guide = "x", framestyle = :zerolines, yminorgrid = true)
 default(titlefont = (20, "times"), legendfontsize = 10, guidefont = (18, :black), tickfont = 12, foreground_color_legend = nothing)

# data_key = :ttmaze_direct_error
# data_key = :ttmaze_uniform_error
# data_key = :ttmaze_round_robin_error
# data_key = :oned_tmaze_dpi_error
# data_key = :oned_tmaze_old_error
# data_key = :oned_tmaze_dmu_error
# data_key = :mc_uniform_error
# data_key = :mc_starts_error
data_key = :twod_grid_world_error_center_dpi

# folder_name = "oned_control"
# folder_name = "oned_rr"
# folder_name = "GPI_Sensitivity"
# folder_name = "NoLegend"
# folder_name = "tabular_control"
# folder_name = "MC_tmp"
# folder_name = "twod_gpi_tb"
# folder_name = "twod_gpi_fine"
# folder_name = "twod_gpi_no_optimism"
# folder_name = "twod_gpi_simple"
# folder_name = "twod_gpi_nstate"
folder_name = "tmp"

# data_home = "../data/Experiment1"
# data_home = "../data/Experiment2_d_pi"
# data_home = "../data/OneDTMaze_Control"
# data_home = "../data/OneDTMaze_RR"
# data_home = "../data/GPI_Sensitivity"
# data_home = "../data/OneDTMaze_GPI_Sensitivity"
# data_home = "../data/MC_Experiments"
# data_home = "../data/MC_Experiments_Final"
# data_home = "../data/TwoDGridWorld_gpi"
# data_home = "../data/TwoDGridWorld_gpi_TB"
# data_home = "../data/TwoDGridWorld_finegrained"
# data_home = "../data/TwoDGridWorld_gpi_no_optimism_Rev_2"
# data_home = "../data/TwoDGridWorld_gpi_simple"
# data_home = "../data/TwoDGridWorld_gpi_nstate"
# data_home = "../data/TwoD"
# data_home = "../data/TwoDGridWorld_GPI_No_Opt_Final_P2"
# data_home = "../data/TwoDFinal"
# data_home = "../data/TwoDGridWorld_GPI_No_Opt_Final"
data_home = "../data/TwoD_All"

function load_data()
    experiment_folders = [data_home]
    # folder_name = "oned_rr_dpi"
    # folder_name = "oned_control_dpi"
    ic = ItemCollection(joinpath(experiment_folders[1], "data"))
    function f(itm)
        itm.parsed_args["seed"] < 31
    end
    ic = search(f, ic)
    return ic
end

function load_best(ic)
# data_key = :oned_tmaze_dpi_error

    algo_divisor_keys = ["behaviour_learner", "demon_learner", "demon_opt", "demon_update","exploration_param"]
    # algo_divisor_keys = ["behaviour_learner", "demon_learner", "demon_opt", "demon_update","behaviour_reward_projector","behaviour_rp_tilings"]
    # sweep_params = ["demon_alpha_init", "demon_eta", "alpha_init"]

    sweep_params = ["behaviour_eta", "demon_eta", "alpha_init", "eta","emphasis_clip_threshold"]

    algo_specs_full = GPU.split_algo(ic, algo_divisor_keys)

    all_algos_ics = [search(ic,algo_spec) for algo_spec in algo_specs_full]
    @show length.(all_algos_ics)

    valid_algos_ind = findall(a -> length(a) != 0, all_algos_ics)

    algo_specs = algo_specs_full[valid_algos_ind]
    algo_ics = all_algos_ics[valid_algos_ind]
    best_per_algo_ics = []
    for (i,algo_ic) in enumerate(algo_ics)
        push!(best_per_algo_ics, GPU.get_best_final_perf(algo_ic,sweep_params, data_key, 0.1))
        # push!(best_per_algo_ics, GPU.get_best(algo_ic,sweep_params, data_key))
        # push!(best_per_algo_ics, GPU.get_most_episodes(algo_ic,sweep_params, data_key))

    end
    @show length.(best_per_algo_ics)
    [GPU.print_params(bic,["behaviour_learner","alpha_init"],["behaviour_eta","demon_eta","demon_learner","demon_update","emphasis_clip_threshold","exploration_param"]) for bic in best_per_algo_ics]

    return best_per_algo_ics
end

function plot_rmse_per_demon(algo_ics, inds_of_interest = nothing)
    best_per_algo_ics = if inds_of_interest isa Nothing
        deepcopy(algo_ics)
    else
        deepcopy(algo_ics)[inds_of_interest]
    end

    data_per_gvf = [GPU.get_stats(GPU.load_results(ic,data_key), per_gvf = true)[1] for ic in best_per_algo_ics]
    std_per_gvf = [GPU.get_stats(GPU.load_results(ic,data_key), per_gvf = true)[2] for ic in best_per_algo_ics]
    num_gvfs = 4
    ps = []
    gvf_labels = ["Distractor" "Constant" "Drifter" "Constant"]
    num_runs = length.(best_per_algo_ics)[1]
    for gvf_ind in 1:num_gvfs
        p = plot()
        for algo_ind in 1:length(data_per_gvf)
            smooth_gvf = GPU.smooth(data_per_gvf[algo_ind][gvf_ind,:],5)
            # label = LU.get_label(best_per_algo_ics[algo_ind])[:label]
            # label = LU.get_params(best_per_algo_ics[algo_ind])[:label]
            # color = LU.get_params(best_per_algo_ics[algo_ind])[:color]
            plot_params = LU.get_params(best_per_algo_ics[algo_ind])
            plot!(p,smooth_gvf, ribbon = std_per_gvf[algo_ind][gvf_ind,:] / sqrt(num_runs), size = (500,500); plot_params...)
            plot!(xlabel="Steps", ylabel = "RMSE",  title = string(gvf_labels[gvf_ind], ""))
        end
        push!(ps,p)
    end
    allp = plot(ps..., layout=(2,2), size = (1200,1200))
    display(allp)
    savefig("./plots/$(folder_name)/RMSE_per_Demon_$(string(data_key)).png")

end

function plot_rmse(algo_ics, inds_of_interest = nothing)
    best_per_algo_ics = if inds_of_interest isa Nothing
        deepcopy(algo_ics)
    else
        deepcopy(algo_ics)[inds_of_interest]
    end

    data = [GPU.smooth(GPU.get_stats(GPU.load_results(ic,data_key))[1],10) for ic in best_per_algo_ics]
    std = [GPU.smooth(GPU.get_stats(GPU.load_results(ic,data_key))[2],10) for ic in best_per_algo_ics]
    # std = [GPU.get_stats(GPU.load_results(ic,data_key))[2] for ic in best_per_algo_ics]

    ylabel = "RMSE"
    # title = "Hall Following in Tabular T-Maze"
    # title = "Control in Tabular T-Maze"
    # title = "Fixed Behaviour in 1D T-Maze"
    # title = "Control in 1D T-Maze"
    title = ""
    tick_scaling = 1000
    xlabel = string("Steps", " (per $(tick_scaling))")

    step_increment = if "logger_interval" in keys(best_per_algo_ics[1].items[1].parsed_args)
        best_per_algo_ics[1].items[1].parsed_args["logger_interval"]
    else
        1000
    end
    num_samples = length(data[1])
    num_runs = length.(best_per_algo_ics)[1]
    plot_params = [LU.get_params(best_per_algo_ics[i]) for i in 1:length(best_per_algo_ics)]

    xticks=collect(1:step_increment:num_samples*step_increment) / tick_scaling
    p = plot(ylabel=ylabel, title=title, xlabel=xlabel, formatter=:plain, figsize=(800,800))
    for i in 1:length(best_per_algo_ics)
        # @show plot_params[i]
        # if best_per_algo_ics[i].items[1].parsed_args["behaviour_reward_projector"] == "tilecoding"
        #     @show "SKETCHY LOOP. ONLY FOR GPI SENSTIVITY"
        #     if best_per_algo_ics[i].items[1].parsed_args["behaviour_rp_tilings"] == 8
        #         plot_params[i][:label] = string(plot_params[i][:label], " Tile Coding")
        #         plot_params[i][:color] =  colorant"#BBBBBB"
        #     elseif best_per_algo_ics[i].items[1].parsed_args["behaviour_rp_tilings"] == 2
        #         plot_params[i][:label] = string(plot_params[i][:label]," Low Tiling")
        #         plot_params[i][:color] =  colorant"#BBBBBB"
        #     end
        # end
        # plot!(p, xticks, data[i], ribbon = std[i]/sqrt(num_runs), legend=:top; plot_params[i]...)
        plot!(p, xticks, data[i], ribbon = std[i]/sqrt(num_runs), legend=:none; plot_params[i]...)
    end
    display(p)
    savefig("./plots/$(folder_name)/RMSE_$(string(data_key)).pdf")
end

function plot_goal_visitation(algo_ics, inds_of_interest = nothing)
    best_per_algo_ics = if inds_of_interest isa Nothing
        deepcopy(algo_ics)
    else
        deepcopy(algo_ics)[inds_of_interest]
    end
    episode_lengths = [GPU.load_results(ic,:episode_length, return_type = "array") for ic in best_per_algo_ics]
    visit_orders = [GPU.load_results(ic,:visit_order, return_type = "array") for ic in best_per_algo_ics]
    max_lengths = [GPU.get_min_length(arrs)-1 for arrs in episode_lengths]
    gvf_labels = ["Distractor" "Constant" "Drifter" "Constant"]

    ps = []
    for i in 1:length(visit_orders)
        @show max_lengths[i]
        @show length(visit_orders[i])
        visit_perc = GPU.goal_visits_per_episode(visit_orders[i], max_lengths[i])
        @show size(visit_perc)
        tmp = [GPU.smooth(visit_perc[i,:],20) for i in 1:4]
        title = LU.get_label(best_per_algo_ics[i])[:label]

        p = plot(tmp, labels = gvf_labels, xlabel="Episode Count", ylabel="Fraction of Goal Visits", ylim=(0.0,1.0), title = title, size = (800,1000), legend=:none)
        push!(ps,p)
    end
    println()
    @show length(visit_orders)
    # plot(ps..., layout=(Int(length(best_per_algo_ics)/2),2))
    plot(ps...)

    savefig("./plots/$(folder_name)/goal_visits_$(string(data_key)).png")
end

function load_and_plot_rmse(inds_of_interest=nothing)
    ic = load_data()
    # ic = search(ic, Dict("demon_learner" => "SR"))
    algo_ic = load_best(ic)
    # plot_rmse_per_demon(algo_ic,inds_of_interest)
    plot_rmse(algo_ic,inds_of_interest)
    # sweep_params = ["demon_alpha_init", "demon_eta", "alpha_init"]
    @show "here"
    [GPU.print_params(bic,["behaviour_learner"],["behaviour_eta","demon_eta","demon_update","emphasis_clip_threshold"]) for bic in algo_ic]
end

function load_and_plot_rmse_per_demon(inds_of_interest=nothing)
    ic = load_data()
    # ic = search(ic, Dict("demon_learner" => "SR"))
    algo_ic = load_best(ic)
    # plot_rmse_per_demon(algo_ic,inds_of_interest)
    plot_rmse_per_demon(algo_ic,inds_of_interest)
    # sweep_params = ["demon_alpha_init", "demon_eta", "alpha_init"]
    [GPU.print_params(bic,["demon_update"],["demon_eta","alpha_init","exploration_param","emphasis_clip_threshold"]) for bic in algo_ic]
end

function load_and_plot_goal_visits(inds_of_interest=nothing)
    ic = load_data()
    algo_ic = load_best(ic)
    plot_goal_visitation(algo_ic,inds_of_interest)
end
function load_and_plot_states(inds_of_interest=nothing)
    ic = load_data()
    algo_ic = load_best(ic)
    plot_states(algo_ic,inds_of_interest)
end
function load_and_plot_episode_lengths(inds_of_interest=nothing)
    ic = load_data()
    algo_ic = load_best(ic)
    # [GPU.print_params(c,["demon_eta", "behaviour_eta"],[]) for c in algo_ic[inds_of_interest]]
    plot_episode_lengths(algo_ic,inds_of_interest)
end

function load_and_box_plot(inds_of_interest=nothing)
    ic = load_data()
    algo_ic = load_best(ic)
    # [GPU.print_params(c,["demon_eta", "behaviour_eta"],[]) for c in algo_ic[inds_of_interest]]
    boxplot(algo_ic;inds_of_interest=inds_of_interest)
end
function load_and_box_plot_back_wall(inds_of_interest=nothing)
    ic = load_data()
    algo_ic = load_best(ic)
    # [GPU.print_params(c,["demon_eta", "behaviour_eta"],[]) for c in algo_ic[inds_of_interest]]
    boxplot_backwall(algo_ic;inds_of_interest=inds_of_interest)
end
function load_and_plot_max_q(inds_of_interest=nothing)
    ic = load_data()
    algo_ic = load_best(ic)
    # [GPU.print_params(c,["demon_eta", "behaviour_eta"],[]) for c in algo_ic[inds_of_interest]]
    plot_max_q(algo_ic,inds_of_interest)
end
function plot_max_q(algo_ic, inds_of_interest)
    best_per_algo_ics = if inds_of_interest isa Nothing
        deepcopy(algo_ic)
    else
        deepcopy(algo_ic)[inds_of_interest]
    end

    max_qs = [GPU.load_results(ic,:max_behaviour_q, return_type = "array") for ic in best_per_algo_ics]
    # #### Generate
    p = plot(xlabel = "Episode Count", ylabel = "Max Q", ylim=(-3,3))
    for (ind,max_q) in enumerate(max_qs)
        l = hcat([q_run for q_run in max_q]...)
        # label = LU.get_label(best_per_algo_ics[i])[:label]
        plot_params = LU.get_params(best_per_algo_ics[ind])
        plot!(p,mean(l,dims=2); plot_params...)
    end
    savefig("./plots/$(folder_name)/max_qs_zoom$(string(data_key)).png")
    # title!("Step Length for Best Performing Algos Last 10 %")
end
function plot_episode_lengths(algo_ic, inds_of_interest)
    best_per_algo_ics = if inds_of_interest isa Nothing
        deepcopy(algo_ic)
    else
        deepcopy(algo_ic)[inds_of_interest]
    end

    episode_lengths = [GPU.load_results(ic,:episode_length, return_type = "array") for ic in best_per_algo_ics]
    for arr in episode_lengths
        for a in arr
            @show length(a)
        end
    end
    # #### Generate
    max_lengths = [GPU.get_min_length(arrs)-1 for arrs in episode_lengths]
    p = plot(xlabel = "Episode Count", ylabel = "Step Length")
    for (ind,episode_lengths) in enumerate(episode_lengths)
        l = hcat([epi[1:max_lengths[ind]] for epi in episode_lengths]...)
        # label = LU.get_label(best_per_algo_ics[i])[:label]
        plot_params = LU.get_params(best_per_algo_ics[ind])
        plot!(p,mean(l,dims=2); plot_params...)
    end
    savefig("./plots/$(folder_name)/episode_lengths_$(string(data_key)).png")
    # title!("Step Length for Best Performing Algos Last 10 %")
end

function load_and_plot_meta_ss(inds_of_interest = nothing; sensitivity_param = "eta")
    other_param = "demon_eta"
    ic = load_data()
    # algo_divisor_keys = ["behaviour_learner", "demon_learner"]
    algo_divisor_keys = ["behaviour_learner", "demon_update"]

    # sweep_params = ["eta"]

    algo_specs_full = GPU.split_algo(ic, algo_divisor_keys)
    all_algos_ics = [search(ic,algo_spec) for algo_spec in algo_specs_full]
    all_algos_ics = if inds_of_interest isa Nothing
        all_algos_ics
    else
        all_algos_ics[inds_of_interest]
    end

    #Plot Sensitivity of eta.
    p = plot(ylabel="RMSE", xlabel = string("Meta Step Size ", sensitivity_param), xaxis=:log, legend=:topleft)
    for algo in all_algos_ics
        eta_vals = diff(algo)[sensitivity_param]
        line = []
        err_bars = []
        for eta in eta_vals
            algo_for_eta = search(algo, Dict(sensitivity_param => eta))
            best_algo_for_eta = GPU.get_best_final_perf(algo_for_eta, other_param, data_key,0.1)

            data = GPU.load_results(best_algo_for_eta,data_key)

            # rmse,rmse_std = GPU.get_stats(data)
            rmse = mean(sum(data, dims = 1),dims = [2,3])[]
            rmse_std = std(mean(data, dims = [1,2]),dims = 3)[]

            num_runs = size(data)[3]

            CI = 1.97 * rmse_std / sqrt(num_runs)

            push!(line, mean(rmse))
            push!(err_bars, CI)
        end
        plot_params = LU.get_params(algo)
        plot!(eta_vals, line, yerror=err_bars,;plot_params...)
    end
    display(p)
    savefig("./plots/$(folder_name)/meta_ss_$(string(data_key)).pdf")
end

function load_and_plot_init_ss(inds_of_interest=nothing; sensitivity_param = "alpha_init")

    ic = load_data()
    algo_divisor_keys = ["behaviour_learner", "demon_learner"]
    sweep_params = ["alpha_init"]

    algo_specs_full = GPU.split_algo(ic, algo_divisor_keys)
    all_algos_ics = [search(ic,algo_spec) for algo_spec in algo_specs_full]
    all_algos_ics = if inds_of_interest isa Nothing
        all_algos_ics
    else
        all_algos_ics[inds_of_interest]
    end    #Plot Sensitivity of eta.
    p = plot(ylabel="RMSE", xlabel = "Initial Step Size", legend=:topright)
    for algo in all_algos_ics
        eta_vals = diff(algo)["alpha_init"]
        line = []
        err_bars = []
        for eta in eta_vals
            algo_for_eta = search(algo, Dict("alpha_init" => eta))
            best_algo_for_eta = GPU.get_best_final_perf(algo_for_eta, ["alpha_init"], data_key,0.1)

            data = GPU.load_results(best_algo_for_eta,data_key)

            # rmse,rmse_std = GPU.get_stats(data)
            rmse = mean(sum(data, dims = 1),dims = [2,3])[]
            rmse_std = std(mean(data, dims = [1,2]),dims = 3)[]

            num_runs = size(data)[3]

            CI = 1.97 * rmse_std / sqrt(num_runs)

            push!(line, mean(rmse))
            push!(err_bars, CI)
        end
        plot_params = LU.get_params(algo)
        plot!(eta_vals, line, yerror=err_bars,;plot_params...)
    end
    display(p)
    savefig("./plots/$(folder_name)/init_ss_$(string(data_key)).pdf")
end

function plot_states(algo_ic, inds_of_interest=nothing)
    best_per_algo_ics = if inds_of_interest isa Nothing
        deepcopy(algo_ic)
    else
        deepcopy(algo_ic)[inds_of_interest]
    end

    ps = []
    for algo in best_per_algo_ics
        states = hcat(GPU.load_results(algo,:states, return_type = "array")...)
        @show size(states)
        #states is 2x length of states
        p = marginalhist(states[1,:],states[2,:])
        push!(ps,p)
    end
    fullp = plot(ps...)
    display(fullp)
    savefig("./plots/$(folder_name)/marginal_hists_$(string(data_key)).png")
end

function boxplot(algo_ic; inds_of_interest = nothing)
    best_per_algo_ics = if inds_of_interest isa Nothing
        deepcopy(algo_ic)
    else
        deepcopy(algo_ic)[inds_of_interest]
    end
    p = plot(ylabel="Right Hill Pseudotermination")
    xticklabels = []
    for (i,algo) in enumerate(best_per_algo_ics)
        # episode_lengths = GPU.load_results(algo,:episode_length, return_type = "array")
        # num_episodes_per_run = [length(e) for e in episode_lengths]
        states_per_run = GPU.load_results(algo,:states, return_type = "array")
        pseudoterm_episodes_per_run = [length(findall(isone,states[1,:])) for states in states_per_run]
        num_episodes_per_run = pseudoterm_episodes_per_run
        plot_params = LU.get_params(algo)
        @show plot_params
        boxplot!(p,num_episodes_per_run; plot_params...)
        push!(xticklabels, plot_params[:label])
    end
    plot!(xticks=([1,2,3],xticklabels),legend=false,formatter=:plain)
    display(p)
    savefig("./plots/$(folder_name)/boxplot_$(string(data_key)).pdf")
end

function boxplot_backwall(algo_ic; inds_of_interest = nothing)
    best_per_algo_ics = if inds_of_interest isa Nothing
        deepcopy(algo_ic)
    else
        deepcopy(algo_ic)[inds_of_interest]
    end
    p = plot(ylabel="Left Hill Pseudotermination")
    xticklabels = []
    for (i,algo) in enumerate(best_per_algo_ics)
        states_per_run = GPU.load_results(algo,:states, return_type = "array")
        pseudoterm_episodes_per_run = [length(findall(iszero,states[1,:])) for states in states_per_run]
        @show pseudoterm_episodes_per_run
        plot_params = LU.get_params(algo)
        boxplot!(p,pseudoterm_episodes_per_run; plot_params...)
        push!(xticklabels, plot_params[:label])
    end
    plot!(xticks=([1,2,3],xticklabels),legend=false,formatter=:plain)
    display(p)
    savefig("./plots/$(folder_name)/boxplot_backwall_$(string(data_key)).pdf")
end
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
function _fileid_reg()
    dig = "\\d"
    num_seps = "[\\.-]"
    id_sep = "[\\b-_]"
    return Regex("^(?<id>(?:$(dig)$(num_seps))*(?:$(dig)$(num_seps))$(id_sep)).*\$")
end

const _FILEID_REG = _fileid_reg()

function _get_file_id(fname::AbstractString)
    fname = basename(fname)
    m = match(_FILEID_REG, fname)
    return isnothing(m) ? "" : m[:id]
end

macro fileid()
    quote
        local file = $(string(__source__.file))
        $(_get_file_id)(file)
    end
end


text:
"""
ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля, на котором могут находиться также внутренние прямоугольные перегородки (все перегородки изолированы друг от друга, прямоугольники могут вырождаться в отрезки)

РЕЗУЛЬТАТ: Робот - в исходном положении, и в 4-х приграничных клетках, две из которых имеют ту же широту, а две - ту же долготу, что и Робот, стоят маркеры.
"""
function coordinat_marker(r::Robot)
    steps=0
    arr_of_steps=[]
    while ((isborder(r,Sud)) && (isborder(r,West)))==false
        push!(arr_of_steps,moves!(r,West))
        push!(arr_of_steps,moves!(r,Sud))
        steps+=2
    end
    vertical=0 #Счётчик - Движение по вертикали
    horizontal=0 #Счётчик - Движение по горизонтали
    for a in 1:steps
        if isodd(a)==true #Проверка на нечётность
            horizontal+=arr_of_steps[a] #Суммирование, значения берутся из массива
        else
            vertical+=arr_of_steps[a]
        end
    end
    side=Nord
    for b in 1:2 #Движение по вертикали и горизонтали
        vertical=coordinater(r,side,vertical)
        side=next(side)
        horizontal=coordinater(r,side,horizontal)
        side=next(side)
    end
    while (steps>0)==true #Аналогичная функция описана в Задаче 5
        side=isodd(steps) ? Ost : Nord
        move_by_points(r,side,arr_of_steps[steps])
        steps-=1
    end
end

function coordinater(r::Robot, side::HorizonSide, num_steps::Int)#Движение и маркировка
    move_by_points(r,side,num_steps)
    putmarker!(r)
    num_steps=moves!(r,side)
    return num_steps
end

function moves!(r::Robot,side::HorizonSide) #Возвращает количество шагов
    num_steps=0
    while isborder(r,side)==false
        move!(r,side)
        num_steps+=1
    end
    return num_steps
end

function move_by_points(r::Robot, side::HorizonSide, num_steps::Int) #Движение в заданную сторону
    for b in 1:num_steps #Количество шагов определено
        move!(r,side)
    end
end

function next(side::HorizonSide) #Возможно использование библиотечного файла
    HorizonSide(mod(Int(side)+3,4))
end
#=Возможна реализация с помощью других функций, в частности, <sum> - суммирование элементов массива=#

max_stars_repo_name: alexgruzberg/mirea-progs

max_issues_repo_name: alexgruzberg/mirea-progs

max_forks_repo_name: alexgruzberg/mirea-progs


text:
module Testfast

export @before_all,
       @setup,
       @teardown,
       @after_all,
       @test,
       run!

type Test
  name::String
  code::Expr
end

type Failure
  exception::ErrorException
  test::Test
end

global before_all = Expr[]
global setup = Expr[]
global teardown = Expr[]
global after_all = Expr[]
global files = String[]
global tests = Test[]
global failures = Failure[]
global total = 0

macro before_all(e::Expr)
  push!(before_all, e)
end

macro setup(e::Expr)
  push!(setup, e)
end

macro teardown(e::Expr)
  push!(teardown, e)
end

macro after_all(e::Expr)
  push!(after_all, e)
end

macro test(name::String, code::Expr)
  push!(tests, Test(name, code))
end

function add_test_files(directory = ".")
  for file in readdir(directory)
    file = "$directory/$file"
    if file == "." || file == ".."
      continue
    elseif oct(stat(file).mode)[1] == '4'
      add_test_files(file)
    elseif split(file, ".")[end-1][end-3:end] == "test"
      push!(files, file)
    end
  end
end

function run_tests()
  global before_all
  global setup
  global teardown
  global after_all
  global tests
  global failures
  global total
  global started

  for file in files
    _before_all = before_all
    _setup = setup
    _teardown = teardown
    _after_all = after_all

    tests = Test[]
    started = time()

    reload(file)

    for expr in before_all
      eval(:(esc($expr)))
    end

    for test in tests
      total += 1

      for expr in setup
        eval(:(esc($expr)))
      end

      try
        eval(:(esc(eval($test.code))))
        print_with_color(:green, ".")
      catch e
        print_with_color(:red, "F")
        push!(failures, Failure(e, test))
      end

      for expr in teardown
        eval(:(esc($expr)))
      end
    end

    for expr in after_all
      eval(:(esc($expr)))
    end

    before_all = _before_all
    setup = _setup
    teardown = _teardown
    after_all = _after_all
  end

end

function print_results()
  println("\n")

  println("$(total) tests run.\n$(total - length(failures)) passed, $(length(failures)) failed in $(round(time() - started, 4)) seconds.")

  if isempty(failures)
    exit(0)
  else
    println("Failures:")
    for failure in failures
      print("$(failure.test.name)\n\t$(failure.exception):\n\t")
      @show failure.test
      println("\n")
    end
    exit(1)
  end

end

function run!()
  # In case we're not in the test directory already...
  try
    cd("test")
  catch
    nothing # We are!
  end

  add_test_files()
  run_tests()
  print_results()
end

end


text:
using Distributions
using LinearAlgebra
using Plots

# Set Parameters
M = 9  # total number of training data
N = 9  # degree of a polynomial
lambda_ = 1.75e-6  # regularization parameter

# Set Function to Generate Data
dist = Normal(0.0, 0.1)
base(x) = sin.(2pi*x)
gen(x) = base(x) + rand(dist, length(x))

# Set Function to Generate Polynomial Matrix
function polynomial_matrix(x, N)
    m = length(x)
    X = fill(1.0, m)
    for j in 1:N
        x_vec = x .^ j
        X = hcat(X, x_vec)
    end
    return X
end

# Set Function to Calculate Coefficient
function weight(x, y, N, lambda_)
    X = polynomial_matrix(x, N)
    w = (X' * X + lambda_ * Matrix{Float64}(I, N+1, N+1)) \ X' * y
    return w
end

# Set Function to Output Value
function f(x, w)
    X = polynomial_matrix(x, length(w)-1)
    y = X * w
    return y
end

# Generate Data
step = 1.0 / (M - 1)
x = collect(0.0:step:1.0)
y = gen(x)
println("x = ", x)
println("y = ", y)

# Calculate Coefficient of Polynomial
w = weight(x, y, N, lambda_)
println("w = ", w)

# Output Graph
gr()
scatter(x, y, xlabel="x", ylabel="y", label="data")
graph_x = collect(0.0:0.01:1.0)
graph_y = f(graph_x, w)
plot!(graph_x, base(graph_x), label="base")
plot!(graph_x, graph_y, label=string("f(x, w)"))

# Save Image
savefig("graph")


text:
@test degree_centrality(g5) == [0.6666666666666666, 0.6666666666666666, 1.0, 0.3333333333333333]
@test indegree_centrality(g5, normalize=false) == [0.0, 1.0, 2.0, 1.0]
@test outdegree_centrality(g5; normalize=false) == [2.0, 1.0, 1.0, 0.0]


text:
using MyExample
using Test

@testset "MyExample.jl" begin
    @test my_f(1, 1) == 5
    @test my_f(2, 1) == 7
    @test my_f(2, 2) == 10
    # Write your tests here.
end


text:
#!/usr/bin/env julia

# Draw a sierpinski tiangle.

using Compose

function sierpinski(n)
    if n == 0
        canvas() << polygon((1,1), (0,1), (1/2, 0))
    else
        t = sierpinski(n - 1)
        compose(canvas(),
                (canvas(1/4,   0, 1/2, 1/2), t),
                (canvas(  0, 1/2, 1/2, 1/2), t),
                (canvas(1/2, 1/2, 1/2, 1/2), t))
    end
end

#img = SVG("sierpinski.svg", 4inch, 4(sqrt(3)/2)inch)
img = D3("sierpinski.js", 4inch, 4(sqrt(3)/2)inch)
draw(img, sierpinski(8) << linewidth(0.1mm) << fill(nothing))


text:
const DEFINED_DIFFRULES = Dict{Tuple{Union{Expr,Symbol},Symbol,Int},Any}()

"""
    @define_diffrule M.f(x) = :(df_dx(\$x))
    @define_diffrule M.f(x, y) = :(df_dx(\$x, \$y)), :(df_dy(\$x, \$y))
    ⋮

Define a new differentiation rule for the function `M.f` and the given arguments, which should
be treated as bindings to Julia expressions. Return the defined rule's key.

The LHS should be a function call with a non-splatted argument list, and the RHS should be
the derivative expression, or in the `n`-ary case, an `n`-tuple of expressions where the
`i`th expression is the derivative of `f` w.r.t the `i`th argument. Arguments should be
interpolated wherever they are used on the RHS.

Note that differentiation rules are purely symbolic, so no type annotations should be used.

Examples:

    @define_diffrule Base.cos(x)          = :(-sin(\$x))
    @define_diffrule Base.:/(x, y)        = :(inv(\$y)), :(-\$x / (\$y^2))
    @define_diffrule Base.polygamma(m, x) = :NaN,       :(polygamma(\$m + 1, \$x))

"""
macro define_diffrule(def)
    @assert isa(def, Expr) && def.head == :(=) "Diff rule expression does not have a left and right side"
    lhs = def.args[1]
    rhs = def.args[2]
    @assert isa(lhs, Expr) && lhs.head == :call "LHS is not a function call"
    qualified_f = lhs.args[1]
    @assert isa(qualified_f, Expr) && qualified_f.head == :(.) "Function is not qualified by module"
    M = qualified_f.args[1]
    f = _get_quoted_symbol(qualified_f.args[2])
    args = lhs.args[2:end]
    rule = Expr(:->, Expr(:tuple, args...), rhs)
    key = Expr(:tuple, Expr(:quote, M), Expr(:quote, f), length(args))
    return esc(quote
        $DiffRules.DEFINED_DIFFRULES[$key] = $rule
        $key
    end)
end

"""
    diffrule(M::Union{Expr,Symbol}, f::Symbol, args...)

Return the derivative expression for `M.f` at the given argument(s), with the argument(s)
interpolated into the returned expression.

In the `n`-ary case, an `n`-tuple of expressions will be returned where the `i`th expression
is the derivative of `f` w.r.t the `i`th argument.

Examples:

    julia> DiffResults.diffrule(:Base, :sin, 1)
    :(cos(1))

    julia> DiffResults.diffrule(:Base, :sin, :x)
    :(cos(x))

    julia> DiffResults.diffrule(:Base, :sin, :(x * y^2))
    :(cos(x * y ^ 2))

    julia> DiffResults.diffrule(:Base, :^, :(x + 2), :c)
    (:(c * (x + 2) ^ (c - 1)), :((x + 2) ^ c * log(x + 2)))
"""
diffrule(M::Union{Expr,Symbol}, f::Symbol, args...) = DEFINED_DIFFRULES[M,f,length(args)](args...)

"""
    hasdiffrule(M::Union{Expr,Symbol}, f::Symbol, arity::Int)

Return `true` if a differentiation rule is defined for `M.f` and `arity`, or return `false`
otherwise.

Here, `arity` refers to the number of arguments accepted by `f`.

Examples:

    julia> DiffResults.hasdiffrule(:Base, :sin, 1)
    true

    julia> DiffResults.hasdiffrule(:Base, :sin, 2)
    false

    julia> DiffResults.hasdiffrule(:Base, :-, 1)
    true

    julia> DiffResults.hasdiffrule(:Base, :-, 2)
    true

    julia> DiffResults.hasdiffrule(:Base, :-, 3)
    false
"""
hasdiffrule(M::Union{Expr,Symbol}, f::Symbol, arity::Int) = haskey(DEFINED_DIFFRULES, (M, f, arity))

"""
    diffrules()

Return a list of keys that can be used to access all defined differentiation rules.

Each key is of the form `(M::Symbol, f::Symbol, arity::Int)`.

Here, `arity` refers to the number of arguments accepted by `f`.

Examples:

    julia> first(DiffRules.diffrules())
    (:Base, :asind, 1)

"""
diffrules() = keys(DEFINED_DIFFRULES)

# For v0.6 and v0.7 compatibility, need to support having the diff rule function enter as a
# `Expr(:quote...)` and a `QuoteNode`. When v0.6 support is dropped, the function will
# always enter in a `QuoteNode` (#23885).
function _get_quoted_symbol(ex::Expr)
    @assert ex.head == :quote
    @assert length(ex.args) == 1 && isa(ex.args[1], Symbol) "Function not a single symbol"
    ex.args[1]
end

function _get_quoted_symbol(ex::QuoteNode)
    @assert isa(ex.value, Symbol) "Function not a single symbol"
    ex.value
end


text:
using DetectAlpha
using Test
using Revise
include("../src/Utils.jl")

@testset "DetectAlpha.jl" begin
    # Write your tests here.
    @test π ≈ 3.14 atol=0.01
    elemens = Elements()
    plutrow = DetectAlpha.get_iso_row(94,elemens)
    @test plutrow !== nothing

    spectra = DetectAlpha.read_alpha_spectrum_file("../Alpha2.csv")
    @test spectra !== nothing

    alphaspectrum = first(spectra)
    @test alphaspectrum.numchannels == 512

    for spectrum in spectra
        find_peaks(spectrum)
    end

end

@testset "Utils.jl" begin
    @test parsebool("false") == false
    @test parsebool("yes") == true
    @test between2and8(5) == true
    @test between2and8(9) == false
    @test between2and8(5//2) == true
    @test between2and8(4//2) == false
end

@testset "AlphaSpectrum" begin
    alphaspectra = DetectAlpha.example_alpha_spectrum()
    @test length(alphaspectra) >= 1
    spectrum = first(alphaspectra)
    hs = DetectAlpha.to_histogram(spectrum)
    @test hs !== nothing

    @test_throws BoundsError find_peak(StepRange(1,1,5128),spectrum)

    @test !valid_peak(find_peak(StepRange(Int32(1),Int32(1),Int32(128)),spectrum))
end


text:
# ---
# title: 975. Odd Even Jump
# id: problem975
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming, Stack, Ordered Map
# link: <https://leetcode.com/problems/odd-even-jump/description/>
# hidden: true
# ---
#
# You are given an integer array `A`. From some starting index, you can make a
# series of jumps. The (1st, 3rd, 5th, ...) jumps in the series are called
# **odd-numbered jumps** , and the (2nd, 4th, 6th, ...) jumps in the series are
# called **even-numbered jumps**. Note that the **jumps** are numbered, not the
# indices.
#
# You may jump forward from index `i` to index `j` (with `i < j`) in the
# following way:
#
#   * During **odd-numbered jumps** (i.e., jumps 1, 3, 5, ...), you jump to the index `j` such that `A[i] <= A[j]` and `A[j]` is the smallest possible value. If there are multiple such indices `j`, you can only jump to the **smallest** such index `j`.
#   * During **even-numbered jumps** (i.e., jumps 2, 4, 6, ...), you jump to the index `j` such that `A[i] >= A[j]` and `A[j]` is the largest possible value. If there are multiple such indices `j`, you can only jump to the **smallest** such index `j`.
#   * It may be the case that for some index `i`, there are no legal jumps.
#
# A starting index is **good** if, starting from that index, you can reach the
# end of the array (index `A.length - 1`) by jumping some number of times
# (possibly 0 or more than once).
#
# Return _the number of **good** starting indices_.
#
#
#
# **Example 1:**
#
#
#
#     Input: A = [10,13,12,14,15]
#     Output: 2
#     Explanation:
#     From starting index i = 0, we can make our 1st jump to i = 2 (since A[2] is the smallest among A[1], A[2], A[3],
#     A[4] that is greater or equal to A[0]), then we cannot jump any more.
#     From starting index i = 1 and i = 2, we can make our 1st jump to i = 3, then we cannot jump any more.
#     From starting index i = 3, we can make our 1st jump to i = 4, so we have reached the end.
#     From starting index i = 4, we have reached the end already.
#     In total, there are 2 different starting indices i = 3 and i = 4, where we can reach the end with some number of
#     jumps.
#
#
# **Example 2:**
#
#
#
#     Input: A = [2,3,1,1,4]
#     Output: 3
#     Explanation:
#     From starting index i = 0, we make jumps to i = 1, i = 2, i = 3:
#
#     During our 1st jump (odd-numbered), we first jump to i = 1 because A[1] is the smallest value in [A[1], A[2],
#     A[3], A[4]] that is greater than or equal to A[0].
#
#     During our 2nd jump (even-numbered), we jump from i = 1 to i = 2 because A[2] is the largest value in [A[2], A[3],
#     A[4]] that is less than or equal to A[1]. A[3] is also the largest value, but 2 is a smaller index, so we can
#     only jump to i = 2 and not i = 3
#
#     During our 3rd jump (odd-numbered), we jump from i = 2 to i = 3 because A[3] is the smallest value in [A[3], A[4]]
#     that is greater than or equal to A[2].
#
#     We can't jump from i = 3 to i = 4, so the starting index i = 0 is not good.
#
#     In a similar manner, we can deduce that:
#     From starting index i = 1, we jump to i = 4, so we reach the end.
#     From starting index i = 2, we jump to i = 3, and then we can't jump anymore.
#     From starting index i = 3, we jump to i = 4, so we reach the end.
#     From starting index i = 4, we are already at the end.
#     In total, there are 3 different starting indices i = 1, i = 3, and i = 4, where we can reach the end with some
#     number of jumps.
#
#
# **Example 3:**
#
#
#
#     Input: A = [5,1,3,4,2]
#     Output: 3
#     Explanation:
#     We can reach the end from starting indices 1, 2, and 4.
#
#
#
#
# **Constraints:**
#
#   * `1 <= A.length <= 2 * 104`
#   * `0 <= A[i] < 105`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module Orbitals

import Base.getindex
import Base.lastindex
import Base.setindex!
import Base.length

abstract type AbstractOrbital end

struct CanonicalOrbital <: AbstractOrbital
    C::Array{Float64,1}
end # module
struct NaturalOrbital <: AbstractOrbital
    C::Array{Float64,1}
    ω::Float64 # occupation number
end
struct LocalOrbital <: AbstractOrbital
    C::Array{Float64,1}
end
struct BruecknerOrbital <: AbstractOrbital
    C::Array{Float64,1}
end

abstract type AbstractOrbitals end

struct CanonicalOrbitals <: AbstractOrbitals
    orbs::Array{CanonicalOrbital}
end
struct NaturalOrbitals <: AbstractOrbitals
    orbs::Array{NaturalOrbitals}
end
struct LocalOrbitals <: AbstractOrbitals
    orbs::Array{LocalOrbitals}
end
struct BruecknerOrbitals <: AbstractOrbitals
    orbs::Array{BruecknerOrbital}
end

function CanonicalOrbitals(C::Array{Float64,2})
    return CanonicalOrbitals([CanonicalOrbital(Array{Float64,1}(C[:,i])) for i in 1:size(C,2)])
end

function toarray(orbs::O) where O <: AbstractOrbitals
    C = hcat([orb.C for orb in orbs.orbs]...)
end

function Base.getindex(orb::O,i) where O <: AbstractOrbitals
    return orb.orbs[i]
end
function Base.getindex(orb::O,R::UnitRange) where O <: AbstractOrbitals
    return O([orb.orbs[i] for i in R])
end

function Base.lastindex(orb::O) where O <: AbstractOrbitals
    return length(orb.orbs)
end

function Base.length(orb::O) where O <: AbstractOrbitals
    return length(orb.orbs)
end

mutable struct OrbDict
    current::String
    ndocc::Int
    nvir::Int
    frozencore::Int
    frozenvir::Int
    cache::Dict{String,O} where O <: AbstractOrbitals
end

function matchphase(orbs,oldC)
    for oi in eachindex(orbs.cache[orbs.current].orbs)
        o = orbs.cache[orbs.current].orbs[oi]
        osorted = reverse(sort(o.C,by=abs))
        csorted = reverse(sort(oldC[:,oi],by=abs))
        if sign(osorted[1]) != sign(csorted[1])
            o.C .*= -1
        end
    end
end
"""
call like readC(ints.orbs,Cmatrix)
"""
function readC(orbs,newC)
    for oi in eachindex(orbs.cache[orbs.current].orbs)
        o = orbs.cache[orbs.current].orbs[oi]
        o.C .= newC[:,oi]
    end
end
function topositive!(orbs::OrbDict)
    for o in orbs.cache[orbs.current].orbs
        sorted = reverse(sort(o.C,by=abs))
        if sorted[1] >= 0
            o.C .*= 1
        else
            o.C .*= -1
        end
    end
end
function rotate!(orbs::OrbDict,U::Array{Float64,2};fc=0,fv=0) where O <: AbstractOrbitals
    C = hcat([orb.C for orb in orbs["OV"].orbs]...)
    C = C*U
    for i=1:length(orbs["OV"])
        _i = i+fc
        orbs.cache[orbs.current].orbs[_i].C[:] .= C[:,i]
    end
end

function OrbDict()
    OrbDict("canonical",0,0,0,0,Dict{String,AbstractOrbitals}())
end

function activate!(OD::OrbDict,entry::String)
    if !(entry in keys(OD.cache))
        error("OrbDict tried to activate a non-existent basis!")
    else
        OD.current = entry
    end
end

function Base.getindex(OD::OrbDict,entry::String)
    gather = occursin("[",entry)
    s = strip(entry,['[',']'])
    all = OD.cache[OD.current]
    type = typeof(all)
    fc = OD.frozencore
    fv = OD.frozenvir
    lookup = Dict{Char,type}(
                             'F' => all[1:OD.ndocc],
                             'O' => all[fc+1:OD.ndocc],
                             'U' => all[OD.ndocc+1:OD.ndocc+OD.nvir],
                             'V' => all[OD.ndocc+1:OD.ndocc+OD.nvir-fv]
                            )
    out = []
    for i in eachindex(s)
        _s = lookup[s[i]]
        if length(_s) != 0
            push!(out,lookup[s[i]])
        end
    end
    out = type(hcat([toarray(o) for o in out]...))
    if gather
        return toarray(out)
    else
        return out
    end
end

function Base.setindex!(OD::OrbDict,val,entry::String)
    OD.cache[entry] = val
end

end #module


text:
push!(LOAD_PATH, "/Users/sean/GitRepos/ParticleMapper/ParticleGridMapper.jl")
push!(LOAD_PATH, "/Users/sean/GitRepos/Chia-Yu_Hu_ParticleMapper/OctreeBH.jl")

using ParticleGridMapper
using StaticArrays
using PyPlot
using HDF5

const N = 3 #spatial dimension
const T = Float64

const BOXSIZE = 1.0  #kpc
const max_depth = 10 #maximum refinement level (RADMC requires < 20)

const nx = ny = 512  #for the column density maps of H2 and CO (npix in RADMC)
const ix, iy = 1, 2  #axes of the column density maps

function read_snap(filename)

    header::Dict = h5readattr(filename, "/Header")
    boxsize::T = header["BoxSize"]
    time::T    = header["Time"]

    N_gas::Int64 = header["NumPart_ThisFile"][1]

    pos_gas::Matrix{T} = h5read(filename, "PartType0/Coordinates");
    vel_gas::Matrix{T} = h5read(filename, "PartType0/Velocities");
    rho::Vector{T}     = h5read(filename, "PartType0/Density");
    u::Vector{T}       = h5read(filename, "PartType0/InternalEnergy");
    m_gas::Vector{T}   = h5read(filename, "PartType0/Masses");
    hsml::Vector{T}    = h5read(filename, "PartType0/SmoothingLength");

    id_gas::Vector{Int64} = h5read(filename, "PartType0/ParticleIDs");

    return N_gas, pos_gas, vel_gas, rho, u, m_gas, hsml, id_gas, boxsize, time
end

function read_chemistry(fname)
    header = h5readattr(fname, "/Header")
    all_species = header["all_species"]
    abund_all   = h5read(fname, "Chemistry/Abundances");
    N_spec = size(abund_all)[1]
    d = Dict(all_species .=> collect(1:N_spec) );
    d[""] = d["PHOTON"] = d["CRP"] = d["CRPHOT"] = 0
    xH2 = abund_all[d["H2"], :]
    xHI = abund_all[d["H"], :]
    xHp = abund_all[d["H+"], :]
    xCO = abund_all[d["CO"], :]
    xCI = abund_all[d["C"], :]
    xCp = abund_all[d["C+"], :]
    xelec = abund_all[d["e-"], :]
    return xH2, xHI, xHp, xCO, xCI, xCp, xelec
end

function vec2svec(vec::Matrix{T}) where {T}
    svec = [SVector{size(vec,1),T}(vec[:,i]) for i in 1:size(vec,2)]
end
function mat2smat(mat::Array{T,3}) where {T}
    smat = [SMatrix{3,3,T}(mat[:,:,i]) for i in 1:size(mat,3)]
end

const XH = 0.71
const BOLTZMANN=1.3806e-16
const PROTONMASS=1.6726e-24
const GRAVCON=6.67e-8
const UnitMass_in_g = 1.989e43
const UnitLength_in_cm    = 3.085678e21
const UnitTime_in_s = 3.08568e+16
const UnitDensity_in_cgs = UnitMass_in_g / UnitLength_in_cm^3
const UnitDensity_in_pccm = UnitDensity_in_cgs/PROTONMASS
const Year_in_s = 31556926.
const fac_col = (UnitMass_in_g/UnitLength_in_cm^2)*(XH/PROTONMASS)

let

i=500
@show i
snap = ""
if i < 10
    snap = "00" * string(i)
elseif i < 100
    snap = "0" * string(i)
else
    snap = string(i)
end

file_path = "/Users/chu/simulations/tallbox/SFSNPI_N1e5_gS10H250dS40H250_soft0p92_from200Myr_SFcut_Z1"
#file_path = "/Users/chu/simulations/tallbox/SFSNPI_N1e7_gS10H250dS40H250_soft0p2_from100Myr_Z1"
#file_path = "/Users/chu/simulations/turbbox/N32"
#file_path = "./"
filename = file_path * "/snap_" * snap * ".hdf5"

Npart, pos, vel, rho, u, mass, hsml, id_gas, boxsize, time = read_snap(filename);

X = vec2svec(pos);

vx = vel[1,:]
vy = vel[2,:]
vz = vel[3,:]

topnode_length = @SVector(ones(N)) * BOXSIZE  #actual length of tree

#center = topnode_length .* 0.5
center = @SVector [0.5,0.5,0.]

part = [DataP2G{N,T}(SVector(X[i]), i, hsml[i], mass[i], 0.0, Int64[]) for i in eachindex(X)]

#build the tree
println("building the tree...")
@time treeNgb = buildtree(part, center, topnode_length);
grid_volumes = get_AMRgrid_volumes(treeNgb)
@assert sum(grid_volumes) ≈ prod(topnode_length) #partition unity check
println( "true maximum depth of tree = ", log2( BOXSIZE / minimum(grid_volumes)^(1/3) ) )

max_depth_true = get_max_tree_depth(treeNgb)
@show max_depth, max_depth_true

treeAMR = deepcopy(treeNgb)
set_max_depth_AMR!(treeAMR, max_depth)
balance_all_level!(treeAMR)

println("getting AMR grid structure...")
@time gridAMR = get_AMRgrid(treeAMR)
num_nodes  = length(gridAMR[gridAMR.==1])
num_leaves = length(gridAMR[gridAMR.==0])
@show num_nodes, num_leaves

boxsizes = @SVector(ones(N)) * BOXSIZE  #for periodic B.C.

PARSEC   = 3.08572e18 #[cm] important to use the exact same value as in radmc to avoid unnecessary sub-pixeling
KILOPARSEC = 1e3 * 3.08572e18

#xmin, ymin, zmin = (center .- 0.5.*boxsizes) .* KILOPARSEC
#xmax, ymax, zmax = (center .+ 0.5.*boxsizes) .* KILOPARSEC
#!!!!!! RADMC default camera is at (0,0), so we let (xmin,xmax) = (-0.5,0.5) instead of (0,1)
#otherwise we will miss 3/4 of the image!!!
xmin, ymin, zmin = -0.5.*boxsizes .* KILOPARSEC
xmax, ymax, zmax = +0.5.*boxsizes .* KILOPARSEC

#the ascii version of amr_grid.inp isn't much bigger from the binary one so we don't bother to go for binary
println("writing to amr_grid.inp...")
@time open("amr_grid.inp", "w") do f
    write(f, "1\n")                       # iformat
    write(f, "1\n")                       # AMR grid style  (0=regular grid, no AMR)
    write(f, "0\n")                       # Coordinate system (Cartesian)
    write(f, "0\n")                       # gridinfo =0 to save disc space
    write(f, "1 1 1\n")                   # included dimension (1 1 1 for 3D)
    write(f, "1 1 1\n")                   # nx ny nz
    write(f, "$(max_depth) $(num_leaves) $(num_leaves+num_nodes)\n")  #depth, # of leaves, # of nodes
    write(f, "$(xmin) $(xmax)\n")                       #base grid coordinate in x
    write(f, "$(ymin) $(ymax)\n")                       #base grid coordinate in y
    write(f, "$(zmin) $(zmax)\n")                       #base grid coordinate in z
    for i in eachindex(gridAMR)
        write(f, "$(gridAMR[i])\n")                       #0 = leaf, 1 = nonleaf
    end
end

println("mapping vx...")
@time map_particle_to_AMRgrid_MFM!(treeAMR, vx, X, hsml, boxsizes, treeNgb=treeNgb)
vxAMR_MFM = get_AMRfield(treeAMR)
println("mapping vy...")
@time map_particle_to_AMRgrid_MFM!(treeAMR, vy, X, hsml, boxsizes)
vyAMR_MFM = get_AMRfield(treeAMR)
println("mapping vz...")
@time map_particle_to_AMRgrid_MFM!(treeAMR, vz, X, hsml, boxsizes)
vzAMR_MFM = get_AMRfield(treeAMR)

open("gas_velocity.binp", "w") do f
    write(f, 1) # iformat
    write(f, 8) # double precision (only for binary input)
    write(f, num_leaves)
    for i in eachindex(vxAMR_MFM)
        write(f, vxAMR_MFM[i]*1e5, vyAMR_MFM[i]*1e5, vzAMR_MFM[i]*1e5) #in cm/sec (instead of km/sec)!!!
        #write(f, vxAMR_MFM[i]*0, vyAMR_MFM[i]*0, vzAMR_MFM[i]*0)
    end
end

println("reading chemistry data...")
fname_base = file_path * "/chem-neqH2Hp-noCOic-TF3-GrRec-" * snap * "-3"
fname = fname_base * ".hdf5"
xH2, xHI, xHp, xCO, xCI, xCp, xelec = read_chemistry(fname);

nHtot = rho.*(UnitDensity_in_pccm*XH)

println("mapping nH2...")
@time map_particle_to_AMRgrid_MFM!(treeAMR, xH2.*nHtot, X, hsml, boxsizes)
nH2_AMR_MFM = get_AMRfield(treeAMR)

NH2 = project_AMRgrid_to_image_thread(nx, ny, ix, iy, treeAMR, center, boxsizes);

println("mapping nCO...")
@time map_particle_to_AMRgrid_MFM!(treeAMR, xCO.*nHtot, X, hsml, boxsizes)
nCO_AMR_MFM = get_AMRfield(treeAMR)

NCO = project_AMRgrid_to_image_thread(nx, ny, ix, iy, treeAMR, center, boxsizes);

XHe = 0.1
mu = @. 1 / (XH * (xHI + xH2 + xHp + XHe + xelec));
temp = u.*mu .* (1e10*PROTONMASS/(1.5*BOLTZMANN));

println("mapping temp...")
@time map_particle_to_AMRgrid_MFM!(treeAMR, temp, X, hsml, boxsizes)
tempAMR_MFM = get_AMRfield(treeAMR)

#ortho-to-para ratio = 3
fac_para = 0.25
fac_ortho = 1 - fac_para
open("numberdens_p-h2.binp", "w") do f
    write(f, 1) # iformat
    write(f, 8) # double precision (only for binary input)
    write(f, num_leaves)
    for i in eachindex(nH2_AMR_MFM)
        write(f, fac_para.*nH2_AMR_MFM[i])
    end
end
open("numberdens_o-h2.binp", "w") do f
    write(f, 1) # iformat
    write(f, 8) # double precision (only for binary input)
    write(f, num_leaves)
    for i in eachindex(nH2_AMR_MFM)
        write(f, fac_ortho.*nH2_AMR_MFM[i])
    end
end

open("numberdens_co.binp", "w") do f
    write(f, 1) # iformat
    write(f, 8) # double precision (only for binary input)
    write(f, num_leaves)
    for i in eachindex(nCO_AMR_MFM)
        write(f, nCO_AMR_MFM[i])
    end
end

open("gas_temperature.binp", "w") do f
    write(f, 1) # iformat
    write(f, 8) # double precision (only for binary input)
    write(f, num_leaves)
    for i in eachindex(tempAMR_MFM)
        write(f, tempAMR_MFM[i] < 1.99e3 ? tempAMR_MFM[i] : 1.99e3) #T_max = 2000 K in RADMC3D
    end
end

L_esc0   = 1e30
L_esc = ones(num_leaves).*L_esc0
open("escprob_lengthscale.binp", "w") do f
    write(f, 1) # iformat
    write(f, 8) # double precision (only for binary input)
    write(f, num_leaves)
    for i in eachindex(L_esc)
        write(f, L_esc[i])
    end
end

vturb0 = 3.0 #km/s
vturb = ones(num_leaves).*vturb0
open("microturbulence.binp", "w") do f
    write(f, 1) # iformat
    write(f, 8) # double precision (only for binary input)
    write(f, num_leaves)
    for i in eachindex(vturb)
        write(f, vturb[i]*1e5) #in cm/sec (instead of km/sec)!!!
    end
end

open("wavelength_micron.inp", "w") do f
    write(f, "2\n")                       # iformat
    write(f, "1e-1\n")
    write(f, "1e4\n")
end

# Write the lines.inp control file
open("lines.inp", "w") do f
    write(f, "2\n")
    write(f, "1\n")
    write(f, "co    leiden    0    0    2\n")
    write(f, "p-h2\n")
    write(f, "o-h2\n")
end

open("radmc3d.inp", "w") do f
    write(f, "lines_mode = 3\n")
    write(f, "lines_slowlvg_as_alternative = 1\n") #more stable LVG (but slower)
    #write(f, "lines_tbg = 0\n")
    write(f, "rto_style = 3\n") #binary output
    write(f, "writeimage_unformatted = 1\n") #binary image output
end

#save NH2 & NCO
fname = file_path * "/data_AMRmaps_H2CO_N" * string(nx) * "_" * snap * ".hdf5"
fid=h5open(fname,"w")
grp_part = create_group(fid,"AMRmaps");
h5write(fname, "AMRmaps/NH2", NH2)
h5write(fname, "AMRmaps/NCO", NCO)
close(fid)

end #let