text:
using QuantLib

## TYPES ##
type ReplicationError
  maturity::Float64
  payoff::PlainVanillaPayoff
  s0::Float64
  sigma::Float64
  r::Float64
  vega::Float64
end

function ReplicationError(optionType::OptionType, maturity::Float64, strike::Float64, s0::Float64, sigma::Float64, r::Float64)
  rDiscount = exp(-r * maturity)
  qDiscount = 1.0
  forward = s0 * rDiscount / rDiscount
  stdDev = sqrt(sigma * sigma * maturity)
  payoff = PlainVanillaPayoff(optionType, strike)

  black = BlackCalculator(payoff, forward, stdDev, rDiscount)

  vega_ = vega(black, maturity)

  return ReplicationError(maturity, payoff, s0, sigma, r, vega_)
end

type ReplicationPathPricer{OT <: OptionType} <: QuantLib.AbstractPathPricer
  optionType::OT
  strike::Float64
  r::Float64
  maturity::Float64
  sigma::Float64
end

## Type Methods ##

function call(rpp::ReplicationPathPricer, path::Path)
  n = length(path)
  # discrete hedging interval
  dt = rpp.maturity / n
  # we assume the stock pays no dividend
  stockDividendYield = 0.0
  # starting...
  t = 0
  stock = path[1] # stock at t=0
  money_account = 0.0 # money_account at t=0

  ### Initial Deal ###
  # Option fair price (black-scholes) at t=0
  rDiscount = exp(-rpp.r * rpp.maturity)
  qDiscount = exp(-stockDividendYield * rpp.maturity)
  forward = stock * qDiscount / rDiscount
  stdDev = sqrt(rpp.sigma * rpp.sigma * rpp.maturity)

  payoff = PlainVanillaPayoff(rpp.optionType, rpp.strike)
  black = BlackCalculator(payoff, forward, stdDev, rDiscount)
  # sell the option, cash in its premium
  money_account += value(black)
  # compute delta
  delta_ = delta(black, stock)
  # delta-hedge the option buying stock
  stockAmount = delta_
  money_account -= stockAmount * stock

  ## Heding during option life ##
  for step = 1:n-1
    # time flows
    t += dt
    # accruing on the money account
    money_account *= exp(rpp.r * dt)

    # stock growth
    stock = path[step + 1]

    # recalculate the option value at the current stock value, and the current time to maturity
    rDiscount = exp(-rpp.r * (rpp.maturity - t))
    qDiscount = exp(-stockDividendYield * (rpp.maturity - t))
    forward = stock * qDiscount / rDiscount
    stdDev = sqrt(rpp.sigma * rpp.sigma * (rpp.maturity - t))
    black = BlackCalculator(payoff, forward, stdDev, rDiscount)
    # recalculate delta
    delta_ = delta(black, stock)
    # re-hedging
    money_account -= (delta_ - stockAmount) * stock
    stockAmount = delta_
  end

  ### Option expiration ###
  # last accrual on my money account
  money_account *= exp(rpp.r * dt)
  # last stock growth
  stock = path[n]
  # the hedger delivers the option payoff to the option holder
  optionPayoff = PlainVanillaPayoff(rpp.optionType, rpp.strike)(stock)
  money_account -= optionPayoff
  # and unwinds the hedge selling his stock position
  money_account += stockAmount * stock
  # final PnL
  return money_account
end

function compute(rp::ReplicationError, nTimeSteps::Int, nSamples::Int)
  cal = QuantLib.Time.TargetCalendar()
  tdy = Dates.today()

  dc = QuantLib.Time.Actual365()

  stateVar = Quote(rp.s0)
  riskFreeRate = FlatForwardTermStructure(tdy, rp.r, dc)
  dividendYield = FlatForwardTermStructure(tdy, 0.0, dc)
  vol = BlackConstantVol(tdy, cal, rp.sigma, dc)
  diffuse = BlackScholesMertonProcess(stateVar, riskFreeRate, dividendYield, vol)

  # Black scholes equation rules the path generator
  # At each step of the log of the stock will ahve drift and sigma^2 variance
  rsg = QuantLib.Math.InverseCumulativeRSG(0, nTimeSteps)
  brownianBridge = false
  myPathGenerator = PathGenerator(diffuse, rp.maturity, nTimeSteps, rsg, brownianBridge)

  myPathPricer = ReplicationPathPricer(rp.payoff.optionType, rp.payoff.strike, rp.r, rp.maturity, rp.sigma)
  statsAccumulator = QuantLib.Math.gen_RiskStatistics()

  mcSim = MonteCarloModel(myPathGenerator, myPathPricer, statsAccumulator, false)

  add_samples!(mcSim, nSamples)

  PLMean = QuantLib.Math.stats_mean(mcSim.sampleAccumulator)
  PLStdDev = QuantLib.Math.stats_std_deviation(mcSim.sampleAccumulator)
  PLSkew = QuantLib.Math.stats_skewness(mcSim.sampleAccumulator)
  PLKurtosis = QuantLib.Math.stats_kurtosis(mcSim.sampleAccumulator)

  println("P&L Mean: ", PLMean)
  println("P&L Standard Deviation: ", PLStdDev)
  println("P&L Skew: ", PLSkew)
  println("P&L Kurtosis: ", PLKurtosis)
end

function main()
  maturity = 1.0 / 12.0 # 1 month
  strike = 100.0
  underlying = 100.0
  vol = 0.20 # 20%
  riskFreeRate = 0.05 # 5%
  rp = ReplicationError(Call(), maturity, strike, underlying, vol, riskFreeRate)

  scenarios = 50000
  hedgesNum = 21
  compute(rp, hedgesNum, scenarios)
  hedgesNum = 84
  compute(rp, hedgesNum, scenarios)
end


text:
#######################
# recursive fibonacci #
#######################

perf_micro_fib(n) = n < 2 ? n : perf_micro_fib(n-1) + perf_micro_fib(n-2)

############
# parseint #
############

function perf_micro_parseint(t)
    local n, m
    for i=1:t
        n = rand(UInt32)
        s = string(n, base=16)
        m = UInt32(parse(Int64, s, base=16))
    end
    @assert m == n
    return n
end

##############
# mandelbrot #
##############

function mandel(z)
    c = z
    maxiter = 80
    for n = 1:maxiter
        if abs(z) > 2
            return n-1
        end
        z = z^2 + c
    end
    return maxiter
end

perf_micro_mandel() = [mandel(complex(r,i)) for i=-1.:.1:1., r=-2.0:.1:0.5]

#############
# quicksort #
#############

function quicksort!(a, lo, hi)
    i, j = lo, hi
    while i < hi
        pivot = a[(lo+hi)>>>1]
        while i <= j
            while a[i] < pivot; i += 1; end
            while a[j] > pivot; j -= 1; end
            if i <= j
                a[i], a[j] = a[j], a[i]
                i, j = i+1, j-1
            end
        end
        if lo < j; quicksort!(a,lo,j); end
        lo, j = i, hi
    end
    return a
end

perf_micro_quicksort(n) = quicksort!(rand(n), 1, n)

########
# πsum #
########

function perf_micro_pisum()
    sum = 0.0
    for j = 1:500
        sum = 0.0
        for k = 1:10000
            sum += 1.0/(k*k)
        end
    end
    return sum
end

###############
# randmatstat #
###############

function stdmean(x)
    n = length(x)
    m = mean(x)
    sqrt(sum(xi->abs2(xi - m), x) / (n - 1)) / m
end

function perf_micro_randmatstat(t)
    n = 5
    v = zeros(t)
    w = zeros(t)
    for i=1:t
        a = randn(n,n)
        b = randn(n,n)
        c = randn(n,n)
        d = randn(n,n)
        P = [a b c d]
        Q = [a b; c d]
        v[i] = tr((P' * P)^4)
        w[i] = tr((Q' * Q)^4)
    end
    return (stdmean(v), stdmean(w))
end

##############
# randmatmul #
##############

perf_micro_randmatmul(t) = rand(t,t)*rand(t,t)

#################
# print_to_file #
#################

using Printf

function perf_printfd(n)
    open("/dev/null", "w") do io
        for i = 1:n
            @printf(io, "%d %d\n", i, i + 1)
        end
    end
end

ext: jl


text:
# Multithreaded multiplication using Threads.@spawn (based on https://github.com/jagot/ThreadedSparseArrays.jl)

## * Threading utilities
struct RangeIterator
    k::Int
    d::Int
    r::Int
end

"""
    RangeIterator(n::Int,k::Int)
Returns an iterator splitting the range `1:n` into `min(k,n)` parts of (almost) equal size.
"""
RangeIterator(n::Int, k::Int) = RangeIterator(min(n,k),divrem(n,k)...)
Base.length(it::RangeIterator) = it.k
endpos(it::RangeIterator, i::Int) = i*it.d+min(i,it.r)
Base.iterate(it::RangeIterator, i::Int=1) = i>it.k ? nothing : (endpos(it,i-1)+1:endpos(it,i), i+1)

"""
    tmul!(y::AbstractVector, A::SparseMatrixCSR, x::AbstractVector, alpha::Number, beta::Number)
    tmul!(y::AbstractVector, A::SparseMatrixCSR, x::AbstractVector)

Evaluates `y = alpha*A*x + beta*y` (`y = A*x`)

In-place multithreaded version of sparse csr matrix - vector multiplication, using the threading provided by Threads.@spawn
"""
function tmul!(y::AbstractVector, A::SparseMatrixCSR, x::AbstractVector, alpha::Number, beta::Number)

    A.n == size(x, 1) || throw(DimensionMismatch())
    A.m == size(y, 1) || throw(DimensionMismatch())

    o = getoffset(A)

    @sync for r in RangeIterator(size(y, 1), matmul_num_threads[])
        Threads.@spawn for row in r
            @inbounds begin
                accu = zero(eltype(y))
                for nz in nzrange(A, row)
                    col = A.colval[nz] + o
                    accu += A.nzval[nz]*x[col]
                end
                y[row] = alpha*accu + beta*y[row]
            end
        end
    end

    return y

end

function tmul!(y::AbstractVector, A::SparseMatrixCSR, x::AbstractVector)

    tmul!(y, A, x, true, false)

end

"""
    tmul(A::SparseMatrixCSR, x::AbstractVector)

Evaluates `A*x`.

Multithreaded version of sparse csr matrix - vector multiplication, using the threading provided by Threads.@spawn
"""
function tmul(A::SparseMatrixCSR, x::AbstractVector)

    y = similar(x, promote_type(eltype(A), eltype(x)), size(A, 1))

    tmul!(y, A, x, true, false)

end


text:
using Flight
using OrdinaryDiffEq
using SciMLBase
using LinearAlgebra
using BenchmarkTools

function ac_benchmark()

    h_trn = AltOrth()
    trn = HorizontalTerrain(altitude = h_trn)
    atm = System(AtmosphereDescriptor())
    ac = System(C172Aircraft());

    kin_init = KinInit( Ob = Geographic(LatLon(), h_trn + 0.9));

    init!(ac, kin_init)

    b = @benchmark f_cont!($ac, $trn, $atm)
    @show median(b)

    ac_mdl = Model(ac, (trn, atm); dt = 0.01, adaptive = false, solver = Heun(), y_saveat = 0.1);

    b = @benchmark step!($ac_mdl) setup=(reinit!($ac_mdl))
    @show median(b)

    b = @benchmark step!($ac_mdl, 1, true) setup=(reinit!($ac_mdl))
    @show median(b)

    ac_mdl = Model(ac, (trn, atm); dt = 0.02, adaptive = false, solver = RK4(), y_saveat = 0.1);

    b = @benchmark step!($ac_mdl) setup=(reinit!($ac_mdl))
    @show median(b)

    b = @benchmark step!($ac_mdl, 1, true) setup=(reinit!($ac_mdl))
    @show median(b)

    ac_mdl = Model(ac, (trn, atm)); #save at all steps

    b = @benchmark step!($ac_mdl) setup=(reinit!($ac_mdl))
    @show median(b)

    b = @benchmark step!($ac_mdl, 1, true) setup=(reinit!($ac_mdl))
    @show median(b)

end

function ac_test02()

    h_trn = AltOrth(811)

    trn = HorizontalTerrain(altitude = h_trn)
    atm = System(AtmosphereDescriptor())
    ac = System(C172Aircraft());

    kin_init = KinInit(v_eOb_b = [10, 0, 0],
                        ω_lb_b = [0, 0, 0],
                        q_nb = REuler(ψ = 0, θ = 0.1, φ = 0),
                        Ob = Geographic(LatLon(ϕ = deg2rad(40.531818), λ = deg2rad(-3.574862)),
                                        h_trn + 0.85));

    init!(ac, kin_init)
    ac.u.throttle = 0
    ac.u.brake_left = 0
    ac.u.brake_right = 0
    atm.u.wind.v_ew_n[1] = 10

    ac_mdl = Model(ac, (trn, atm); dt = 0.01, adaptive = false, solver = Heun());
    # reinit!(ac_mdl, tf = 20)
    step!(ac_mdl, 5, true)
    ac.u.pedals = 1
    step!(ac_mdl, 5, true)

    plot_settings = (linewidth=2, margin = 10mm, guidefontsize = 12)
    plots(ac_mdl; save_path = joinpath("tmp", "plots_test02"), plot_settings...)

    return nothing

end

function ac_test03()

end


text:
module FieldOperationsTests

using Gridap.Arrays
using Gridap.Fields
using Gridap.TensorValues
using FillArrays

# using BenchmarkTools
using Test

# Operations

d = 2
p = Point(1.0,2.0)

np = 3
x = fill(p,np)

v = VectorValue{d}(1.0,1.0)
f = MockField(v)
fp = evaluate(f,p)

df = f+f
∇fp = 2*evaluate(gradient(f),p)
test_field(df,p,2*fp,grad=2.0*∇fp)

c = return_cache(df,p)
# @btime evaluate!($c,$df,$p)

∇df = ∇(df)
c = return_cache(∇df,p)
# @btime evaluate!($c,$∇df,$p)

df = f-f
test_field(df,p,fp-fp,grad=0.0*∇fp)

df = ConstantField(2.0)*f
test_field(df,p,fp*2.0,grad=2.0*∇fp)

c = return_cache(df,p)
# @btime evaluate!($c,$df,$p)

q(x) = 2*x
∇q = gradient(q)

f = GenericField(q)
fp = evaluate(f,p)
∇f = ∇(f)
∇fp = evaluate(∇f,p)
@test ∇fp == TensorValue(2.0, 0.0, 0.0, 2.0)
test_field(f,p,fp,grad=∇fp)

c = return_cache(f,p)
# @btime evaluate!($c,$f,$p)

df = f+f
test_field(df,p,2*fp,grad=2.0*∇fp)

c = return_cache(df,p)
# @btime evaluate!($c,$df,$p)

df = f-f
test_field(df,p,fp-fp,grad=0.0*∇fp)

df = ConstantField(2.0)*f
evaluate(df,p)
∇df = ∇(df)
∇dfp = 2*∇fp
evaluate(∇df,p)
test_field(df,p,fp*2.0,grad=∇dfp)

c = return_cache(df,p)
# @btime evaluate!($c,$df,$p)
c = return_cache(∇df,p)
# @btime evaluate!($c,$∇df,$p)

df = f⋅f
∇df = ∇(df)
dfp = (∇fp⋅fp)*2
test_field(df,p,fp⋅fp,grad=dfp)

c = return_cache(df,p)
# @btime evaluate!($c,$df,$p)
c = return_cache(∇df,p)
# @btime evaluate!($c,$∇df,$p)

bdf = df
evaluate(bdf,x)
∇bdf = gradient(bdf)
evaluate(∇bdf,x)
bdfp = fill(fp⋅fp,np)
∇bdfp = fill(dfp,np)
test_field(bdf,x,bdfp,grad=∇bdfp)

c = return_cache(bdf,x)
# @btime evaluate!($c,$bdf,$x)
c = return_cache(∇bdf,x)
# @btime evaluate!($c,$∇bdf,$x)

# Composition

q(x) = VectorValue(x[1]^2,x[2]^2)

f = GenericField(q)
bf = f
∇bf = ∇(bf)

evaluate(bf,x)
bdfx = evaluate(bf,evaluate(bf,x))
bdf = evaluate(Operation(bf),bf)
@test evaluate(bdf,x) == bdfx
test_field(bdf,x,bdfx)#,grad=∇bdfp)

bdfx = evaluate(bf+bf,x)
bdf = evaluate(Operation(+),bf,bf)
@test evaluate(bdf,x) == bdfx
test_field(bdf,x,bdfx)#,grad=∇bdfp)

bdfp = evaluate(bf⋅bf,p)
bdfx = evaluate(bf⋅bf,x)
bdf = bf⋅bf
@test evaluate(bdf,x) == bdfx
test_field(bdf,x,bdfx)#,grad=∇bdfp)

evaluate(bf,x)
evaluate(∇bf,x)
evaluate(bdf,x)
evaluate(∇bdf,x)
∇bdfx = evaluate(∇bf,x).⋅evaluate(bf,x)*2
evaluate(∇bdf,x) == ∇bdfx
test_field(bdf,x,bdfx,grad=∇bdfx)

# Composition

bdf = bf∘bf
∇bdf = ∇(bdf)
bdfx = evaluate(bdf,x)
∇bdfx = evaluate(∇bf,evaluate(bf,x)).⋅evaluate(∇bf,x)
test_field(bdf,x,bdfx,grad=∇bdfx)

c = return_cache(∇bdf,x)
# @btime evaluate!(c,∇bdf,x)

end # module


text:
module SpringCollab2020NoRefillField

using ..Ahorn, Maple

@mapdef Trigger "SpringCollab2020/NoRefillField" NoRefillField(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight)

const placements = Ahorn.PlacementDict(
    "No Refill Field (Spring Collab 2020)" => Ahorn.EntityPlacement(
        NoRefillField,
        "rectangle",
    ),
)

end


text:
module TestRegress
	using Test
	using kNN
	using StatsBase
	using Statistics
	using Random

	Random.seed!(1)
	n = 1_000
	x = 10 * randn(n)
	y = sin.(x) + 0.5 * randn(n)

	fit = kernelregression(x, y)
	grid = minimum(x):0.1:maximum(x)
	predictions = predict(fit, grid)

	fit = kernelregression(x, y, kernel = :gaussian)
	grid = minimum(x):0.1:maximum(x)
	predictions = predict(fit, grid)

	fit = kernelregression(x, y, kernel = :gaussian, bandwidth = 1.0)
	grid = minimum(x):0.1:maximum(x)
	predictions = predict(fit, grid)
end


text:
module Simulate_ContactBoxOnTable

using Modia3D
import Modia3D.ModiaMath

vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)    # material of SolidFileMesh
vmat2 = Modia3D.Material(color="Green" , transparency=0.5)

#=
workpiece_Lx = 0.04
workpiece_Ly = 0.04
workpiece_Lz = 0.08
=#

zDim = 0.02
boxDim = 0.04
Lx = 0.04
Ly = 0.04
Lz = 0.08

rsmall = 0.001   # 0.03    # 0.03
@assembly ContactBoxOnTable begin
  world = Modia3D.Object3D(visualizeFrame=true)

  box = Modia3D.ContactBox(world, scale=[Lx, Ly, Lz], rsmall=rsmall,
                           massProperties="DryWood", material=vmat1,
                           contactMaterial = "DryWood", r=[0.1, 0.1, Lz/2], fixed=false)
  table = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidBox(0.6,0.6,zDim) , "DryWood", vmat2; contactMaterial = "DryWood"); r=[0.3, 0.3, -zDim/2], fixed=true )
end

gravField = Modia3D.UniformGravityField(g=9.81, n=[0,0,-1])
contactBoxOnTable = ContactBoxOnTable(sceneOptions=Modia3D.SceneOptions(gravityField=gravField,visualizeFrames=false, defaultFrameLength=0.2,nz_max = 100,
                                      enableContactDetection=true, elasticContactReductionFactor=1e-4,
                                      visualizeContactPoints=false, visualizeSupportPoints=false))

# Modia3D.visualizeAssembly!( threeD )

model = Modia3D.SimulationModel( contactBoxOnTable )
#ModiaMath.print_ModelVariables(model)
result = ModiaMath.simulate!(model; stopTime=2.0, tolerance=1e-5, interval=0.001, log=false)

ModiaMath.plot(result, ["box.box.r[3]","box.box.v[3]"])

println("... success of Simulate_ContactBoxOnTable.jl!")

end


text:
######################################
############################  PREAMBLE
######################################

### Set the environment
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))

### precompile and add packages to the namespace
using Revise, Random, Mimi, MimiIWG

######################################
######################  SET PARAMETERS
######################################

seed = 42
N = 100 # set thee number of trials
years = collect(2020:5:2060)
discount_rates = [0.02,0.025,0.03,0.05]

######################################
################################  DICE
######################################

# :HFC134a
Random.seed!(seed)
MimiIWG.run_scc_mcs(DICE, gas=:HFC134a, trials=N,
                    perturbation_years=years,
                    discount_rates=discount_rates,
                    domestic = true,
                    tables= true,
                    output_dir=joinpath(@__DIR__, "..", "data/hfc134a/dice"))

######################################
################################  PAGE
######################################

# :HFC134a
Random.seed!(seed)
MimiIWG.run_scc_mcs(PAGE, gas=:HFC134a, trials=N,
                    perturbation_years=years,
                    discount_rates=discount_rates,
                    domestic = true,
                    tables= true,
                    output_dir=joinpath(@__DIR__, "..", "data/hfc134a/page"))

######################################
################################  FUND
######################################

# :HFC134a
Random.seed!(seed)
MimiIWG.run_scc_mcs(FUND, gas=:HFC134a, trials=N,
                    perturbation_years=years,
                    discount_rates=discount_rates,
                    domestic = true,
                    tables= true,
                    output_dir=joinpath(@__DIR__, "..", "data/hfc134a/fund"))

# ## END OF SCRIPT. Have a great day!

max_stars_repo_name: bryanparthum/MimiIWG-HFC

max_issues_repo_name: bryanparthum/MimiIWG-HFC

max_forks_repo_name: bryanparthum/MimiIWG-HFC


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

using Markdown
using InteractiveUtils

# ╔═╡ 9396fcac-b0b6-11eb-3a60-9f2ce25df953
begin
	using Pkg
	Pkg.activate(".")
	Pkg.instantiate()
end

# ╔═╡ 4ecf1267-aac9-4e8f-a522-5595b1205f43
begin
	using CSV, RCall, DataFrames, MixedModels
	using LinearAlgebra, Statistics, StatsBase
	using CategoricalArrays, Arrow
	using AlgebraOfGraphics
	using AlgebraOfGraphics: linear
	using CairoMakie  # for Scatter
end

# ╔═╡ eebef932-1d63-41e5-998f-9748379c43af
md"""
# Mixed Models Tutorial: Conditional Modes

Ths script uses a subset of data reported in Fühner, Golle, Granacher, & Kliegl (2021).
Physical fitness in third grade of primary school:
A mixed model analysis of 108,295 children and 515 schools.

To circumvent delays associated with model fitting we work with a reduced data set and less complex models than those in the reference publication. All the data to reproduce the models in the publication are used here, too; the script requires only a few changes to specify the more complex models in the paper.

The script is structured in three main sections:

1. **Setup** with reading and examing the data, plotting the main results, and specifying the contrasts for the fixed factor `Test`.
2. **Extraction of conditional modes**
3. **Caterpillar plot**

## 1. Setup
### 1.0 Packages and functions
"""

# ╔═╡ faafd359-fb38-4f49-9dfd-19cb4f0c5a54
begin
	function viewdf(x)
		DataFrame(Arrow.Table(take!(Arrow.write(IOBuffer(), x))))
	end
end

# ╔═╡ ee85abd9-0172-44d3-b03a-c6a780d33c72
md"""
### 1.1 Readme for 'EmotikonSubset.rds'

1. Cohort: 9 levels; 2011-2019
2. School: 46 levels
3. Child: 11566 levels; all children were between 6.0 and 6.99 years at legal keydate (30 September) of school enrollement.
4. Sex: 5893 girls, 5673 boys
5. age: test date - middle of month of birthdate (ranges between 7.9 and 9.2)
6. Test: 5 levels
     + Endurance (`Run`):  6 minute endurance run [m]; to nearest 9m in 9x18m field
     + Coordination (`Star_r`): star coordination run [m/s]; 9x9m field, 4 x diagonal = 50.912 m
     + Speed(`S20_r`): 20-meters sprint [m/s]
     + Muscle power low (`SLJ`): standing long jump [cm]
     + Muscle power up (`BPT`): 1-kg medicine ball push test [m]
 7. score - see units

### 1.2 Preprocessing
#### Read data
"""

# ╔═╡ 39443fb0-64ec-4a87-b072-bc8ad6fa9cf4
dat = rcopy(R"readRDS('./data/EmotikonSubset.rds')");

# ╔═╡ 0051ec83-7c30-4d28-9dfa-4c6f5d0259fa
md"""
#### Transformations
"""

# ╔═╡ 7440c439-9d4c-494f-9ac6-18c9fb2fe144
begin
	transform!(dat, :age, :age => (x -> x .- 8.5) => :a1); # centered age (linear)
	transform!(dat,  :a1, :a1  => (x -> x.^2) => :a2);     # centered age (quadr.)
	select!(groupby(dat,  :Test), :, :score => zscore => :zScore); # z-score
end;

# ╔═╡ 16c3cdaa-49fa-46d5-a844-03094329fe4c
viewdf(dat)

# ╔═╡ 64cc1f8e-f831-4a53-976f-dc7600b5634d
# ... by Test and Sex
begin
	dat2 = combine(groupby(dat, [:Test, :Sex]),
                             :score => mean, :score  => std,
                             :zScore => mean, :zScore => std)
	viewdf(dat2)
end

# ╔═╡ c6c6f056-b8b9-4190-ac14-b900bafa04df
md"""
## 2 LMMs and extraction of conditional modes
### 2.1 _SeqDiffCoding_ of `Test`

_SeqDiffCoding_ was used in the publication. This specification tests pairwise
differences between the five neighboring levels of `Test`, that is:

+ H1: `Star_r` - `Run` (2-1)
+ H2: `S20_r` - `Star_r` (3-2)
+ H3: `SLJ` - `S20_r` (4-3)
+ H4: `BPT` - `SLJ` (5-4)

Various options for contrast coding are the topic of the *MixedModelsTutorial_contrasts.jl*
notbebook.
"""

# ╔═╡ c5326753-a03b-4739-a82e-90ffa7c1ebdb
begin
	contr = merge(
        Dict(nm => Grouping() for nm in (:School, :Child, :Cohort)),
		Dict(:Sex => EffectsCoding(; levels=["Girls", "Boys"])),
	    Dict(:Test => SeqDiffCoding(; levels=["Run", "Star_r", "S20_r", "SLJ", "BPT"])),
        Dict(:TestHC => HelmertCoding(; levels=["S20_r", "SLJ", "Star_r", "Run", "BPT"])),
	   )
end;

# ╔═╡ e3cf1ee6-b64d-4356-af23-1dd5c7a0fec6
md"""
#### 2.2 LMMs  `m11` for the reduced data

We fit LMM `m11_Cohort` (only nine levels, six VCs for Test; no CPs) and `m11_School` (46 levels; VCs and CPs for test and linear effect of age)
"""

# ╔═╡ 0dd0d060-81b4-4cc0-b306-dda7589adbd2
f11_Cohort =  @formula zScore ~ 1 + Test * a1 * Sex +  zerocorr(1 + Test | Cohort)

# ╔═╡ f68f6884-7c37-4868-b0b7-01ad450ee383
m11_Cohort = fit(MixedModel, f11_Cohort, dat)

# ╔═╡ ef1dfbf4-700c-4093-9824-01368c511531
VarCorr(m11_Cohort)

# ╔═╡ ace9a2f8-075a-49b5-94eb-12be4ffb6568
f11_School =  @formula zScore ~ 1 + Test * a1 * Sex +  (1 + Test + a1 | School)

# ╔═╡ 167fbd71-2a7b-4afa-8ffb-815dcda63129
m11_School = fit(MixedModel, f11_School, dat)

# ╔═╡ 68ee8abb-04e6-45ee-8f0a-033e7d2ef57b
VarCorr(m11_School)

# ╔═╡ dd8b5c05-c1c8-4778-b179-1ad39f50f8d3
f12 = @formula zScore ~ 1 + Test * a1 * Sex +
                       (1 + Test + a1 | School) + zerocorr(1 + Test | Cohort)

# ╔═╡ 25376bb3-4673-4c89-8600-6e3da85971b4
m12 = fit(MixedModel, f12, dat)

# ╔═╡ e7cdd374-57c3-4d3c-9b1b-b6a211733b81
VarCorr(m12)

# ╔═╡ 5ffbc2ca-02df-4dd7-8e3b-f692d7d12adb

# ╔═╡ 1a2586a6-ff9d-4067-a5af-6c0ee6929e81
md"""
### 2.2 Conditional modes
"""

# ╔═╡ d2c0979e-ed52-4b56-b3ec-58db1609843e
cms = raneftables(m11_School)

# ╔═╡ 776cec2c-c566-4f40-81f5-12546c6a4da7
cms_Chrt = DataFrame(cms.School);

# ╔═╡ cf3f5312-6702-4069-95aa-862dc54d9a34
cvs = condVar(m11_School)

# ╔═╡ 10ef55d6-6048-4353-92e5-e4621bfb130e
md"""
### 3 MixedModelsMakie functions

Information on random effects conditional modes/means, variances, etc.
"""

# ╔═╡ cbe2197e-12c7-4d55-8817-d25d7b7cc5bd
struct RanefInfo{T<:AbstractFloat}
    cnames::Vector{String}
    levels::Vector
    ranef::Matrix{T}
    stddev::Matrix{T}
end

# ╔═╡ 32acfdf0-ecb6-4acb-904f-d7b963771d6d
md"""

ranefinfo(m::LinearMixedModel)

Return a `NamedTuple{fnames(m), NTuple(k, RanefInfo)}` from model `m`
"""

# ╔═╡ 90ce17cc-d497-4605-aa62-be69343e03eb
function ranefinfo(m::LinearMixedModel{T}) where {T}
    fn = fnames(m)
    val = sizehint!(RanefInfo[], length(fn))
    for (re, eff, cv) in zip(m.reterms, ranef(m), condVar(m))
        push!(
            val,
            RanefInfo(
                re.cnames,
                re.levels,
                Matrix(adjoint(eff)),
                Matrix(adjoint(dropdims(sqrt.(mapslices(diag, cv, dims=1:2)), dims=2))),
            )
        )
    end
    NamedTuple{fn}((val...,))
end

# ╔═╡ cd4d8409-cb5e-4f38-bd63-c54b8f09e5ff
function caterpillar!(f::Figure, r::RanefInfo; orderby=1)
    rr = r.ranef
    vv = view(rr, :, orderby)
    ord = sortperm(vv)
    y = axes(rr, 1)
    cn = r.cnames
    axs = [Axis(f[1, j]) for j in axes(rr, 2)]
    linkyaxes!(axs...)
    for (j, ax) in enumerate(axs)
        xvals = view(rr, ord, j)
        scatter!(ax, xvals, y, color=(:red, 0.2))
        errorbars!(ax, xvals, y, view(r.stddev, ord, j), direction=:x)
        ax.xlabel = cn[j]
        ax.yticks = y
        j > 1 && hideydecorations!(ax, grid=false)
    end
    axs[1].yticks = (y, r.levels[ord])
    f
end

# ╔═╡ bd10bfa8-9e22-4927-b83e-eb80f328c6a3
md"""
    caterpillar(m::LinearMixedModel, gf::Symbol)

Returns a `Figure` of a "caterpillar plot" of the random-effects means and prediction intervals

A "caterpillar plot" is a horizontal error-bar plot of conditional means and standard deviations
of the random effects.
"""

# ╔═╡ 727059df-d7a5-4714-8dcb-6568d9119510
function caterpillar(m::LinearMixedModel, gf::Symbol=first(fnames(m)))
    caterpillar!(Figure(resolution=(1000,800)), ranefinfo(m)[gf])
end

# ╔═╡ 5a62529e-cc10-4e2c-8d75-125951f2dc34
md"""
  caterpillar!(f::Figure, r::RanefInfo; orderby=1)

Add Axes of a caterpillar plot from `r` to `f`.

The order of the levels on the vertical axes is increasing `orderby` column
of `r.ranef`, usually the `(Intercept)` random effects.
"""

# ╔═╡ a9043efa-17e5-498f-a744-7c0690446263
md"""
## 4. Caterpillar plots
### 4.1 LMM `m11_Cohort`
"""

# ╔═╡ 07a00ce2-4658-49ff-a384-4b5230691a7b
caterpillar(m11_Cohort)

# ╔═╡ 53ed436c-82b9-4331-bef4-0e034c233e7f
md"""
### 4.2 LMM `m11_School`
"""

# ╔═╡ 44b78423-e28c-4787-921e-e4ba4173a31c
caterpillar(m11_School)

# ╔═╡ 8e15a0b8-b445-4bec-8016-7b6eacf3c889
md"""
### 4.3 LMM `m12`
"""

# ╔═╡ 90dc96b9-3cc6-47d3-84d0-6c3eb8a05a5d
last(fnames(m12))

# ╔═╡ 12518738-2269-4019-9572-bfe6c9dc224f
caterpillar(m12)

# ╔═╡ Cell order:
# ╠═9396fcac-b0b6-11eb-3a60-9f2ce25df953
# ╟─eebef932-1d63-41e5-998f-9748379c43af
# ╠═4ecf1267-aac9-4e8f-a522-5595b1205f43
# ╠═faafd359-fb38-4f49-9dfd-19cb4f0c5a54
# ╟─ee85abd9-0172-44d3-b03a-c6a780d33c72
# ╠═39443fb0-64ec-4a87-b072-bc8ad6fa9cf4
# ╟─0051ec83-7c30-4d28-9dfa-4c6f5d0259fa
# ╠═7440c439-9d4c-494f-9ac6-18c9fb2fe144
# ╠═16c3cdaa-49fa-46d5-a844-03094329fe4c
# ╠═64cc1f8e-f831-4a53-976f-dc7600b5634d
# ╟─c6c6f056-b8b9-4190-ac14-b900bafa04df
# ╠═c5326753-a03b-4739-a82e-90ffa7c1ebdb
# ╟─e3cf1ee6-b64d-4356-af23-1dd5c7a0fec6
# ╠═0dd0d060-81b4-4cc0-b306-dda7589adbd2
# ╠═f68f6884-7c37-4868-b0b7-01ad450ee383
# ╠═ef1dfbf4-700c-4093-9824-01368c511531
# ╠═ace9a2f8-075a-49b5-94eb-12be4ffb6568
# ╠═167fbd71-2a7b-4afa-8ffb-815dcda63129
# ╠═68ee8abb-04e6-45ee-8f0a-033e7d2ef57b
# ╠═dd8b5c05-c1c8-4778-b179-1ad39f50f8d3
# ╠═25376bb3-4673-4c89-8600-6e3da85971b4
# ╠═e7cdd374-57c3-4d3c-9b1b-b6a211733b81
# ╠═5ffbc2ca-02df-4dd7-8e3b-f692d7d12adb
# ╠═1a2586a6-ff9d-4067-a5af-6c0ee6929e81
# ╠═d2c0979e-ed52-4b56-b3ec-58db1609843e
# ╠═776cec2c-c566-4f40-81f5-12546c6a4da7
# ╠═cf3f5312-6702-4069-95aa-862dc54d9a34
# ╟─10ef55d6-6048-4353-92e5-e4621bfb130e
# ╠═cbe2197e-12c7-4d55-8817-d25d7b7cc5bd
# ╟─32acfdf0-ecb6-4acb-904f-d7b963771d6d
# ╠═90ce17cc-d497-4605-aa62-be69343e03eb
# ╠═cd4d8409-cb5e-4f38-bd63-c54b8f09e5ff
# ╟─bd10bfa8-9e22-4927-b83e-eb80f328c6a3
# ╠═727059df-d7a5-4714-8dcb-6568d9119510
# ╟─5a62529e-cc10-4e2c-8d75-125951f2dc34
# ╟─a9043efa-17e5-498f-a744-7c0690446263
# ╠═07a00ce2-4658-49ff-a384-4b5230691a7b
# ╟─53ed436c-82b9-4331-bef4-0e034c233e7f
# ╠═44b78423-e28c-4787-921e-e4ba4173a31c
# ╟─8e15a0b8-b445-4bec-8016-7b6eacf3c889
# ╠═90dc96b9-3cc6-47d3-84d0-6c3eb8a05a5d
# ╠═12518738-2269-4019-9572-bfe6c9dc224f


text:
using Documenter, RiskAdjustedLinearizations

makedocs(
    modules = [RiskAdjustedLinearizations],
    doctest = false,
    strict = false,
    clean = false,
    format = Documenter.HTML(; # prettyurls = get(ENV, "CI", "false") == "true",
                             canonical = "https://juliadocs.github.io/Documenter.jl/stable/",
                             assets=String[],
                             ),
    sitename = "RiskAdjustedLinearizations.jl",
    authors = "William Chen",
    linkcheck = false,
    pages = ["Home" => "index.md",
             "Risk-Adjusted Linearizations" => "risk_adjusted_linearization.md",
             "Sparse Arrays and Jacobians" => "sparse_arrays_jacs.md",
             "Numerical Algorithms" => "numerical_algorithms.md",
             "Example" => "example.md",
             "Caching" => "caching.md",
             "Diagnostics" => "diagnostics.md",
             "Tips" => "tips.md",
             ]
)

deploydocs(;
           repo = "github.com/chenwilliam77/RiskAdjustedLinearizations.jl.git",
           target = "build",
           devbranch = "main",
           versions = ["stable" => "v^", "v#.#"],
)


text:
# Automatically generated test file
function super_exponential4()

g = [
54.598150033144236,
21.380942759123343,
9.487735836358526,
4.770733181967603,
2.718281828459045,
1.755054656960298,
1.284025416687741,
1.000000000000000,
1.064494458917860,
1.284025416687741,
1.755054656960298,
2.718281828459045,
4.770733181967603,
9.487735836358526,
21.380942759123343,
54.598150033144236]

ζ_vec = [1751.97, 403.647, 92.9988, 21.4266, 4.93661, 1.13738, 0.262047, 0.0603748, 0.0139101, 0.00320484]

I_sols = [
[1, 16],
[1, 5, 16],
[1, 3, 14, 16],
[1, 2, 4, 14, 16],
[1, 2, 4, 13, 15, 16],
[1, 2, 3, 6, 13, 15, 16],
[1, 2, 3, 5, 12, 14, 15, 16],
[1, 2, 3, 4, 7, 12, 14, 15, 16],
[1, 2, 3, 4, 6, 11, 13, 14, 15, 16],
[1, 2, 3, 4, 5, 8, 11, 13, 14, 15, 16],
[1, 2, 3, 4, 5, 7, 10, 12, 13, 14, 15, 16],
[1, 2, 3, 4, 5, 6, 8, 11, 12, 13, 14, 15, 16],
[1, 2, 3, 4, 5, 6, 8, 10, 11, 12, 13, 14, 15, 16],
[1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12, 13, 14, 15, 16],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]]

f_sols = [
4774.560611981282818,
1751.967784733572444,
195.390035748150694,
107.549941104629397,
26.991994220677043,
14.014950973534724,
4.348921243498808,
2.031256405141903,
0.719701827511926,
0.319354306106106,
0.114089286055787,
0.043210036683377,
0.008458732570034,
0.004006052663499,
0.000000000000909]

return g, ζ_vec, I_sols, f_sols

end


text:
module TestSpawnBackground

using Test
using ThreadPools
using Statistics

# Put in the function to test this in compiled form
# to make sure there is no yield etc introduced from running interpreted
function uses_all_threads()
    bg_nthreads = Threads.nthreads() - 1
    bg_threads = zeros(bg_nthreads)
    futures = map(1:bg_nthreads) do i
        return spawnbg() do
            id = Threads.threadid()
            bg_threads[id - 1] = id
            return id
        end
    end
    foreach(wait, futures)
    return sum(bg_threads) == sum(2:(bg_nthreads + 1))
end

function busy_wait(time_s)
    t = time()
    while time() - t < time_s
    end
    return
end

function count_occurence(list)
    occurences = Dict{Int,Int}()
    for elem in list
        i = get!(occurences, elem, 0)
        occurences[elem] = i + 1
    end
    return occurences
end

function spam_threads(f, spam_factor)
    bg_nthreads = Threads.nthreads() - 1
    n_executions = bg_nthreads * spam_factor
    thread_ids = []
    time_spent = @elapsed begin
        futures = map(1:n_executions) do i
            return spawnbg() do
                f()
                return Threads.threadid()
            end
        end
        thread_ids = map(fetch, futures)
    end
    return time_spent, thread_ids
end

function spam_threads_busy(time_waiting, spam_factor)
    return spam_threads(spam_factor) do
        return busy_wait(time_waiting)
    end
end

@testset "threading" begin
    nthreads = Threads.nthreads()
    bg_nthreads = nthreads - 1
    if bg_nthreads == 0
        @test fetch(spawnbg(()-> Threads.threadid())) == 1
    else
        @testset "scheduling" begin
            # When we quickly schedule nthreads work items, the implementation should use all threads
            @test uses_all_threads()

            spam_factor = 5
            time_spent, thread_ids = spam_threads(() -> nothing, spam_factor)
            occurences = count_occurence(thread_ids)
            # We should spread out work to all threads when spamming lots of tasks
            @test all(x -> x in keys(occurences), 2:bg_nthreads)
            # a few threads may get more work items, but the mean should be equal to the spamfactor
            @test spam_factor == mean(values(occurences))

            time_spent, thread_ids = spam_threads_busy(0.5, spam_factor)
            occurences = count_occurence(thread_ids)
            @test all(x -> x in keys(occurences), 2:bg_nthreads)
            @test spam_factor == mean(values(occurences))
            # I'm not sure how stable this will be on the CI, we may need to tweak the atol
            @test time_spent ≈ 0.5 * spam_factor atol = 0.1
        end
        @testset "Queue contains all threads, after work is done" begin
            @test length(unique(ThreadPools.AVAILABLE_THREADS[].data)) == bg_nthreads
        end
    end

    @testset "error handling" begin
        @test_throws CapturedException checked_fetch(spawnbg(() -> error("hey")))
    end

end

end


text:
using ControlSimTools
using Test

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


text:
# 1D Burgers equation: Eq. (6.3) from [1]
# Space discretization: Eq. (6.4) from [1]

# Reference
# [1] Liu, J. P., Kolden, H. Ø., Krovi, H. K., Loureiro, N. F., Trivisa, K., &
# Childs, A. M. (2020). Efficient quantum algorithm for dissipative nonlinear
# differential equations. arXiv preprint arXiv:2011.03185.

# Burgers equation:
@taylorize function burgers!(du, u, p, t)

    np = 4 # discretization points
    nu = 0.05 # kinematic viscosity
    Δx = 1.0/(np-1)
    c1 = nu/Δx^2 # c1 = ν/Δx²
    c2 = 0.25/Δx # c2 = 1/(4Δx)

    du[1] = zero(u[1])
    du[np] = zero(u[np])
    for i in 2:np-1
        du[i] = c1*(u[i-1]-2*u[i]+u[i+1]) - c2*(u[i+1]^2 - u[i-1]^2)
    end

end


text:
using JuMP, Ipopt
import Base:show,gamma
export DemonCooling, t, gamma, spec, steps, tol, hamiltonian, timer, probability, cooling!
const MAX_COOLING_FAILURE_COUNT = 200
cooling_module_id = 1

type DemonCooling{N}<:AQCprocessor{N}
    t::Real
    gamma::Real
    H::AbstractMatrix
    k::Int # k steps
    tol::Int # tolerance for error correction, at least $tol zeros
    optimize::Int
    failure_count::Int
end

t(dc::DemonCooling) = dc.t
gamma(dc::DemonCooling) = dc.gamma
spec(dc::DemonCooling) = dc.spectrum
steps(dc::DemonCooling) = dc.k
tol(dc::DemonCooling) = dc.tol
timer(dc::DemonCooling) = dc.t

function DemonCooling(H::AbstractSparseMatrix; k::Int=5, tol::Int=2, optimize=:off)
    if optimize == :static
        t,gamma_ = static_opt_parameters(H)
        return DemonCooling{Int(log2(size(H)[1]))}(t,gamma_,H,k,tol,0,0)
    elseif optimize == :off
        maxEigen = eigs(H;nev=1,which=:LR)[1][1]|>real
        t = π/(2*maxEigen)
        gamma_ = 0
        return DemonCooling{Int(log2(size(H)[1]))}(t,gamma_,H,k,tol,1,0)
    elseif optimize == :runtime
        return DemonCooling{Int(log2(size(H)[1]))}(0.0,0.0,H,k,tol,2,0)
    end
end

function static_opt_parameters(H::AbstractSparseMatrix)
    eigens = eigfact(full(H))
    spec = real(eigens[:values])
    val,idx = findmin(spec)
    m = Model()
    @variable(m, t)
    @variable(m, gamma_)
    @NLobjective(m, Min, exp(sum(1-2*sin(e*t-gamma_) for e in spec)/(1-2*sin(val*t-gamma_))))
    solve(m)
    return getvalue(t),getvalue(gamma_)
end

function multi_cooling_obj(t,gamma_,e,k,j,p)
    return binomial(k,j)*p*(0.5-sin(e*t-gamma_)^j*(0.5+sin(e*t-gamma_))^(k-j))
end

function ∇multi_cooling_obj(g,t,gamma_,e,k,j,p)
    phi = e*t-gamma_
    g[1] = e*(k-j)*cos(phi)*(0.5-sin(phi))^j*(0.5+sin(phi))^(k-j-1)-e*j*cos(phi)*(0.5-sin(phi))^(j-1)*(0.5+sin(phi))^(k-j)

    g[2] = (j-k)*cos(phi)*(0.5-sin(phi))^j*(0.5+sin(phi))^(k-j-1)+j*cos(phi)*(0.5-sin(phi))^(j-1)*(0.5+sin(phi))^(-j+k)
end

function cooling_obj_register(H::AbstractSparseMatrix,k::Int, tol::Int, state::AbstractVector)
    eigens = eigfact(full(H))
    spec = real(eigens[:values])
    val,idx = findmin(spec)
    vecs = eigens[:vectors]
    p = abs2(state.'*eigens[:vectors])
    ground = vecs[:,idx]
    obj(x,y) = sum(sum(
            multi_cooling_obj(x,y,spec[i],k,j,p[i])*multi_cooling_obj(x,y,spec[idx],k,j,p[idx])
            for i in eachindex(spec)) for j =tol:k)
    ∇obj(g,x,y) = sum(sum(
                        ∇multi_cooling_obj(g,x,y,spec[i],k,j,p[i])*multi_cooling_obj(x,y,spec[idx],k,j,p[idx])+∇multi_cooling_obj(g,x,y,spec[idx],k,j,p[idx])*multi_cooling_obj(x,y,spec[i],k,j,p[i])
                        for i in eachindex(spec)) for j =tol:k)
    cooling_obj_name = Symbol("cool_obj$(cooling_module_id)")
    global cooling_module_id += 1
    JuMP.register(cooling_obj_name,2,obj,∇obj)
    return cooling_obj_name
end

function runtime_opt_parameters(dc::DemonCooling,aqc::AQC)
    cooling_obj_name = cooling_obj_register(dc.H,dc.k,dc.tol,coeff(aqc.state))
    m = Model()
    spec = real(eigvals(full(dc.H)))
    minEigen, maxEigen = minimum(spec), maximum(spec)
    @variable(m, t_)
    @variable(m, gamma_)
    @constraint(m, -pi<=minEigen*t_-gamma_<=pi)
    @constraint(m, -pi<=maxEigen*t_-gamma_<=pi)
    @NLobjective(m, Max, $(cooling_obj_name)(t_,gamma_))
    solve(m)
    return getvalue(t_),getvalue(gamma_)
end

DemonCooling(H::AbstractMatrix;k::Int=5, tol::Int=2) = DemonCooling(sparse(H);k=k,tol=tol)

function show{N}(io::IO,dc::DemonCooling{N})
    println(io,"----------------------------")
    println(io,"$N qubits demon-like cooling")
    println(io,"----------------------------")
    println(io,"t = $(dc.t)")
    println(io,"gamma = $(dc.gamma)")
    println(io,"for $(dc.k) steps")
    println(io,"tolerence: larger than $(dc.tol) zeros")
end

function run!{N}(dc::DemonCooling{N}, aqc::AQC)
    if dc.optimize == 2
        dc.t,dc.gamma = runtime_opt_parameters(dc,aqc)
    end
    ret=cooling!(aqc.state, dc)
    run_t = 0.0
    if ret==false
        dc.failure_count += 1
        aqc.state = BellState(N)
        run_t = run!(aqc)
        # aqc.process_list = append!(aqc.process_list[1:aqc.status],aqc.process_list)
        # aqc.state = BellState(N)
        return run_t+t(dc)
    elseif ret
        return t(dc)
    else
    end
end

# return a list of probability for each
# reference: README
function probability{V<:AbstractVector,T,N}(
    state::ComputState{V,T,N},
    dc::DemonCooling{N}
    )
    ret = zeros(Float64,dc.k+1)
    eigens = eigfact(full(dc.H))
    spec = real(eigens[:values])
    A = spec
    map!(x->0.5*sin(x*dc.t-dc.gamma), A)
    p = coeff(state).'*eigens[:vectors]
    for j = 0:dc.k
        for i = 1:2^N
            ret[j+1] += binomial(dc.k,j)*abs2(p[i])*(0.5-A[i])^j*(0.5+A[i])^(dc.k-j)
        end
    end
    return ret
end

# return a cooled state
function cooling!{V<:AbstractVector,T,N}(
    state::ComputState{V,T,N},
    dc::DemonCooling{N}
    )
    p = probability(state, dc)
    k = steps(dc)
    eigens = eigfact(full(dc.H))
    spec = real(eigens[:values])
    A = spec
    map!(x->0.5*sin(x*dc.t-dc.gamma), A)

    dice = rand()
    if dice<sum(p[1:tol(dc)])
        return false
    end
    for j=tol(dc)+1:k
        if dice<sum(p[1:j+1])
            state.coeffs = 0.5^k*binomial(k,j)*(I-im*exp(im*dc.gamma)*expm(-im*full(dc.H)*dc.t))^j*(I+im*exp(im*dc.gamma)*expm(-im*full(dc.H)*dc.t))^(k-j)*state.coeffs
            normalize!(state.coeffs)
            return true
        end
    end
    return false
end

max_stars_repo_name: Roger-luo/openQC

max_issues_repo_name: Roger-luo/openQC

max_forks_repo_name: Roger-luo/openQC


text:
using Plots
using DataFrames
using Distributions
using ShiftedArrays
using CategoricalArrays

function simulate_queue(
    until::Real,
    interval_length::Real,
    # job_size_dist::Distribution,
    p1_arrival_dist::Distribution,
    p2_arrival_dist::Distribution,
    service_dist_pub::Distribution,
    service_dist_priv::Distribution,
    p1_prob_public::Real,
    p2_prob_public::Real,)
    now = 0.0
    p1_next_arrival = rand(p1_arrival_dist)
    p2_next_arrival = rand(p2_arrival_dist)
    p1_priv_next_departure = Inf
    p2_priv_next_departure = Inf
    public_next_departure = Inf
    events = [
        p1_next_arrival,
        p2_next_arrival,
        p1_priv_next_departure,
        p2_priv_next_departure,
        public_next_departure,
    ]
    p1_queue = Float64[]
    p2_queue = Float64[]
    public_queue = Float64[]
    public_queue_who = Int8[]
    event_data = DataFrame(time = Union{Float64,Missing}[], eventtype = String[])
    interval_data = DataFrame(interval = Int[], p1_pp = Float64[], p2_pp = Float64[])
    interval = 1
    while now < until
        push!(interval_data, [interval, p1_prob_public, p2_prob_public])
        while now < interval_length * (interval)
            nextevent = argmin(events)
    # p1_next_arrival
            if nextevent == 1
                now = p1_next_arrival
                where_to = rand(Bernoulli(p1_prob_public))
                if where_to == 1
                    if isempty(public_queue)
                        public_next_departure = p1_next_arrival + rand(service_dist_pub)
                    end
                    push!(public_queue, p1_next_arrival)
                    push!(public_queue_who, 1)
                    push!(event_data, [now "p1-public_arrival"])
                else
                    if isempty(p1_queue)
                        p1_priv_next_departure = p1_next_arrival + rand(service_dist_priv)
                    end
                    push!(p1_queue, p1_next_arrival)
                    push!(event_data, [now "p1-private_arrival"])
                end
                p1_next_arrival += rand(p1_arrival_dist)

    # p2_next_arrival
            elseif nextevent == 2
                now = p2_next_arrival
                where_to = rand(Bernoulli(p2_prob_public))
                if where_to == 1
                    if isempty(public_queue)
                        public_next_departure = p2_next_arrival + rand(service_dist_pub)
                    end
                    push!(public_queue, p2_next_arrival)
                    push!(public_queue_who, 2)
                    push!(event_data, [now "p2-public_arrival"])
                else
                    if isempty(p2_queue)
                        p2_priv_next_departure = p2_next_arrival + rand(service_dist_priv)
                    end
                    push!(p2_queue, p2_next_arrival)
                    push!(event_data, [now "p2-private_arrival"])
                end
                p2_next_arrival += rand(p2_arrival_dist)

                # p1_priv_next_departure
                elseif nextevent == 3
                    now = p1_priv_next_departure
                    push!(event_data, [now "p1-private_departure"])
                    popfirst!(p1_queue)
                    p1_priv_next_departure += isempty(p1_queue) ? Inf : rand(service_dist_priv)

                # p2_priv_next_departure
                elseif nextevent == 4
                    now = p2_priv_next_departure
                    push!(event_data, [now "p2-private_departure"])
                    popfirst!(p2_queue)
                    p2_priv_next_departure += isempty(p2_queue) ? Inf : rand(service_dist_priv)

                # public_next_departure
                else
                    now = public_next_departure
                    push!(event_data, [now "p$(popfirst!(public_queue_who))-public_departure"])
                    popfirst!(public_queue)
                    public_next_departure += isempty(public_queue) ? Inf : rand(service_dist_pub)
            end
            events = [
                p1_next_arrival,
                p2_next_arrival,
                p1_priv_next_departure,
                p2_priv_next_departure,
                public_next_departure,
            ]
        end
        interval += 1
    # choose new probabilities (i.e. strategies) here
        # p1_prob_public = p1_prob_public
        # p2_prob_public = p2_prob_public
    end
    return event_data[1:end - 1, :], interval_data
end

"""
simulate_queue

Alternative version where you input only numbers
"""
function simulate_queue(sim_length::Real,
                        int_length::Real,
                        p1_AR::Real,
                        p2_AR::Real,
                        pu_SR::Real,
                        pr_SR::Real,
                        js::Real,)
    processing_rate(base_rate, job_size) = base_rate * (1 / job_size)
    pst, qst = equilibrium_probability(
        p1_AR,
        p2_AR,
        js,
        α -> processing_rate(pr_SR, α),
        α -> processing_rate(pu_SR, α)
    )

    simulate_queue(
        sim_length,
        int_length,
        Exponential(1 / p1_AR),
        Exponential(1 / p2_AR),
        Exponential(processing_rate(pu_SR, js)),
        Exponential(processing_rate(pr_SR, js)),
        pst,
        qst
    )
end

function prepare_event_data_for_plot(event_data::DataFrame, colours::AbstractVector = [:red, :orange, :blue, :green])
    @assert length(colours) == 4
    arr_dept_to_effect(x) = x == "arrival" ? 1 : -1
    event_plot_data = copy(event_data)
    event_plot_data[:, :effect] = arr_dept_to_effect.(event_plot_data[:, :event])
    sort!(event_plot_data, :time)
    event_plot_data = transform!(groupby(event_plot_data, :source),
        :effect => cumsum,
        :time => Base.Fix2(lag, -1) => :changetime,
        ungroup = false
    )
    event_plot_data = transform!(event_plot_data, :effect_cumsum => Base.Fix2(lag, -1) => :next_effect_cumsum, ungroup = true)

    colour_dict = Dict("p1-private" => colours[1], "p1-public" => colours[2], "p2-private" => colours[3], "p2-public" => colours[4])
    translate(x) = colour_dict[x]
    event_plot_data[:, :groupcolour] = translate.(event_plot_data[!, :source])
    return event_plot_data
end

function plot_event_data(plottable_event_data::DataFrame)
    p = plot()
    for g in groupby(plottable_event_data, :source)
        plot!(p, g.time, g.effect_cumsum, colour = first(g.groupcolour), label = first(g.source))
    end
    return p
end

# plot_event_data(prepare_event_data_for_plot(event_data))

function plot_path(df)
    p = plot(title = "Simulation of a queue", xlabel = "Time", ylabel = "Queue length")
    @inbounds for row in eachrow(df)
        plot!(p,
            [row.time, row.changetime],
            [row.effect_cumsum, row.effect_cumsum],
            colour = row.groupcolour,
            label = false
        )
        plot!(p,
            [row.changetime, row.changetime],
            [row.effect_cumsum, row.next_effect_cumsum],
            colour = row.groupcolour,
            label = false)
    end
    firstlinesdf = combine(groupby(df, :source)) do sdf
        sdf[argmin(skipmissing(sdf.time)), :]
    end
    @inbounds for row in eachrow(firstlinesdf)
        plot!(p,
            [0, row.time],
            [0, 0],
            colour = row.groupcolour,
            label = row.source
        )
    end

    tempdf = transform(df[:, [:time, :effect]], :effect => cumsum, :time => Base.Fix2(lag, -1) => :changetime)
    transform!(tempdf, :effect_cumsum => Base.Fix2(lag, -1) => :next_effect_cumsum)
    @inbounds for row in eachrow(tempdf)
        plot!(p,
            [row.time, row.changetime],
            [row.effect_cumsum, row.effect_cumsum],
            colour = :black,
            label = false
        )
        plot!(p,
            [row.changetime, row.changetime],
            [row.effect_cumsum, row.next_effect_cumsum],
            colour = :black,
            label = false,
            legend = :topleft)
    end
    return p
end

# p = plot_path(prepare_event_data_for_plot(event_data))

# savefig(p, "outputs/simulated_queue_plot.pdf")

function match_events(df::AbstractDataFrame)
    source = df[1, :source]
    arrivals, departures = groupby(df, :event)
    arrivals = DataFrame(arrivals)
    departures = DataFrame(departures)
    # Fill rows in departures so that it has the same length a arrivals
    for _ in 1:(nrow(arrivals) - nrow(departures))
        push!(departures, [missing source "departures"])
    end
    rename!(departures, :time => :departure_time)
    rename!(arrivals, :time => :arrival_time)
    select!(departures, Not([:source, :event]))
    # Concatenate arrivals and departures so that the n-th arrival is on the same row as the n-th departure
    hcat(departures, arrivals, makeunique=true)
end

"""
Warning: Will give missing if there are no arrivals of either private or public for a certain player in a certain time
"""
function evaluate_cost(df)
    if nrow(df) != 2
        return missing
    end
    private_wait_time = df[2, :mean_response_time]
    public_wait_time = df[1, :mean_response_time]
    player = df[1, :player]
    p = first(df[1, Regex("$(player)")])
    p * public_wait_time + (1 - p) * private_wait_time
end


text:
using ModelingToolkit, SparseArrays, Test, GalacticOptim, Optim

@variables x y
@parameters a b
loss = (a - x)^2 + b * (y - x^2)^2
sys1 = OptimizationSystem(loss,[x,y],[a,b],name=:sys1)
sys2 = OptimizationSystem(loss,[x,y],[a,b],name=:sys2)

@variables z
@parameters β
loss2 = sys1.x - sys2.y + z*β
combinedsys = OptimizationSystem(loss2,[z],[β],systems=[sys1,sys2],name=:combinedsys)

equations(combinedsys)
states(combinedsys)
parameters(combinedsys)

calculate_gradient(combinedsys)
calculate_hessian(combinedsys)
generate_function(combinedsys)
generate_gradient(combinedsys)
generate_hessian(combinedsys)
ModelingToolkit.hessian_sparsity(combinedsys)

u0 = [
    sys1.x=>1.0
    sys1.y=>2.0
    sys2.x=>3.0
    sys2.y=>4.0
    z=>5.0
]
p = [
    sys1.a => 6.0
    sys1.b => 7.0
    sys2.a => 8.0
    sys2.b => 9.0
    β => 10.0
]

prob = OptimizationProblem(combinedsys,u0,p,grad=true)
sol = solve(prob,NelderMead())
@test sol.minimum < -1e5

prob2 = remake(prob,u0=sol.minimizer)
sol = solve(prob,BFGS(initial_stepnorm=0.0001),allow_f_increases=true)
@test sol.minimum < -1e8
sol = solve(prob2,BFGS(initial_stepnorm=0.0001),allow_f_increases=true)
@test sol.minimum < -1e9

rosenbrock(x, p) =  (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
x0 = zeros(2)
_p  = [1.0, 100.0]

f = OptimizationFunction(rosenbrock,ModelingToolkit.AutoModelingToolkit(),x0,_p,grad=true,hess=true)
prob = OptimizationProblem(f,x0,_p)
sol = solve(prob,Optim.Newton()


text:
using CSV
using CategoricalArrays
using DataFrames
using Distributions
using Flux
using Gadfly
using HypothesisTests
using Omega
using PlotlyJS
using Plots
using RDatasets
using StatsBase
using StatsPlots
using Turing
using UnicodePlots

max_stars_repo_name: lykmapipo/data-science-learning

max_issues_repo_name: lykmapipo/data-science-learning

max_forks_repo_name: lykmapipo/data-science-learning


text:
# ------------------------------------------------------------------------------------------
# # Functions
#
# Topics:
# 1. How to declare a function
# 2. Duck-typing in Julia
# 3. Mutating vs. non-mutating functions
# 4. Some higher order functions
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## How to declare a function
# Julia gives us a few different ways to write a function. The first requires the `function`
# and `end` keywords
# ------------------------------------------------------------------------------------------

function sayhi(name)
    println("Hi $name, it's great to see you!")
end

function f(x)
    x^2
end

# ------------------------------------------------------------------------------------------
# We can call either of these functions like this:
# ------------------------------------------------------------------------------------------

sayhi("C-3PO")

f(42)

# ------------------------------------------------------------------------------------------
# Alternatively, we could have declared either of these functions in a single line
# ------------------------------------------------------------------------------------------

sayhi2(name) = println("Hi $name, it's great to see you!")

f2(x) = x^2

sayhi2("R2D2")

f2(42)

# ------------------------------------------------------------------------------------------
# Finally, we could have declared these as "anonymous" functions
# ------------------------------------------------------------------------------------------

sayhi3 = name -> println("Hi $name, it's great to see you!")

f3 = x -> x^2

sayhi3("Chewbacca")

f3(42)

# ------------------------------------------------------------------------------------------
# ## Duck-typing in Julia
# *"If it quacks like a duck, it's a duck."*

# Julia functions will just work on whatever inputs make sense.

# For example, `sayhi` works on the name of this minor tv character, written as an
# integer...
# ------------------------------------------------------------------------------------------

sayhi(55595472)

# ------------------------------------------------------------------------------------------
# And `f` will work on a matrix.
# ------------------------------------------------------------------------------------------

A = rand(3, 3)
A

f(A)

# ------------------------------------------------------------------------------------------
# `f` will also work on a string like "hi" because `*` is defined for string inputs as
# string concatenation.
# ------------------------------------------------------------------------------------------

f("hi")

# ------------------------------------------------------------------------------------------
# On the other hand, `f` will not work on a vector. Unlike `A^2`, which is well-defined, the
# meaning of `v^2` for a vector, `v`, is not a well-defined algebraic operation.
# ------------------------------------------------------------------------------------------

v = rand(3)

f(v)

# ------------------------------------------------------------------------------------------
# ## Mutating vs. non-mutating functions
#
# As we mentioned before, the convention is that functions followed by `!` alter their
# contents and functions lacking `!` do not.
#
# For example, let's look at the difference between `sort` and `sort!`.
#
# ------------------------------------------------------------------------------------------

v = [3, 5, 2]

sort(v)

v

# ------------------------------------------------------------------------------------------
# `sort(v)` returns a sorted array that contains the same elements as `v`, but `v` is left
# unchanged.

#
# On the other hand, when we run `sort!(v)`, the contents of v are sorted within the array
# `v`.
# ------------------------------------------------------------------------------------------

sort!(v)

v

# ------------------------------------------------------------------------------------------
# ## Changing names vs. mutation
#
# Recall that variable assignments is just a naming game. Now consider this attempt at a
# mutating function:
# ------------------------------------------------------------------------------------------

function setzero!(x)
    x = 0
    return x
end
y = 1

# ------------------------------------------------------------------------------------------
# Do you think this will be able to set `y` to zero? Try it:
# ------------------------------------------------------------------------------------------

setzero!(y)

y

# ------------------------------------------------------------------------------------------
# `y` is just a name for the number `1`, and so was `x` inside the `setzero!` function until
# we decided it would be better served as a name for `0`! So this function doesn't do any
# mutation at all — it just returns a `0`. If we want to actually mutate the argument, we
# need to use indexed assignment (or use a real `!` function that actually does what it
# claims to do). We also need to use a _mutable_ value as the argument — like an array.
# ------------------------------------------------------------------------------------------

function setzero!(x)
    x[1] = 0
    return x
end
y = [1]

setzero!(y)
y

# ------------------------------------------------------------------------------------------
# ## Some higher order functions
#
# ### map
#
# `map` is a "higher-order" function in Julia that *takes a function* as one of its input
# arguments.
# `map` then applies that function to every element of the data structure you pass it. For
# example, executing
#
# ```julia
# map(f, [1, 2, 3])
# ```
# will give you an output array where the function `f` has been applied to all elements of
# `[1, 2, 3]`
# ```julia
# [f(1), f(2), f(3)]
# ```
# ------------------------------------------------------------------------------------------

map(f, [1, 2, 3])

# ------------------------------------------------------------------------------------------
# Here we've squared all the elements of the vector `[1, 2, 3]`, rather than squaring the
# vector `[1, 2, 3]`.
#
# To do this, we could have passed to `map` an anonymous function rather than a named
# function, such as
# ------------------------------------------------------------------------------------------

x -> x^3

# ------------------------------------------------------------------------------------------
# via
# ------------------------------------------------------------------------------------------

map(x -> x^3, [1, 2, 3])

# ------------------------------------------------------------------------------------------
# and now we've cubed all the elements of `[1, 2, 3]`!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### broadcast
#
# `broadcast` is another higher-order function like `map`. `broadcast` is a generalization
# of `map`, so it can do every thing `map` can do and more. The syntax for calling
# `broadcast` is the same as for calling `map`
# ------------------------------------------------------------------------------------------

broadcast(f, [1, 2, 3])

# ------------------------------------------------------------------------------------------
# and again, we've applied `f` (squared) to all the elements of `[1, 2, 3]` - this time by
# "broadcasting" `f`!
#
# Some syntactic sugar for calling `broadcast` is to place a `.` between the name of the
# function you want to `broadcast` and its input arguments. For example,
#
# ```julia
# broadcast(f, [1, 2, 3])
# ```
# is the same as
# ```julia
# f.([1, 2, 3])
# ```
# ------------------------------------------------------------------------------------------

f.([1, 2, 3])

# ------------------------------------------------------------------------------------------
# Notice again how different this is from calling
# ```julia
# f([1, 2, 3])
# ```
# We can square every element of a vector, but we can't square a vector!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# To drive home the point, let's look at the difference between
#
# ```julia
# f(A)
# ```
# and
# ```julia
# f.(A)
# ```
# for a matrix `A`:
# ------------------------------------------------------------------------------------------

A = [i + 3*j for j in 0:2, i in 1:3]

f(A)

# ------------------------------------------------------------------------------------------
# As before we see that for a matrix, `A`,
# ```
# f(A) = A^2 = A * A
# ```
#
# On the other hand,
# ------------------------------------------------------------------------------------------

B = f.(A)

# ------------------------------------------------------------------------------------------
# contains the squares of all the entries of `A`.
#
# This dot syntax for broadcasting allows us to write relatively complex compound
# elementwise expressions in a way that looks natural/closer to mathematical notation. For
# example, we can write
# ------------------------------------------------------------------------------------------

A .+ 2 .* f.(A) ./ A

# ------------------------------------------------------------------------------------------
# instead of
# ------------------------------------------------------------------------------------------

broadcast(x -> x + 2 * f(x) / x, A)

# ------------------------------------------------------------------------------------------
# and the two will perform exactly the same.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Exercises
#
# #### 6.1
# Write a function `add_one` that adds 1 to its input.
# ------------------------------------------------------------------------------------------

@assert add_one(1) == 2

@assert add_one(11) == 12

# ------------------------------------------------------------------------------------------
# #### 6.2
# Use `map` or `broadcast` to increment every element of matrix `A` by `1` and assign it to
# a variable `A1`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### 6.3
# Use the broadcast dot syntax to increment every element of matrix `A1` by `1` and store it
# in variable `A2`
# ------------------------------------------------------------------------------------------

@assert A2 == [3 4 5; 6 7 8;9 10 11]

# ------------------------------------------------------------------------------------------
# Please click on `Validate` on the top, once you are done with the exercises.
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/intro-to-julia/long-version/030 Functions and Broadcasting.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using Pkg
# Pkg.develop(path=joinpath(@__DIR__, "../../DojoEnvironments"))
Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()

# ## visualizer
vis = Visualizer()
open(vis)

# ## setup
using Dojo
using IterativeLQR
using LinearAlgebra
using FiniteDiff
using DojoEnvironments

using BenchmarkTools

mech = get_mechanism(:pendulum)
initialize!(mech, :pendulum, angle=0.1)
storage = simulate!(mech, 1.0)
Main.@profiler storage = simulate!(mech, 1000.0)
@benchmark storage = simulate!(mech, 10.0)

visualize(mech, storage, vis=vis)

function impulse_map(mechanism, joint::JointConstraint, body::Body)
    relative = body.id == joint.parent_id ? :parent : :child
    pbody = get_body(mechanism, joint.parent_id)
    cbody = get_body(mechanism, joint.child_id)
    tra = impulse_map(relative, joint.translational,
        pbody, cbody,
        joint.impulses[2][joint_impulse_index(joint, 1)])
    rot = impulse_map(relative, joint.rotational,
        pbody, cbody,
        joint.impulses[2][joint_impulse_index(joint, 2)])
    return hcat(tra, rot)
end

joint = mech.joints[1]
body = mech.bodies[1]
impulse_map(mech, joint, body)
@benchmark $impulse_map($mech, $joint, $body)

function constraint_jacobian_configuration(relative::Symbol, joint::Joint{T,Nλ,Nb},
        xa::AbstractVector, qa::Quaternion,
        xb::AbstractVector, qb::Quaternion,
        η) where {T,Nλ,Nb}

    ∇comp = szeros(T,Nb,6)
    ∇mincoord = minimal_coordinates_jacobian_configuration(relative, joint, xa, qa, xb, qb, attjac=true)
    ∇unlim = joint_constraint_jacobian_configuration(relative, joint, xa, qa, xb, qb, η)

    return [∇comp; ∇mincoord; -∇mincoord; ∇unlim]
end

function joint_constraint_jacobian_configuration(relative::Symbol, joint::Joint{T},
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion, η) where {T}
    X, Q = displacement_jacobian_configuration(relative, joint, xa, qa, xb, qb, attjac=true)
    return constraint_mask(joint) * [X Q]
end

function impulse_map(relative::Symbol, joint::Joint{T,Nλ,0},
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion,
    η) where {T,Nλ}
    J = constraint_jacobian_configuration(relative, joint, xa, qa, xb, qb, η)
    G = [[Diagonal(sones(T,3)) szeros(T,3,3)]; [szeros(4,3) LVᵀmat(relative == :parent ? qa : qb)]]
    return Diagonal([sones(3); 0.5 * sones(3)]) * transpose(J * G)
end

function constraint_jacobian_configuration(relative::Symbol, joint::Joint{T,Nλ,0},
        xa::AbstractVector, qa::Quaternion,
        xb::AbstractVector, qb::Quaternion, η) where {T,Nλ}
    joint_constraint_jacobian_configuration(relative, joint, xa, qa, xb, qb, η)
end

relative = :parent
joint_tra = mech.joints[1].translational
joint_rot = mech.joints[1].rotational
xa = srand(3)
xb = srand(3)
qa = Quaternion(normalize(rand(4))...)
qb = Quaternion(normalize(rand(4))...)
η = 0.0

pbody = mech.origin
cbody = mech.bodies[1]
impulse_map(relative, joint_rot, pbody, cbody, η)
Main.@code_warntype impulse_map(relative, joint_rot, pbody, cbody, η)
@benchmark $impulse_map($relative, $joint_rot, $pbody, $cbody, $η)

impulse_map(relative, joint_rot, xa, qa, xb, qb, η)
Main.@code_warntype impulse_map(relative, joint_rot, xa, qa, xb, qb, η)
@benchmark $impulse_map($relative, $joint_rot, $xa, $qa, $xb, $qb, $η)

# function displacement_jacobian_configuration_unstable(relative::Symbol, joint::Rotational{T},
#         xa::AbstractVector, qa::Quaternion,
#         xb::AbstractVector, qb::Quaternion;
#         attjac::Bool=true, vmat::Bool=true) where T
#     X = szeros(T, 3, 3)
#     if relative == :parent
# 		Q = Lᵀmat(joint.orientation_offset) * Rmat(qb) * Tmat()
# 		attjac && (Q *= LVᵀmat(qa))
#     elseif relative == :child
# 		Q = Lᵀmat(joint.orientation_offset) * Lᵀmat(qa)
# 		attjac && (Q *= LVᵀmat(qb))
# 	end
# 	vmat && (Q = Vmat() * Q)
# 	return X, Q
# end
#
#
# function displacement_jacobian_configuration3(relative::Symbol, joint::Rotational{T},
#         xa::AbstractVector, qa::Quaternion,
#         xb::AbstractVector, qb::Quaternion;
# 		) where T
#     X = szeros(T, 3, 3)
#     if relative == :parent
# 		Q = Lᵀmat(joint.orientation_offset) * Rmat(qb) * Tmat()
#     elseif relative == :child
# 		Q = Lᵀmat(joint.orientation_offset) * Lᵀmat(qa)
# 	end
# 	return X, Vmat() * Q
# end
#
# function displacement_jacobian_configuration3(relative::Symbol, joint::Translational{T},
#     xa::AbstractVector, qa::Quaternion,
#     xb::AbstractVector, qb::Quaternion) where T
#
#     vertices = joint.vertices
#
#     if relative == :parent
#         d = xb + vector_rotate(vertices[2], qb) - (xa + vector_rotate(vertices[1], qa)) # in the world frame
#         X = -rotation_matrix(inv(qa))
#         Q = -rotation_matrix(inv(qa)) * ∂rotation_matrix∂q(qa, vertices[1])
#         Q += ∂rotation_matrix_inv∂q(qa, d)
#     elseif relative == :child
#         X = rotation_matrix(inv(qa))
#         Q = rotation_matrix(inv(qa)) * ∂rotation_matrix∂q(qb, vertices[2])
#     end
# 	return X, Vmat() * Q
# end
#
# function displacement_jacobian_configuration3(relative::Symbol, joint::Joint{T},
#         xa::AbstractVector, qa::Quaternion,
#         xb::AbstractVector, qb::Quaternion;
#         attjac::Bool=true) where T
#     X, Q = displacement_jacobian_configuration3(relative, joint, xa, qa, xb, qb)
# 	G = relative == :parent ? LVᵀmat(qa) : LVᵀmat(qb)
# 	if attjac
# 		Q = Q * G
# 		return X, Q::SMatrix{3,3,T,9}
# 	else
# 		return X, Q::SMatrix{3,4,T,12}
# 	end
# end

function displacement_jacobian_configuration(relative::Symbol, joint::Joint{T},
        xa::AbstractVector, qa::Quaternion,
        xb::AbstractVector, qb::Quaternion;
        attjac::Bool=true) where T
    X, Q = displacement_jacobian_configuration(relative, joint, xa, qa, xb, qb)
	G = relative == :parent ? LVᵀmat(qa) : LVᵀmat(qb)
	if attjac
		Q = Q * G
		return X, Q::SMatrix{3,3,T,9}
	else
		return X, Q::SMatrix{3,4,T,12}
	end
end

function displacement_jacobian_configuration(relative::Symbol, joint::Translational{T},
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion) where T

    vertices = joint.vertices

    if relative == :parent
        d = xb + vector_rotate(vertices[2], qb) - (xa + vector_rotate(vertices[1], qa)) # in the world frame
        X = -rotation_matrix(inv(qa))
        Q = -rotation_matrix(inv(qa)) * ∂rotation_matrix∂q(qa, vertices[1])
        # Q += ∂rotation_matrix_inv∂q(qa, d)
    elseif relative == :child
        X = rotation_matrix(inv(qa))
        Q = rotation_matrix(inv(qa)) * ∂rotation_matrix∂q(qb, vertices[2])
    end
	return X, Q
end

relative = :parent
child = :child
joint_tra = mech.joints[1].translational
joint_rot = mech.joints[1].rotational
xa = srand(3)
xb = srand(3)
qa = Quaternion(normalize(rand(4))...)
qb = Quaternion(normalize(rand(4))...)
η = 0.0
displacement_jacobian_configuration(
	relative, joint_tra, xa, qa, xb, qb)
displacement_jacobian_configuration(
	relative, joint_rot, xa, qa, xb, qb)
# Main.@code_warntype displacement_jacobian_configuration(
# 	relative, joint_rot, xa, qa, xb, qb, attjac=false)
# Main.@code_warntype displacement_jacobian_configuration(
# 	relative, joint_rot, xa, qa, xb, qb, attjac=true)
Main.@code_warntype displacement_jacobian_configuration(
	relative, joint_tra, xa, qa, xb, qb)
Main.@code_warntype displacement_jacobian_configuration(
	relative, joint_tra, xa, qa, xb, qb)
# @benchmark $displacement_jacobian_configuration(
	# $relative, $joint_rot, $xa, $qa, $xb, $qb, attjac=false)
# @benchmark $displacement_jacobian_configuration(
	# $relative, $joint_rot, $xa, $qa, $xb, $qb, attjac=true)
@benchmark $displacement_jacobian_configuration(
	$relative, $joint_tra, $xa, $qa, $xb, $qb)
@benchmark $displacement_jacobian_configuration(
	$relative, $joint_tra, $xa, $qa, $xb, $qb)

@benchmark $displacement_jacobian_configuration(
	$child, $joint_tra, $xa, $qa, $xb, $qb)
@benchmark $displacement_jacobian_configuration(
	$child, $joint_tra, $xa, $qa, $xb, $qb)

qa = Quaternion(normalize(rand(4))...)
v = srand(3)
∂rotation_matrix∂q(qa, v)
@benchmark $∂rotation_matrix∂q($qa, $v)

function joint_constraint_jacobian_configuration(relative::Symbol, joint::Joint{T},
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion, η) where {T}
    X, Q = displacement_jacobian_configuration(relative, joint, xa, qa, xb, qb, attjac=false)
    return constraint_mask(joint) * [X Q]
end

relative = :parent
joint_tra = mech.joints[1].translational
joint_rot = mech.joints[1].rotational
xa = srand(3)
xb = srand(3)
qa = Quaternion(normalize(rand(4))...)
qb = Quaternion(normalize(rand(4))...)
η = 0.0
constraint_jacobian_configuration(relative, joint_rot, xa, qa, xb, qb, η)
Main.@code_warntype constraint_jacobian_configuration(relative, joint_rot, xa, qa, xb, qb, η)
@benchmark $constraint_jacobian_configuration($relative, $joint_rot, $xa, $qa, $xb, $qb, $η)

relative = :parent
joint_tra = mech.joints[1].translational
joint_rot = mech.joints[1].rotational
xa = srand(3)
xb = srand(3)
qa = Quaternion(normalize(rand(4))...)
qb = Quaternion(normalize(rand(4))...)
η = 0.0
joint_constraint_jacobian_configuration(relative, joint_rot, xa, qa, xb, qb, η)
Main.@code_warntype joint_constraint_jacobian_configuration(relative, joint_rot, xa, qa, xb, qb, η)
@benchmark $joint_constraint_jacobian_configuration($relative, $joint_rot, $xa, $qa, $xb, $qb, $η)

constraint_jacobian_configuration(relative, joint_rot, xa, qa, xb, qb, η)
Main.@code_warntype constraint_jacobian_configuration(relative, joint_rot, xa, qa, xb, qb, η)
@benchmark $constraint_jacobian_configuration($relative, $joint_rot, $xa, $qa, $xb, $qb, $η)

impulse_map(relative, joint_rot, xa, qa, xb, qb, η)
Main.@code_warntype impulse_map(relative, joint_rot, xa, qa, xb, qb, η)
@benchmark $impulse_map($relative, $joint_rot, $xa, $qa, $xb, $qb, $η)

displacement_jacobian_configuration(relative, joint, xa, qa, xb, qb)

aa = szeros(10,10)
bb = szeros(5,5)
diagonal_cat(aa, bb)
@benchmark $diagonal_cat($aa, $bb)

aa = DenseMatrix(Diagonal(sones(4)))
aaa = SMatrix{4,4,Float64,16}(aa)

joint = mech.joints[1]
joint.impulses

constraint_jacobian(joint)

@benchmark $constraint_jacobian($joint)

joint_impulse_index(joint, 1)
joint_impulse_index(joint, 2)

joint = mech.joints[1]
joint.impulses
joint_tra = joint.translational
joint_rot = joint.rotational

imp_tra = szeros(3)
imp_rot = szeros(length(joint_impulse_index(joint, 2)))
constraint_jacobian(joint_tra, imp_tra)
constraint_jacobian(joint_rot, imp_rot)

# @benchmark $constraint_jacobian($joint_tra, $imp_tra)
@benchmark $constraint_jacobian($joint_rot, $imp_rot)

joint_impulse_index(joint, 1)
joint_impulse_index(joint, 2)

a = 10
a = 10
a = 10
a = 10
a = 10
a = 10
a = 10

function reset2!(joint::JointConstraint{T,N,Nc}; scale::T=1.0) where {T,N,Nc}
    Nλ_tra = joint_length(joint.translational)
    Nb_tra = limits_length(joint.translational)
    Nλ_rot = joint_length(joint.rotational)
    Nb_rot = limits_length(joint.rotational)
    joint.impulses[1] = [scale * sones(2Nb_tra); szeros(Nλ_tra); scale * sones(2Nb_rot); szeros(Nλ_rot)]
    joint.impulses[2] = [scale * sones(2Nb_tra); szeros(Nλ_tra); scale * sones(2Nb_rot); szeros(Nλ_rot)]
    return
end

joint = mech.joints[1]

reset!(joint)
Main.@code_warntype reset!(joint)
@benchmark $reset!($joint)

a = 10
a = 10
a = 10
a = 10
a = 10

function constraint_jacobian_configuration(mechanism::Mechanism{T,Nn,Ne,Nb}, body::Body{T}; reg::T=Dojo.REG) where {T,Nn,Ne,Nb}
    state = body.state
    timestep = mechanism.timestep
    mass = body.mass
    inertia = body.inertia

    x2, q2 = current_configuration(state)
    x3, q3 = next_configuration(state, timestep)

    I3 = SMatrix{3,3,T,9}(Diagonal(sones(T,3)))
    Z33 = szeros(T, 3, 3)
    Z34 = szeros(T, 3, 4)

    # dynamics
    dynT = I3 * mass / timestep
    dynR = -2.0 / timestep * LVᵀmat(q2)' * Tmat() * (∂Rᵀmat∂q(Vᵀmat() * inertia * Vmat() * Lmat(q2)' * vector(q3)) + Rmat(q3)' * Vᵀmat() * inertia * Vmat() * Lmat(q2)')

    state.D = [[dynT; Z33] [Z34; dynR]] * integrator_jacobian_velocity(body, timestep)
    state.D += [[reg * I3; Z33] [Z33; reg * I3]]

    # inputs
    nothing

    # impulses
    for id in connections(mechanism.system, body.id)
        Ne < id <= Ne + Nb && continue # body
        impulses_jacobian_velocity!(mechanism, body, get_node(mechanism, id))
    end

    return state.D
end

body = mech.bodies[1]
constraint_jacobian_configuration(mech, body)
Main.@code_warntype constraint_jacobian_configuration(mech, body)
@benchmark $constraint_jacobian_configuration($mech, $body)

a = 10
a = 10
a = 10
a = 10
a = 10
a = 10

mech = get_mechanism(:quadruped)
initialize!(mech, :quadruped)
# Main.@profiler
storage = simulate!(mech, 1.0)
storage = simulate!(mech, 1.0)

visualize(mech, storage, vis=vis)
#

using BenchmarkTools

# function constraint_jacobian_configuration(relative::Symbol, joint::Joint{T,Nλ,Nb},
#         xa::AbstractVector, qa::Quaternion,
#         xb::AbstractVector, qb::Quaternion,
#         η) where {T,Nλ,Nb}
#
#     ∇comp = szeros(T,Nb,7)
#     ∇mincoord = minimal_coordinates_jacobian_configuration(relative, joint, xa, qa, xb, qb, attjac=false)
#     ∇unlim = joint_constraint_jacobian_configuration(relative, joint, xa, qa, xb, qb, η)
#
#     return [∇comp; ∇mincoord; -∇mincoord; ∇unlim]
# end

relative = :parent
joint = mech.joints[5].rotational
xa = srand(3)
xb = srand(3)
qa = Quaternion(normalize(rand(4))...)
qb = Quaternion(normalize(rand(4))...)
η = 0.0

constraint_jacobian_configuration(relative, joint, xa, qa, xb, qb, η)
@benchmark $constraint_jacobian_configuration($relative, $joint, $xa, $qa, $xb, $qb, $η)

a = 10
a = 10
a = 10
a = 10

# function joint_impulse_index(joint::JointConstraint{T,N,Nc}, i::Int) where {T,N,Nc}
#     s = 0
#     for j = 1:i-1
#         element = (joint.translational, joint.rotational)[j]
#         s += impulses_length(element)
#     end
#     joint_impulse_index((joint.translational, joint.rotational)[i], s) # to be allocation free
#     # joint_impulse_index((joint.translational, joint.rotational)[i], s) # to be allocation free
#     # joint_impulse_index((joint.translational, joint.rotational)[i], s) # to be allocation free
# end
#
# function joint_impulse_index(joint::Joint{T,Nλ,Nb,N}, s::Int) where {T,Nλ,Nb,N}
#     ind = SVector{N,Int}(s+1:s+N)
#     return ind
# end

joint_con0 = mech.joints[1]
joint0 = joint_con0.translational
joint0 = joint_con0.rotational
i0 = 1
s0 = 1
@benchmark $joint_impulse_index($joint_con0, $i0)
@benchmark $joint_impulse_index($joint0, $s0)

for (i,el) in enumerate(3,4,5)
    @show el
end

rosrun unitree_controller unitree_servo # let the robot stretch legs
rosrun unitree_controller unitree_move_kinetic # place the robot back to origin

roslaunch unitree_gazebo normal.launch rname:=a1 wname:=stairs_new

rotation_matrix0(q::Quaternion) = VRᵀmat(q) * LVᵀmat(q)
function ∂rotation_matrix∂q0(q::Quaternion, p::AbstractVector; attjac::Bool=false)
    M = ∂VRᵀmat∂q(LVᵀmat(q) * p) + VRᵀmat(q) * ∂LVᵀmat∂q(p)
    attjac && (M *= LVᵀmat(q))
    return M
end
function ∂rotation_matrix_inv∂q0(q::Quaternion, p::AbstractVector; attjac::Bool=false)
    M = ∂rotation_matrix∂q0(inv(q), p) * Tmat()
    attjac && (M *= LVᵀmat(q))
    return M
end

using FiniteDiff

p = srand(3)
q = Quaternion(normalize(srand(4))...)
J0 = FiniteDiff.finite_difference_jacobian(q -> rotation_matrix0(Quaternion(q...)) * p, vector(q)) * LVᵀmat(q)
J1 = ∂rotation_matrix∂q0(q, p, attjac=true)
J0 - J1

p = srand(3)
q = Quaternion(normalize(srand(4))...)
J0 = FiniteDiff.finite_difference_jacobian(q -> rotation_matrix0(inv(Quaternion(q...))) * p, vector(q)) * LVᵀmat(q)
J1 = ∂rotation_matrix_inv∂q0(q, p, attjac=true)
J0 - J1

p = srand(3)
q = Quaternion(normalize(srand(4))..., true)
inv(q)
# q = Quaternion(1,0,0,0.0)
# J0 = FiniteDiff.finite_difference_jacobian(q -> rotation_matrix(inv(Quaternion(q..., true))) * p, vector(q)) * LVᵀmat(q)
J0 = FiniteDiff.finite_difference_jacobian(q -> rotation_matrix(inv(Quaternion(q..., true))) * p, vector(q))# * LVᵀmat(q)
J1 = ∂rotation_matrix_inv∂q(q, p)
J0 - J1

rotation_matrix(q) = VRᵀmat(q) * LVᵀmat(q)


text:
c = ContinuousClock(10)

x = recorded(DynamicEntry, c, 1.0)
y = recorded(DynamicEntry, c, 1.0 + 1.0im)

@test x == recorded(DynamicEntry, c, 1.0) == 1.0
@test y == recorded(DynamicEntry, c, 1.0 + 1.0im) == 1.0 + 1.0im
@test hash(x) == hash(1.0)
@test hash(x, UInt(1)) == hash(1.0, UInt(1))
@test hash(y) == hash(1.0 + 1.0im)
@test hash(y, UInt(1)) == hash(1.0 + 1.0im, UInt(1))

test_show(x, "1.0")
test_show(y, "1.0 + 1.0im")
@test repr(x) == "1.0"
@test repr(y) == "1.0 + 1.0im"
test_recipe(TimeSeries((x,)), [([0], [1.0])])
test_recipe(TimeSeries((y,)), [([0], [1.0 + 1.0im])])

@test isnum(Real, x)
@test isnum(Complex, y)
@test isnum(Complex, x + y)
@test issubtype(typeof(x), Real)
@test issubtype(typeof(y), Complex)
@test issubtype(typeof(x + y), Complex)

@test promote_type(RReal{Int}, Int) == Int
@test promote_type(RNumber{Int}, Int) == Int
@test promote_type(RReal{Float64}, Int) == Float64
@test promote_type(RNumber{Float64}, Int) == Float64
@test promote_type(RReal{Float64}, RReal{Int}) == Float64
@test promote_type(RReal{Float64}, RNumber{Int}) == Float64
@test promote_type(RNumber{Float64}, RReal{Int}) == Float64
@test promote_type(RNumber{Float64}, RNumber{Int}) == Float64

@test x + x == 2
@test y + y == 2 + 2im
@test x - x == 0
@test y - y == 0
@test x * x == 1
@test y * y == 2im
@test x / x == 1
@test y / y == 1
@test x \ x == 1
@test y \ y == 1
@test x^x == 1
@test y^y == (1 + im)^(1 + im)

@test convert(Int, x) == 1
@test convert(RNumber{Int}, x) == 1

xs = [x[1]]
ys = [y[1]]

increase!(c, 1) # t = 1
x[CartesianIndex(1)] = x + x
y[1] += 1

@test x == 2
@test y == 2 + 1im
push!(xs, x)
push!(ys, y)

increase!(c, 1) # t = 2
x[CartesianIndex()] += UInt(1)
y[1, 1] = y + 2im

@test x == 3
@test y == 2 + 3im
push!(xs, x)
push!(ys, y)

increase!(c, 1) # t = 3

x[CartesianIndex(1)] = x - UInt(1)
y[CartesianIndex(1, 1)] = y' - 1

@test x == 2
@test y == 1 - 3im
push!(xs, x)
push!(ys, y)

increase!(c, 1) # t = 4
x[] = x * x
y[] = y / y

@test x == 4
@test y == 1
push!(xs, x)
push!(ys, y)

increase!(c, 1) # t = 5
x[] = x^2
y[] = (y^2) * im

@test x == 16
@test y == 1im
push!(xs, x)
push!(ys, y)

increase!(c, 1) # t = 6
x[] += real(y)
y[CartesianIndex(1, 1)] += imag(y)

@test x == 16
@test y == 1 + 1im
push!(xs, x)
push!(ys, y)

increase!(c, 1) # t = 7
x[] = 2 \ x
y[] = y^y

@test x == 8
@test y == (1 + im)^(1 + im)
push!(xs, x)
push!(ys, y)

increase!(c, 1) # t = 8

x[] = exp(x)
y[] = float(y)

@test x == exp(8)
@test y == (1 + 1im)^(1 + 1im)
push!(xs, x)
push!(ys, y)

increase!(c, 1) # t = 9

x[] = x / x
y[] = y[1]' + y[1, 1]

@test x == 1
@test y == 2 * real((1 + 1im)^(1 + 1im))
push!(xs, x)
push!(ys, y)

e1 = getentries(x)
e2 = getentries(y)

@test getts(e1) == 0:9
@test getts(e2) == 0:9
@test getvs(e1) == xs
@test getvs(e2) == ys


text:
# AbstractSurrogateWrapper
Base.broadcastable( sw :: AbstractSurrogateWrapper ) = Ref(sw)

# constructor
init_wrapper( :: Type{<:AbstractSurrogateWrapper}, objf :: AbstractObjective, mod :: SurrogateModel, meta :: SurrogateMeta ) :: AbstractSurrogateWrapper = nothing

num_outputs( sw :: AbstractSurrogateWrapper ) :: Int = 0
get_objf( sw :: AbstractSurrogateWrapper ) :: AbstractObjective = nothing
get_model( sw :: AbstractSurrogateWrapper ) :: SurrogateModel = nothing
get_meta( sw :: AbstractSurrogateWrapper ) :: SurrogateMeta = nothing

# AbstractSurrogateContainer
Base.broadcastable( sc :: AbstractSurrogateContainer ) = Ref(sc)

# constructor
function init_surrogates( :: Type{<:AbstractSurrogateContainer}, :: AbstractMOP, :: AbstractIterData,
	:: AbstractDB, :: AbstractConfig )
	return nothing
end

# update methods
function update_surrogates!( :: AbstractSurrogateContainer, :: AbstractMOP, :: AbstractIterData,
	:: AbstractDB; ensure_fully_linear = false )
	return nothing
end

function improve_surrogates!( :: AbstractSurrogateContainer, :: AbstractMOP, :: AbstractIterData,
	:: AbstractDB; ensure_fully_linear = false, kwargs... )
	return nothing
end

# mandatory functions

list_of_wrappers( sc :: AbstractSurrogateContainer ) :: AbstractVector{<:AbstractSurrogateWrapper } = nothing

"""
	_output_model_mapping(sc, l)

Return a `Tuple{Int,Int}` with the index `i` of the surrogate model and its output `k`
describing output `l` of an MOP with scalarized outputs `1:l`.

# Example

Suppose an MOP is setup with the 3 scalar objectives f1, f2 and f3 and
there is one surrogate model modelling f1 and f3 (in that order) and
a second model for f2. `sc` holds the surrogates in that order. Then:
```julia_repl
julia> _output_model_mapping(sc,1)
(1,1)
julia> _output_model_mapping(sc,2)
(2,1)
julia> _output_model_mapping(sc,3)
(1,2)
```
"""
_output_model_mapping( sc :: AbstractSurrogateContainer, l :: Int ) :: Tuple{Int, Int} = nothing

"""
	_sc_output_model_mapping(sc, ℓ)

Return a `Tuple{Int,Int}` with the index `i` of the surrogate model and its output `k`
describing output `ℓ` of `sc`, i.e., as returned by `eval_models(sc, x)`.
"""
_sc_output_model_mapping( sc :: AbstractSurrogateContainer, ℓ :: Int ) :: Tuple{Int,Int} = nothing

## Derived

# get_wrapper can be improved
get_wrapper(sc :: AbstractSurrogateContainer, i :: Int) = list_of_wrappers(sc)[i]

function num_outputs( sc :: AbstractSurrogateContainer ) :: Int
    return sum( num_outputs(sw) for sw ∈ list_of_wrappers(sc) )
end

function fully_linear( sc :: AbstractSurrogateContainer )
    return all( fully_linear(get_model(sw)) for sw ∈ list_of_wrappers(sc) )
end

"Evaluate all surrogate models stored in `sc` at scaled site `x̂`."
function eval_models( sc :: AbstractSurrogateContainer, x̂ :: Vec )
	return vcat( (eval_models( get_model(sw) , x̂) for sw ∈ list_of_wrappers(sc) )...)
end

"Evaluate Jacobian of surrogate models stored in `sc` at scaled site `x̂`."
function get_jacobian( sc :: AbstractSurrogateContainer, x̂ :: Vec )
    model_jacobians = [ get_jacobian( get_model(sw), x̂) for sw ∈ list_of_wrappers(sc) ]
    return vcat( model_jacobians... )
end

@doc """
Return a function handle to be used with `NLopt` for output `l` of `sc`.
Index `l` is assumed to be an *internal* index in the range of `1, …, n_objfs`,
where `n_objfs` is the total number of (scalarized) objectives modelled by `sc`.
"""
function get_optim_handle( sc :: AbstractSurrogateContainer, l :: Int )
    i, ℓ = sc_output_model_mapping(sc, l)
    sw = get_wrapper( sc, i )
    return _get_optim_handle( get_model(sw), ℓ )
end

function get_saveable_type( sc :: AbstractSurrogateContainer )
	return Tuple{ (get_saveable_type( get_meta(sw) ) for sw in list_of_wrappers(sc))... }
end

get_saveable( :: Nothing ) = nothing
get_saveable_type( :: Nothing ) = Nothing

function get_saveable( sc :: AbstractSurrogateContainer )
	return Tuple( get_saveable( get_meta(sw) ) for sw in list_of_wrappers(sc) )
end

# These are used for PS constraints
# one can in theory provide vector constraints but most solvers fail then
@doc """
    eval_models( sc, x̂, l )

Return model value for output `l` of `sc` at `x̂`.
Index `l` is assumed to be an *internal* index in the range of 1,…,n_objfs,
where n_objfs is the total number of (scalarized) objectives stored in `sc`.
"""
function eval_models( sc :: AbstractSurrogateContainer, x̂ :: Vec, l :: Int )
    i, ℓ = sc_output_model_mapping(sc,l)
    sw = get_wrapper(sc.surrogates, i )
    return eval_models( get_model(sw), x̂, ℓ )
end

@doc """
    get_gradient( sc, x̂, l )

Return a gradient for output `l` of `sc` at `x̂`.
Index `l` is assumed to be an internal index in the range of 1,…,n_objfs,
where n_objfs is the total number of (scalarized) objectives stored in `sc`.
"""
function get_gradient( sc :: AbstractSurrogateContainer, x̂ :: Vec, l :: Int )
    i, ℓ = sc_output_model_mapping( sc, l )
    sw = get_wrapper( sc, i )
    return get_gradient( get_model(sw), x̂, ℓ );
end


text:
module MatchCore
using MLStyle
using MLStyle.Err

export @sswitch, ellipsis_split, backend, P_partial_struct_decons
using MLStyle.AbstractPatterns
using MLStyle.AbstractPatterns.BasicPatterns

"""
[a, b..., c] -> :vec3 => [a], b, [c]
[a, b, c]    -> :vec => [a, b, c]
"""
function ellipsis_split(args::AbstractArray{T, 1}) where {T}
    ellipsis_index = findfirst(args) do arg
        Meta.isexpr(arg, :...)
    end
    if ellipsis_index === nothing
        Val(:vec) => args
    else
        Val(:vec3) => (
            args[1:ellipsis_index-1],
            args[ellipsis_index].args[1],
            args[ellipsis_index+1:end],
        )
    end
end

@nospecialize
function qt2ex(ex::Any)
    if ex isa Expr
        Meta.isexpr(ex, :$) && return ex.args[1]
        Expr(
            :call,
            Expr,
            QuoteNode(ex.head),
            Expr(:vect, (qt2ex(e) for e in ex.args if !(e isa LineNumberNode))...),
        )
    elseif ex isa Symbol
        QuoteNode(ex)
    else
        ex
    end
end

const backend = RedyFlavoured.backend

function P_partial_struct_decons(t, partial_fields, ps, prepr::AbstractString = "$t")
    function tcons(_...)
        t
    end
    comp = PComp(prepr, tcons;)
    function extract(sub, i::Int, ::Any, ::Any)
        :($sub.$(partial_fields[i]))
    end
    decons(comp, extract, ps)
end

basic_ex2tf(eval::Function, a) =
    isprimitivetype(typeof(a)) ? literal(a) : error("invalid literal $a")
basic_ex2tf(eval::Function, l::LineNumberNode) = wildcard
basic_ex2tf(eval::Function, q::QuoteNode) = literal(q.value)
basic_ex2tf(eval::Function, s::String) = literal(s)
basic_ex2tf(eval::Function, n::Symbol) = n === :_ ? wildcard : P_capture(n)

function basic_ex2tf(eval::Function, ex::Expr)
    !(x) = basic_ex2tf(eval, x)
    hd = ex.head
    args = ex.args
    n_args = length(args)
    if hd === :||
        @assert n_args === 2
        l, r = args
        or(!l, !r)
    elseif hd === :&&
        @assert n_args === 2
        l, r = args
        and(!l, !r)

    elseif hd === :if
        @assert n_args === 2
        let cond = args[1]
            guard() do _, scope, _
                see_captured_vars(cond, scope)
            end
        end
    elseif hd === :&
        @assert n_args === 1
        val = args[1]
        guard() do target, scope, _
            see_captured_vars(:($target == $val), scope)
        end
    elseif hd === :let
        bind = args[1]
        @assert bind isa Expr
        if bind.head === :(=)
            @assert bind.args[1] isa Symbol
            P_bind(bind.args[1], bind.args[2], see_capture = true)
        else
            @assert bind.head === :block
            binds = Function[
                P_bind(arg.args[1], arg.args[2], see_capture = true) for arg in bind.args
            ]
            push!(binds, wildcard)
            and(binds)
        end
    elseif hd === :(::)
        if n_args === 2
            p, ty = args
            ty = eval(ty)::TypeObject
            and(P_type_of(ty), !p)
        else
            @assert n_args === 1
            ty = args[1]
            ty = eval(ty)::TypeObject
            P_type_of(ty)
        end
    elseif hd === :vect
        tag, split = ellipsis_split(args)
        return tag isa Val{:vec} ? P_vector([!e for e in split]) :
               let (init, mid, tail) = split
            P_vector3([!e for e in init], !mid, [!e for e in tail])
        end
    elseif hd === :tuple
        P_tuple([!e for e in args])
    elseif hd === :call
        let f = args[1], args′ = view(args, 2:length(args))
            n_args′ = n_args - 1
            t = eval(f)
            if t === Core.svec
                tag, split = ellipsis_split(args′)
                return tag isa Val{:vec} ? P_svec([!e for e in split]) :
                       let (init, mid, tail) = split
                    P_svec3([!e for e in init], !mid, [!e for e in tail])
                end
            end
            all_field_ns = fieldnames(t)
            partial_ns = Symbol[]
            patterns = Function[]
            if n_args′ >= 1 && Meta.isexpr(args′[1], :parameters)
                kwargs = args′[1].args
                args′ = view(args′, 2:length(args′))
            else
                kwargs = []
            end
            if length(all_field_ns) === length(args′)
                append!(patterns, [!e for e in args′])
                append!(partial_ns, all_field_ns)
            elseif length(partial_ns) !== 0
                error(
                    "count of positional fields should be 0 or the same as the fields($all_field_ns)",
                )
            end
            for e in kwargs
                if e isa Symbol
                    e in all_field_ns ||
                        error("unknown field name $e for $t when field punnning.")
                    push!(partial_ns, e)
                    push!(patterns, P_capture(e))
                elseif Meta.isexpr(e, :kw)
                    key, value = e.args
                    key in all_field_ns ||
                        error("unknown field name $key for $t when field punnning.")
                    @assert key isa Symbol
                    push!(partial_ns, key)
                    push!(patterns, and(P_capture(key), !value))
                end
            end
            P_partial_struct_decons(t, partial_ns, patterns)
        end
    elseif hd === :quote
        !qt2ex(args[1])
    else
        error("not implemented expr=>pattern rule for '($hd)' Expr.")
    end
end

const case_sym = Symbol("@case")
"""a minimal implementation of sswitch
this is incomplete and only for bootstrapping, do not use it.
"""
macro sswitch(val, ex)
    @assert Meta.isexpr(ex, :block)
    branches = Pair{Function, Tuple{LineNumberNode, Int}}[]
    k = 0
    ln = __source__
    terminal = Dict{Int, Any}()
    body = nothing
    for i in eachindex(ex.args)
        stmt = ex.args[i]
        if Meta.isexpr(stmt, :macrocall) &&
           stmt.args[1] === case_sym &&
           length(stmt.args) == 3
            k += 1
            pattern = try
                basic_ex2tf(__module__.eval, stmt.args[3])
            catch e
                throw(PatternCompilationError(ln, e))
            end
            push!(branches, (pattern => (ln, k)))
            body = terminal[k] = Expr(:block)
        else
            stmt isa LineNumberNode && (ln = stmt)
            k === 0 || push!(body.args, stmt)
        end
    end

    block = backend(val, branches, terminal, __source__; hygienic = false)
    esc(init_cfg(block))
end
@specialize

end # module end


text:
using RedPitayaDAQServer
using PyPlot

# obtain the URL of the RedPitaya
include("config.jl")

rp = RedPitayaCluster(URLs[1:2:3])

dec = 8
modulus = 4800
base_frequency = 125000000
periods_per_step = 10
samples_per_period = div(modulus, dec)#*periods_per_step
periods_per_frame = 100 # about 0.5 s frame length
frame_period = dec*samples_per_period*periods_per_frame / base_frequency

slow_dac_steps_per_frame = div(periods_per_frame, periods_per_step)

stopADC(rp)
masterTrigger(master(rp), false)

@info samples_per_period, periods_per_frame, slow_dac_steps_per_frame
decimation(rp, dec)
samplesPerPeriod(rp, samples_per_period)
periodsPerFrame(rp, periods_per_frame)
passPDMToFastDAC(master(rp), true)
slowDACStepsPerFrame(rp, slow_dac_steps_per_frame)

numSlowDACChan(master(rp), 2)
lutA = collect(range(0,0.3,length=slow_dac_steps_per_frame))
lutB = collect(ones(slow_dac_steps_per_frame))

lutEnableDACA = ones(Bool, slow_dac_steps_per_frame)
lutEnableDACA[1:2:end] .= false
lutEnableDACB = ones(Bool, slow_dac_steps_per_frame)
lutEnableDACB[2:2:end] .= false

modeDAC(rp, "STANDARD")
frequencyDAC(rp,1,1, base_frequency / modulus)

freq = frequencyDAC(rp,1,1)
println(" frequency = $(freq)")
jumpSharpnessDAC(rp, 1, 0.01) # controls the sharpness of the jump for the square
signalTypeDAC(rp, 1 , "SQUARE")
phaseDAC(rp, 1, 1, 0.0 ) # Phase has to be given in between 0 and 1

ramWriterMode(rp, "TRIGGERED")
triggerMode(rp, "EXTERNAL")

amplitudeDAC(rp, 1, 1, 0.1)
lut = collect(cat(lutA,lutB*0.1,dims=2)')
enableLUT = collect( cat(lutEnableDACA,lutEnableDACB,dims=2)' )
seq = ArbitrarySequence(lut, enableLUT, slow_dac_steps_per_frame, 1, 0.0, 0.0)
appendSequence(master(rp), seq)
success = prepareSequence(master(rp))

startADC(rp)
masterTrigger(rp, true)

uCurrentPeriod = readData(rp, 0, 1)
stopADC(rp)
masterTrigger(rp, false)

lut = collect(cat(lutA,lutB*0.4,dims=2)')
amplitudeDAC(rp, 1, 1, 0.1)
lut = collect(cat(lutA,lutB*0.4,dims=2)')
enableLUT = collect( cat(lutEnableDACA,lutEnableDACB,dims=2)' )
seq = ArbitrarySequence(lut, enableLUT, slow_dac_periods_per_frame, 1, 0.0, 0.0)
appendSequence(master(rp), seq)
success = prepareSequence(master(rp))

startADC(rp)
masterTrigger(rp, true)

uCurrentPeriod2 = readData(rp, 0, 1)

fig = figure(1)
clf()
subplot(2,1,1)
plot(vec(uCurrentPeriod[:,1,:,:]))
plot(vec(uCurrentPeriod[:,2,:,:]))
legend(("Rx1", "Rx2"))
subplot(2,1,2)
plot(vec(uCurrentPeriod2[:,1,:,:]))
plot(vec(uCurrentPeriod2[:,2,:,:]))
legend(("Rx1", "Rx2"))

savefig("images/slowDACMultiChannelWithOffset.png")

stopADC(rp)
masterTrigger(rp, false)

max_stars_repo_name: tknopp/RedPitayaDAQServer

max_issues_repo_name: tknopp/RedPitayaDAQServer

max_forks_repo_name: tknopp/RedPitayaDAQServer


text:
include("intervals.jl")
include("functions.jl")


text:
module Resampling

using Juno
using PyCall
using StatsBase
using DSB2018
using DSB2018.Data
using DSB2018.Model
using ColorTypes
using Images

import StatsBase: predict
@pyimport skimage.restoration as skimage_restoration
@pyimport skimage.morphology as skimage_morphology
@pyimport scipy.ndimage as ndi
@pyimport skimage.filters as skimage_filters

function re_sample(img::Matrix{Float64}, n::Int64)
    jx = collect(1 : length(img))
    w = StatsBase.ProbabilityWeights(img[:])

    out = zeros(length(jx))
    Juno.@progress "Sampling" for k = 1 : n
        x = sample(jx, w, (length(jx), ))
        out[x] .+= 1.
    end

    reshape(out, size(img))/n
end

function re_sample(img::DSB2018.Data.Image, n::Int64)
    return re_sample(Float64.(Gray.(img.image)), n)
end

function is_bi_modal(x::X) where {X <: AbstractArray}
    dx = diff(x)
    return sum(sign.(dx[1:(end-1)]) .!= sign.(dx[2:end])) == 3
end

function bi_modal_picks(x::X) where {X <: AbstractArray}
    pick1 = findfirst(diff(x) .< 0.0)
    pick2 = findlast(diff(x) .> 0.0)
    return pick1, pick2
end

end