text:
# QuEST_jl/src/base/data_structure_functions.jl
#

import .._QUEST_LIB

function createQuESTEnv() :: QuEST_Types.QuESTEnv

    if _QUEST_LIB[]==Ptr{Cvoid}(0)
        QuEST_init()
    else
        @assert 4*ccall(:getQuEST_PREC,Cint,(),) == sizeof(Qreal)
    end

    return ccall(:createQuESTEnv, QuEST_Types.QuESTEnv, () )

end

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

function createCloneQureg(qureg ::QuEST_Types.Qureg, env ::QuEST_Types.QuESTEnv) ::QuEST_Types.Qureg
    return ccall(:createCloneQureg, QuEST_Types.Qureg, (QuEST_Types.Qureg,QuEST_Types.QuESTEnv), qureg, env)
end

function createComplexMatrixN(numQubits ::Integer) ::QuEST_Types.ComplexMatrixN
    @assert 1 ≤ numQubits ≤ 50

    return ccall(:createComplexMatrixN, QuEST_Types.ComplexMatrixN,
                 (Cint,),
                 numQubits)
end

function createDensityQureg(numQubits ::Integer,
                            env       ::QuEST_Types.QuESTEnv) ::QuEST_Types.Qureg
    @assert 1 ≤ numQubits ≤ 50
    return ccall(:createDensityQureg, QuEST_Types.Qureg,
                 (Cint,     QuEST_Types.QuESTEnv),
                 numQubits, env)
end

function createDiagonalOp(numQubits     :: Integer,
                          env           :: QuEST_Types.QuESTEnv)    :: QuEST_Types.DiagonalOp
    return ccall(:createDiagonalOp, QuEST_Types.DiagonalOp,
                 (Cint,     QuEST_Types.QuESTEnv),
                 numQubits, env)
end

function createPauliHamil(numQubits     :: Integer,
                          numSumTerms   :: Integer)     :: QuEST_Types.PauliHamil
    return ccall(:createPauliHamil, QuEST_Types.PauliHamil,
                 (Cint, Cint),
                 numQubits, numSumTerms)
end

function createPauliHamilFromFile(filename ::String)   :: QuEST_Types.PauliHamil
    return ccall(:createPauliHamilFromFile, QuEST_Types.PauliHamil,
                 (Cstring,),
                 filename)
end

function createQureg(numQubits ::Integer,
                     env       ::QuEST_Types.QuESTEnv) ::QuEST_Types.Qureg
    @assert 1 ≤ numQubits ≤ 50

    return ccall(:createQureg, QuEST_Types.Qureg,
                 (Cint,     QuEST_Types.QuESTEnv),
                 numQubits, env)
end

function destroyComplexMatrixN(M ::QuEST_Types.ComplexMatrixN) ::Nothing
    ccall(:destroyComplexMatrixN, Cvoid,
          (QuEST_Types.ComplexMatrixN,),
          M)
    nothing
end

function destroyDiagonalOp(op       :: QuEST_Types.DiagonalOp,
                           env      :: QuEST_Types.QuESTEnv)    :: Nothing
    ccall(:destroyDiagonalOp, Cvoid,
          (QuEST_Types.DiagonalOp, QuEST_Types.QuESTEnv),
          op,                      env)
    return nothing
end

function destroyPauliHamil(hamil        ::QuEST_Types.PauliHamil)   :: Nothing
    ccall(:destroyPauliHamil, Cvoid, (QuEST_Types.PauliHamil,), hamil)
    return nothing
end

function destroyQuESTEnv(env ::QuEST_Types.QuESTEnv) ::Nothing
    ccall(:destroyQuESTEnv, Cvoid, (QuEST_Types.QuESTEnv,), env)
    return nothing
end

function destroyQureg(qureg ::QuEST_Types.Qureg, env ::QuEST_Types.QuESTEnv) ::Nothing
    ccall(:destroyQureg, Cvoid, (QuEST_Types.Qureg,QuEST_Types.QuESTEnv), qureg, env)
    return nothing
end

@doc raw"
Function `make_QuEST_matrix(M ::Matrix{Qreal}) ::ComplexMatrixN`
Convenience function for creating (`createComplexMatrixN()`) and filling a QuEST matrix data structure (`struct ComplexMatrixN`).
### Input
* `M` must be a 2ⁿ×2ⁿ matrix for 1 ≤ 𝑛 ≤ 50 (haha) [`@assert`]
### Return value
* An 𝑛 qubit matrix.
"
function make_QuEST_matrix(M ::Matrix{Base.Complex{Qreal}}) ::QuEST_Types.ComplexMatrixN
    (R,C) = size(M)
    @assert R==C
    @assert R ≥ 2
    numQubits = Int(round( log2(R) ))
    @assert 2^numQubits == R
    @assert 1 ≤ numQubits ≤ 50

    MQ = createComplexMatrixN(numQubits)
    for c = 1:R                       # Julia is column major
        for r = 1:R
            re_MQ_r = unsafe_load(MQ.real,r)
            im_MQ_r = unsafe_load(MQ.imag,r)
            unsafe_store!(re_MQ_r, real(M[r,c]), c)
            unsafe_store!(im_MQ_r, imag(M[r,c]), c)
        end
    end
    return MQ
end

@doc raw"
Function
     fill_ComplexMatrix!(M ::ComplexMatrixN, F ::Function) ::Nothing
Fills the entries of `M`: ∀(k,ℓ): 𝑀[k,ℓ] = 𝐹(k,ℓ).
(Indices are in 1, ..., 2ⁿ.)
### Input
* `M` must have been ''created'' (`createComplexMatrixN()`)
### Output
* Entries of `M` are overwritten.
"
function fill_ComplexMatrix!(M ::QuEST_Types.ComplexMatrixN, M_ ::Function) ::Nothing
    N = 2^M.numQubits
    for k = 1:N
        re_M_k = unsafe_load(M.real,k)
        im_M_k = unsafe_load(M.imag,k)
        for ℓ = 1:N
            unsafe_store!(re_M_k, real(M_(k,ℓ)), ℓ)
            unsafe_store!(im_M_k, imag(M_(k,ℓ)), ℓ)
        end
    end
    nothing
end

function initComplexMatrixN(m       ::QuEST_Types.ComplexMatrixN,
                            real_   ::Vector{Qreal},
                            imag_   ::Vector{Qreal})        :: Nothing

    @assert 1<<(2*m.numQubits) == length(real_) == length(imag_)
    ccall(:initComplexMatrixN, Cvoid, (QuEST_Types.ComplexMatrixN, Ptr{Qreal}, Ptr{Qreal}), m, real_, imag_)
    return nothing
end

function initDiagonalOp(op      :: QuEST_Types.DiagonalOp,
                        real_   :: Vector{Qreal},
                        imag_   :: Vector{Qreal})  ::Nothing

    @assert 1<<op.numQubits == length(real_) == length(imag_)
    ccall(:initDiagonalOp, Cvoid, (QuEST_Types.DiagonalOp, Ptr{Qreal}, Ptr{Qreal}), op, real_, imag_)
    return nothing
end

function initPauliHamil(hamil       :: QuEST_Types.PauliHamil,
                        coeffs      :: Vector{Qreal},
                        codes       :: Vector{QuEST_Types.pauliOpType})    ::Nothing

    @assert length(codes) == hamil.numSumTerms*hamil.numQubits
    ccall(:initPauliHamil, Cvoid, (QuEST_Types.PauliHamil, Ptr{Qreal}, Ptr{QuEST_Types.pauliOpType}), hamil, coeffs, codes)
    return nothing
end

function setDiagonalOpElems(op          ::QuEST_Types.DiagonalOp,
    startInd    ::Integer,
    real_       ::Vector{Qreal},
    imag_       ::Vector{Qreal},
    numElems    ::Integer)  ::Nothing
ccall(:setDiagonalOpElems,
Cvoid,
(QuEST_Types.DiagonalOp, Clonglong, Ptr{Qreal}, Ptr{Qreal}, Clonglong),
op,
Clonglong(startInd),
real_,
imag_,
Clonglong(numElems))
return nothing
end

function syncDiagonalOp(op        ::QuEST_Types.DiagonalOp)     ::Nothing
    ccall(:syncDiagonalOp, Cvoid, (QuEST_Types.DiagonalOp,), op)
    return nothing
end


text:
# Population pyramids

# Vegalite

using VegaLite, VegaDatasets

dataset("population") |>
@vlplot(
    width=300,
    height=200,
    :bar,
    transform=[
        {filter="datum.year==2000"},
        {calculate="datum.sex==2 ? 'Female' : 'Male'", as="gender"},
        {calculate="datum.sex==2 ? -datum.people : datum.people", "as"="signed_people" }
    ],
    x={"sum(signed_people)", axis={title="population", format=:s} },
    y={"age:o", axis=nothing, sort="descending" },
    color={"gender:n", scale={range=["#675193", "#ca8861"]},legend={orient=:top,title=nothing}},
    config={view={stroke=nothing},axis={grid=false}}
)

ext: jl

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
function _is_up_to_day(repodir, url)
    rhash = _remote_HEAD_hash(url)
    chash = _HEAD_hash(repodir)
    return rhash == chash
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Description
#
#   Miscellaneous functions used by all back-ends.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    function _str_escaped(str::AbstractString)

Return the escaped string representation of `str`.

"""
_str_escaped(str::AbstractString) =
    # NOTE: Here we cannot use `escape_string(str)` because it also adds the
    # character `"` to the list of characters to be escaped.
    sprint(escape_string, str, "", sizehint = lastindex(str))


text:
function euler!(t1,t2,dt,tm,c,daymax,n::Int64,rk,dc::Function)

#=**************************************************************************
This is a function in an ozone NOx formaldehyde mechanism, adapted from Fortran
by Obin Sturm (UC Davis) in October 2019.
The original comments and file description are below.

    created by: Mike Kleeman (April 1999)
                ECI 298 Air Quality Modeling
                University of California, Davis

    The purpose of this subroutine is to use an explicit Euler method to
    integrate a set of coupled ordinary differential equations between t1
    and t2.

    Inputs:
     t1   - start time
     t2   - stop time
     dt   - fixed integration time step
     tm   - concentration recording timestep
     c    - initial active species concentrations
     n    - number of active species
     diff - name of the subroutine to evaluate the derivative dc/dt = diff

    Outputs:
     C    -  species concentrations reported at specific time intervals
**************************************************************************=#
    #include("modlspc.jl")
    f = zeros(1,n)
    C = zeros(ceil(Int,t2/tm)+1,n)
    C[1,:] = copy(c[1,:])
    S = zeros(ceil(Int,t2/tm),length(r))
    s_aggregate = zeros(Float64,1,maxrxn)
    J = zeros(ceil(Int,t2/tm)+1,1)
    J[1,1] = copy(constn[2])

    #return C
    #c = zeros(1,n)
# preliminary steps
    if n > maxrxn
        error("too many reactions")
    end
    Time = t1
# loop over the time step until we come to t2
    icount=Int64(1)

    while (Time < t2)
        # oscillating HV value during the day
        constn[2] = daymax*max(sin(Time/720*pi - pi/2),0.0)

        f = dc(c,f,n)
        #       print*,'report from euler:c,f'
        #       do i=1,n
        #        print*,i,c(i),f(i)
        #       enddo
        #       stop 'euler marker1'

        c .= c .+ f*dt
        if any(c.<0.0)
            c[c.<0] .= 0
            #s_aggregate[1] = NaN  # Added March 2020: marker to show possible mass balance violation
            #println("Mass balance prob")
        end

        Time = Time + dt
        s_aggregate = s_aggregate .+ r.*dt

        if mod(icount,floor(Int,1000*tm)) == 0
            C[icount ÷ floor(Int,1000*tm)+1,:] = copy(c[1,:])
            S[icount ÷ floor(Int,1000*tm),:]   = copy(s_aggregate)
            J[icount ÷ floor(Int,1000*tm)+1,1] = copy(constn[2])
            s_aggregate = zeros(Float64,1,maxrxn)
        end
        icount=icount+1

    end

    return C,S[:,1:10],J
end

function diff!(c,f,n)

#=**************************************************************************
    This is a function in an ozone NOx formaldehyde mechanism, adapted from Fortran
    by Obin Sturm (UC Davis) in October 2019.
    The original comments and file description are below.
    created by: Mike Kleeman (April 1999)
             ECI 298 Air Quality Modeling
             University of California, Davis

 The purpose of this subroutine is to evaluate the right hand side of
 the equation
                    dc/dt = f
 Inputs:
  c    - active species concentrations
  n    - number of active species

 Outputs:
  f    - net rate of change of c
**************************************************************************=#
    #include("modlspc.jl")

# call the subroutine to evaluate the rates for the chemical reaction sys
    global s, r, fr, rlr = difun!(constn,c,rpssa,rk,r,fr,rlr)
    f[maxact+1:(maxbo+maxact)] = bldup!(r,f[maxact+1:(maxact+maxbo)])

# calculate the net rate of reaction for active species
    f[1:maxact] .= fr[1:maxact] .- rlr[1:maxact].*c[1:maxact]
    return f
end

function jcalc(c,n,dc,rj)
#=**************************************************************************
This is a function in an ozone NOx formaldehyde mechanism, adapted from Fortran
by Obin Sturm (UC Davis) in October 2019.

Note that this subroutine was never called in the original mechanism.

The original comments and file description are below.

    created by: Mike Kleeman (April 1999)
                ECI 298 Air Quality Modeling
                University of California, Davis

   The purpose of this subroutine is to numerically evaluate the Jacobian
   of the right hand side of the equation
                      dc/dt = f

   Inputs:
    c   - active species concentration
    n   - number of active species
    dc  - name for the function f

   Outputs:
    rj  - jacobian matrix for function f
**************************************************************************=#
end


text:
using Pkg
Pkg.activate(".")

using BenchmarkTools
using Distributions
using Flux.Optimise:update!
using Revise
using Flux.Optimise
using Random
using BSON: @save, @load
using LaTeXStrings
using ColorSchemes
using Match
using Plots

Random.seed!(123);

includet("./models/robust_rnn.jl")
includet("./models/solvers.jl")
includet("./models/REN.jl")
includet("./models/ffREN.jl")
includet("./models/models.jl")
includet("utils.jl")
includet("./data_handler.jl")
includet("train.jl")

plotlyjs()

function save_results(file_name, model, tloss, val_res, gamma,
                      model_options, data_options, train_options)

    set = data_options.set
    path = string("./results/sys_id/", set, "/")

    try
        print("Making directory: ", path, "...")
        mkdir(path)
        println("Done!")
    catch
        println("Directory already exists")
    end

    print("Saving model:", path,  file_name, ".params...")
    @save string(path, file_name, ".params") model
    println("Done")

    save_data = (model_options = model_options,
                 data_options = data_options,
                 train_options = train_options,
                 gamma_lower = gamma,
                 tloss = tloss,
                 val = val_res)

    print("Saving experimental data:", path,  file_name, ".bson...")
    @save string(path, file_name, ".bson") save_data
    println("Done")
end

function run_sys_id_test(file_name, model_options, data_options, train_options)

    Random.seed!(123);

    download_and_extract(data_options.set)
    train, val, stats    = @match data_options.set begin
        "cascaded_tanks" => load_cascaded_tanks()
        "silverbox" => load_silverbox()
        "wing_flutter" => load_wing_flutter()
        "f16" => load_f16()
        "wiener_hammerstein" => load_wiener_hammerstein()
        "heat_exchanger" => load_liquid_saturated_heat_exchanger()
    end

    println("Making model: ", model_options...)
    model = @match model_options.model begin
        "bounded_ren" => bounded_ren(model_options.model_args...)
        "stable_ren" => stable_ren(model_options.model_args...)
        "bounded_ffren" => bounded_ff_rnn(model_options.model_args...)
        "stable_ffren" => stable_ff_rnn(model_options.model_args...)
        "lstm" => lstm(model_options.model_args...)
        "rnn" => rnn(model_options.model_args...)
        "robust_rnn_stable" => robust_rnn_stable(model_options.model_args...)
        "robust_rnn_lipschitz" => robust_rnn_lipschitz(model_options.model_args...)
    end

    # Make training and validation datasets
    ut = Iterators.partition(train[1], data_options.seq_len)
    yt = Iterators.partition(train[2], data_options.seq_len)
    train_data = zip(Iterators.partition(train[2], data_options.seq_len))

    train_data = zip(ut, yt)
    val_data = [val]

    if model_options.model == "robust_rnn_stable" || model_options.model == "robust_rnn_lipschitz"
        opt = Flux.Optimise.ADAM(train_options.schedule.η)
        tloss, vloss = snlsdp!(train_data, train_data, model, opt; Epochs=train_options.Epochs, β0=1000)
    else

        # Account for number of steps in epoch
        sched = train_options.schedule
        opt = Flux.Optimiser(Flux.Optimise.ExpDecay(sched.η, sched.decay_rate,
                                                        length(train_data) * sched.decay_steps,
                                                        sched.min_lr),
                                    ClipValue(train_options.clip_grad),
                                    Flux.Optimise.ADAM(train_options.schedule.η))

        tloss = my_train!(train_data, val_data, model, opt; Epochs=train_options.Epochs)
    end
    val_res = validate(model, val_data, stats)
    println(val_res[1]["nrmse"])
    gamma = estimate_lipschitz_lower(model; seq_len=200, maxIter=10,step_size=1E-1)

    save_results(file_name, model, tloss, val_res, gamma, model_options, data_options, train_options)
    return model, tloss
end

function load_sys_id_test(file_path)
    @load string(file_path, ".params") model
    @load string(file_path, ".bson") save_data
    return model, save_data
end

function recalculate_lipschitz(file_path; step_size=1E-2, clip_at=1E-2, init_var=1E-5)
    # Recalulates Lipschtiz lower bound and saves.
    model, results = load_sys_id_test(file_path)
    gamma_lower = estimate_lipschitz_lower(model; step_size=step_size, clip_at=clip_at, init_var=init_var, maxIter=300, seq_len=3000)

    save_data = (model_options = results.model_options,
                 data_options = results.data_options,
                 train_options = results.train_options,
                 gamma_lower = gamma_lower,
                 tloss = results.tloss,
                 val = results.val)

    @save string(file_path, ".bson") save_data
    return gamma_lower
end

# # Silverbox
nx = 10
nv = 40
ϵ = 1E-2
solver = PeacemanRachford(tol=1E-4, α=1.0, maxIter=2000, verbose=false, cg=false)
(model::implicit_rnn)(x0, u) = model(x0, u, solver)  # capture solver

# Silver box experiments
data_options = (seq_len = 1024, set = "silverbox")
train_options = (Epochs = 50, clip_grad = 1E-1, seed = 123,
                    schedule = (η = 1E-3, decay_rate = 0.1, decay_steps = 10, min_lr = 1E-6))

model_options = (model = "stable",
                model_args = (nu = 1, ny = 1, nx = nx, nv = nv))

model_name = string("stable_rren_", nx, "_", nv)
M = run_sys_id_test(model_name, model_options, data_options, train_options)

for γ in [1.0, 2.5, 3.5, 5.0, 10.0]
    model_options = (model = "bounded",
                    model_args = (nu = 1, ny = 1, nx = nx, nv = nv, γ = γ))

    model_name = string("bounded_rren_", nx, "_", nv, "_", γ)
    run_sys_id_test(model_name, model_options, data_options, train_options)
end

model_name = "lstm_30"
model_options = (model = "lstm", model_args = (nu = 1, nv = 30, ny = 1))
run_sys_id_test(model_name, model_options, data_options, train_options)

model_name = "rnn_60"
model_options = (model = "rnn", model_args = (nu = 1, nv = 60, ny = 1))
run_sys_id_test(model_name, model_options, data_options, train_options)

# # F16
nx = 75
nv = 150
solver = PeacemanRachford(tol=1E-4, α=0.5, maxIter=10000, verbose=false, cg=false)
(model::implicit_param)(x0, u) = model(x0, u, solver)  # capture solver

data_options = (seq_len = 1024, set = "f16")
train_options = (Epochs = 70, clip_grad = 1E-1, seed = 123,
                    schedule = (η = 1E-3, decay_rate = 0.1, decay_steps = 20, min_lr = 1E-6))

model_options = (model = "stable_ffren",
                    model_args = (nu = 2, nx = nx, nv = nv, ny = 3))
model_name = string("stable_ffren_", nx, "_", nv)
run_sys_id_test(model_name, model_options, data_options, train_options)

model_options = (model = "stable_ren",
                model_args = (nu = 2, nx = nx, nv = nv, ny = 3))
model_name = string("stable_ren_", nx, "_", nv)
run_sys_id_test(model_name, model_options, data_options, train_options)

for γ in [4.0, 6.0, 10.0, 20.0, 20.0, 40.0, 60.0, 100.0]
    model_options = (model = "bounded_ren",
                    model_args = (nu = 2, nx = nx, nv = nv, ny = 3, γ = γ))
    model_name = string("bounded_ren_", nx, "_", nv, "_", γ)
    run_sys_id_test(model_name, model_options, data_options, train_options)

    model_options = (model = "bounded_ffren",
    model_args = (nu = 2, nx = nx, nv = nv, ny = 3, γ = γ))
    model_name = string("bounded_ffren_", nx, "_", nv, "_", γ)
    run_sys_id_test(model_name, model_options, data_options, train_options)
end

model_options = (model = "lstm", model_args = (nu = 2, nv = 170, ny = 3))
model_name = string("lstm_", model_options.model_args.nv)
run_sys_id_test(model_name, model_options, data_options, train_options)

train_options = (Epochs = 70, clip_grad = 1E-1, seed = 123,
                    schedule = (η = 1E-4, decay_rate = 0.1, decay_steps = 20, min_lr = 1E-6))

model_options = (model = "rnn", model_args = (nu = 2, nv = 340, ny = 3))
model_name = string("rnn_", model_options.model_args.nv)
run_sys_id_test(model_name, model_options, data_options, train_options)

train_options = (Epochs = 70, clip_grad = 1E-1, seed = 123,
                    schedule = (η = 1E-4, decay_rate = 0.1, decay_steps = 20, min_lr = 1E-6))

model_options = (model = "robust_rnn_stable",
                model_args = (nu = 2, nx = nx, nv = nv, ny = 3))
model_name = string("robust_rnn_stable", nx, "_", nv)
run_sys_id_test(model_name, model_options, data_options, train_options)

for γ in [4.0, 6.0, 10.0, 20.0, 20.0, 40.0, 60.0, 100.0]
    model_options = (model = "robust_rnn_lipschitz",
                    model_args = (nu = 2, nx = nx, nv = nv, ny = 3, γ = γ))
    model_name = string("robust_rnn_lipschitz", nx, "_", nv, "_", γ)
    run_sys_id_test(model_name, model_options, data_options, train_options)
end

# Wiener Hammerstein
nx = 40
nv = 100
solver = PeacemanRachford(tol=1E-6, α=1.0, maxIter=2000, verbose=false, cg=false)
(model::implicit_rnn)(x0, u) = model(x0, u, solver)  # capture solver

data_options = (seq_len = 512, set = "wiener_hammerstein")
train_options = (Epochs = 60, clip_grad = 1E-1, seed = 123,
                    schedule = (η = 1E-3, decay_rate = 0.1, decay_steps = 40, min_lr = 1E-6))

model_options = (model = "stable_ffren",
                model_args = (nu = 2, nx = nx, nv = nv, ny = 1))
model_name = string("stable_ffren_", nx, "_", nv)
M, tloss = run_sys_id_test(model_name, model_options, data_options, train_options)

for γ in [1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 5.0, 8.0]
    model_options = (model = "bounded_ffren",
                    model_args = (nu = 2, nx = nx, nv = nv, ny = 1, γ = γ))

    model_name = string("bounded_ffren_", nx, "_", nv, "_", γ)
    run_sys_id_test(model_name, model_options, data_options, train_options)
end

model_options = (model = "lstm", model_args = (nu = 2, nv = 100, ny = 1))
model_name = string("lstm_", model_options.model_args.nv)
run_sys_id_test(model_name, model_options, data_options, train_options)

model_options = (model = "rnn", model_args = (nu = 2, nv = 200, ny = 1))
model_name = string("rnn_", model_options.model_args.nv)
run_sys_id_test(model_name, model_options, data_options, train_options)

train_options = (Epochs = 250, clip_grad = 1E-1, seed = 123,
                    schedule = (η = 1E-3, decay_rate = 0.1, decay_steps = 40, min_lr = 1E-6))

model_options = (model = "robust_rnn_stable",
                model_args = (nu = 2, nx = nx, nv = nv, ny = 1))
model_name = string("robust_rnn_stable", nx, "_", nv)
run_sys_id_test(model_name, model_options, data_options, train_options)

for γ in [1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 5.0, 8.0]
    model_options = (model = "robust_rnn_lipschitz",
                    model_args = (nu = 2, nx = nx, nv = nv, ny = 1, γ = γ))
    model_name = string("robust_rnn_lipschitz", nx, "_", nv, "_", γ)
    run_sys_id_test(model_name, model_options, data_options, train_options)
end

recalculate_lipschitz("./results/sys_id/heat_exchanger/regularized_rren_40_40")
recalculate_lipschitz("./results/sys_id/heat_exchanger/stable_rren_40_40")
recalculate_lipschitz("./results/sys_id/heat_exchanger/lstm_65")
recalculate_lipschitz("./results/sys_id/heat_exchanger/rnn_130")

# Heat exchanger
nx = 40
nv = 40
solver = PeacemanRachford(tol=1E-4, α=1.0, maxIter=10000, verbose=false, cg=false)
(model::implicit_rnn)(x0, u) = model(x0, u, solver)  # capture solver

data_options = (seq_len = 50, set = "heat_exchanger")
train_options = (Epochs = 250, clip_grad = 1E-1, seed = 123,
                    schedule = (η = 1E-3, decay_rate = 0.5, decay_steps = 200, min_lr = 1E-4))

model_options = (model = "stable", model_args = (nu = 1, ny = 1, nx = nx, nv = nv))
model_name = string("stable_rren_", nx, "_", nv)
run_sys_id_test(model_name, model_options, data_options, train_options)

model_options = (model = "regularized", model_args = (nu = 1, ny = 1, nx = nx, nv = nv, γ = 2.0))
model_name = string("regularized_rren_", nx, "_", nv)
run_sys_id_test(model_name, model_options, data_options, train_options)

model_options = (model = "lstm", model_args = (nu = 1, nv = 65, ny = 1))
model_name = string("lstm_", model_options.model_args.nv)
run_sys_id_test(model_name, model_options, data_options, train_options)

model_options = (model = "rnn", model_args = (nu = 1, nv = 130, ny = 1))
model_name = string("rnn_", model_options.model_args.nv)
run_sys_id_test(model_name, model_options, data_options, train_options)

#  Run these to recalculate empirical lipschitz lower bounds
recalculate_lipschitz("./results/sys_id/f16/rnn_340")
recalculate_lipschitz("./results/sys_id/f16/lstm_170")

recalculate_lipschitz("./results/sys_id/f16/bounded_rren_75_150_10.0")
recalculate_lipschitz("./results/sys_id/f16/bounded_rren_75_150_20.0")
recalculate_lipschitz("./results/sys_id/f16/bounded_rren_75_150_40.0")
recalculate_lipschitz("./results/sys_id/f16/bounded_rren_75_150_100.0")
recalculate_lipschitz("./results/sys_id/f16/bounded_rren_75_150_60.0")
recalculate_lipschitz("./results/sys_id/f16/stable_rren_75_150")

recalculate_lipschitz("./results/sys_id/f16/bounded_ffren_75_150_10.0")
recalculate_lipschitz("./results/sys_id/f16/bounded_ffren_75_150_20.0")
recalculate_lipschitz("./results/sys_id/f16/bounded_ffren_75_150_40.0")
recalculate_lipschitz("./results/sys_id/f16/bounded_ffren_75_150_100.0")
recalculate_lipschitz("./results/sys_id/f16/bounded_ffren_75_150_60.0")
recalculate_lipschitz("./results/sys_id/f16/stable_ffren_75_150")

recalculate_lipschitz("./results/sys_id/f16/robust_rnn_lipschitz75_150_4.0")
recalculate_lipschitz("./results/sys_id/f16/robust_rnn_lipschitz75_150_6.0")
recalculate_lipschitz("./results/sys_id/f16/robust_rnn_lipschitz75_150_10.0")
recalculate_lipschitz("./results/sys_id/f16/robust_rnn_lipschitz75_150_20.0")
recalculate_lipschitz("./results/sys_id/f16/robust_rnn_lipschitz75_150_40.0")
recalculate_lipschitz("./results/sys_id/f16/robust_rnn_lipschitz75_150_60.0")
recalculate_lipschitz("./results/sys_id/f16/robust_rnn_lipschitz75_150_100.0")
recalculate_lipschitz("./results/sys_id/f16/robust_rnn_stable75_150")

recalculate_lipschitz("./results/sys_id/wiener_hammerstein/rnn_200")
recalculate_lipschitz("./results/sys_id/wiener_hammerstein/lstm_100")
recalculate_lipschitz("./results/sys_id/wiener_hammerstein/bounded_ffren_40_100_1.5")
recalculate_lipschitz("./results/sys_id/wiener_hammerstein/bounded_ffren_40_100_2.5")
recalculate_lipschitz("./results/sys_id/wiener_hammerstein/bounded_ffren_40_100_3.5")
recalculate_lipschitz("./results/sys_id/wiener_hammerstein/bounded_ffren_40_100_8.0")
recalculate_lipschitz("./results/sys_id/wiener_hammerstein/bounded_ffren_40_100_15.0")
recalculate_lipschitz("./results/sys_id/wiener_hammerstein/stable_ffren_40_100")

max_stars_repo_name: imanchester/REN

max_issues_repo_name: imanchester/REN

max_forks_repo_name: imanchester/REN


text:
#module NaquadahEvents
# using Cairo, Gtk, Gtk.ShortNames

export  splotch, ScrollEvent, # not sure all this needs exported (?)
        DragEvent, MouseDownEvent, ClickEvent #, Event,

# ======================================================================================
#
# ======================================================================================
function DragEvent(document, widget, event)
    pressed, released = document.event.pressed[1], document.event.pressed[end]
    splotch(widget, pressed,1.0,0.0,0.0)

    ctx = getgc(widget)
    move_to(ctx,pressed.x, pressed.y)
    line_to(ctx,released.x, released.y)
    stroke(ctx)

    splotch(widget, released,1.0,0.0,0.0)
    println("drag from X $(pressed.x), Y $(pressed.y) to X $(released.x), Y $(released.y)")
    document.event.pressed = []
end

# ======================================================================================
#
# ======================================================================================
onArea(x, y,  shape ) = onArea(x, y,  getBorderBox(shape)... ) # WAS: getContentBox

function onArea(x, y,    l, t, w, h)
        return x > l && x < l+w && y > t && y < t+h ? true : false
end
# ======================================================================================
function insideParentClip(x, y, node)
    return onArea( x, y,  lastClipParent(node).shape )
end
# ======================================================================================
function order(a,b)
    a<b ? (a,b) : (b,a)
end
# ======================================================================================
function order2(l,r,t,b)
    l, r = order(l,r)
    t, b = order(t,b)
    return (l,r,t,b)
end
# ======================================================================================
function boxesOverlap(shape, l,r,t,b)
    return onArea(l, t, shape) || onArea(l, b, shape) || onArea(r, t, shape) || onArea(r, b, shape)
end
# ======================================================================================
# Also create: function findNodesAtPoint(node, x, y)
# ======================================================================================
function findNodesInBox(node, left, top, right, bottom)
    items = []
    shape = node.shape
    box = getContentBox(shape)
    started = false

    # TODO: also include some test to exclude "unselectable" nodes.
    function looper(rows)
        for r in rows
            for n in r.nodes
                        if isdefined(n, :rows)
                            looper(n.rows)
                        else
                                s = n.shape
                                if onArea(left, top, s)         # ¡Selection Start
                                    push!(items, n)
                                    started = true
                                    if onArea(right, bottom, s)
                                        started = false
                                        return items
                                    end
                                elseif onArea(right, bottom, s) # Selection End!
                                    push!(items, n)
                                    started = false
                                    return items
                                elseif started
                                    push!(items, n)
                                end

                        end
            end
        end
        return items
    end

    looper(node.rows)

    return items
end
# ======================================================================================
#
# a,b,c,d = [1],[2],[3],[4]
# l = [a b b c d a a]
# same(list, object) = list .=== object
# ======================================================================================
function selectText(ctx, document, pressed, released)

    left, right, top, bottom = order2(pressed.x,released.x, pressed.y,released.y)
    selected = findNodesInBox(document.children[1].children[3], left, top, right, bottom)

    clipPath = getContentBox(document.children[1].children[3].shape)
    # Objective here is to remove duplicates and also repaint tainted rows
    selectedNodes = document.event.selectedNodes

    found = false
for oldnode in selectedNodes
    for node in selected
        found = false
         if node ===  oldnode
             found = true
         end
    end
    if !found
        c = inheritColor(oldnode.parent)

        setcolor( ctx, c...)
        rectangle(ctx,  oldnode.shape.left-1, oldnode.shape.top-1, oldnode.shape.width+2, oldnode.shape.height+2 )
        fill(ctx)
        DrawText(ctx, oldnode, clipPath)
    end
end

    # buffer = []
    # for shape in selected
    #     if shape.top > bottom || (shape.top+shape.height) < top # outside range and needs removed
    #         setcolor( ctx, shape.parent.shape.color...)
    #         rectangle(ctx,  shape.left, shape.top, shape.width, shape.height )
    #         fill(ctx)
    #         DrawText(ctx, shape, clipPath)
    #     else
    #         push!(buffer,shape)
    #     end
    # end
    # selected = buffer

    if pressed.y < released.y # just in case the selection is made from bottom to top
        DrawSelectedText(ctx, selected, pressed.x, released.x, pressed.y, released.y, clipPath)
    else
        DrawSelectedText(ctx, selected, released.x, pressed.x, released.y, pressed.y, clipPath)
    end
    document.event.selectedNodes = selected

end
# ======================================================================================
#
# ======================================================================================
function DrawEvent(document, widget, event)
    println("move x$(event.x) y$(event.y)")
end
# ======================================================================================
#
# ======================================================================================
function MouseDownEvent(document, widget, event)
    down = document.eventsList.mousedown
    #println("clicked: ", event.x, event.y)
    for node in down

        if onArea(event.x, event.y, node.shape ) && insideParentClip( event.x, event.y, node)
            ctx = getgc(widget)
            page = document.children[1].children[3]

            if node.shape.color[2] > 0.9
                node.temp = deepcopy(node.shape)
                node.shape.color = [0.996094, 0.5, 0.5]
            else
                node.shape = deepcopy(node.temp)
            end
            drawNode(ctx, document, node)
            reveal(widget)
        end
    end
end
#setAttribute(node, "color", "red")
#AtributesToLayout(document, node, false)
# ======================================================================================

# ======================================================================================
function setAttribute(d, key::String, attribute::String)
    println(d.DOM)
    d.DOM[key] = attribute
end
# ======================================================================================
#
# ======================================================================================
function ClickEvent(document, widget, event)
    clicked = document.eventsList.click
    splotch(widget, event, 1.0,0.0,0.0)
    document.event.pressed = []
    for node in clicked
        if onArea(event.x, event.y, node.shape ) && insideParentClip( event.x, event.y, node)
            if haskey(node.DOM, "click")
                click = node.DOM["click"]
                println(click)
                # if haskey(click, "code")
                    code = Meta.parse(click)
                    ex = :($(Expr(:toplevel, :(document = $document), :($code))))
                    eval(ex)
                    # WAS: eval(parse(click["code"]))
                # end
            end
            #"click":{"code":"println(\"This is some code!\")", "preventDefault":true}

        end
    end

end
# ======================================================================================
#
# ======================================================================================
function OnObject(document, x, y, all=true)

end
# ======================================================================================
#
# ======================================================================================
function ScrollEvent(document, widget, event)
  ctx = getgc(widget)
  node = document.children[1].children[3] # TODO: fix! ..test to see if mouse is over object.

  Unit = 70.0

  # I am scrolling(jumping) by 30px here but Opera scrolls by about 50px
  # Opera lacks smoothness too but it seems to transition-scroll by the 50px
  # ...so using the mouse wheel it is impossible to move less than that increment.
  function drawTwice()
      setcolor(ctx, node.shape.color...)
      rectangle(ctx,  node.shape.left,  node.shape.top, node.shape.width,  node.shape.height )
      fill(ctx);

      DrawViewport(ctx, document, node)
      reveal(widget)
  end

  # SCROLL UP!
  if event.direction == 0 && node.scroll.y < node.shape.top # was: 0
    diff = abs(node.scroll.y)
    if diff < Unit
      Unit = diff
    end
    node.scroll.y += (Unit*.5)
    VmoveAllChildren(node, (Unit*.5), false)
    drawTwice()
    node.scroll.y += (Unit*.5)
    VmoveAllChildren(node, (Unit*.5), false)
    drawTwice()
    # SCROLL DOWN!
elseif event.direction == 1 && (node.scroll.contentHeight + node.scroll.y + node.shape.top) > node.shape.height
      diff = (node.scroll.contentHeight + node.scroll.y + node.shape.top) - node.shape.height
      if diff < Unit
        Unit = diff
      end
      node.scroll.y -= (Unit*.5)
      VmoveAllChildren(node, -(Unit*.5), false) # -30
      drawTwice()
      node.scroll.y -= (Unit*.5)
      VmoveAllChildren(node, -(Unit*.5), false) # -30
      drawTwice()
  end

  # #set_source_rgb(ctx, 1,1,1)
  # setcolor( ctx, node.shape.color...)
  # rectangle(ctx,  node.shape.left,  node.shape.top, node.shape.width,  node.shape.height )
  # fill(ctx);
  #
  # DrawViewport(ctx, document, node)
  # reveal(widget)

end

function drawTwice()
    setcolor( ctx, node.shape.color...)
    rectangle(ctx,  node.shape.left,  node.shape.top, node.shape.width,  node.shape.height )
    fill(ctx);

    DrawViewport(ctx, document, node)
    reveal(widget)
end
# ======================================================================================
# Draw a splotch
# CALLED FROM:
# ======================================================================================
function splotch(widget,event,r,g,b)
      ctx = getgc(widget)
            deg = (pi/180.0)
        s1, e1 = 1 * deg, 100 * deg
        s2, e2 = 120 * deg, 220 * deg
        s3, e3 = 240 * deg, 340 * deg
        move_to(ctx, event.x, event.y)
        set_line_width(ctx, 2.56);
        set_antialias(ctx,4)
        set_source_rgb(ctx, r, g, b)
        arc(ctx, event.x, event.y, 2, 0, 2pi) # 0, 2pi
        stroke(ctx)

        #set_line_width(ctx, 3);
        set_source_rgba(ctx, r, g, b, 0.7)
        arc(ctx, event.x, event.y, 5, s1, e1) # 0, 2pi
        stroke(ctx)
        arc(ctx, event.x, event.y, 5, s2, e2) # 0, 2pi
        stroke(ctx)
        arc(ctx, event.x, event.y, 5, s3, e3) # 0, 2pi
        stroke(ctx)
        # arc_negative(cr, xc, yc, node.shape.radius,
        #     node.shape.angle[1] * (pi/180.0),
        #     node.shape.angle[2] * (pi/180.0));

        #set_line_width(ctx, 2);
        set_source_rgba(ctx, r, g, b, 0.5)
        arc(ctx, event.x, event.y, 8, s1, e1)
        stroke(ctx)
        arc(ctx, event.x, event.y, 8, s2, e2)
        stroke(ctx)
        arc(ctx, event.x, event.y, 8, s3, e3)
        stroke(ctx)

        #set_line_width(ctx, 1);
        set_source_rgba(ctx, r, g, b, 0.3)
        arc(ctx, event.x, event.y, 11, s1, e1)
        stroke(ctx)
        arc(ctx, event.x, event.y, 11, s2, e2)
        stroke(ctx)
        arc(ctx, event.x, event.y, 11, s3, e3)
        stroke(ctx)
        reveal(widget)
    end
#end # module


text:
# This file provides an introductory tour for Legolas.jl. Run lines in the REPL
# to inspect output at each step! Tests are littered throughout to demonstrate
# functionality in a concrete manner, and so that we can ensure examples stay
# current as the package evolves.

using Legolas, Arrow, Tables, Test

using Legolas: @row, Schema

#####
##### Introduction to the `@row` Macro and `Legolas.Row` constructors
#####
# The most interesting bit of Legolas.jl functionality is the package's `@row` macro, which callers can use
# define new [`Tables.AbstractRow`-compliant](https://tables.juliadata.org/stable/#Tables.AbstractRow-1)
# row types that exhibit some opinionated (but desirable!) properties w.r.t. composability/extensibility.
#
# These row type properties are fundamental to the wider data curation patterns that Legolas.jl seeks to
# facilitate, so let's explore them before we dig deeper into Legolas.jl's other table-centric utilities.

# Declare a `Legolas.Schema` with name `"my-schema"` at version `1` whose *required fields*
# `a`, `b`, `c`, `d`, and `e` are specified via the provided assignment expressions, then return
# a corresponding row type that matches the declared schema:
const MyRow = @row("my-schema@1",
                   a::Real = sin(a),
                   b::String = string(a, b, c),
                   c = [a, b, c],
                   d::Int,
                   e)

# `MyRow` is thus an alias for the type returned by the `@row` macro:
@test MyRow == Legolas.Row{typeof(Schema("my-schema@1"))}

# The `MyRow` type has several useful constructors. Let's start with the constructor that
# accepts all required fields as keyword arguments:
row = MyRow(a=1.5, b="hello", c="goodbye", d=2, e=["anything"])

# By examining `row`'s fields, we can see how the assignment expressions from `MyRow`'s
# `@row` declaration were applied in a simple linear fashion to the input arguments:
@test row.a == sin(1.5)
@test row.b == string(sin(1.5), "hello", "goodbye")
@test row.c == [sin(1.5), string(sin(1.5), "hello", "goodbye"), "goodbye"]
@test row.d == 2
@test row.e == ["anything"]

# In fact, the field assignment expressions provided to the `@row` macro are interpolated
# nearly as-is into the underlying code generated by `@row`. For example, the relevant code
# generated by `MyRow`'s `@row` declaration looks roughly like:
#
#     function Legolas._transform(::typeof(Legolas.Schema("my-schema", 1));
#                                 a=missing, b=missing, c=missing, d=missing, e=missing,
#                                 other...)
#         a::Real = sin(a)
#         b::String = string(a, b, c)
#         c::Any = [a, b, c]
#         d::Int = d
#         e::Any = e
#         return (; a, b, c, d, e, other...)
#     end
#
# This `Legolas._transform` method is invoked at the core of the `MyRow` constructor. As you
# might have noticed, this method has two interesting properties we haven't yet demonstrated.
#
# Here we demonstrate the first property - required fields have a `missing` value by default:
@test isequal(MyRow(a=1.5, b="hello", c="goodbye", d=2), MyRow((; a=1.5, b="hello", c="goodbye", d=2, e=missing)))
@test_throws MethodError MyRow(a=1.5, b="hello", c="goodbye") # correctly throws a `MethodError` when evaluating `d::Int = missing`

# And here's a demonstration of the second property - callers may pass in any other fields in
# addition to the required fields:
row = MyRow(a=1.5, b="hello", c="goodbye", my_other_field=":)", d=2, e=["anything"])
@test row.my_other_field == ":)"

# Finally, there's also a `MyRow` constructor that accepts any `Tables.AbstractRow`-compliant value,
# and extracts all input fields from that value. Here, we demonstrate with a `NamedTuple`:
@test row == MyRow((a=1.5, b="hello", c="goodbye", d=2, e=["anything"], my_other_field=":)"))

# To summarize:
#
# - Inputs to `Legolas.Row` constructors...
#     - ...may be any `Tables.AbstractRow`-compliant value
#     - ...may contain required fields in any order
#     - ...may elide required fields, in which case the constructor will assume them to be `missing`
#     - ...may contain any other fields in addition to the required fields
# - Outputs of `Legolas.Row` constructors...
#     - ...will contain all required fields ("missing" fields are explicitly presented with `missing` values)
#     - ...will contain all provided non-required fields

#####
##### Extending Existing Rows/Schemas
#####
# Row types declared via `@row` can inherit the fields and transformations specified by other `@row`-declared types.
# Niftily, the properties of `Legolas.Row` that we demonstrated above enable this extension mechanism to be
# implemented via composition under the hood.

# Declare a row type whose schema is named `"my-child-schema"` at version `1` that inherits the fields of the
# `my-schema@1` schema that we defined in the previous section.
const MyChildRow = @row("my-child-schema@1" > "my-schema@1",
                        d, # "declaring" the underlying `my-schema@1` field here so that it
                           # can be referenced in our definition for the `f` field.
                        f::Int = f + d,
                        g::String,
                        c = last(c))

# The constructor for `MyChildRow` will first apply its parent's transformation before applying its own. The
# effect of this behavior can be seen clearly in the `c` field value in the following example:
input = (a=1.5, b="hello", c="goodbye", d=2, e=["anything"], f=3, g="foo")
child = MyChildRow(input)
@test child.a == sin(1.5)
@test child.b == string(sin(1.5), "hello", "goodbye")
@test child.c == "goodbye"
@test child.d == 2
@test child.e == ["anything"]
@test child.f == 5
@test child.g == "foo"

# Note that even though we didn't write down any constraints on `d` in our `my-child-schema@1` definition,
# that field still undergoes the parent transformation (defined by `my-schema@1`) where it is constrained
# to `d::Int`.
@test_throws InexactError MyChildRow(Tables.rowmerge(child; d=1.5))

# A note on syntax: You might ask "Why use `>` as the inheritance operator instead of `<:`?" There are actually three reasons. Firstly,
# `<:` is canonically a subtyping operator that implies the Liskov substition principle, but because Legolas allow arbitrary RHS code in
# required field declarations, a "child" row is not de facto substitutable for its parent. Secondly, `>` implies the actual ordering that
# row transformations are applied in; the parent transformation comes before the child transformation. Thirdly, the child row will usually
# (though technically not always) have a greater total number of required fields than the parent row.

#####
##### Validating/Writing/Reading `Arrow.Table`s with Legolas.jl
#####

# Legolas provides special methods for reading/writing/validating Arrow tables that utilize `Legolas.Schema`s. To
# start exploring these methods, we'll first construct a dummy table using `row::MyRow` from the previous section:
rows = [row,
        Tables.rowmerge(row; b="a different one"),
        Tables.rowmerge(row; e=:anything)]

# We can validate that `rows` is compliant w.r.t. `schema` via `Legolas.validate`, which will throw a descriptive
# error if the `Tables.schema(rows)` mismatches with the required columns/types specified by the `schema`.
schema = Schema("my-schema", 1)
@test schema === Schema("my-schema@1") # `Schema("my-schema", 1)` is an alternative to `Schema("my-schema@1")`
@test (Legolas.validate(rows, schema); true)
invalid = vcat(rows, Tables.rowmerge(row; a="this violates the schema's `a::Real` requirement"))
@test_throws ArgumentError("field `a` has unexpected type; expected <:Real, found Any") Legolas.validate(invalid, schema)

# This highlights two important properties regarding `Legolas.Schema` validation:
#
# - First, it's okay that the `e` field isn't present in this `Tables.Schema` because `my-schema` permits `e::Missing`.
# - Second, field ordering is unimportant and is not considered when determining whether a give `Tables.`
@test (Legolas.validate(Tables.Schema((:c, :d, :a, :b), Tuple{Vector,Int,Float64,String}), schema); true)

# Legolas also provides `Legolas.write(path_or_io, table, schema; kwargs...)`, which wraps `Arrow.write(path_or_io, table; kwargs...)`
# and performs two additional operations atop the usual operations performed by that function:
#
# - By default, the provided Tables.jl-compliant `table` is validated against `schema` via `Legolas.validate` before
#   it is actually written out. Note that this can be disabled by passing `validate=false` to `Legolas.write`.
#
# - `Legolas.write` ensures that the written-out Arrow table's metadata contains a `"legolas_schema_qualified"`
#   key whose value is `Legolas.schema_qualified_string(schema)`. This field enables consumers of the table to
#   perform automated (or manual) schema discovery/evolution/validation.

schema = Schema("my-child-schema", 1) # For this example, let's use a schema that has a parent
rows = [child,
        Tables.rowmerge(child; b="a different one"),
        Tables.rowmerge(child; e=:anything)]
io = IOBuffer()
Legolas.write(io, rows, schema)
t = Arrow.Table(seekstart(io))

table_isequal(a, b) = isequal(Legolas.materialize(a), Legolas.materialize(b))

@test Arrow.getmetadata(t) == Dict("legolas_schema_qualified" => "my-child-schema@1>my-schema@1")
@test table_isequal(t, Arrow.Table(Arrow.tobuffer(rows)))
@test table_isequal(t, Arrow.Table(Legolas.tobuffer(rows, schema))) # `Legolas.tobuffer` is analogous to `Arrow.tobuffer`

invalid = vcat(rows, Tables.rowmerge(child; a="this violates the schema's `a::Real` requirement"))
@test_throws ArgumentError("field `a` has unexpected type; expected <:Real, found Any") Legolas.tobuffer(invalid, schema)

# Similarly, Legolas provides `Legolas.read(path_or_io)`, which wraps `Arrow.Table(path_or_io)`
# and performs `Legolas.validate` on the resulting `Arrow.Table` before returning it.
@test table_isequal(Legolas.read(Legolas.tobuffer(rows, schema)), t)
msg = """
      could not extract valid `Legolas.Schema` from provided Arrow table;
      is it missing the expected custom metadata and/or the expected
      \"legolas_schema_qualified\" field?
      """
@test_throws ArgumentError(msg) Legolas.read(Arrow.tobuffer(rows))
invalid_but_has_schema_metadata = Arrow.tobuffer(invalid;
                                                 metadata = ("legolas_schema_qualified" => Legolas.schema_qualified_string(schema),))
@test_throws ArgumentError("field `a` has unexpected type; expected <:Real, found Union{Missing, Float64, String}") Legolas.read(invalid_but_has_schema_metadata)

# A note about one additional benefit of `Legolas.read`/`Legolas.write`: Unlike their Arrow.jl counterparts,
# these functions are relatively agnostic to the types of provided path arguments. Generally, as long as a
# given `path` supports `Base.read(path)::Vector{UInt8}` and `Base.write(path, bytes::Vector{UInt8})` then
# `path` will work as an argument to `Legolas.read`/`Legolas.write`. At some point, we'd like to make similar
# upstream improvements to Arrow.jl to render its API more path-type-agnostic.


text:
# `tree_base(pd)` returns the same results as Profile.print().
# It's equivalent to `tree`, but ~3X faster, so I'm keeping it here.

function tree_base(bt::Vector{Vector{UInt64}}, counts::Vector{Int},
                   lidict::LineInfoFlatDict, level::Int, fmt::ProfileFormat, noisefloor::Int)
    if level > fmt.maxdepth
        return []::Any
    end
    # Organize backtraces into groups that are identical up to this level
    if fmt.combine
        # Combine based on the line information
        d = Dict{StackFrame,Vector{Int}}()
        for i = 1:length(bt)
            ip = bt[i][level + 1]
            key = lidict[ip]
            indx = Base.ht_keyindex(d, key)
            if haskey(d, key)
                push!(d[key], i)
            else
                d[key] = [i]
            end
        end
        # Generate counts
        dlen = length(d)
        lilist = Vector{StackFrame}(dlen)
        group = Vector{Vector{Int}}(dlen)
        n = Vector{Int}(dlen)
        i = 1
        for (key, v) in d
            lilist[i] = key
            group[i] = v
            n[i] = sum(counts[v])
            i += 1
        end
    else
        # Combine based on the instruction pointer
        d = Dict{UInt64,Vector{Int}}()
        for i = 1:length(bt)
            key = bt[i][level+1]
            if haskey(d, key)
                push!(d[key], i)
            else
                d[key] = [i]
            end
        end
        # Generate counts, and do the code lookup
        dlen = length(d)
        lilist = Vector{StackFrame}(dlen)
        group = Vector{Vector{Int}}(dlen)
        n = Vector{Int}(dlen)
        i = 1
        for (key, v) in d
            lilist[i] = lidict[key]
            group[i] = v
            n[i] = sum(counts[v])
            i += 1
        end
    end
    # Order the line information
    if length(lilist) > 1
        p = Profile.liperm(lilist)
        lilist = lilist[p]
        group = group[p]
        n = n[p]
    end
    # Recurse to the next level
    len = Int[length(x) for x in bt]
    out = Node[]
    for i = 1:length(lilist)
        n[i] < fmt.mincount && continue
        n[i] < noisefloor && continue
        idx = group[i]
        keep = len[idx] .> level+1
        if any(keep)
            idx = idx[keep]
            children = tree_base(bt[idx], counts[idx], lidict, level + 1, fmt,
                                 fmt.noisefloor > 0 ?
                                 floor(Int, fmt.noisefloor * sqrt(n[i])) : 0)
        else
            children = []
        end
        push!(out, Node(lilist[i], n[i], children))
    end
    return out
end

function tree_base(data::Vector{UInt64}, lidict::LineInfoFlatDict, fmt::ProfileFormat)
    if !fmt.C
        data = purgeC(data, lidict)
    end
    bt, counts = tree_aggregate(data)
    level = 0
    len = Int[length(x) for x in bt]
    keep = len .> 0
    # Using UNKNOWN as the root works, but it should ideally be a different flag value...
    # Or use a Nullable. - @cstjean
    return Node(UNKNOWN, -1, tree_base(bt[keep], counts[keep], lidict, level, fmt, 0))
end

function tree_base(data::Vector, lidict::LineInfoDict, fmt::ProfileFormat)
    newdata, newdict = flatten(data, lidict)
    return tree_base(newdata, newdict, fmt)
end

"""
    tree(pd::ProfileData; C = false, combine = true, maxdepth::Int = typemax(Int),
         mincount::Int = 0, noisefloor = 0)

Returns a tree view of the profiling data `pd`

The keyword arguments can be any combination of:

 - `C` -- If `true`, backtraces from C and Fortran code are shown (normally they are excluded).

 - `combine` -- If `true` (default), instruction pointers are merged that correspond to the same line of code.

 - `maxdepth` -- Limits the depth higher than `maxdepth` in the `:tree` format.

 - `sortedby` -- Controls the order in `:flat` format. `:filefuncline` (default) sorts by the source
    line, whereas `:count` sorts in order of number of collected samples.

 - `noisefloor` -- Limits frames that exceed the heuristic noise floor of the sample (only applies to format `:tree`).
    A suggested value to try for this is 2.0 (the default is 0). This parameter hides samples for which `n <= noisefloor * √N`,
    where `n` is the number of samples on this line, and `N` is the number of samples for the callee.

 - `mincount` -- Limits the printout to only those lines with at least `mincount` occurrences.
"""
function tree_base(pd::ProfileData;
                   C = false,
                   combine = true,
                   maxdepth::Int = typemax(Int),
                   mincount::Int = 0,
                   noisefloor = 0)
    tree_base(pd.data, pd.lidict, ProfileFormat(C = C,
            combine = combine,
            maxdepth = maxdepth,
            mincount = mincount,
            noisefloor = noisefloor))
end


text:
using Statistics
using Plots
using Flux
using OceanParameterizations
using Oceananigans.Grids
using BSON
using OrdinaryDiffEq, DiffEqSensitivity
using WindMixing
using JLD2
using FileIO
using OceanTurb
using CairoMakie
include("modified_pacalowski_philander_model.jl")

PATH = pwd()

DATA_PATH = joinpath(PATH, "extracted_training_output", "NDE_training_modified_pacalowski_philander_1sim_-1e-3_2_extracted.jld2")
FILE_PATH = joinpath(pwd(), "Output", "profiles_fluxes_modified_pacalowski_philander_1sim_-1e-3_test")

file = jldopen(DATA_PATH, "r")

losses = file["losses"]

minimum(losses)
size = length(losses)

train_files = file["training_info/train_files"]

Plots.plot(1:1:size, losses, yscale=:log10)
Plots.xlabel!("Iteration")
Plots.ylabel!("Loss mse")

𝒟train = WindMixing.data(train_files, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=false)

test_files = ["-1e-3"]
𝒟test = WindMixing.data(test_files, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=false)
uw_NN = file["neural_network/uw"]
vw_NN = file["neural_network/vw"]
wT_NN = file["neural_network/wT"]

trange = 1:1:1153
plot_data = NDE_profile(uw_NN, vw_NN, wT_NN, test_files[1], 𝒟test, 𝒟train, trange, unscale=true, modified_pacalowski_philander=true)

ds = jldopen(joinpath(pwd(), "Data", "three_layer_constant_fluxes_hr192_Qu1.0e-03_Qb0.0e+00_f1.0e-04_Nh256_Nz128__statistics.jld2"))

Nz = ds["grid/Nz"]
Lz = ds["grid/Lz"]
u₀ = ds["timeseries/u/0"][:]
v₀ = ds["timeseries/v/0"][:]
T₀ = ds["timeseries/T/0"][:]
Fu = ds["boundary_conditions/u_top"]
Fθ = ds["boundary_conditions/θ_top"]
f₀ = ds["parameters/coriolis_parameter"]

ΔRi = 1.
ν₀ = 1e-4
ν₋ = 1e-1
Riᶜ = 0.25
Pr = 1.

constants = OceanTurb.Constants(Float64, f=f₀)
parameter = ModifiedPacanowskiPhilanderParameters(Cν₀=ν₀, Cν₋=ν₋, Riᶜ=Riᶜ, ΔRi=ΔRi)
model = ModifiedPacanowskiPhilanderModel(N=Nz, L=Lz, stepper=:BackwardEuler, constants=constants, parameters=parameter)

model.bcs[1].top = OceanTurb.FluxBoundaryCondition(Fu)
model.bcs[3].top = OceanTurb.FluxBoundaryCondition(Fθ)

model.solution[1].data[1:Nz] .= u₀
model.solution[2].data[1:Nz] .= v₀
model.solution[3].data[1:Nz] .= T₀

Δt = 10.0
times = [ds["timeseries/t/$i"] for i in keys(ds["timeseries/t"])]
Nt = length(times)

U_solution = zeros(Nz, Nt)
V_solution = zeros(Nz, Nt)
T_solution = zeros(Nz, Nt)

U′W′_solution = zeros(Nz+1, Nt)
V′W′_solution = zeros(Nz+1, Nt)
W′T′_solution = zeros(Nz+1, Nt)

Ri_solution = zeros(Nz+1, Nt)

function get_diffusive_flux(field_index, model)
    flux = FaceField(model.grid)
    field = model.solution[field_index]
    K = model.timestepper.eqn.K[field_index]
    for i in interiorindices(flux)
        @inbounds flux[i] = - K(model, i) * ∂z(field, i)
    end
    return flux
end

function get_richardson_number_profile(model)
    Ri = FaceField(model.grid)
    for i in interiorindices(Ri)
        @inbounds Ri[i] = local_richardson(model, i)
    end
    return Ri
end

for n in 1:Nt
    OceanTurb.run_until!(model, Δt, times[n])
    @info "Time = $(times[n])"
    U_solution[:, n] .= model.solution[1][1:Nz]
    V_solution[:, n] .= model.solution[2][1:Nz]
    T_solution[:, n] .= model.solution[3][1:Nz]

    U′W′_solution[:, n] .= get_diffusive_flux(1, model)[1:Nz+1]
    V′W′_solution[:, n] .= get_diffusive_flux(2, model)[1:Nz+1]
    W′T′_solution[:, n] .= get_diffusive_flux(3, model)[1:Nz+1]

    U′W′_solution[Nz+1, n] = Fu

    Ri_solution[:, n] = get_richardson_number_profile(model)[1:Nz+1]
end

Ri_solution = clamp.(Ri_solution, -1, 2)

frame = Node(1)

truth_u = @lift plot_data["truth_u"][:,$frame]
truth_v = @lift plot_data["truth_v"][:,$frame]
truth_T = @lift plot_data["truth_T"][:,$frame]

test_u = @lift plot_data["test_u"][:,$frame]
test_v = @lift plot_data["test_v"][:,$frame]
test_T = @lift plot_data["test_T"][:,$frame]

truth_uw = @lift plot_data["truth_uw"][:,$frame]
truth_vw = @lift plot_data["truth_vw"][:,$frame]
truth_wT = @lift plot_data["truth_wT"][:,$frame]

test_uw = @lift plot_data["test_uw"][:,$frame]
test_vw = @lift plot_data["test_vw"][:,$frame]
test_wT = @lift plot_data["test_wT"][:,$frame]

truth_Ri = @lift clamp.(plot_data["truth_Ri"][:,$frame], -1, 2)
test_Ri = @lift clamp.(plot_data["test_Ri"][:,$frame], -1, 2)

u_mpp = @lift U_solution[:,$frame]
v_mpp = @lift V_solution[:,$frame]
T_mpp = @lift T_solution[:,$frame]

uw_mpp = @lift U′W′_solution[:,$frame]
vw_mpp = @lift V′W′_solution[:,$frame]
wT_mpp = @lift W′T′_solution[:,$frame]

Ri_mpp = @lift Ri_solution[:,$frame]

u_max = maximum([maximum(plot_data["truth_u"]), maximum(plot_data["test_u"]), maximum(U_solution)])
u_min = minimum([minimum(plot_data["truth_u"]), minimum(plot_data["test_u"]), minimum(U_solution)])

v_max = maximum([maximum(plot_data["truth_v"]), maximum(plot_data["test_v"]), maximum(V_solution)])
v_min = minimum([minimum(plot_data["truth_v"]), minimum(plot_data["test_v"]), minimum(V_solution)])

T_max = maximum([maximum(plot_data["truth_T"]), maximum(plot_data["test_T"]), maximum(T_solution)])
T_min = minimum([minimum(plot_data["truth_T"]), minimum(plot_data["test_T"]), minimum(T_solution)])

uw_max = maximum([maximum(plot_data["truth_uw"]), maximum(plot_data["test_uw"]), maximum(U′W′_solution)])
uw_min = minimum([minimum(plot_data["truth_uw"]), minimum(plot_data["test_uw"]), minimum(U′W′_solution)])

vw_max = maximum([maximum(plot_data["truth_vw"]), maximum(plot_data["test_vw"]), maximum(V′W′_solution)])
vw_min = minimum([minimum(plot_data["truth_vw"]), minimum(plot_data["test_vw"]), minimum(V′W′_solution)])

wT_max = maximum([maximum(plot_data["truth_wT"]), maximum(plot_data["test_wT"]), maximum(W′T′_solution)])
wT_min = minimum([minimum(plot_data["truth_wT"]), minimum(plot_data["test_wT"]), minimum(W′T′_solution)])

SIMULATION_NAME="Modified Pacalowski Wind-Mixing, Training Data"
plot_title = @lift "$SIMULATION_NAME: time = $(round(times[$frame]/86400, digits=2)) days"
fig = Figure(resolution=(1920, 1080))

colors=["navyblue", "hotpink2", "forestgreen"]

u_str = "u / m s⁻¹"
v_str = "v / m s⁻¹"
T_str = "T / °C"
uw_str = "uw / m² s⁻²"
vw_str = "vw / m² s⁻²"
wT_str = "wT / m s⁻¹ °C"

zc = plot_data["depth_profile"]
zf = plot_data["depth_flux"]
z_str = "z / m"

zc_mpp = model.grid.zc
zf_mpp = model.grid.zf

ax_u = fig[1, 1] = Axis(fig, xlabel=u_str, ylabel=z_str)
u_lines = [lines!(ax_u, truth_u, zc, linewidth=3, color=colors[1]), lines!(ax_u, test_u, zc, linewidth=3, color=colors[2]), lines!(ax_u, u_mpp, zc_mpp, linewidth=3, color=colors[3])]
CairoMakie.xlims!(ax_u, u_min, u_max)
CairoMakie.ylims!(ax_u, minimum(zc), 0)

ax_v = fig[1, 2] = Axis(fig, xlabel=v_str, ylabel=z_str)
v_lines = [lines!(ax_v, truth_v, zc, linewidth=3, color=colors[1]), lines!(ax_v, test_v, zc, linewidth=3, color=colors[2]), lines!(ax_v, v_mpp, zc_mpp, linewidth=3, color=colors[3])]
CairoMakie.xlims!(ax_v, v_min, v_max)
CairoMakie.ylims!(ax_v, minimum(zc), 0)

ax_T = fig[1, 3] = Axis(fig, xlabel=T_str, ylabel=z_str)
T_lines = [lines!(ax_T, truth_T, zc, linewidth=3, color=colors[1]), lines!(ax_T, test_T, zc, linewidth=3, color=colors[2]), lines!(ax_T, T_mpp, zc_mpp, linewidth=3, color=colors[3])]
CairoMakie.xlims!(ax_T, T_min, T_max)
CairoMakie.ylims!(ax_T, minimum(zc), 0)

ax_uw = fig[2, 1] = Axis(fig, xlabel=uw_str, ylabel=z_str)
uw_lines = [lines!(ax_uw, truth_uw, zf, linewidth=3, color=colors[1]), lines!(ax_uw, test_uw, zf, linewidth=3, color=colors[2]), lines!(ax_uw, uw_mpp, zf_mpp, linewidth=3, color=colors[3])]
CairoMakie.xlims!(ax_uw, uw_min, uw_max)
CairoMakie.ylims!(ax_uw, minimum(zf), 0)

ax_vw = fig[2, 2] = Axis(fig, xlabel=vw_str, ylabel=z_str)
vw_lines = [lines!(ax_vw, truth_vw, zf, linewidth=3, color=colors[1]), lines!(ax_vw, test_vw, zf, linewidth=3, color=colors[2]), lines!(ax_vw, vw_mpp, zf_mpp, linewidth=3, color=colors[3])]
CairoMakie.xlims!(ax_vw, vw_min, vw_max)
CairoMakie.ylims!(ax_vw, minimum(zf), 0)

ax_wT = fig[2, 3] = Axis(fig, xlabel=wT_str, ylabel=z_str)
wT_lines = [lines!(ax_wT, truth_wT, zf, linewidth=3, color=colors[1]), lines!(ax_wT, test_wT, zf, linewidth=3, color=colors[2]), lines!(ax_wT, wT_mpp, zf_mpp, linewidth=3, color=colors[3])]
CairoMakie.xlims!(ax_wT, wT_min, wT_max)
CairoMakie.ylims!(ax_wT, minimum(zf), 0)

ax_Ri = fig[2, 4] = Axis(fig, xlabel="Ri", ylabel=z_str)
Ri_lines = [lines!(ax_Ri, truth_Ri, zf, linewidth=3, color=colors[1]), lines!(ax_Ri, test_Ri, zf, linewidth=3, color=colors[2]), lines!(ax_Ri, Ri_mpp, zf_mpp, linewidth=3, color=colors[3])]
CairoMakie.xlims!(ax_Ri, -1, 2)
CairoMakie.ylims!(ax_Ri, minimum(zf), 0)

legend = fig[1, 4] = Legend(fig, u_lines, ["Oceananigans.jl LES", "NDE Prediction", "OceanTurb.jl Modified Pacalowski-Philander"])
supertitle = fig[0, :] = Label(fig, plot_title, textsize=30)
trim!(fig.layout)

record(fig, "$FILE_PATH.mp4", 1:length(times), framerate=30) do n
    @info "Animating mp4 frame $n/$(length(times))..."
    frame[] = n
end


text:
using ChainRulesTestUtils
using Documenter

makedocs(
    modules=[ChainRulesTestUtils],
    format=Documenter.HTML(prettyurls=false, assets=["assets/chainrules.css"]),
    sitename="ChainRulesTestUtils",
    authors="JuliaDiff contributors",
    strict=true,
    checkdocs=:exports,
)

const repo = "github.com/JuliaDiff/ChainRulesTestUtils.jl.git"
const PR = get(ENV, "TRAVIS_PULL_REQUEST", "false")
if PR == "false"
    # Normal case, only deploy docs if merging to master or release tagged
    deploydocs(repo=repo)
else
    @info "Deploying review docs for PR #$PR"
    # TODO: remove most of this once https://github.com/JuliaDocs/Documenter.jl/issues/1131 is resolved

    # Overwrite Documenter's function for generating the versions.js file
    foreach(Base.delete_method, methods(Documenter.Writers.HTMLWriter.generate_version_file))
    Documenter.Writers.HTMLWriter.generate_version_file(_, _) = nothing
    # Overwrite necessary environment variables to trick Documenter to deploy
    ENV["TRAVIS_PULL_REQUEST"] = "false"
    ENV["TRAVIS_BRANCH"] = "master"

    deploydocs(
        devurl="preview-PR$(PR)",
        repo=repo,
    )
end


text:
# various utility functions

"""
Return a vector of action=>Q pairs ranked by Q.

For debugging/ checking purposes - this may be slow
"""
function ranked_actions(planner::AbstractMCTSPlanner, state)
    @assert planner.tree != nothing "in ranked_actions, planner did not have a tree; use keep_tree=true when constructing the solver"
    ranked_actions(planner.tree, state)
end

function ranked_actions(tree::MCTSTree, state)
    sanodes = children(StateNode(tree, tree.state_map[state]))
    sorted_sanodes = sort(collect(sanodes), by=q, rev=true)
    return [action(n)=>q(n) for n in sorted_sanodes]
end

function ranked_actions(tree::DPWTree, state)
    snode = tree.s_lookup[state]
    sanodes = tree.children[snode]
    tuples = [(tree.a_labels[n], tree.q[n]) for n in sanodes]
    return sort!(tuples, by=last, rev=true)
end


text:
const JULIA_LOGO = begin
    R = RGB565(0.796,0.235,0.2)
    r = RGB565(0.835,0.388,0.361)
    G = RGB565(0.22,0.596,0.149)
    g = RGB565(0.376,0.678,0.318)
    P = RGB565(0.584,0.345,0.698)
    p = RGB565(0.667,0.475,0.757)
    z = RGB565(0.0,0.0,0.0)

    permutedims([z z z z z z z z;
                 z z z G G z z z;
                 z z G g g G z z;
                 z z G g g G z z;
                 z R R G G P P z;
                 R r r R P p p P;
                 R r r R P p p P;
                 z R R z z P P z], (2,1)) # transpose to correct orientation
end


text:
using LinearAlgebra

A = [3 1 2;
     1 3 1;
     1 0 1]

B = [1 0 2;
     1 3 6;
     2 2 0]

C = [1 0 0;
     0 1 1;
     1 1 2]

D = [0 0 2;
     1 0 1;
     2 2 1]

println("(LHS)=", inv(A+B*C*D))
println("(RHS))=", inv(A)-inv(A)*B*inv(inv(C)+D*inv(A)*B)*D*inv(A))


text:
using IICGP
using Dates

rootdir = joinpath(homedir(), "Documents/git/ICGP-results/")
# rootdir = dirname(@__DIR__)
resdir = joinpath(rootdir, "results/")
fetch_backup(rootdir, clean=true)

# Filters parameters
min_date = DateTime(2021, 09, 01)
max_date = DateTime(2021, 10, 02)
savedir_index = 1
labels = Vector{String}()
colors = Vector{Symbol}()
games_12 = ["boxing" "assault" "freeway" "solaris" "defender" "gravitar" "space_invaders" "private_eye" "asteroids" "breakout" "frostbite" "riverraid"]
games = games_12# ["gravitar"]
reducers = ["pooling"]# Vector{String}()
dotime = false
dosave = true

baselines = true

#=
labels = ["1 + λ (\"deterministic\" game)", "GA (stochastic games)"] # Vector{String}()
colors = [:skyblue3, :limegreen] # Vector{Symbol}()
# https://juliagraphics.github.io/Colors.jl/stable/namedcolors/
savedir_index = 2
=#

for g in games
    exp_dirs, games = get_exp_dir(
        resdir, min_date=min_date, max_date=max_date, games=[g],
        reducers=reducers
    )
    process_results(
        exp_dirs, games, dotime, dosave, ma=1, baselines=baselines,
        labels=labels, colors=colors, savedir_index=savedir_index
    )
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 2940cb92-4d76-11eb-3628-493a998f3033
using BenchmarkTools

# ╔═╡ bdb2792e-4d71-11eb-1c5b-03e16bad7067
md"""
# Analysis of Algorithms
Theoretical Study of Program Performance and Resources Usage

**What is more important that Performance ?**
"""

# ╔═╡ 0a158d10-4d72-11eb-3bf5-0ff33f2b9b9d
begin
	hint(text) = Markdown.MD(Markdown.Admonition("hint", "Hint", [text]))
	hint(md"""
	- Correctness
	- Simplicity
	- Maintainability: Modularity,
	- Rapid Prototyping (Programming Time)
	- Robustness and Stability of the software
	- Wider Collection of features and Functionality
	- Security
	- User Friendliness
	- Scalability => related to performance
	""")
end

# ╔═╡ dc8d89be-4d72-11eb-13e8-718f21bc1362
md"""
**Why Study Analysis of Algorithms and Performance though its in bottom of the above list ?**
- Performance related to User Friendliness (too long lag makes it less user friendly), also related to Security
- Performance draws the line between feasible and not feasible. (If its not fast enough or uses too much memory it is not useable) => bring you to cutting edge of entrepreneurship
- Algorithms give language to talk about program behaviour and Quantify the resource measure of the algorithm
- Efficiency is a concern for large scale problems (scalability): Solving Large Problems with Large Inputs
- Cases where the solution might not be really clear on how to do (or even to state clearly) such as Ml, CV, NLP, Playing Games

Examples
- Shortage Path between two points given map of all highways
- Best Assignment of Resources to People or Jobs to machines, Scheduling
- Measuring Similarity of two Documents
- Human Genome
- Social Network in Facebook
"""

# ╔═╡ f3972612-4d73-11eb-3d7a-ff39ee463047
md"""
#### Two functions might produce the same output but have different efficiency?
"""

# ╔═╡ e8a90d58-4d73-11eb-3d34-b93907fef7bc
begin
	function f1(n)
	    for i in 1:2
	        n += 100
		end
	    return n
	end

	function f2(n)
	    for i in 1:100
	        n += 2
		end
	    return n
	end

	f1(10),f2(10)
end

# ╔═╡ 9a9926ca-4d73-11eb-35fb-8529f7e050d4
md"""
## How should we measure the running time of the algorithms?

Experimental Approach
1. Write a program that implements the program
2. Run the program with the datasets of varying size and composition
3. Measure system time to get an accurate measure of the actual running time.

we can use `@time` macro to measure runtime.

System time might not be accurate enough.
- There are multiple process running in the background which leads to random noise in runtime measurements => Take Measurement Multiple Times.
- Compilation time and Runtime might need to be separated.

Popular tool for Julia [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl) creates a `@benchmarkable` object and uses `tune!` to figure out the params (such as number of samples to evaluate) and even set params by oneself before executing the benchmark trials and collecting samples. One can also run GC before starting a benchmark with `gctrial=true` and before taking samples with `gcsample=true`. Also set time to run with `seconds=300` to prolong the  benchmark time.
"""

# ╔═╡ 9a1431b8-4d73-11eb-3381-8bb9213a52cd
b = @benchmark sum(sin.(cos.(rand(100_000))))

# ╔═╡ 99e369c0-4d73-11eb-1d56-53a3da9bedda
begin
	using Plots
	histogram(b.times)
end

# ╔═╡ f201e03c-4d73-11eb-26d4-3150337daec6
md"""
#### Limitations of experimental study
- need to implement the algorithm before we can test the running time.
- we can only perform experiment on a limited set of inputs, this may not be indicative of running time on different set of data
- for comparison of two algorithms, the same hardware and software environment as well as operating environment needs to be maintained.

Benchmarking on different systems will result in different speeds need a more generalized method to analyze and communicate the performance of the algorithm.
- use high-level description of algorithm (pseudocodes) instead of testing implemented program
- take into account all possible inputs
- evaluate efficiency independent of hardware and software environment as well as tiny implementation details.

"""

# ╔═╡ c32d52c6-4d7a-11eb-3a20-019cdb3c889b
md"""
Pseudocode Convention for low-level Primitive operations.
- Data assignment: use = for assignment
- Arthimetic and logical operations: == for equality, (=,-,*,/) use standard mathematical symbols
- Programming Constructs
    - Conditionals: if ... then ... [else if ... ] [else ...]
    - Loops
        - While ... do ...
        - repeat ... until ...
        - for ... do ...
    - array indexing: A[i], A[i,j]
    - Methods
        - Declaration: **Algorithm** methodname(arg1, arg2)
        - calls: methodname(args)
        - returns: return value

Julia Code is very similar to Pseudocode
"""

# ╔═╡ eb203a78-4d7a-11eb-2e5a-13c68c97a58c
md"""
## Asymptotic Analysis: How does runtime scales with input size ?
Ignore Machine dependent constants, look at the growth of the running $T(n)$ as $n -> \infty$ rather than actually running time.

Three Types of Analysis
- Worst-case (usually): $T(n)$ = max time on any input of size n
- Average-case (sometimes): $T(n)$ = expected time over all inputs of size n. (Expectation needs assumption of statistical distribution of time)
- Best-case (Not useful: eg sorting a sorted array)

"""

# ╔═╡ a3d6e466-4d80-11eb-1d4e-933988883a5b
md"""
## Notations
Bounded From Above

> Big-O : $f(n) = O(g(n)) \text{ or } f \preceq g, \quad \exists N,c \ : \ \forall n \geq N, f(n) \leq c \cdot g(n)$
- eg: $3n^2+5n+2 \leq 3n^2+5n^2+2n^2$ both are $O(n^2)$ as they have the same growth rate. Smaller items are general omited
- eg: $O(\log(n))$ is same for $\log_2(n)$, $\log_{10}(n)$, $\log_x(n)$  base doesn't matter as they only differ by a constant multiple
- since $n^a \prec n^b \quad \forall 0 <a<b : n=O(n^2),\  \sqrt(n) = O(n)$ O(n^2) grows faster
- also $n^a \prec b^n \quad \forall a> 0 b>0 : n^5=O((\sqrt2)^n),\  n^{100} = O(1.1^n)$
- also $(\log n)^a \prec n^b \quad \forall (a,b>0): (\log n)^3 = O(\sqrt(n)), n \log n = O(n^2)$

"""

# ╔═╡ f273f432-4d8d-11eb-0d4b-e330161da74b
md"""

log $(@bind plotlog html"<input type=checkbox >")
sqrt $(@bind plotsqrt html"<input type=checkbox >")
n $(@bind plotn html"<input type=checkbox >")
n logn $(@bind plotnlogn html"<input type=checkbox >")
n² $(@bind plotnsq html"<input type=checkbox >")
2ⁿ $(@bind plot2n html"<input type=checkbox >")

`plotrange = ` $(@bind plotrange html"<input type=range min=2 max=50 step=1 value=10>")
"""

# ╔═╡ 1a95c1d8-4d7b-11eb-2c23-abfd4e673b3c
begin
	#using LaTeXStrings
	plot([n->log(n), n->sqrt(n),n->n, n->n*log(n),n->n^2,n->2^n], 1:0.1:plotrange,label=(["log(n)" "sqrt(n)" "n" "n*log(n)" "n^2" "2^n"]))
end

# ╔═╡ 82abcbc4-4d84-11eb-05e7-b7fdd509052b
md"""
$$\log{n} \prec \sqrt{n} \prec n \prec n \log{n} \prec n^2 \prec 2^n$$
"""

# ╔═╡ 399e9e2a-4d8e-11eb-315a-319497f05e8b
md"""
For function f,g : $\mathbb N \rightarrow \mathbb R^+$ we say that
- Big-Ω: $f(n) = \Omega(g(n)) \text{ or } f \succeq g \text{ if for some }c, f(n) \geq c \cdot g(n)$ f grows slower than g
- Big-Θ:$f(n) = \Theta(g(n)) \text{ or } f \asymp g \text{ if }, f = O(g) \text{ and } f= \Omega(g)$ f grows at the same rate as g
- little-o : $ f(n) = o(g(n)) \text{ or } f \prec g \text{ if } f(n)/g(n) \rightarrow 0 \text{ as } n \rightarrow \infty $ f grows slower than g
"""

# ╔═╡ 0522e1fe-4d90-11eb-32fe-dd3a37e7d498
md"""
Big-O looks only asymptotic detail (useful for very large data) and loses important information about constant multiples. But change by a factor of 100 or 2 is still a big deal. Once we have a good asymptotic run time we can then look into other details such as benchmarking, interaction with memory hierarchy,...
"""

# ╔═╡ 256fd75e-4d73-11eb-0ba3-cd5276a8640c
md"""
## References
- [MIT 6.006 Introduction to Algorithms](https://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-006-introduction-to-algorithms-fall-2011/)
- [MIT 6.046 Design and Analysis of Algorithms](https://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-046j-design-and-analysis-of-algorithms-spring-2015/)
- [UCSandiegoX DSA](https://www.edx.org/micromasters/ucsandiegox-algorithms-and-data-structures)
- [GTX DSA](https://www.edx.org/professional-certificate/gtx-data-structures-and-algorithms)
- [Algorithms](https://en.wikipedia.org/wiki/Algorithm)

"""

# ╔═╡ Cell order:
# ╟─bdb2792e-4d71-11eb-1c5b-03e16bad7067
# ╟─0a158d10-4d72-11eb-3bf5-0ff33f2b9b9d
# ╟─dc8d89be-4d72-11eb-13e8-718f21bc1362
# ╟─f3972612-4d73-11eb-3d7a-ff39ee463047
# ╠═e8a90d58-4d73-11eb-3d34-b93907fef7bc
# ╟─9a9926ca-4d73-11eb-35fb-8529f7e050d4
# ╠═2940cb92-4d76-11eb-3628-493a998f3033
# ╠═9a1431b8-4d73-11eb-3381-8bb9213a52cd
# ╠═99e369c0-4d73-11eb-1d56-53a3da9bedda
# ╟─f201e03c-4d73-11eb-26d4-3150337daec6
# ╟─c32d52c6-4d7a-11eb-3a20-019cdb3c889b
# ╟─eb203a78-4d7a-11eb-2e5a-13c68c97a58c
# ╠═a3d6e466-4d80-11eb-1d4e-933988883a5b
# ╟─f273f432-4d8d-11eb-0d4b-e330161da74b
# ╟─1a95c1d8-4d7b-11eb-2c23-abfd4e673b3c
# ╟─82abcbc4-4d84-11eb-05e7-b7fdd509052b
# ╠═399e9e2a-4d8e-11eb-315a-319497f05e8b
# ╟─0522e1fe-4d90-11eb-32fe-dd3a37e7d498
# ╟─256fd75e-4d73-11eb-0ba3-cd5276a8640c


text:
function utility (cons)

#This function takes consumption as an argument and returns utility. The
#utility functin is CRRA except that we add a very small number (eps) to
#consumption so that the computer can deal wi

#global gamma

if cons<=0
   error("Error in utility. Consumption is <=0");
end

if gamma == 1
    utils = log(cons);
else
    utils = ((cons)^(1-gamma)  )/(1-gamma);
end

return utils

end

max_stars_repo_name: floswald/ucl-econ-julia

max_issues_repo_name: floswald/ucl-econ-julia

max_forks_repo_name: floswald/ucl-econ-julia


text:
#function testSWGalewskiSphere

using CGDycore

# Physical parameters
Param=CGDycore.PhysParameters();

# Grid
nz=1;
Param.nPanel=32;
Param.H=2;
Param.Grid=CGDycore.CubedGrid(Param.nPanel,CGDycore.OrientFaceSphere,Param);

Param.Grid.nz=nz;
Param.Grid.zP=zeros(nz,1);
Param.Grid.z=zeros(nz+1,1);
Param.Grid.dz=Param.H/nz;
Param.Grid.zP[1]=Param.Grid.dz/2;
for i=2:nz
  Param.Grid.zP[i]=Param.Grid.zP[i-1]+Param.Grid.dz;
end
for i=2:nz+1
  Param.Grid.z[i]=Param.Grid.z[i-1]+Param.Grid.dz;
end

# Model
Param.ModelType="Curl";
Param.Deep=false;
Param.HeightLimit=30000.0;
Param.Damping=false;
Param.Coriolis=true;
Param.CoriolisType="Sphere";
Param.Buoyancy=false;
Param.Source=false;
Param.Equation="Shallow";
Param.TopoS="";
Param.ProfVel="Galewsky";
Param.ProfRho="Galewsky";
Param.ProfTheta="Galewsky";
Param.NumV=5;
Param.RhoPos=1;
Param.uPos=2;
Param.vPos=3;
Param.wPos=4;
Param.ThPos=5;
Param.Thermo="";#"Energy"
Param.alphaG=1.0/3.0;
Param.betaG=1.0/15.0;
Param.hH=120.0;
Param.H0G=10000.0;
Param.uM=80.0;
Param.lat0G=pi/7.0;
Param.lat1G=pi/2.0-Param.lat0G;
Param.eN=exp(-4.0/(Param.lat1G-Param.lat0G)^2.0);
Param.RefProfile=false

# Discretization
OrdPoly=4;
OrdPolyZ=1;
(CG,Param)=CGDycore.Discretization(OrdPoly,OrdPolyZ,CGDycore.JacobiSphere3,Param);
LRef=11*1.e5;
dx=2*pi*Param.RadEarth/4/Param.nPanel/OrdPoly;
Param.HyperVisc=true;
Param.HyperDCurl=2.e14;
Param.HyperDGrad=2.e14;
Param.HyperDDiv=0; # Scalars
Param.Upwind=false;

# Output
Param.Flat=true;
Param.vtk=0;
Param.vtkFileName="Galewsky";
Param.RadPrint=Param.H;
vtkGrid=CGDycore.vtkCGGrid(CG,CGDycore.TransSphere,CGDycore.Topo,Param);
Param.cNames = [
  "Rho",
  "u",
  "v",
  "Th",
  "Vort"
]

# Initial conditions
Param.NumV=5;
U=zeros(CG.NumG,nz,Param.NumV);
U[:,:,Param.RhoPos]=CGDycore.Project(CGDycore.fRho,CG,Param);
(U[:,:,Param.uPos],U[:,:,Param.vPos])=CGDycore.ProjectVec(CGDycore.fVel,CG,Param);
U[:,:,Param.ThPos]=CGDycore.Project(CGDycore.fTheta,CG,Param).*U[:,:,Param.RhoPos];

# Integration
time=0;
IntMethod="RungeKutta";
dtau=100;
Param.RK=CGDycore.RungeKuttaMethod("RK3");
SimDays=6;
PrintDay=.5;
nIter=24*3600*SimDays/dtau;
PrintInt=24*3600*PrintDay/dtau;
# Print initial conditions
Param.vtk=CGDycore.vtkOutput(U,vtkGrid,CG,Param);
#
str = IntMethod
if str == "RungeKutta"
    @time begin
      for i=1:nIter
        @info "Iteration: $i"
        U .= CGDycore.RungeKuttaExplicit(U,dtau,CGDycore.FcnNHCurlVec,CG,Param);
        if mod(i,PrintInt)==0
          Param.vtk=CGDycore.vtkOutput(U,vtkGrid,CG,Param);
        end
      end
end
end


text:
function varnames(graph)
    result = Set{VarName}()

    for stmt in values(graph)
        if stmt isa AutoGibbs.Tilde
            vn = AutoGibbs.getvn(stmt)
            if !isnothing(vn)
                push!(result, vn)
            end
        end
    end

    return result
end

track_with_context(model, ctx) = trackdependencies(model, VarInfo(), SampleFromPrior(), ctx)

macro testdependencies(model, varname_exprs...)
    result = quote
        let m = $(esc(model)),
            vns = tuple($(DynamicPPL.varname.(varname_exprs)...))

            graph_default = trackdependencies(m)
            @test varnames(graph_default) ⊇ Set{VarName}(vns)
            # @test check_dependencies(graph_default)

            # graph_likelihood = track_with_context(m, LikelihoodContext())
            # $(test_dependency_invariants(:graph_likelihood))
        end
    end

    Base.remove_linenums!(result)
    return result
end

macro test_nothrow(ex)
    orig_ex = Expr(:inert, ex)
    truthy_ex = quote
        $(esc(ex))
        true
    end

    result = quote
        try
            Test.Returned($truthy_ex, nothing, $(QuoteNode(__source__)))
        catch _e
            _e isa Test.InterruptException && rethrow()
            Test.Returned(false, nothing, $(QuoteNode(__source__)))
        end
    end

    Base.remove_linenums!(result)
    return :(Test.do_test($result, $orig_ex))
end

function issimilar(d1::DiscreteNonParametric, d2::DiscreteNonParametric; atol::Real=0)
    return (support(d1) == support(d2)) &&
        all(isapprox.(probs(d1), probs(d2); atol=atol))
end

function issimilar(p1::Product, p2::Product; atol::Real=0)
    return all(issimilar(d1, d2) for (d1, d2) in zip(p1.v, p2.v))
end


text:
using UnPack, OrderedCollections, LinearAlgebra, JLD2, SparseArrays

# Load guesses
sssout = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "..", "test", "reference", "crw_sss.jld2"), "r")

zguess = vec(sssout["z_rss"])
yguess = vec(sssout["y_rss"])
Psiguess = sssout["Psi_rss"]

mutable struct CoeurdacierReyWinant{T <: Real}
    σr::T # SD of interest rate shocks
    σy::T # SD of endowment shocks
    β::T  # intertemporal discount rate
    γ::T  # risk aversion coefficient
    θ::T
    ρr::T # persistence of interest rate
    ρy::T # persistence of endowment
    rr::T # long-run interest rate
    yy::T # long-run endowment
end

function CoeurdacierReyWinant(; σr::T = .025, σy::T = .025, β::T = .96, γ::T = 2.,
                              θ::T = 1., ρr::T = .9, ρy::T = .9, rr::T = .01996, yy::T = log(θ)) where {T <: Real}
    return CoeurdacierReyWinant{T}(σr, σy, β, γ, θ, ρr, ρy, rr, yy)
end

function crw(m::CoeurdacierReyWinant{T}; Ψ = nothing, sparse_jacobian::Vector{Symbol} = Symbol[],
             sparse_arrays::Bool = false) where {T <: Real}
    @unpack σr, σy, β, γ, θ, ρr, ρy, rr, yy = m

    # Nₜ = exp(rₜ) * Aₜ₋₁ + Yₜ, where Aₜ is foreign assets and Yₜ is the endowment
    # The jump variables are consumption, expected return on assets Xₜ = 𝔼ₜ[Rₜ₊₁], and
    # Wₜ = 𝔼ₜ[Yₜ₊₁]
    S  = OrderedDict{Symbol, Int}(:N => 1, :r => 2, :y => 3) # State variables
    J  = OrderedDict{Symbol, Int}(:c => 1, :x => 2, :w => 3) # Jump variables
    SH = OrderedDict{Symbol, Int}(:εr => 1, :εy => 2)        # Exogenous shocks
    Nz = length(S)
    Ny = length(J)
    Nε = length(SH)

    function μ(F, z, y) # note that y here refers to jump variables
        F[S[:N]] = exp(y[J[:w]]) + exp(y[J[:x]]) * (z[S[:N]] - exp(y[J[:c]]))
        F[S[:r]] = (1 - ρr) * rr + ρr * z[S[:r]]
        F[S[:y]] = (1 - ρy) * yy + ρy * z[S[:y]]
    end

    function ξ(F, z, y)
        F[J[:c]] = log(β) + γ * y[J[:c]]    # Euler equation
        F[J[:x]] = -y[J[:x]]                # rₜ₊₁ - xₜ, rational expectations
        F[J[:w]] = -y[J[:w]]                # yₜ₊₁ - wₜ
    end

    # The cache is initialized as zeros so we only need to fill non-zero elements
    function Λ(F, z, y)
        F[S[:N], J[:x]] = z[S[:N]] - exp(y[J[:c]])
        F[S[:N], J[:w]] = 1.
    end

    # The cache is initialized as zeros so we only need to fill non-zero elements
    function Σ(F, z, y)
        F[S[:r], SH[:εr]] = σr
        F[S[:y], SH[:εy]] = σy
    end

    Γ₅ = zeros(T, Ny, Nz)
    Γ₅[J[:c], S[:r]] = 1.
    Γ₅[J[:x], S[:r]] = 1.
    Γ₅[J[:w], S[:y]] = 1.

    Γ₆ = zeros(T, Ny, Ny)
    Γ₆[J[:c], J[:c]] = -γ

    if sparse_arrays
        Γ₅ = sparse(Γ₅)
        Γ₆ = sparse(Γ₆)
    end

    z = zguess
    y = yguess
    if isnothing(Ψ)
        Ψ = Psiguess
    end

    if sparse_arrays
        return RiskAdjustedLinearization(μ, Λ, Σ, ξ, Γ₅, Γ₆, crw_ccgf, z, y, Ψ, Nε; sparse_jacobian = sparse_jacobian,
                                         Λ_cache_init = dims -> spzeros(dims...),
                                         Σ_cache_init = dims -> spzeros(dims...),
                                         jump_dependent_shock_matrices = true)
    else
        return RiskAdjustedLinearization(μ, Λ, Σ, ξ, Γ₅, Γ₆, crw_ccgf, z, y, Ψ, Nε; sparse_jacobian = sparse_jacobian,
                                         jump_dependent_shock_matrices = true)
    end
end

crw_cₜ(m, zₜ) = exp(m.y[1] + (m.Ψ * (zₜ - m.z))[1])

# Evaluates m_{t + 1} + r_{t + 1}
function crw_logSDFxR(m, zₜ, εₜ₊₁, Cₜ)
    zₜ₊₁, yₜ₊₁ = simulate(m, εₜ₊₁, zₜ)

    return log(m_crw.β) - m_crw.γ * (yₜ₊₁[1] - log(Cₜ)) + zₜ₊₁[2]
end

# Calculate 𝔼ₜ[exp(mₜ₊₁ + rₜ₊₁)] via quadrature
std_norm_mean = zeros(2)
std_norm_sig  = ones(2)
crw_𝔼_quadrature(f::Function) = gausshermite_expectation(f, std_norm_mean, std_norm_sig, 10)

# Calculate implied state variable(s)
function crw_endo_states(m, zₜ, zₜ₋₁, c_impl)
    # rₜ, yₜ are exogenous while Nₜ = exp(rₜ) * Aₜ₋₁ + Yₜ is entirely pre-determined.
    # Thus, our implied state variable will be foreign asset Aₜ = Nₜ - Cₜ.

    # zₜ₋₁ may be the previous period's implied state, so we start from there
    # to calculate Aₜ₋₁.
    yₜ₋₁ = m.y + m.Ψ * (zₜ₋₁ - m.z) # Calculate implied jump variables last period
    Cₜ₋₁ = exp(yₜ₋₁[1])             # to get the implied consumption last period.
    Aₜ₋₁ = zₜ₋₁[1] - Cₜ₋₁           # Given that consumption, we compute implied foreign assets yesterday.
    Nₜ   = exp(zₜ[2]) * Aₜ₋₁ + exp(zₜ[3]) # Now we can get implied resources available today.

    return vcat(zₜ, Nₜ - exp(c_impl)) # This gives us implied foreign assets today, along with other state variables
end

function crw_ccgf(F, α, z)
    # F .= .5 * diag(α * α') # slower but this is the underlying math
    sum!(F, α.^2) # faster implementation
    F .*= .5
end


text:
# Autogenerated wrapper script for wolfSSL_jll for i686-linux-gnu
export libwolfssl

JLLWrappers.@generate_wrapper_header("wolfSSL")
JLLWrappers.@declare_library_product(libwolfssl, "libwolfssl.so.24")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libwolfssl,
        "lib/libwolfssl.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using BenchmarkTools, Combinatorics

function missingperm(arr::Vector)
    allperms = String.(permutations(arr[1]))  # revised for type safety
    for perm in allperms
        if perm ∉ arr return perm end
    end
end

arr = ["ABCD", "CABD", "ACDB", "DACB", "BCDA", "ACBD", "ADCB", "CDAB", "DABC", "BCAD",
       "CADB", "CDBA", "CBAD", "ABDC", "ADBC", "BDCA", "DCBA", "BACD", "BADC", "BDAC",
       "CBDA", "DBCA", "DCAB"]
@show missingperm(arr)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode