text:
module ds_DistributedBoosting

using RCall

export Unibeta, Covarmat, DSLogin, Boostscratch
export ds_loadPkg, ds_login, ds_check, ds_start, ds_logout
export calc_covarmat, calc_covarmat!, calc_unibeta
export boost!, reboost!
export ds_boosting
export selectionofcovs, getselections, getcovarfromtriangular

"""
	Unibeta(unibeta, labels)

Stucture Unibeta, containing effect estimates and corresponding lables.

#Arguments:
- `unibeta::Array{Float64,1}`: univariable effect estimates
- `labels::Array{String,1}`: labels of univariable effect estimates (should be of same length as unibeta and the same order)
"""
mutable struct Unibeta
	unibeta::Array{Float64,1}
	labels::Array{String,1}
end

"""
	Covarmat(covarmat, labels)

Stucture Covarmat, containing 2-dimensional covariance matrix and corresponding lables.

#Arguments:
- `covarmat::Array{Float64,2}`: 2-dimensional covariance matrix (should be symmetrical)
- `labels::Array{String,1}`: labels of 2-dimensional covairance matrix (should be of same length as one-dim of covarmat and of the same order)
"""
mutable struct Covarmat
	covarmat::Array{Float64,2}
	labels::Array{String,1}
end

"""
	DSLogin(url, user, password, table, servernames)

Stucture DSLogin, containing information to loginto DataShield.

#Arguments:
- `url::Array{String,1}`: Vector with urls to DataSHIELD Servers
- `user::Array{String,1}`: Username
- `password::Array{String,1}`: Password or name of key
- `table::Array{String,1}`: Table name to be loaded
- `servernames::Array{String,1}`: Vector with servernames (personal choice)
"""
mutable struct DSLogin
	url::Array{String,1}
	user::String
	password::String
	table::Array{String,1}
	servernames::Array{String,1}
end

"""
	Boostscratch(actualstepno, actualbeta, actualscore,
					pooledunibeta, pooledcovarmat, actualupdates,
					selections, wantedlabels, usedlabels, nu,
					stepno,	actualsel, actualnom, a,
					x, maxvar, y, labels)

Stucture Boostscratch, containing all relevant information of the ds_DistributedBoosting algorithm.

#Arguments:
- `actualstepno::Int`: Actual stepno in boosting algorithm
- `actualbeta::Array{Float64,1}`: Actual effect estimate vector (same size as labels)
- `actualscore::Array{Float64,1}`: Actual score function value vector (same size as labels)
- `pooledunibeta::Unibeta`: Pooled univariable effect estiamtes (same size as labels)
- `pooledcovarmat::Covarmat`: Actual pooled covariance matrix (expanded in the course of the algorithm)
- `actualupdates::Array{Float64,1}`: Actual update size (before multiplication with nu)
- `selections::Array{String,1}`: Vector of selected variables (same size as stepno, variables can be in the vector more than once)
- `wantedlabels::Array{String,1}`: Vector containing variable names of the ones for which new covariances need to be called
- `usedlabels::Array{String,1}`: Vector containing variable names of the ones where already covariances are called
- `nu::Float64`: Shrinkage parameter (Experience: 0.1)
- `stepno::Int`: Number of boosting steps to perform maximally
- `actualsel::Int`: Position of last selected variable
- `actualnom::Array{Float64,1}`: Actual denominator for boosting
- `a::Int`: Number of variables for which covariances should be called before the first boosting step (Experience: 20)
- `x::Int`: Number of variables for which covariances should be called additionally to the ones selected by the heuristic (Experience: 20)
- `maxvar::Int`: Number of variables to be selected maximally
- `y::String`: Name of endpoint variable
- `labels::Array{String,1}`: Vector containing names of potential predictors to be considered
"""
mutable struct Boostscratch
	actualstepno::Int
	actualbeta::Array{Float64,1}
	actualscore::Array{Float64,1}
	pooledunibeta::Unibeta
	pooledcovarmat::Covarmat
	actualupdates::Array{Float64,1}
	selections::Array{String,1}
	wantedlabels::Array{String,1}
	usedlabels::Array{String,1}
	nu::Float64
	stepno::Int
	actualsel::Int
	actualnom::Array{Float64,1}
	a::Int
	x::Int
	maxvar::Int
	y::String
	labels::Array{String,1}
	login::DSLogin
end

include("Session.jl")
include("calc.jl")
include("getselections.jl")
include("boost.jl")
include("boosting.jl")

end


text:
# activation structs
struct ActivationRamp
    t0::Float64     # onset time relative to end of last ramp
    r::Float64      # ramp time
    a1::Float64     # new activation level
end

struct ActivationProfile
    a0::Float64     # starting activation level
    ramps::Vector{ActivationRamp}
end
ActivationProfile(a0, r::ActivationRamp) = ActivationProfile(a0, [r])
ActivationProfile(a0, ramps...) = ActivationProfile(a0, [ActivationRamp(ramps[i:i+2]...) for i = 1:3:length(ramps)])
ActivationProfile(x::AbstractArray) = ActivationProfile(x[1], x[2:end]...)

# ramp function - time relative to onset time (t-t0)
ramp(t, a0, r, a1) = a0 + (a1 - a0) * (t / r)^3 * (6 * (t / r)^2 - 15 * (t / r) + 10)
ramp(t, a0, p::ActivationRamp) = ramp(t, a0, p.r, p.a1)

# calculate activation from time - t relative to end of last ramp
activation(t, a0, p::ActivationRamp) = t < p.t0 ? a0 : t ≤ p.t0 + p.r ? ramp(t - p.t0, a0, p) : p.a1

function activation(t, p::ActivationProfile)
    a0 = p.a0   # previous activation level
    _t = 0.0    # end of last ramp time
    for ramp in p.ramps
        if t < _t + ramp.t0 + ramp.r
            return activation(t - _t, a0, ramp)
        else
            _t += ramp.t0 + ramp.r
            a0 = ramp.a1
        end
    end

    return p.ramps[end].a1
end

# find time activation goes on for
time(p::ActivationProfile) = [r.t0 + r.r for r in p.ramps] |> sum

# custom plotting
Plots.plot(t, p::ActivationProfile) = plot(x -> activation(x, p), t)
Plots.plot(p::ActivationProfile) = plot(t -> activation(t, p), 0.0, time(p))

# convert ActivationProfile to vec
Base.Vector(a::ActivationProfile) = [a.a0, vcat([[b.t0, b.r, b.a1] for b in a.ramps]...)...]

ext: jl


text:
using Includeonce
using Test


text:
module PathTracing

# The following started out as a port of smallpt
# https://www.kevinbeason.com/smallpt
#
# Unlike the original, we're not restricting ourselves to 99 lines - we're
# permitting ourselves to expand and generalize a bit for the sake of clarity.

using StaticArrays, LinearAlgebra, PNGFiles, Random, Colors

const Vec = SVector{3,Float64}
const V = SA{Float64}
const C = SA{Float32}
const CVec = SVector{3,Float32}

struct Ray
    o::Vec # Origin
    d::Vec # Direction (normalized)
end

@enum Refl_t DIFF SPEC REFR

struct Sphere
    rad::Float64
    p::Vec  # Position
    e::CVec # Emission
    c::CVec # Color
    refl::Refl_t
end

# Returns distance, 0 if nohit
function intersect(s::Sphere, r::Ray)::Float64
    # Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0
    op = s.p - r.o
    b = dot(op, r.d)
    det = b^2 - dot(op,op) + s.rad^2
    det < 0 && return 0
    eps = 1e-4
    det = sqrt(det)
    t = b - det
    t > eps && return t
    t = b + det
    t > eps && return t
    return 0
end

function intersect(scene, r::Ray)
    inf = 1e20
    t = inf
    id = 0
    for i = 1:length(scene)
        d = intersect(scene[i], r)
        if d != 0 && d < t
            t = d
            id = i
        end
    end
    (t < inf, t, id)
end

clamp01(x) = clamp(x, zero(x), one(x))

function radiance(scene, r::Ray, depth::Int, Xi)
    (did_hit, t, id) = intersect(scene, r)
    did_hit || return C[0,0,0]      # if miss, return black
    obj = scene[id]                 # the hit object
    x = r.o + r.d*t
    n = normalize(x-obj.p)
    nl = dot(n, r.d) < 0 ? n : -n
    f = obj.c
    p = maximum(f) # max refl
    depth += 1
    if depth > 5
        if rand(Xi) < p
            f = f*(1/p)
        else
            return obj.e # R.R.
        end
    end
    if obj.refl == DIFF     # Ideal DIFFUSE reflection
        r1 = 2*π*rand(Xi)
        r2 = rand(Xi)
        r2s = sqrt(r2)
        w = nl
        u = normalize((abs(w[1]) > 0.1 ? V[0,1,0] : V[1,0,0]) × w)
        v = w × u
        d = normalize(u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2))
        return obj.e + f .* radiance(scene, Ray(x,d), depth, Xi)
    elseif obj.refl == SPEC # Ideal SPECULAR reflection
        return obj.e + f .* radiance(scene, Ray(x, r.d - n*2*dot(n,r.d)), depth, Xi)
    else # obj.refl == REFR # Ideal dielectric REFRACTION
        reflRay = Ray(x, r.d - n*2*dot(n,r.d))
        into = dot(n, nl) > 0  # Ray from outside going in?
        nc = 1
        nt = 1.5
        nnt = into ? nc/nt : nt/nc
        ddn = dot(r.d, nl)
        cos2t = 1 - nnt^2 * (1 - ddn^2)
        if cos2t < 0 # Total internal reflection
            return obj.e + f .* radiance(scene, reflRay, depth, Xi)
        end
        tdir = normalize(r.d*nnt - n*((into ? 1 : -1)*(ddn*nnt+sqrt(cos2t))))
        a = nt - nc
        b = nt + nc
        R0 = a^2/b^2
        c = 1 - (into ? -ddn : dot(tdir, n))
        Re = R0 + (1-R0)*c^5
        Tr = 1-Re
        P = 0.25 + 0.5*Re
        RP = Re/P
        TP = Tr/(1-P)
        return obj.e + f .* (depth>2 ?
            (rand(Xi)<P ? radiance(scene, reflRay,depth,Xi)*RP : radiance(scene, Ray(x,tdir),depth,Xi)*TP) :
            radiance(scene, reflRay,depth,Xi)*Re + radiance(scene, Ray(x,tdir),depth,Xi)*Tr
        )
    end
end

function render(scene, w, h, samps)
    subpix_samps = samps ÷ 4
    cam = Ray(V[50,52,295.6], normalize(V[0,-0.042612,-1])) # cam pos, dir
    cx = 0.5135 * V[w/h, 0, 0]
    cy = 0.5135 * normalize(cx × cam.d)
    c = zeros(CVec, w, h)
    Xi = Random.GLOBAL_RNG # TODO: row-based seeding for multithreading
    rowcount = Threads.Atomic{Int}(0)
    Threads.@threads for y = 1:h # Loop over image rows
        rowcount[] += 1
        #@info "Rendering ($samps spp)" progress=(rowcount[])/h
        for x = 1:w # Loop cols
            # Xi = Random.MersenneTwister(y^3)
            for sy = 0:1      # 2x2 subpixel rows
                for sx = 0:1  # 2x2 subpixel cols
                    r = V[0,0,0]
                    for _ = 1:subpix_samps
                        # Bilinear (tent) filter for subpixels
                        r1 = 2*rand(Xi)
                        r2 = 2*rand(Xi)
                        dx = r1 < 1 ? sqrt(r1)-1 : 1-sqrt(2-r1)
                        dy = r2 < 1 ? sqrt(r2)-1 : 1-sqrt(2-r2)
                        d = normalize(cx*(((sx + 0.5 + dx)/2 + x-1)/w - 0.5) +
                                      cy*(((sy + 0.5 + dy)/2 + y-1)/h - 0.5) + cam.d)
                        r += radiance(scene, Ray(cam.o + d*140, d), 0, Xi)
                        # Camera rays are pushed ^^^^^ forward to start in interior
                    end
                    # The original includes clamp01 here (bias↑, noise↓)
                    c[x,y] += 0.25 * (1.00/subpix_samps) .* r
                end
            end
        end
    end
    c
end

to_sRGB(v) = RGB(clamp01.(v).^(1/2.2)...)

function main()
    scene = cornell_box_scene()
    c = render(scene, 128, 96, 128)
    #c = render(scene, 512, 384, 4)
    #c = render(scene, 1024, 768, 4)
    #c = render(scene, 2, 2, 4)
    img = to_sRGB.(permutedims(reverse(c, dims=2)))
    # PNGFiles.save("image.png", img, filters=0)
    img
end

function cornell_box_scene()
    # radius, position, emission, color, material
    [
      Sphere(1e5,  V[ 1e5+1,40.8,81.6],  C[0,0,0],    C[.75,.25,.25], DIFF), #Left
      Sphere(1e5,  V[-1e5+99,40.8,81.6], C[0,0,0],    C[.25,.25,.75], DIFF), #Rght
      Sphere(1e5,  V[50,40.8, 1e5],      C[0,0,0],    C[.75,.75,.75], DIFF), #Back
      Sphere(1e5,  V[50,40.8,-1e5+170],  C[0,0,0],    C[0,0,0],       DIFF), #Frnt
      Sphere(1e5,  V[50, 1e5, 81.6],     C[0,0,0],    C[.75,.75,.75], DIFF), #Botm
      Sphere(1e5,  V[50,-1e5+81.6,81.6], C[0,0,0],    C[.75,.75,.75], DIFF), #Top
      Sphere(16.5, V[27,16.5,47],        C[0,0,0],    C[1,1,1]*.999,  SPEC), #Mirr
      Sphere(16.5, V[73,16.5,78],        C[0,0,0],    C[1,1,1]*.999,  REFR), #Glas
      Sphere(600,  V[50,681.6-.27,81.6], C[12,12,12], C[0,0,0],       DIFF), #Lite
      # Sphere( 4.5, V[50,30,62],          C[12,12,12],    C[0,0,0],       DIFF), #Light
    ]
end

end


text:
module BlockTrees

import FMMTrees

struct BlockTree{T}
    test_cluster::T
    trial_cluster::T
end

testcluster(blocktree::BlockTree) = blocktree.test_cluster
trialcluster(blocktree::BlockTree) = blocktree.trial_cluster

FMMTrees.root(t::BlockTree) = (FMMTrees.root(testcluster(t)), FMMTrees.root(trialcluster(t)))
function FMMTrees.data(t::BlockTree, n)
    (
        FMMTrees.data(testcluster(t), n[1]),
        FMMTrees.data(trialcluster(t), n[2]))
end

function FMMTrees.children(b::BlockTree, node)
    test_chds = FMMTrees.children(testcluster(b), node[1])
    trial_chds = FMMTrees.children(testcluster(b), node[2])
    ((ch[1],ch[2]) for ch in Iterators.product(test_chds, trial_chds))
end

function FMMTrees.haschildren(b::BlockTree, node)
    !FMMTrees.haschildren(testcluster(b), node[1]) && return false
    !FMMTrees.haschildren(trialcluster(b), node[2]) && return false
    return true
end

function FMMTrees.LevelledTrees.numlevels(bt)
    num_test_levels = FMMTrees.LevelledTrees.numlevels(testcluster(bt))
    num_trial_levels = FMMTrees.LevelledTrees.numlevels(trialcluster(bt))
    @assert num_test_levels == num_trial_levels
    return num_test_levels
end

end # module BlockTrees


text:
@testset "preconditioned Crank Nicolson tests" begin

using DiffEqNoiseProcess, Test, Random
using Statistics
using DiffEqBase
using DiffEqBase.EnsembleAnalysis

##
# Tests of pCN
##

W = WienerProcess(0.0,0.0,0.0)

dt = 0.1
calculate_step!(W,dt,nothing,nothing)

for i in 1:10
  accept_step!(W,dt,nothing,nothing)
end

_W = deepcopy(W)

# test with ρ=1
W2 = pCN!(_W, 1.0)
W2a = pCN(W, 1.0)
WWrapper = NoiseWrapper(_W)

# test source
@test W2.source.W == WWrapper.source.W
@test W2.source.Z == WWrapper.source.Z
@test W2.source.W == W.W
@test W2.source.u == W.u
@test W2.source.Z == W.Z
@test W2.source.W == W2a.source.W
@test W2.source.Z == W2a.source.Z

# multi dimensional
W = WienerProcess(0.0,zeros(4),zeros(4))

dt = 0.1
calculate_step!(W,dt,nothing,nothing)

for i in 1:10
  accept_step!(W,dt,nothing,nothing)
end

_W = deepcopy(W)

W2 = pCN!(_W, 1.0)
W2a = pCN(W, 1.0)
WWrapper = NoiseWrapper(_W)

# test source
@test W2.source.W == WWrapper.source.W
@test W2.source.Z == WWrapper.source.Z
@test W2.source.W == _W.W
@test W2.source.u == _W.u
@test W2.source.Z == _W.Z
@test W2.source.W == W.W
@test W2.source.Z == W.Z
@test W2.source.W == W2a.source.W
@test W2.source.Z == W2a.source.Z

# test ρ!=0 and ρ!=1
_W = deepcopy(W)
W3 = pCN!(_W, 0.2)
W3a = pCN(W, 0.2)

# test source

@test W3.source.W != W.W
@test W3.source.u != W.u
@test W3.source.Z == W.Z # no action on auxilary process
@test W3.source.W == W3a.source.W
@test W3.source.u == W3a.source.u
@test W3.source.Z == W3a.source.Z

# inplace
W = WienerProcess!(0.0,zeros(4),zeros(4))

dt = 0.1
calculate_step!(W,dt,nothing,nothing)

for i in 1:10
  accept_step!(W,dt,nothing,nothing)
end

_W = deepcopy(W)

# test with ρ=0
W2 = pCN!(_W, 0.0)
W2a = pCN(W, 0.0)
WWrapper = NoiseWrapper(_W)

# test source
@test W2.source.W != W.W
@test W2.source.Z == W.Z
@test W2.source.W == W2a.source.W
@test W2.source.Z == W2a.source.Z

# statistics test
W = WienerProcess(0.0,0.0,0.0)
prob = NoiseProblem(W,(0.0,1.0))
sol = solve(prob,dt=0.1)

function prob_func(prob,i,repeat)
  _sol = deepcopy(sol)
  Wtmp = pCN!(_sol,1.0)
  remake(prob, noise=Wtmp)
end

ensemble_prob = EnsembleProblem(prob, prob_func=prob_func)
@time sim = solve(ensemble_prob,dt=0.1,trajectories=100)

# Spot check the mean and the variance
qs = 0:0.1:1
for i in 1:11
  q = qs[i]
  @test ≈(timestep_mean(sim,i),sol.W[i],atol=1e-2)
  @test ≈(timestep_meanvar(sim,i)[2],zero(sol.W[i]),atol=1e-2)
end

# ρ = 0.5 test
ρ = 0.5
W = WienerProcess(0.0,0.0,0.0)
prob = NoiseProblem(W,(0.0,1000.0))
sol = solve(prob,dt=0.001)

prob2 = NoiseProblem(pCN(sol, ρ),(0.0,1000.0))
sol2 = solve(prob2,dt=0.001)

computedW(sim, indx) = (sim.u[indx+1]-sim.u[indx])/sqrt(sim.t[indx+1]-sim.t[indx])
dWnew = []
dWold = []
for i in 1:length(sol.t[1:end-1])
  push!(dWnew,computedW(sol,i))
  push!(dWold,computedW(sol2,i))
end
@show cor(dWnew,dWold)

@test ≈(cor(dWnew,dWold),ρ,atol=1e-2)

# Noise Grid tests

trange = 0.0:0.001:1000.0

# scalar
Ws = cumsum([0.0;[sqrt(trange[i+1]-ti)*randn()
        for (i,ti) in enumerate(trange[1:end-1])]])
NG = NoiseGrid(trange,Ws)
NG2 = pCN(NG,1.0) # ρ = 1.0
NG3 = pCN(NG,0.5) # ρ = 0.5

@test ≈(cor(diff(NG.W)./sqrt.(diff(NG.t)), diff(NG2.W)./sqrt.(diff(NG2.t))),1,atol=1e-2)
@test ≈(cor(diff(NG.W)./sqrt.(diff(NG.t)), diff(NG3.W)./sqrt.(diff(NG3.t))),0.5,atol=1e-2)

# array
dim = 4
Ws = cumsum([[zeros(dim)];[sqrt(trange[i+1]-ti)*randn(dim)
        for (i,ti) in enumerate(trange[1:end-1])]])
NG = NoiseGrid(trange,Ws)
NG1 = pCN(NG,1.0) # ρ = 1.0
NG2 = pCN(NG,0.5) # ρ = 0.5

for i=1:dim
  @test ≈(cor(getindex.(diff(NG.W)./sqrt.(diff(NG.t)),i),getindex.(diff(NG1.W)./sqrt.(diff(NG1.t)),i)),1.0,atol=1e-2)
  @test ≈(cor(getindex.(diff(NG.W)./sqrt.(diff(NG.t)),i),getindex.(diff(NG2.W)./sqrt.(diff(NG2.t)),i)),0.5,atol=1e-2)
end

end


text:
"""
    AbstractJuliaType
"""
abstract type AbstractJuliaType end

struct JuliaUnknown{T<:CLType} <: AbstractJuliaType
    x::T
end

struct JuliaUnsupported <: AbstractJuliaType end

struct JuliaCpointer <: AbstractJuliaType
    ref::CLPointer
end
struct JuliaCconstarray <: AbstractJuliaType
    ref::CLConstantArray
end
struct JuliaCincompletearray <: AbstractJuliaType
    ref::CLIncompleteArray
end

struct JuliaCrecord{T<:CLCursor} <: AbstractJuliaType
    sym::Symbol
    cursor::T
end
JuliaCrecord() = JuliaCrecord(Symbol())
JuliaCrecord(x::AbstractString) = JuliaCrecord(Symbol(x), getNullCursor())

struct JuliaCenum{T<:CLCursor} <: AbstractJuliaType
    sym::Symbol
    cursor::T
end
JuliaCenum() = JuliaCenum(Symbol())
JuliaCenum(x::AbstractString) = JuliaCenum(Symbol(x), getNullCursor())

struct JuliaCtypedef <: AbstractJuliaType
    sym::Symbol
end
JuliaCtypedef() = JuliaCtypedef(Symbol())
JuliaCtypedef(x::AbstractString) = JuliaCtypedef(Symbol(x))

struct JuliaCfunction <: AbstractJuliaType
    sym::Symbol
end
JuliaCfunction() = JuliaCfunction(Symbol())
JuliaCfunction(x::AbstractString) = JuliaCfunction(Symbol(x))

"""
    System Independent Types
"""
abstract type AbstractJuliaSIT <: AbstractJuliaType end

struct JuliaCbool <: AbstractJuliaSIT end
struct JuliaCuchar <: AbstractJuliaSIT end
struct JuliaCshort <: AbstractJuliaSIT end
struct JuliaCushort <: AbstractJuliaSIT end
struct JuliaCint <: AbstractJuliaSIT end
struct JuliaCuint <: AbstractJuliaSIT end
struct JuliaClonglong <: AbstractJuliaSIT end
struct JuliaCulonglong <: AbstractJuliaSIT end
struct JuliaCintmax_t <: AbstractJuliaSIT end
struct JuliaCuintmax_t <: AbstractJuliaSIT end
struct JuliaCfloat <: AbstractJuliaSIT end
struct JuliaCdouble <: AbstractJuliaSIT end
struct JuliaComplexF32 <: AbstractJuliaSIT end
struct JuliaComplexF64 <: AbstractJuliaSIT end
struct JuliaCptrdiff_t <: AbstractJuliaSIT end
struct JuliaCssize_t <: AbstractJuliaSIT end
struct JuliaCsize_t <: AbstractJuliaSIT end
struct JuliaCvoid <: AbstractJuliaSIT end
struct JuliaNoReturn <: AbstractJuliaSIT end
struct JuliaPtrCvoid <: AbstractJuliaSIT end
struct JuliaCstring <: AbstractJuliaSIT end
struct JuliaPtrUInt8 <: AbstractJuliaSIT end
struct JuliaPtrPtrUInt8 <: AbstractJuliaSIT end
struct JuliaAny <: AbstractJuliaSIT end
struct JuliaRefAny <: AbstractJuliaSIT end

struct JuliaCuint128 <: AbstractJuliaSIT end
struct JuliaCint128 <: AbstractJuliaSIT end
struct JuliaCschar <: AbstractJuliaSIT end
struct JuliaClongdouble <: AbstractJuliaSIT end
struct JuliaComplex <: AbstractJuliaSIT end
struct JuliaChalf <: AbstractJuliaSIT end
struct JuliaCfloat16 <: AbstractJuliaSIT end

struct JuliaCuint64_t <: AbstractJuliaSIT end
struct JuliaCuint32_t <: AbstractJuliaSIT end
struct JuliaCuint16_t <: AbstractJuliaSIT end
struct JuliaCuint8_t <: AbstractJuliaSIT end
struct JuliaCint64_t <: AbstractJuliaSIT end
struct JuliaCint32_t <: AbstractJuliaSIT end
struct JuliaCint16_t <: AbstractJuliaSIT end
struct JuliaCint8_t <: AbstractJuliaSIT end

struct JuliaCuintptr_t <: AbstractJuliaSIT end  # Csize_t
struct JuliaCtm <: AbstractJuliaSIT end  # Libc.TmStruct
struct JuliaCFILE <: AbstractJuliaSIT end  # Libc.FILE

"""
    System Dependent Types
"""
abstract type AbstractJuliaSDT <: AbstractJuliaType end

struct JuliaCchar <: AbstractJuliaSDT end
struct JuliaClong <: AbstractJuliaSDT end
struct JuliaCulong <: AbstractJuliaSDT end
struct JuliaCwchar_t <: AbstractJuliaSDT end

"""
    tojulia
Map a `CLType` type to a `AbstractJuliaType` type.
"""
function tojulia end

tojulia(x::CLType) = JuliaUnknown(x)
tojulia(x::CLVoid) = JuliaCvoid()
tojulia(x::CLBool) = JuliaCbool()
tojulia(x::CLChar_U) = JuliaCchar()
tojulia(x::CLUChar) = JuliaCuchar()
tojulia(x::CLChar16) = JuliaUnknown(x)  # C++11 char16_t
tojulia(x::CLChar32) = JuliaUnknown(x)  # C++11 char32_t
tojulia(x::CLUShort) = JuliaCushort()
tojulia(x::CLUInt) = JuliaCuint()
tojulia(x::CLULong) = JuliaCulong()
tojulia(x::CLULongLong) = JuliaCulonglong()
tojulia(x::CLUInt128) = JuliaCuint128()
tojulia(x::CLChar_S) = JuliaCchar()
tojulia(x::CLSChar) = JuliaCschar()  # Int8
tojulia(x::CLWChar) = JuliaCwchar_t()
tojulia(x::CLShort) = JuliaCshort()
tojulia(x::CLInt) = JuliaCint()
tojulia(x::CLLong) = JuliaClong()
tojulia(x::CLLongLong) = JuliaClonglong()
tojulia(x::CLInt128) = JuliaCint128()
tojulia(x::CLFloat) = JuliaCfloat()
tojulia(x::CLDouble) = JuliaCdouble()
tojulia(x::CLLongDouble) = JuliaClongdouble()  # depends on compiler
tojulia(x::CLComplex) = JuliaComplex()  # ComplexF32 or ComplexF64
tojulia(x::CLHalf) = JuliaChalf()  # Float16
tojulia(x::CLFloat16) = JuliaCfloat16()  # Float16
# tojulia(x::CLFloat128) = JuliaCfloat128() # see Quadmath.jl
tojulia(x::CLNullPtr) = JuliaUnknown(x)  # C++11 nullptr
tojulia(x::CLPointer) = JuliaCpointer(x)
tojulia(x::CLBlockPointer) = JuliaUnknown(x)  # ObjectveC's block pointer
function tojulia(x::CLElaborated)
    jlty = tojulia(getNamedType(x))
    @assert jlty isa JuliaCenum || jlty isa JuliaCrecord
    return jlty
end
tojulia(x::CLInvalid) = JuliaUnknown(x)
tojulia(x::CLFunctionProto) = JuliaCfunction(spelling(getTypeDeclaration(x)))
tojulia(x::CLFunctionNoProto) = JuliaCfunction(spelling(getTypeDeclaration(x)))
function tojulia(x::CLEnum)
    c = getTypeDeclaration(x)
    JuliaCenum(Symbol(spelling(c)), c)
end
function tojulia(x::CLRecord)
    c = getTypeDeclaration(x)
    JuliaCrecord(Symbol(spelling(c)), c)
end
tojulia(x::CLTypedef) = JuliaCtypedef(spelling(getTypeDeclaration(x)))
tojulia(x::CLUnexposed) = JuliaUnknown(x)
tojulia(x::CLConstantArray) = JuliaCconstarray(x)
tojulia(x::CLIncompleteArray) = JuliaCincompletearray(x)

# helper functions
is_jl_basic(x::AbstractJuliaType) = false
is_jl_basic(x::AbstractJuliaSIT) = true
is_jl_basic(x::AbstractJuliaSDT) = true

is_jl_unknown(x::AbstractJuliaType) = false
is_jl_unknown(x::JuliaUnknown) = true

is_jl_unsupported(x::AbstractJuliaType) = false
is_jl_unsupported(x::JuliaUnsupported) = true

is_jl_char(x::AbstractJuliaType) = false
is_jl_char(x::JuliaCchar) = true

is_jl_wchar(x::AbstractJuliaType) = false
is_jl_wchar(x::JuliaCwchar_t) = true

is_jl_pointer(x::AbstractJuliaType) = false
is_jl_pointer(x::JuliaCpointer) = true

is_jl_array(x::AbstractJuliaType) = false
is_jl_array(x::JuliaCconstarray) = true
is_jl_array(x::JuliaCincompletearray) = true

is_jl_funcptr(x::AbstractJuliaType) = false
is_jl_funcptr(x::CLPointer) = Clang.is_function(getPointeeType(getCanonicalType(x)))
is_jl_funcptr(x::CLConstantArray) = Clang.is_function(getElementType(getCanonicalType(x)))
is_jl_funcptr(x::CLIncompleteArray) = Clang.is_function(getElementType(getCanonicalType(x)))
is_jl_funcptr(x::JuliaCpointer) = is_jl_funcptr(x.ref)
is_jl_funcptr(x::JuliaCconstarray) = is_jl_funcptr(x.ref)
is_jl_funcptr(x::JuliaCincompletearray) = is_jl_funcptr(x.ref)

function get_jl_leaf_pointee_type(x::JuliaCpointer)
    is_jl_funcptr(x) && return JuliaPtrCvoid()
    jlptree = tojulia(getPointeeType(x.ref))
    return get_jl_leaf_type(jlptree)
end

function get_jl_leaf_eltype(x::Union{JuliaCconstarray,JuliaCincompletearray})
    is_jl_funcptr(x) && return JuliaPtrCvoid()
    jlelty = tojulia(getElementType(x.ref))
    return get_jl_leaf_type(jlelty)
end

get_jl_leaf_type(x::AbstractJuliaType) = x
get_jl_leaf_type(x::JuliaCpointer) = get_jl_leaf_pointee_type(x)
get_jl_leaf_type(x::JuliaCconstarray) = get_jl_leaf_eltype(x)
get_jl_leaf_type(x::JuliaCincompletearray) = get_jl_leaf_eltype(x)

is_jl_integer(x::AbstractJuliaType) = false
is_jl_integer(x::JuliaCshort) = true
is_jl_integer(x::JuliaCushort) = true
is_jl_integer(x::JuliaCint) = true
is_jl_integer(x::JuliaCuint) = true
is_jl_integer(x::JuliaClonglong) = true
is_jl_integer(x::JuliaCulonglong) = true
is_jl_integer(x::JuliaCintmax_t) = true
is_jl_integer(x::JuliaCuintmax_t) = true
is_jl_integer(x::JuliaCssize_t) = true
is_jl_integer(x::JuliaCsize_t) = true
is_jl_integer(x::JuliaCuint128) = true
is_jl_integer(x::JuliaCint128) = true
is_jl_integer(x::JuliaCuint64_t) = true
is_jl_integer(x::JuliaCuint32_t) = true
is_jl_integer(x::JuliaCuint16_t) = true
is_jl_integer(x::JuliaCuint8_t) = true
is_jl_integer(x::JuliaCint64_t) = true
is_jl_integer(x::JuliaCint32_t) = true
is_jl_integer(x::JuliaCint16_t) = true
is_jl_integer(x::JuliaCint8_t) = true
is_jl_integer(x::JuliaClong) = true
is_jl_integer(x::JuliaCulong) = true


text:
@testset ExtendedTestSet "squared_euclidean" begin
    @testset ExtendedTestSet "squared_euclidean" begin
        x = boolean_indicator([0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0])
        tfm = SquaredEuclidean(x)
        answer = [1, 0, 1, 4, 1, 0, 0, 0, 0, 0, 1]
        @test transform(x, tfm) == answer
    end

    @testset ExtendedTestSet "squared_euclidean" begin
        x = boolean_indicator(
            [
                0 1 1 1 0
                1 1 1 1 1
                1 0 0 0 1
                1 0 0 0 1
                1 0 0 0 1
                1 1 1 1 1
                0 1 1 1 0
            ]
        )
        tfm = SquaredEuclidean(x)
        answer = [
            1.0 0.0 0.0 0.0 1.0
            0.0 0.0 0.0 0.0 0.0
            0.0 1.0 1.0 1.0 0.0
            0.0 4.0 4.0 4.0 0.0
            0.0 1.0 1.0 1.0 0.0
            0.0 0.0 0.0 0.0 0.0
            1.0 0.0 0.0 0.0 1.0
        ]
        @test transform(x, tfm) == answer
    end

    @testset ExtendedTestSet "squared_euclidean" begin
        x = boolean_indicator(
            [
                0 1 1 1 0
                1 1 1 1 1
                1 0 0 0 1
                1 0 0 0 1
                1 0 0 0 1
                1 1 1 1 1
                0 1 1 1 0
            ]
        )
        tfm = SquaredEuclidean(x)
        answer = [
            1.0 0.0 0.0 0.0 1.0
            0.0 0.0 0.0 0.0 0.0
            0.0 1.0 1.0 1.0 0.0
            0.0 4.0 4.0 4.0 0.0
            0.0 1.0 1.0 1.0 0.0
            0.0 0.0 0.0 0.0 0.0
            1.0 0.0 0.0 0.0 1.0
        ]
        @test transform!(x, tfm, Threads.nthreads()) == answer
    end
end


text:
function writesol(filename,ifile,xdg)

if size(xdg,1)>0
    ncx = size(xdg,2);
else
    ncx = 0;
end
# if size(udg,1)>0
#     ncu = size(udg,2);
# else
#     ncu = 0;
# end
# if size(odg,1)>0
#     nco = size(odg,2);
# else
#     nco = 0;
# end
# if size(wdg,1)>0
#     ncw = size(wdg,2);
# else
#     ncw = 0;
# end

ndims = zeros(12,1);
ndims[1] = ncx;
# ndims[2] = ncu;
# ndims[3] = nco;
# ndims[4] = ncw;

nsize = zeros(20,1);
nsize[1] = length(ndims[:]);
nsize[2] = length(xdg[:]);
# nsize[3] = length(udg[:]);
# nsize[4] = length(odg[:]);
# nsize[5] = length(wdg[:]);

if (ifile>0)
    print("Writing solution into file " * string(ifile) * "\n");
    fileID = open(string(filename, string(string(ifile), ".bin")),"w");
else
    print("Writing solution into file\n");
    fileID = open(string(filename, ".bin"),"w");
end

write(fileID,Float64(length(nsize[:])));
write(fileID,Float64.(nsize[:]));

if nsize[1]>0
    write(fileID,Float64.(ndims[:]));
end
if nsize[2]>0
    write(fileID,xdg[:]);
end
# if nsize[3]>0
#     write(fileID,udg[:]);
# end
# if nsize[4]>0
#     write(fileID,odg[:]);
# end
# if nsize[5]>0
#     write(fileID,wdg[:]);
# end

close(fileID);

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
""" Diagrams in a category and their morphisms.
"""
module Diagrams
export Diagram, DiagramHom, id, op, co, shape, diagram, shape_map, diagram_map

using ...GAT
import ...Theories: dom, codom, id, compose, ⋅, ∘, munit
using ...Theories: Category, composeH
import ..Categories: ob_map, hom_map
using ..FinCats, ..FreeDiagrams
using ..FinCats: mapvals
import ..FinCats: force, collect_ob, collect_hom
import ..Limits: limit, colimit, universal

# TODO: Implement these functions more generally, and move elsewhere.

""" Opposite of a category or, more generally, 1-cell dual of a 2-category.
"""
function op end

""" 2-cell dual of a 2-category.
"""
function co end

# Data types
############

""" Diagram in a category.

Recall that a *diagram* in a category ``C`` is a functor ``D: J → C``, where for
us the *shape category* ``J`` is finitely presented. Although such a diagram is
captured perfectly well by a `FinDomFunctor`, there are several different
notions of morphism between diagrams. This simple wrapper type exists to
distinguish them. See [`DiagramHom`](@ref) for more about the morphisms.
"""
struct Diagram{T,C<:Cat,D<:Functor{<:FinCat,C}}
  diagram::D
end
Diagram{T}(F::D) where {T,C<:Cat,D<:Functor{<:FinCat,C}} = Diagram{T,C,D}(F)

Diagram{T}(d::Diagram) where T = Diagram{T}(d.diagram)
Diagram(args...) = Diagram{id}(args...)

""" Functor underlying a diagram object.
"""
diagram(d::Diagram) = d.diagram

""" The *shape* or *indexing category* of a diagram.

This is the domain of the underlying functor.
"""
shape(d::Diagram) = dom(diagram(d))

Base.:(==)(d1::Diagram{T}, d2::Diagram{S}) where {T,S} =
  T == S && diagram(d1) == diagram(d2)

ob_map(d::Diagram, x) = ob_map(diagram(d), x)
hom_map(d::Diagram, f) = hom_map(diagram(d), f)

force(d::Diagram{T}) where T = Diagram{T}(force(diagram(d)))

function Base.show(io::IO, d::Diagram{T}) where T
  print(io, "Diagram{$T}(")
  show(io, diagram(d))
  print(io, ")")
end

""" Morphism of diagrams in a category.

In fact, this type encompasses several different kinds of morphisms from a
diagram ``D: J → C`` to another diagram ``D′: J′ → C``:

1. `DiagramHom{id}`: a functor ``F: J → J′`` together with a natural
   transformation ``ϕ: D ⇒ F⋅D′``
2. `DiagramHom{op}`: a functor ``F: J′ → J`` together with a natural
   transformation ``ϕ: F⋅D ⇒ D′``
3. `DiagramHom{co}`: a functor ``F: J → J′`` together with a natural
   transformation ``ϕ: F⋅D′ ⇒ D``.

Note that `Diagram{op}` is *not* the opposite category of `Diagram{id}`, but
`Diagram{op}` and `Diagram{co}` are opposites of each other. Explicit support is
included for both because they are useful for different purposes: morphisms of
type `DiagramHom{op}` induce morphisms between the limits of the diagrams,
whereas morphisms of type `DiagramHom{co}` generalize morphisms of polynomial
functors.
"""
struct DiagramHom{T,C<:Cat,F<:FinFunctor,Φ<:FinTransformation,D<:Functor{<:FinCat,C}}
  shape_map::F
  diagram_map::Φ
  precomposed_diagram::D
end
DiagramHom{T}(shape_map::F, diagram_map::Φ, precomposed_diagram::D) where
    {T,C,F<:FinFunctor,Φ<:FinTransformation,D<:Functor{<:FinCat,C}} =
  DiagramHom{T,C,F,Φ,D}(shape_map, diagram_map, precomposed_diagram)

DiagramHom{T}(f::DiagramHom) where T =
  DiagramHom{T}(f.shape_map, f.diagram_map, f.precomposed_diagram)
DiagramHom(args...) = DiagramHom{id}(args...)

DiagramHom{T}(ob_maps, hom_map, D::Diagram{T}, D′::Diagram{T}) where T =
  DiagramHom{T}(ob_maps, hom_map, diagram(D), diagram(D′))
DiagramHom{T}(ob_maps, D::Union{Diagram{T},FinDomFunctor},
              D′::Union{Diagram{T},FinDomFunctor}) where T =
  DiagramHom{T}(ob_maps, nothing, D, D′)

function DiagramHom{id}(ob_maps, hom_map, D::FinDomFunctor, D′::FinDomFunctor)
  f = FinFunctor(mapvals(cell1, ob_maps), hom_map, dom(D), dom(D′))
  DiagramHom{id}(f, mapvals(x -> cell2(D′,x), ob_maps), D, D′)
end
function DiagramHom{op}(ob_maps, hom_map, D::FinDomFunctor, D′::FinDomFunctor)
  f = FinDomFunctor(mapvals(cell1, ob_maps), hom_map, dom(D′), dom(D))
  DiagramHom{op}(f, mapvals(x -> cell2(D,x), ob_maps), D, D′)
end
function DiagramHom{co}(ob_maps, hom_map, D::FinDomFunctor, D′::FinDomFunctor)
  f = FinDomFunctor(mapvals(cell1, ob_maps), hom_map, dom(D), dom(D′))
  DiagramHom{co}(f, mapvals(x -> cell2(D′,x), ob_maps), D, D′)
end

function DiagramHom{id}(f::FinFunctor, components, D::FinDomFunctor, D′::FinDomFunctor)
  ϕ = FinTransformation(components, D, f⋅D′)
  DiagramHom{id}(f, ϕ, D′)
end
function DiagramHom{op}(f::FinFunctor, components, D::FinDomFunctor, D′::FinDomFunctor)
  ϕ = FinTransformation(components, f⋅D, D′)
  DiagramHom{op}(f, ϕ, D)
end
function DiagramHom{co}(f::FinFunctor, components, D::FinDomFunctor, D′::FinDomFunctor)
  ϕ = FinTransformation(components, f⋅D′, D)
  DiagramHom{co}(f, ϕ, D′)
end

cell1(pair::Union{Pair,Tuple{Any,Any}}) = first(pair)
cell1(x) = x
cell2(D::FinDomFunctor, pair::Union{Pair,Tuple{Any,Any}}) = last(pair)
cell2(D::FinDomFunctor, x) = id(codom(D), ob_map(D, x))

shape_map(f::DiagramHom) = f.shape_map
diagram_map(f::DiagramHom) = f.diagram_map

Base.:(==)(f::DiagramHom{T}, g::DiagramHom{S}) where {T,S} =
  T == S && shape_map(f) == shape_map(g) && diagram_map(f) == diagram_map(g) &&
  f.precomposed_diagram == g.precomposed_diagram

ob_map(f::DiagramHom, x) = (ob_map(f.shape_map, x), component(f.diagram_map, x))
hom_map(f::DiagramHom, g) = hom_map(f.shape_map, g)

collect_ob(f::DiagramHom) =
  collect(zip(collect_ob(f.shape_map), components(f.diagram_map)))
collect_hom(f::DiagramHom) = collect_hom(f.shape_map)

function Base.show(io::IO, f::DiagramHom{T}) where T
  J = dom(shape_map(f))
  print(io, "DiagramHom{$T}([")
  join(io, mapvals(x -> ob_map(f,x), ob_generators(J), iter=true), ", ")
  print(io, "], [")
  join(io, mapvals(g -> hom_map(f,g), hom_generators(J), iter=true), ", ")
  print(io, "], ")
  show(IOContext(io, :compact=>true, :hide_domains=>true), diagram(dom(f)))
  print(io, ", ")
  show(IOContext(io, :compact=>true, :hide_domains=>true), diagram(codom(f)))
  print(io, ")")
end

# Categories of diagrams
########################

dom_diagram(f::DiagramHom{id}) = dom(diagram_map(f))
dom_diagram(f::DiagramHom{op}) = f.precomposed_diagram
dom_diagram(f::DiagramHom{co}) = codom(diagram_map(f))
codom_diagram(f::DiagramHom{id}) = f.precomposed_diagram
codom_diagram(f::DiagramHom{op}) = codom(diagram_map(f))
codom_diagram(f::DiagramHom{co}) = f.precomposed_diagram

dom(f::DiagramHom{T}) where T = Diagram{T}(dom_diagram(f))
codom(f::DiagramHom{T}) where T = Diagram{T}(codom_diagram(f))

function id(d::Diagram{T}) where T
  F = diagram(d)
  DiagramHom{T}(id(dom(F)), id(F), F)
end

function compose(f::DiagramHom{id}, g::DiagramHom{id})
  DiagramHom{id}(
    shape_map(f) ⋅ shape_map(g),
    diagram_map(f) ⋅ (shape_map(f) * diagram_map(g)),
    codom_diagram(g))
end
function compose(f::DiagramHom{op}, g::DiagramHom{op})
  DiagramHom{op}(
    shape_map(g) ⋅ shape_map(f),
    (shape_map(g) * diagram_map(f)) ⋅ diagram_map(g),
    dom_diagram(f))
end
function compose(f::DiagramHom{co}, g::DiagramHom{co})
  DiagramHom{co}(
    shape_map(f) ⋅ shape_map(g),
    (shape_map(f) * diagram_map(g)) ⋅ diagram_map(f),
    codom_diagram(g))
end

# TODO: The diagrams in a category naturally form a 2-category, but for now we
# just implement the category struture.

@instance Category{Diagram,DiagramHom} begin
  @import dom, codom, compose, id
end

op(d::Diagram{op}) = Diagram{co}(d)
op(d::Diagram{co}) = Diagram{op}(d)
op(f::DiagramHom{op}) = DiagramHom{co}(f)
op(f::DiagramHom{co}) = DiagramHom{op}(f)

# Any functor ``F: C → D`` induces a functor ``Diag(F): Diag(C) → Diag(D)`` by
# post-composition and post-whiskering.

function compose(d::Diagram{T}, F::Functor; kw...) where T
  Diagram{T}(compose(diagram(d), F; kw...))
end
function compose(f::DiagramHom{T}, F::Functor; kw...) where T
  DiagramHom{T}(shape_map(f), composeH(diagram_map(f), F; kw...),
                compose(f.precomposed_diagram, F; kw...))
end

# Limits and colimits
#####################

# In a cocomplete category `C`, colimits define a functor `Diag{id,C} → C`.
# Dually, in a complete category `C`, limits define functor `Diag{op,C} → C`.

limit(d::Diagram{op}; alg=nothing) = limit(diagram(d), alg)
colimit(d::Diagram{id}; alg=nothing) = colimit(diagram(d), alg)

function universal(f::DiagramHom{op}, dom_lim, codom_lim)
  J′ = shape(codom(f))
  cone = Multispan(apex(dom_lim), map(ob_generators(J′)) do j′
    j, g = ob_map(f, j′)
    πⱼ = legs(dom_lim)[j]
    compose(πⱼ, g)
  end)
  universal(codom_lim, cone)
end

function universal(f::DiagramHom{id}, dom_colim, codom_colim)
  J = shape(dom(f))
  cocone = Multicospan(apex(codom_colim), map(ob_generators(J)) do j
    j′, g = ob_map(f, j)
    ιⱼ′ = legs(codom_colim)[j′]
    compose(g, ιⱼ′)
  end)
  universal(dom_colim, cocone)
end

# Monads of diagrams
####################

# TODO: Define monad multiplications that go along with the units.

function munit(::Type{Diagram{T}}, C::Cat, x; shape=nothing) where T
  if isnothing(shape)
    shape = FinCat(1)
  else
    @assert is_discrete(shape) && length(ob_generators(shape)) == 1
  end
  Diagram{T}(FinDomFunctor([x], shape, C))
end

function munit(::Type{DiagramHom{T}}, C::Cat, f;
               dom_shape=nothing, codom_shape=nothing) where T
  f = hom(C, f)
  d = munit(Diagram{T}, C, dom(C, f), shape=dom_shape)
  d′= munit(Diagram{T}, C, codom(C, f), shape=codom_shape)
  j = only(ob_generators(shape(d′)))
  DiagramHom{T}([Pair(j, f)], d, d′)
end

end


text:
# Make sure we're running from the examples directory.
# cd(@__DIR__())

# reload("Eratosthenes")
# sleep(0.1)

module Temp

# using BenchmarkTools
using Eratosthenes # Bring in the top-level sim environment.
using Plots # For plotting (Pkg.add("Plots"); Pkg.add("PyPlot");)
using Printf
import HDF5 # For loading the log file (Pkg.add("HDF5");)
# import Juno # For the progress bar

# Include a module that's not part of Eratosthenes itself. Anything we "use"
# here will be visible during the setup portion of the scenario. So we can
# "setup" things that Eratosthenes doesn't even explicitly know about it!
# include("EratosthenesExampleSoftware.jl")
# using .EratosthenesExampleSoftware

# Create the hierarchy of objects used by the sim.
scenario = setup("scenarios/reaction-wheel.yaml")
# scenario.sim.log_file = ""

# # Just run it.
# simulate(scenario)

# # Profile it.
# scenario.sim.log_file = ""
# @profile simulate(scenario)
# open("profiler.txt", "w") do s
#     Profile.print(IOContext(s, :displaysize => (24, 500)))
# end
# Profile.clear()

# Run it with text progress (e.g., from VS Code).
simulate(scenario) do k, n
    if (k-1) % 100 == 0
        @printf("Progress: % 5.1f%% done.\n", round(100*(k-1)/(n-1), digits=1))
    end
    return true
end

# Create a Juno waitbar. This do-block syntax makes sure it gets "closed" even
# if there's an error.
# Juno.progress(name="Simulating...") do progress_bar
#
#     # Run the simulation. Anything inside this do-block will be run at the
#     # beginning of every major time step. This is useful for updating a progress
#     # bar or a plot.
#     simulate(scenario) do k, n
#         if k % 30 == 0
#             Juno.progress(progress_bar, (k-1)/n) # Update Juno's progress bar.
#         end
#         return true # Return false to end the sim.
#     end
#
# end

# Open the logs and plot some things.
if !isempty(scenario.sim.log_file)

    # Read some data from the logs. This do block takes care of
    # error-handling, so we can just read whatever we want and return it.
    t, q_BI, ω_BI_B, rw1, rw2, rw3 = # , t_st, q_BI_meas, t_gyro, ω_BI_B_meas
        HDF5.h5open(scenario.sim.log_file, "r") do logs
            (read(logs, "/smallsat1/body/state/time"),
             read(logs, "/smallsat1/body/state/data/q_BI"),
             read(logs, "/smallsat1/body/state/data/ω_BI_B"),
            #  read(logs, "/cube1/rw1/state/data"),
            #  read(logs, "/cube1/rw2/state/data"),
            #  read(logs, "/cube1/rw3/state/data"),
             read(logs, "/smallsat1/rw1/state/data/ω_RI_G"),
             read(logs, "/smallsat1/rw2/state/data/ω_RI_G"),
             read(logs, "/smallsat1/rw3/state/data/ω_RI_G"),
             #read(logs, "/cube1/star_tracker/outputs/time"),
             #read(logs, "/cube1/star_tracker/outputs/data"),
             #read(logs, "/cube1/gyro/outputs/time"),
             #read(logs, "/cube1/gyro/outputs/data")
             )
        end

    # Scalars are logged as 1-by-n (don't know why). Convert to just n for
    # the sake of plotting.
    t      = dropdims(t, dims=1)
    #t_st   = dropdims(t_st, dims=1)
    #t_gyro = dropdims(t_gyro, dims=1)

    # Choose a friendly plotting package.
    # pyplot()
    # plotlyjs()
    gr()

    # Define each plot that we'll need.
    display(plot(t, transpose(q_BI),
                 label  = ["q1" "q2" "q3" "q4"],
                 xlabel = "Time (s)",
                 ylabel = "Attitude Quaternion"))
    display(plot(t, 180/π * transpose(ω_BI_B),
                 label  = ["ω1" "ω2" "ω3"],
                 xlabel = "Time (s)",
                 ylabel = "Rotation Rate (deg/s)"))

    # display(plot(randn(5,5))) # Dummy. This doesn't show up, but the next will?

    display(plot(t, 30/π*[rw1[1,:] rw2[1,:] rw3[1,:]],
                 label  = ["rw1" "rw2" "rw3"],
                 xlabel = "Time (s)",
                 ylabel = "Wheel Rate (RPM)"))

    # # For coupled reaction wheels:
    # display(plot(t, 30/π*(rw1 - transpose(ω_BI_B)),
    #              xlabel = "Time (s)",
    #              ylabel = "RW1 Rate (RPM)"))
    # display(plot(t, 30/π*(rw1 - transpose(ω_BI_B)),
    #              xlabel = "Time (s)",
    #              ylabel = "RW1 Rate (RPM)")) # I need two of this plot to get it plotted once??
    # display(plot(t, 30/π*([0. -1. 0.; 1. 0. 0.; 0. 0. 1.] * rw2 - transpose(ω_BI_B)),
    #              xlabel = "Time (s)",
    #              ylabel = "RW2 Rate (RPM)"))
    # display(plot(t, 30/π*([0. 0. -1.; 0. 1. 0.; 1. 0. 0.] * rw3 - transpose(ω_BI_B)),
    #              xlabel = "Time (s)",
    #              ylabel = "RW3 Rate (RPM)"))

    # # For simple reaction wheels:
    # display(plot(t, 30/π*dropdims(rw1,dims=1),
    #              xlabel = "Time (s)",
    #              ylabel = "RW1 Rate (RPM)"))
    # display(plot(t, 30/π*dropdims(rw2,dims=1),
    #              xlabel = "Time (s)",
    #              ylabel = "RW2 Rate (RPM)"))
    # display(plot(t, 30/π*dropdims(rw3,dims=1),
    #              xlabel = "Time (s)",
    #              ylabel = "RW3 Rate (RPM)"))

    # Save the individual plots before showing them together.
    # savefig(p1, "out/sim-results-q.png")
    # savefig(p2, "out/sim-results-ω.png")

    # Display those plots as subplots using the default "2" layout. Put any
    # common descriptors here.
    # display(plot(p1, p2, p3,
    #              layout = 3,
    #              xlabel = "Time (s)"))

end # plotting

end # Temp


text:
import Base: ==, -, *, +, append!, convert, copy, delete!, deleteat!,
  firstindex, getindex, hash, in, isempty, isequal, lastindex, length,
  push!, merge, merge!, read, setindex!, show, size, sizeof, sort!,
  sort, summary, write
import DSP: filtfilt


text:
################################################################################
# PACKAGES USED BY ALL SUBMODULES
################################################################################
using Cascadia
using Gumbo
using HTTP
using DataFrames

# Used by shared.jl
using Dates

################################################################################
# Data Structures
################################################################################

const Nullable{T} = Union{T, Nothing}

################################################################################
# Convenience Methods
################################################################################

function convert_to_seconds(minutes::AbstractString)::Integer
    # A small workaround to get the total seconds
    mins = Dates.DateTime(minutes, "MM:SS")
    # Round down to the nearest date
    dfloor = DateTime(Date(mins))
    secs_obj = convert(Dates.Second, Dates.Period(mins - dfloor))
    return Dates.value(secs_obj)
end

function convert_dict_to_dataframe(array_of_dict::Array{Dict{Symbol, Any}})::DataFrame
    return vcat(DataFrame.(array_of_dict)...)
end

function add_missing_columns!(plays::Array{Dict{Symbol, Any}})
    all_keys = Set{Symbol}()
    [union!(all_keys, keys(play)) for play in plays]
    for play in plays
        missing_keys = setdiff(all_keys, keys(play))
        new_dict = Dict(key=>nothing for key in missing_keys)
        merge!(play, new_dict)
    end
end

function fix_name(name::String)::String
    return replace(uppercase(replace(name, r"\.|\_" => "")), r"\s|\-"=>".")
end


text:
# implementation of the GPUCompiler interfaces for generating SPIR-V code

# https://github.com/llvm/llvm-project/blob/master/clang/lib/Basic/Targets/SPIR.h
# https://github.com/KhronosGroup/LLVM-SPIRV-Backend/blob/master/llvm/docs/SPIR-V-Backend.rst
# https://github.com/KhronosGroup/SPIRV-LLVM-Translator/blob/master/docs/SPIRVRepresentationInLLVM.rst

const SPIRV_LLVM_Translator_unified_jll = LazyModule("SPIRV_LLVM_Translator_unified_jll", UUID("85f0d8ed-5b39-5caa-b1ae-7472de402361"))
const SPIRV_Tools_jll = LazyModule("SPIRV_Tools_jll", UUID("6ac6d60f-d740-5983-97d7-a4482c0689f4"))

## target

export SPIRVCompilerTarget

Base.@kwdef struct SPIRVCompilerTarget <: AbstractCompilerTarget
end

llvm_triple(::SPIRVCompilerTarget) = Int===Int64 ? "spir64-unknown-unknown" : "spirv-unknown-unknown"

# SPIRV is not supported by our LLVM builds, so we can't get a target machine
llvm_machine(::SPIRVCompilerTarget) = nothing

llvm_datalayout(::SPIRVCompilerTarget) = Int===Int64 ?
    "e-i64:64-v16:16-v24:32-v32:32-v48:64-v96:128-v192:256-v256:256-v512:512-v1024:1024" :
    "e-p:32:32-i64:64-v16:16-v24:32-v32:32-v48:64-v96:128-v192:256-v256:256-v512:512-v1024:1024"

## job

# TODO: encode debug build or not in the compiler job
#       https://github.com/JuliaGPU/CUDAnative.jl/issues/368
runtime_slug(job::CompilerJob{SPIRVCompilerTarget}) = "spirv"

function process_module!(job::CompilerJob{SPIRVCompilerTarget}, mod::LLVM.Module)
    # calling convention
    for f in functions(mod)
        # JuliaGPU/GPUCompiler.jl#97
        #callconv!(f, LLVM.API.LLVMSPIRFUNCCallConv)
    end
end

function process_entry!(job::CompilerJob{SPIRVCompilerTarget}, mod::LLVM.Module, entry::LLVM.Function)
    entry = invoke(process_entry!, Tuple{CompilerJob, LLVM.Module, LLVM.Function}, job, mod, entry)

    if job.source.kernel
        # calling convention
        callconv!(entry, LLVM.API.LLVMSPIRKERNELCallConv)
    end

    return entry
end

function finish_module!(job::CompilerJob{SPIRVCompilerTarget}, mod::LLVM.Module, entry::LLVM.Function)
    ctx = context(mod)
    entry = invoke(finish_module!, Tuple{CompilerJob, LLVM.Module, LLVM.Function}, job, mod, entry)

    if job.source.kernel
        # HACK: Intel's compute runtime doesn't properly support SPIR-V's byval attribute.
        #       they do support struct byval, for OpenCL, so wrap byval parameters in a struct.
        entry = wrap_byval(job, mod, entry)
    end

    # add module metadata
    ## OpenCL 2.0
    push!(metadata(mod)["opencl.ocl.version"],
          MDNode([ConstantInt(Int32(2); ctx),
                  ConstantInt(Int32(0); ctx)]; ctx))
    ## SPIR-V 1.5
    push!(metadata(mod)["opencl.spirv.version"],
          MDNode([ConstantInt(Int32(1); ctx),
                  ConstantInt(Int32(5); ctx)]; ctx))

    return entry
end

@unlocked function mcgen(job::CompilerJob{SPIRVCompilerTarget}, mod::LLVM.Module,
                         format=LLVM.API.LLVMAssemblyFile)
    # The SPIRV Tools don't handle Julia's debug info, rejecting DW_LANG_Julia...
    strip_debuginfo!(mod)

    ModulePassManager() do pm
        # SPIR-V does not support trap, and has no mechanism to abort compute kernels
        # (OpKill is only available in fragment execution mode)
        add!(pm, ModulePass("RemoveTrap", rm_trap!))

        # the LLVM to SPIR-V translator does not support the freeze instruction
        # (SPIRV-LLVM-Translator#1140)
        add!(pm, ModulePass("RemoveFreeze", rm_freeze!))

        run!(pm, mod)
    end

    # translate to SPIR-V
    input = tempname(cleanup=false) * ".bc"
    translated = tempname(cleanup=false) * ".spv"
    write(input, mod)
    SPIRV_LLVM_Translator_unified_jll.llvm_spirv() do translator
        proc = run(ignorestatus(`$translator --spirv-debug-info-version=ocl-100 -o $translated $input`))
        if !success(proc)
            error("""Failed to translate LLVM code to SPIR-V.
                     If you think this is a bug, please file an issue and attach $(input).""")
        end
    end

    # validate
    # XXX: parameterize this on the `validate` driver argument
    # XXX: our code currently doesn't pass the validator
    if Base.JLOptions().debug_level >= 2 && false
        SPIRV_Tools_jll.spirv_val() do validator
            proc = run(ignorestatus(`$validator $translated`))
            if !success(proc)
                error("""Failed to validate generated SPIR-V.
                         If you think this is a bug, please file an issue and attach $(input) and $(translated).""")
            end
        end
    end

    # optimize
    # XXX: parameterize this on the `optimize` driver argument
    # XXX: the optimizer segfaults on some of our code
    optimized = tempname(cleanup=false) * ".spv"
    if false
        SPIRV_Tools_jll.spirv_opt() do optimizer
            proc = run(ignorestatus(`$optimizer -O --skip-validation $translated -o $optimized`))
            if !success(proc)
                error("""Failed to optimize generated SPIR-V.
                         If you think this is a bug, please file an issue and attach $(input) and $(translated).""")
            end
        end
    end

    output = if format == LLVM.API.LLVMObjectFile
        read(translated)
    else
        # disassemble
        SPIRV_Tools_jll.spirv_dis() do disassembler
            read(`$disassembler $translated`, String)
        end
    end

    rm(input)
    rm(translated)
    #rm(optimized)

    return output
end

# reimplementation that uses `spirv-dis`, giving much more pleasant output
function code_native(io::IO, job::CompilerJob{SPIRVCompilerTarget}; raw::Bool=false, dump_module::Bool=false)
    obj, _ = codegen(:obj, job; strip=!raw, only_entry=!dump_module, validate=false)
    mktemp() do input_path, input_io
        write(input_io, obj)
        flush(input_io)

        SPIRV_Tools_jll.spirv_dis() do disassembler
            if io == stdout
                run(`$disassembler $input_path`)
            else
                mktemp() do output_path, output_io
                    run(`$disassembler $input_path -o $output_path`)
                    asm = read(output_io, String)
                    print(io, asm)
                end
            end
        end
    end
end

## LLVM passes

# remove llvm.trap and its uses from a module
function rm_trap!(mod::LLVM.Module)
    job = current_job::CompilerJob
    changed = false
    @timeit_debug to "remove trap" begin

    if haskey(functions(mod), "llvm.trap")
        trap = functions(mod)["llvm.trap"]

        for use in uses(trap)
            val = user(use)
            if isa(val, LLVM.CallInst)
                unsafe_delete!(LLVM.parent(val), val)
                changed = true
            end
        end

        @compiler_assert isempty(uses(trap)) job
        unsafe_delete!(mod, trap)
    end

    end
    return changed
end

# remove freeze and replace uses by the original value
# (KhronosGroup/SPIRV-LLVM-Translator#1140)
function rm_freeze!(mod::LLVM.Module)
    job = current_job::CompilerJob
    changed = false
    @timeit_debug to "remove freeze" begin

    for f in functions(mod), bb in blocks(f), inst in instructions(bb)
        if inst isa LLVM.FreezeInst
            orig = first(operands(inst))
            replace_uses!(inst, orig)
            @compiler_assert isempty(uses(inst)) job
            unsafe_delete!(bb, inst)
            changed = true
        end
    end

    end
    return changed
end

# wrap byval pointers in a single-value struct
function wrap_byval(@nospecialize(job::CompilerJob), mod::LLVM.Module, f::LLVM.Function)
    ctx = context(mod)
    ft = eltype(llvmtype(f)::LLVM.PointerType)::LLVM.FunctionType
    @compiler_assert LLVM.return_type(ft) == LLVM.VoidType(ctx) job

    # find the byval parameters
    byval = BitVector(undef, length(parameters(ft)))
    if LLVM.version() >= v"12"
        for i in 1:length(byval)
            attrs = collect(parameter_attributes(f, i))
            byval[i] = any(attrs) do attr
                kind(attr) == kind(EnumAttribute("byval", 0; ctx))
            end
        end
    else
        # XXX: byval is not round-trippable on LLVM < 12 (see maleadt/LLVM.jl#186)
        args = classify_arguments(job, ft)
        filter!(args) do arg
            arg.cc != GHOST
        end
        for arg in args
            if arg.cc == BITS_REF
                byval[arg.codegen.i] = true
            end
        end
    end

    # generate the wrapper function type & definition
    new_types = LLVM.LLVMType[]
    for (i, param) in enumerate(parameters(ft))
        typ = if byval[i]
            st = LLVM.StructType([eltype(param)]; ctx)
            LLVM.PointerType(st, addrspace(param))
        else
            param
        end
        push!(new_types, typ)
    end
    new_ft = LLVM.FunctionType(LLVM.VoidType(ctx), new_types)
    new_f = LLVM.Function(mod, "", new_ft)
    linkage!(new_f, linkage(f))
    for (arg, new_arg) in zip(parameters(f), parameters(new_f))
        LLVM.name!(new_arg, LLVM.name(arg))
    end

    # emit IR performing the "conversions"
    new_args = Vector{LLVM.Value}()
    Builder(ctx) do builder
        entry = BasicBlock(new_f, "entry"; ctx)
        position!(builder, entry)

        # perform argument conversions
        for (i, param) in enumerate(parameters(new_f))
            if byval[i]
                ptr = struct_gep!(builder, param, 0)
                push!(new_args, ptr)
            else
                push!(new_args, param)
            end
        end

        # map the arguments
        value_map = Dict{LLVM.Value, LLVM.Value}(
            param => new_args[i] for (i,param) in enumerate(parameters(f))
        )

        value_map[f] = new_f
        clone_into!(new_f, f; value_map,
                    changes=LLVM.API.LLVMCloneFunctionChangeTypeGlobalChanges)

        # apply byval attributes again (`clone_into!` didn't due to the type mismatch)
        for i in 1:length(byval)
            attrs = parameter_attributes(new_f, i)
            if byval[i]
                if LLVM.version() >= v"12"
                    push!(attrs, TypeAttribute("byval", eltype(new_types[i]); ctx))
                else
                    push!(attrs, EnumAttribute("byval", 0; ctx))
                end
            end
        end

        # fall through
        br!(builder, collect(blocks(new_f))[2])
    end

    # remove the old function
    # NOTE: if we ever have legitimate uses of the old function, create a shim instead
    fn = LLVM.name(f)
    @assert isempty(uses(f))
    replace_metadata_uses!(f, new_f)
    unsafe_delete!(mod, f)
    LLVM.name!(new_f, fn)

    # clean-up
    # NOTE: byval wrapping happens very late, after optimization
    ModulePassManager() do pm
        # merge GEPs
        instruction_combining!(pm)

        # fold the entry bb into the rest of the function
        cfgsimplification!(pm)

        run!(pm, mod)
    end

    return new_f
end


text:
function check_args(y_true::Vector, y_pred::Vector)
    @assert length(y_true) == length(y_pred)
    return(true)
end


text:
"""
Data managers for saving and loading data. Stolen from [Config.jl](https://github.com/ajjacobs/Config.jl)

"""

import FileIO
import HDF5, BSON, JLD2
import Pkg.TOML

"""
    DataManger

Abstract type for various filetype managers.
"""
abstract type SaveManager end

# get manager filetype extension
extension(data_manager::SaveManager) = extension(typeof(data_manager))

# save data, overwriting existing
save(data_manager::SaveManager, path, data) = save(typeof(data_manager), path, data)

# save data, adding to existing
save!(data_manager::SaveManager, path, data) = save!(typeof(data_manager), path, data)

# load data
load(data_manager::SaveManager, path) = load(typeof(data_manager), path)
# load(data_manager::SaveManager, path::String) = raise(DomainError("loand() not defined!"))

# ===============
# --- H D F 5 ---
# ===============

struct HDF5Manager <: SaveManager end

extension(data_manager::Type{HDF5Manager}) = ".h5"

# Saving/Loading data
function _save(data_manager::Type{HDF5Manager}, path, data, writeMode)
    HDF5.h5open(path, writeMode) do f
      for (k,v) in data
          write(f, k, v)
      end
    end
end

save!(data_manager::Type{HDF5Manager}, path, data) = _save(data_manager, path, data, "cw")
save(data_manager::Type{HDF5Manager}, path, data) = _save(data_manager, path, data, "w")

function load(data_manager::Type{HDF5Manager}, path)
    data = Dict()
    HDF5.h5open(path) do f
        keys = names(f)
        for k in keys
            data[k] = read(f[k])
        end
    end
    return data
end

# ===============
# --- B S O N ---
# ===============

struct BSONManager <: SaveManager end

extension(data_manager::Type{BSONManager}) = ".bson"

function save(data_manager::Type{BSONManager}, path, data)
    BSON.bson(path, data)
end

function save!(data_manager::Type{BSONManager}, path, data)
    try
        priorData = BSON.load(path)
        newData = merge(data, BSON.load(path))
        BSON.bson(path, newData)
    catch
        BSON.bson(path, data)
    end
end

function load(data_manager::Type{BSONManager}, path)
    return BSON.load(path)
end

# ===============
# --- J L D 2 ---
# ===============

struct JLD2Manager <: SaveManager end

extension(manager::Type{JLD2Manager}) = ".jld2"

save(manager::Type{JLD2Manager}, path, data) = FileIO.save(path, data) # JLD2.@save
function save!(manager::Type{JLD2Manager}, path, data)
    try
        priorData = load(manager, path)
        merge!(priorData, data)
        save(manager, path, priorData)
    catch
        save(manager, path, data)
    end
end

load(manager::Type{JLD2Manager}, path) = FileIO.load(path)

# ===============
# --- T O M L ---
# ===============

struct TOMLManager <: SaveManager end

extension(manager::Type{TOMLManager}) = ".toml"

function save(manager::Type{TOMLManager}, path, data)
    open(path) do io
        TOML.print(io, data)
    end
end

load(manager::Type{TOMLManager}, path) = TOML.parsefile(path)


text:
# mcmc_inference.jl
# 2019-11-11
# David Merrell
#
# Contains functions relevant for inference
# in our network reconstruction setting.

"""
Generic MCMC inference wrapper.
Arguments:
    * gen_model: a Gen probabilistic program
    * model_args: a tuple of arguments for the "gen_model"
    * update_loop_fn: a function which updates the variables in "gen_model"
    * update_loop_args: a tuple of arguments for "update_loop_fn"
    * update_results_fn: a function which updates the data structure holding MCMC results
    * update_results_args: arguments for "update_results_fn"
    * timeout: maximum execution time (in seconds) (default 1 hour)
    * n_steps: maximum number of iterations (default no limit)
    * thinning: only store results once every "thinning"th iteration
"""
function mcmc_inference(gen_model,
                        model_args,
                        observations,
                        update_loop_fn::Function,
                        update_loop_args::Vector,
                        update_results_fn::Function,
                        update_results_args::Vector;
                        timeout::Float64=3600.0,
                        n_steps::Int64=-1,
                        thinning::Int64=1,
                        chain_id::String="")

    # start the timer
    t_start = time()
    t_elapsed = 0.0

    # Generate an initial trace
    tr, _ = Gen.generate(gen_model, model_args,
			 observations)

    # The results we care about
    results = update_results_fn(nothing, tr, update_results_args)

    # Sampling loop
    print_str = string("Sampling for ",string(timeout), " seconds.")
    if n_steps > 0
        print_str = string(print_str, " (Or ", n_steps, " iterations.)")
    else
        n_steps = typemax(Int64)
        print_str = string(print_str, " (No iteration limit.)")
    end
    println(print_str)

    prop_count = 0
    t_print = 0.0
    while t_elapsed < timeout && prop_count <= n_steps

        # thinning loop (if applicable)
        for i=1:thinning

            # Print progress
            if (prop_count > 0) && (t_elapsed - t_print >= 20.0)
                println("[", chain_id, "]\t", prop_count," iterations;\t", round(t_elapsed), " seconds")
                t_print = t_elapsed
            end

            # update the variables
            tr, acc = update_loop_fn(tr, update_loop_args)
            prop_count += 1

            t_elapsed = time() - t_start
        end

        # update the results
        results = update_results_fn(results, tr, update_results_args)

    end

    return results
end

max_stars_repo_name: gitter-lab/ssps

max_issues_repo_name: gitter-lab/ssps

max_forks_repo_name: gitter-lab/ssps


text:
@testset "Lorenz 96" begin

    using DifferentialEquations, Random, LinearAlgebra
    using AnalogDataAssimilation

    rng = MersenneTwister(123)
    F = 8
    J = 40::Int64
    parameters = [F, J]
    dt_integration = 0.05
    dt_states = 1
    dt_obs = 4
    var_obs = randperm(rng, J)[1:20]
    nb_loop_train = 20
    nb_loop_test = 5
    sigma2_catalog = 0.0
    sigma2_obs = 2.0

    ssm = StateSpaceModel(
        lorenz96,
        dt_integration,
        dt_states,
        dt_obs,
        parameters,
        var_obs,
        nb_loop_train,
        nb_loop_test,
        sigma2_catalog,
        sigma2_obs,
    )

    # 5 time steps (to be in the attractor space)

    u0 = F .* ones(Float64, J)
    u0[J÷2] = u0[J÷2] + 0.01

    tspan = (0.0, 5.0)
    p = [F, J]
    prob = ODEProblem(lorenz96, u0, tspan, p)
    sol = solve(prob, reltol = 1e-6, saveat = dt_integration)

    u0 = last(sol.u)

    xt, yo, catalog = generate_data(ssm, u0)

    local_analog_matrix = BitArray{2}(
        diagm(
            -2 => trues(xt.nv - 2),
            -1 => trues(xt.nv - 1),
            0 => trues(xt.nv),
            1 => trues(xt.nv - 1),
            2 => trues(xt.nv - 2),
            J - 2 => trues(xt.nv - (J - 2)),
            J - 1 => trues(xt.nv - (J - 1)),
        ) + transpose(
            diagm(J - 2 => trues(xt.nv - (J - 2)), J - 1 => trues(xt.nv - (J - 1))),
        ),
    )

    neighborhood = local_analog_matrix
    regression = :local_linear
    sampling = :gaussian

    @testset "Global analog with local linear regression " begin

        f = AnalogForecasting(
            100,
            xt,
            catalog,
            regression = regression,
            sampling = sampling,
        )
        DA = DataAssimilation(f, xt, ssm.sigma2_obs)
        @time x̂ = forecast(DA, yo, AnEnKS(500), progress = false)
        rmse = RMSE(xt, x̂)
        println("RMSE(global analog  DA) = $rmse ")

        @test rmse < 5.0

    end

    @testset "Global analog with classic computation " begin

        DA = DataAssimilation(ssm, xt)
        @time x̂ = forecast(DA, yo, AnEnKS(500), progress = false)
        rmse = RMSE(xt, x̂)
        println("RMSE(global classic DA) = $rmse ")

        @test rmse < 5.0

    end

    @testset "Local analog with local linear regression " begin

        f = AnalogForecasting(100, xt, catalog, neighborhood, regression, sampling)
        DA = DataAssimilation(f, xt, ssm.sigma2_obs)
        @time x̂ = forecast(DA, yo, AnEnKS(500), progress = false)
        rmse = RMSE(xt, x̂)
        println("RMSE(local analog DA) = $rmse ")

        @test rmse < 5.0

    end

    @test true

end


text:
using Morsel
using FactCheck
import Requests
req = Requests

app = Morsel.app()

get(app, "/") do req, res
    "root_get"
end

@async start(app, 8000)

# make sure the app has run
response = req.get("http://localhost:8000/")
while (response.status != 200)
    response = req.get("http://localhost:8000/")
end

facts("route") do
    facts("get, post, post, put, delete, request") do
        # patch doesn't support...
        route(app, POST | PUT | DELETE , "/") do req, res
            "root"
        end

        context("get") do
            response = req.get("http://localhost:8000/")
            @fact response.status => 200
            @fact response.data => "root_get"
        end

        context("post, put, delete") do
            @fact req.post("http://localhost:8000/").data => "root"
            @fact req.put("http://localhost:8000/").data => "root"
            @fact req.delete("http://localhost:8000/").data => "root"
        end
    end

    facts("Add route when the app is running") do
        get(app, "/running") do req, res
            "running"
        end
        @fact req.get("http://localhost:8000/running").status => 200
    end

    facts("dynamic route") do
        context("the request will succeed if the type is right") do
            route(app, GET, "/users/<id::Int>" ) do req, res
                string("User id is:", req.params[:id])
            end

            @fact req.get("http://localhost:8000/users/10").status => 200
            @fact req.get("http://localhost:8000/users/abc").status => 404 # wrong type
        end
        context("can get the value by req.params[:id]") do
            @fact req.get("http://localhost:8000/users/10").data => "User id is:10"
        end
        context("regex route") do
            route(app, GET, "/name/<name::%[a-z]{3}[0-9]{2}>") do req, res
                req.params[:name]
            end

            @fact req.get("http://localhost:8000/name/abc33").data => "abc33"
            @fact req.get("http://localhost:8000/name/abc333").status => 404 # not match the regex
        end
    end

    facts("namespace") do
        namespace(app, "/namespace") do app
            get(app, "/hello/") do req, res
                "hello namespace"
            end
        end

        @fact req.get("http://localhost:8000/namespace/hello").data => "hello namespace"
    end
end

facts("redirect") do
    get(app, "/redirect302") do req, res
        redirect(res, "/") # the default status is 302

    end

    response = req.get("http://localhost:8000/redirect302")
    @fact response.status => 302
    @fact response.headers["Location"] => "/"
end


text:
# Test type relations

using Distributions
using ForwardDiff: Dual

@assert UnivariateDistribution <: Distribution
@assert MultivariateDistribution <: Distribution
@assert MatrixDistribution <: Distribution

@assert DiscreteDistribution <: Distribution
@assert ContinuousDistribution <: Distribution

@assert DiscreteUnivariateDistribution <: DiscreteDistribution
@assert DiscreteUnivariateDistribution <: UnivariateDistribution
@assert ContinuousUnivariateDistribution <: ContinuousDistribution
@assert ContinuousUnivariateDistribution <: UnivariateDistribution
@assert DiscreteMultivariateDistribution <: DiscreteDistribution
@assert DiscreteMultivariateDistribution <: MultivariateDistribution
@assert ContinuousMultivariateDistribution <: ContinuousDistribution
@assert ContinuousMultivariateDistribution <: MultivariateDistribution
@assert DiscreteMatrixDistribution <: DiscreteDistribution
@assert DiscreteMatrixDistribution <: MatrixDistribution
@assert ContinuousMatrixDistribution <: ContinuousDistribution
@assert ContinuousMatrixDistribution <: MatrixDistribution

@testset "Test Sample Type" begin
    for T in (Float64,Float32,Dual{Nothing,Float64,0})
        @testset "Type $T" begin
            for d in (MvNormal,MvLogNormal,MvNormalCanon,Dirichlet)
                dist = d(map(T,ones(2)))
                @test eltype(dist) == T
                @test eltype(rand(dist)) == eltype(dist)
            end
            dist = Distributions.mvtdist(map(T,1.0),map(T,[1.0 0.0; 0.0 1.0]))
            @test eltype(dist) == T
            @test eltype(rand(dist)) == eltype(dist)
        end
    end
end


text:
import Base: LineEdit, REPL

#==
Some of this code is derived from Cxx.jl.

> Copyright (c) 2013-2016: Keno Fischer and other contributors
> Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
> The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
> THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
==#

function process_line(line)
    try
        QuoteNode(eval(Main, insane(line*"\n")))
    catch e
        print_with_color(:red, STDERR, "ERROR: ", sprint(showerror, e), "\n")
    end
end

function create_insane_REPL(prompt, name, repl=Base.active_repl, main_mode=repl.interface.modes[1])
    mirepl = isdefined(repl,:mi) ? repl.mi : repl
    # Setup insane panel
    panel = LineEdit.Prompt(prompt;
        on_enter = s -> try
            insane(String(LineEdit.buffer(s).data)*"\n")
            true
        catch e
            isa(e, BoundsError) ? false : rethrow(e)
        end)

    panel.on_done = REPL.respond(process_line,repl,panel)

    main_mode == mirepl.interface.modes[1] &&
        push!(mirepl.interface.modes,panel)

    hp = main_mode.hist
    hp.mode_mapping[name] = panel
    panel.hist = hp

    search_prompt, skeymap = LineEdit.setup_search_keymap(hp)
    mk = REPL.mode_keymap(main_mode)

    b = Dict{Any,Any}[skeymap, mk, LineEdit.history_keymap, LineEdit.default_keymap, LineEdit.escape_defaults]
    panel.keymap_dict = LineEdit.keymap(b)

    panel
end

function run_insane_REPL(prompt="λ> ", name=:insane, key='<')
    repl = Base.active_repl
    mirepl = isdefined(repl,:mi) ? repl.mi : repl
    main_mode = mirepl.interface.modes[1]

    panel = create_insane_REPL(prompt, name)

    # Install this mode into the main mode
    const insane_keymap = Dict{Any,Any}(
        key => function (s,args...)
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                buf = copy(LineEdit.buffer(s))
                LineEdit.transition(s, panel) do
                    LineEdit.state(s, panel).input_buffer = buf
                end
            else
                LineEdit.edit_insert(s,key)
            end
        end
    )
    main_mode.keymap_dict = LineEdit.keymap_merge(main_mode.keymap_dict, insane_keymap);
    nothing
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 9f60e65c-904e-11eb-3f83-639b5c45bb88
begin
    import Pkg
    Pkg.activate(mktempdir())

    Pkg.add([
        Pkg.PackageSpec(name="DataFrames")
    ])
	using DataFrames
end

# ╔═╡ Cell order:
# ╠═9f60e65c-904e-11eb-3f83-639b5c45bb88

max_stars_repo_name: likanzhan/statistics

max_issues_repo_name: likanzhan/statistics

max_forks_repo_name: likanzhan/statistics


text:
"""
ArithmeticBlocks
"""

export ProductBlock

mutable struct ProductBlock <: AbstractArithmeticBlock
    inport::Vector{InPort}
    outport::Vector{OutPort}

    function ProductBlock()
        @createblock new(Vector{InPort}(), Vector{OutPort}()) 2 1
    end
end

"""
IO
"""

function Base.show(io::IO, blk::ProductBlock)
    println(io, "ProductBlock()")
end

"""
to expr
"""

function _toexpr(blk::ProductBlock)
    Expr(:call, :*, _toexpr(blk.inport[1]), _toexpr(blk.inport[2]))
end


text:
using AdFem
using PyPlot
mmesh = Mesh(10, 10, 1/10)
close("all")
visualize_mesh(mmesh)
savefig("mesh1.png")

mmesh = Mesh_v2(10, 10, 1/10)
close("all")
visualize_mesh(mmesh)
savefig("mesh2.png")

mmesh = Mesh_v3(10, 10, 1/10)
close("all")
visualize_mesh(mmesh)
savefig("mesh3.png")


text:
using Random
using Statistics
using StatsBase
using Plots
import Base.+
import Base.-

struct Point{T}
	x::T
	y::T
end

Base.:+(a::Point{Float64}, b::Point{Float64}) = Point{Float64}(a.x+b.x , a.y+b.y)
Base.:-(a::Point{Float64}, b::Point{Float64}) = Point{Float64}(a.x-b.x , a.y-b.y)
Base.:*(a::Point{Float64}, b::Point{Float64}) = Point{Float64}(a.x*b.x , a.y*b.y)

function gauss_3D(p::Point{Float64})
	return exp( -(p.x)^2 - (p.y)^2 )
end

function gradino(p::Point{Float64})
	if (p.x > -1 && p.x < 1 && p.y > -1 && p.y < 1)
		return 1.0
	else
		return 0.0
	end
end

function cilindro(p::Point{Float64})
	if (p.x^2 + p.y^2 < 1 )
		return 1.0
	else
		return 0.0
	end
end

function Integrate(fun::Function, xl::Point{Float64}, xu::Point{Float64}, cal::Int64)

	#; Primo step, cal/3. valutazioni della funzione
	rng = MersenneTwister(0)
	pt = [Point{Float64}(rand(rng), rand(rng)) for n = 1:Int64(floor(calls/2.))]

	#; Cambio coordinate, per mettermi nel range di integrazione scelto
	p = ( pt .* Ref(xu-xl) ) .+ Ref(xl)
	density_unif = 1 / Int64(floor(calls/2.))

	#; Valuto l'integrale
	f = [ fun(point) * (xu.x-xl.x) * (xu.y-xl.y) * density_unif  for point in p ]
	f_abs = abs.(f)

	Integral      = sum(f)
	Integral_abs  = sum(f_abs)

	#; Valore dell'integrale di prova
	Int_Vero = 2.23099
	println("0. Integrale = ",(Integral - Int_Vero)," ", Integral)

	#; Calcolo la nuova densita` di probabilita`
	g = (f_abs ./ Integral_abs)

	#for i = 1:length(g)
	#	println(pt[i].x," ",pt[i].y," ",g[i])
	#end

	#println("1. Controllo integrale pdf = ",sum(g))

	#println(1+log(2,cal/2.) + log(2, 1+skewness(g)/((6*(cal/2. -2)/((cal/2. + 1)*(cal/2. + 2)))^0.5)  ) )

	NSIDE = 1000

	#; posso fare qui la separazione dei punti nelle varie BOX

	w_m	 = zeros(Int64(NSIDE^2))
	Int_box  = zeros(Int64(NSIDE^2))
	Int_abs_box = zeros(Int64(NSIDE^2))
	dix = diy = diz = 1.0/NSIDE
	xx = Array{Float64}(undef, 0)
	yy = Array{Float64}(undef, 0)

	Box_X = Array{Float64}(undef, 0)
	Box_Y = Array{Float64}(undef, 0)

	for iter = 1:20

		#; Organizzo la pdf nel BOX
		for idx_g = 1:length(g)
			box_id_x = Int64(ceil(pt[idx_g].x / dix)) -1
			box_id_y = Int64(ceil(pt[idx_g].y / diy)) -1
			ID = (box_id_x) + (box_id_y)*NSIDE
			w_m[ID+1] += g[idx_g]
		end

		#println("Densita` in BOX")
		#for i in 1:length(w_m)
		#	println(w_m[i])
		#end
		#println(iter+1,". ","Controllo Integrale PDF = ", sum(w_m))
		control = 0
		pt_2_all = Array{Point{Float64}}(undef, 0)
		f_2_all  = Array{Float64}(undef, 0)
		p_2 = Array{Point{Float64}}(undef, 0)

		for ix = (1:NSIDE)
			for iy = (1:NSIDE)
				p_start = Point{Float64}((ix-1) * dix , (iy-1) * diy)
				p_final = Point{Float64}(ix * dix, iy * diy)
				ID_BOX = ( (ix-1) + ( (iy-1) * NSIDE ) ) + 1

				#; Densita` Uniforme
				#density    = 1/(NSIDE*NSIDE)

				#; Densita` Pesata
				density    = w_m[ID_BOX]

				hit_in_box = Int64(floor(density * (cal/2.)))
				control    += hit_in_box
				#println("Hit = ", hit_in_box)
				if hit_in_box != 0
					#; Qui va costruita la nuova densita` di punti
					pt_2 = [Point{Float64}((p_final.x-p_start.x)*rand(rng)+p_start.x, (p_final.y-p_start.y)*rand(rng)+p_start.y) for n = 1:hit_in_box]

					#; Cambio Coordinate
					p_2 = ( pt_2 .* Ref(xu-xl) ) .+ Ref(xl)

					#; Valuto l'integrale
					f_2 = [fun(point) * (xu.x-xl.x) * (xu.y-xl.y) / density   for point in p_2 ]

					#println(p_start.x, " ", p_final.x, " ", p_start.y," ", p_final.y, " ", sum(f_2)/hit_in_box)

					append!(pt_2_all, pt_2)
					append!(f_2_all, f_2)

					append!(Box_X, ix)
					append!(Box_Y, iy)

					#println("Box = ", ix, " ", iy)

				else
					continue
				end
			end
		end

		# xu = Point{Float64}(x_max, y_max)
		# xl = Point{Float64}(x_min, y_min)

		# println( xmin, " ", xmax, " ", ymin, " ", ymax )

		pt = pt_2_all

		Integral += (sum(f_2_all) / (control*NSIDE*NSIDE) )
		Integral_abs = sum(abs.(f_2_all))
		g = (abs.(f_2_all) ./ Integral_abs)
		w_m     = zeros(Int64(NSIDE^2))

		println(iter,". Integrale = ", (Integral / (iter + 1) ) - Int_Vero, " ", Integral/(iter + 1) )

		prova = (pt_2_all .* Ref(xu-xl) ) .+ Ref(xl)

		for i = 1:length(g)
			xx=append!(xx, prova[i].x)
			yy=append!(yy, prova[i].y)
		end

		# println("Per le X")
		# println(      ((findmin(Box_X)[1]-1)/100 * (xu.x-xl.x)) + xl.x   )
		# println(      ((findmax(Box_X)[1])/100 * (xu.x-xl.x)) + xl.x  )
		#
		# println("Per le Y")
		# println(((findmin(Box_Y)[1]-1)/100 * (xu.x-xl.x)) + xl.x)
		# println(((findmax(Box_Y)[1])/100 * (xu.x-xl.x)) + xl.x )
		#
		# # xu_n = Point{Float64}(((findmax(Box_X)[1])/100 * (xu.x-xl.x)) + xl.x, ((findmax(Box_Y)[1])/100 * (xu.x-xl.x)) + xl.x     )
		# # xl_n = Point{Float64}(((findmin(Box_X)[1])/100 * (xu.x-xl.x)) + xl.x, ((findmin(Box_Y)[1])/100 * (xu.x-xl.x)) + xl.x     )
		# #
		# # xu = xu_n
		# # xl = xl_n
		#
		# println(xu)
		# Box_X = Array{Float64}(undef, 0)
		# Box_Y = Array{Float64}(undef, 0)

	end #end iteration

	return Integral,xx,yy
end

#################### Integration #####################

calls = 100000

xu = Point{Float64}(1., 1.)
xl = Point{Float64}(-1.,-1.)

result,x,y = Integrate(cilindro, xl, xu, calls)
plot(x,y, seriestype=:scatter)

#println("Valore dell'integrale = ", result)

# il punto e` che gli estremi di integrazione devono essere rimodellati
# dove la funzione non e` diversa da zero!

max_stars_repo_name: gruffoc/AFO

max_issues_repo_name: gruffoc/AFO

max_forks_repo_name: gruffoc/AFO


text:
using stoosc
rc("font", family="arial")

MODEL = "model01"

if MODEL == "model01"
    I = "0.02"  # 0.02, 0.05, 0.1
    σ_string_arr = ["0.001", "0.003", "0.01", "0.03"]
    title_string = "Reentrainment after jet lag (I = $I)"
elseif MODEL == "model02"
    I = "0.5"
    σ_string_arr = ["0.02", "0.1", "0.15", "0.2"]
    title_string = "Limit cycle osillator (I = $I)"
elseif MODEL == "model02b"
    I = "0.5"
    σ_string_arr = ["0.02", "0.1", "0.15", "0.2"]
    title_string = "Noise-induced osillator (I = $I)"
end

phase_diff_mean_arr = []
phase_diff_std_arr = []
pre_phase_diff_mean_arr = []
pre_phase_diff_std_arr = []
for σ_string in σ_string_arr
    global phase_diff_mean_arr, phase_diff_std_arr
    local filename

    xarr = []
    pre_xarr = []

    for i = 1:10
        filename = "sde_I=$(I)_σ=$(σ_string)_r=$(i).csv"
        fullfilename = joinpath("jet_lags", MODEL, filename)
        df = load_data(fullfilename)

        x = (mean(df[:, "pre_phases"]) .- df[:, "post_phases"])
        push!(xarr, x)

        pre_x = (mean(df[:, "pre_phases"]) .- df[:, "pre_phases"])
        push!(pre_xarr, pre_x)
    end

    phase_diff_mean = mean(xarr)
    phase_diff_std = std(xarr)

    push!(phase_diff_mean_arr, phase_diff_mean)
    push!(phase_diff_std_arr, phase_diff_std)

    pre_phase_diff_mean = mean(pre_xarr)
    pre_phase_diff_std = std(pre_xarr)

    push!(pre_phase_diff_mean_arr, pre_phase_diff_mean)
    push!(pre_phase_diff_std_arr, pre_phase_diff_std)

end

colors = ["gray", "black", "red", "blue"]
fig, ax = subplots(figsize=(6.5, 2.5))
pre_max = 5
["σ = 0.03", "σ = 0.01", "σ = 0.003", "σ = 0.001"]
for i = 1:length(σ_string_arr)
    ax.errorbar(-pre_max:-1, pre_phase_diff_mean_arr[i][21-pre_max:end], yerr=pre_phase_diff_std_arr[i][21-pre_max:end], fmt="o", color=colors[i], label="σ = $(σ_string_arr[i])")
    ax.errorbar(1:20, phase_diff_mean_arr[i], yerr=phase_diff_std_arr[i], fmt="o", color=colors[i])
end
ax.legend(edgecolor="black", framealpha=1.0, ncol=1, loc=2)
ax.plot([-pre_max-.5, 20+0.5], [0, 0], color="black", linestyle="--")
ax.plot([0, 0], [-0.1, 0.5], color="black", linestyle="--")
ax.set_xlabel("Days")
ax.set_ylabel("Phase difference (days)")
ax.set_title(title_string, pad=0.1)
ax.set_xlim(-5.5, 20.5)
ax.set_ylim(-0.1, 0.5)
fig.tight_layout()

filename = joinpath("jet_lags", MODEL, "jet_lag_I=$(I).svg")
save_figure(fig, filename)

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
##
# For Chebyshev T. Note the shift in indexing is fine due to the AbstractFill
##
Base.@propagate_inbounds _forwardrecurrence_next(n, A::Vcat{<:Any,1,<:Tuple{<:Number,<:AbstractFill}}, B::Zeros, C::Ones, x, p0, p1) =
    _forwardrecurrence_next(n, A.args[2], B, C, x, p0, p1)

Base.@propagate_inbounds _clenshaw_next(n, A::Vcat{<:Any,1,<:Tuple{<:Number,<:AbstractFill}}, B::Zeros, C::Ones, x, c, bn1, bn2) =
    _clenshaw_next(n, A.args[2], B, C, x, c, bn1, bn2)

# Assume 1 normalization
_p0(A) = one(eltype(A))

function initiateforwardrecurrence(N, A, B, C, x, μ)
    T = promote_type(eltype(A), eltype(B), eltype(C), typeof(x))
    p0 = convert(T, μ)
    p1 = convert(T, muladd(A[1],x,B[1])*p0)
    @inbounds for n = 2:N
        p1,p0 = _forwardrecurrence_next(n, A, B, C, x, p0, p1),p1
    end
    p0,p1
end

for (get, vie) in ((:getindex, :view), (:(Base.unsafe_getindex), :(Base.unsafe_view)))
    @eval begin
        Base.@propagate_inbounds @inline $get(P::OrthogonalPolynomial{T}, x::Number, n::OneTo) where T =
            copyto!(Vector{T}(undef,length(n)), $vie(P, x, n))

        $get(P::OrthogonalPolynomial{T}, x::AbstractVector, n::AbstractUnitRange{Int}) where T =
            copyto!(Matrix{T}(undef,length(x),length(n)), $vie(P, x, n))
    end
end

function copyto!(dest::AbstractVector, V::SubArray{<:Any,1,<:OrthogonalPolynomial,<:Tuple{<:Number,<:OneTo}})
    P = parent(V)
    x,n = parentindices(V)
    A,B,C = recurrencecoefficients(P)
    forwardrecurrence!(dest, A, B, C, x, _p0(P))
end

function forwardrecurrence_copyto!(dest::AbstractMatrix, V)
    checkbounds(dest, axes(V)...)
    P = parent(V)
    xr,jr = parentindices(V)
    A,B,C = recurrencecoefficients(P)
    shift = first(jr)
    Ã,B̃,C̃ = A[shift:∞],B[shift:∞],C[shift:∞]
    for (k,x) = enumerate(xr)
        p0, p1 = initiateforwardrecurrence(shift, A, B, C, x, _p0(P))
        _forwardrecurrence!(view(dest,k,:), Ã, B̃, C̃, x, p0, p1)
    end
    dest
end
copyto!(dest::AbstractMatrix, V::SubArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{<:AbstractVector,<:AbstractUnitRange}}) = forwardrecurrence_copyto!(dest, V)
copyto!(dest::LayoutMatrix, V::SubArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{<:AbstractVector,<:AbstractUnitRange}}) = forwardrecurrence_copyto!(dest, V)

function copyto!(dest::AbstractVector, V::SubArray{<:Any,1,<:OrthogonalPolynomial,<:Tuple{<:Number,<:UnitRange}})
    checkbounds(dest, axes(V)...)
    P = parent(V)
    x,jr = parentindices(V)
    A,B,C = recurrencecoefficients(P)
    shift = first(jr)
    Ã,B̃,C̃ = A[shift:∞],B[shift:∞],C[shift:∞]
    p0, p1 = initiateforwardrecurrence(shift, A, B, C, x, _p0(P))
    _forwardrecurrence!(dest, Ã, B̃, C̃, x, p0, p1)
    dest
end

getindex(P::OrthogonalPolynomial, x::Number, n::AbstractVector) = layout_getindex(P, x, n)
getindex(P::OrthogonalPolynomial, x::AbstractVector, n::AbstractVector) = layout_getindex(P, x, n)
getindex(P::SubArray{<:Any,1,<:OrthogonalPolynomial}, x::AbstractVector) = layout_getindex(P, x)
getindex(P::OrthogonalPolynomial, x::Number, n::Number) = P[x,oneto(n)][end]

unsafe_layout_getindex(A...) = sub_materialize(Base.unsafe_view(A...))

Base.unsafe_getindex(P::OrthogonalPolynomial, x::Number, n::AbstractUnitRange) = unsafe_layout_getindex(P, x, n)
Base.unsafe_getindex(P::OrthogonalPolynomial, x::AbstractVector, n::AbstractUnitRange) = unsafe_layout_getindex(P, x, n)
Base.unsafe_getindex(P::OrthogonalPolynomial, x::Number, n::AbstractVector) = Base.unsafe_getindex(P,x,oneto(maximum(n)))[n]
Base.unsafe_getindex(P::OrthogonalPolynomial, x::AbstractVector, n::AbstractVector) = Base.unsafe_getindex(P,x,oneto(maximum(n)))[:,n]
Base.unsafe_getindex(P::OrthogonalPolynomial, x::AbstractVector, n::Number) = Base.unsafe_getindex(P, x, 1:n)[:,end]
Base.unsafe_getindex(P::OrthogonalPolynomial, x::Number, ::Colon) = Base.unsafe_getindex(P, x, axes(P,2))
Base.unsafe_getindex(P::OrthogonalPolynomial, x::Number, n::Number) = Base.unsafe_getindex(P,x,oneto(n))[end]

getindex(P::OrthogonalPolynomial, x::Number, jr::AbstractInfUnitRange{Int}) = view(P, x, jr)
getindex(P::OrthogonalPolynomial, x::AbstractVector, jr::AbstractInfUnitRange{Int}) = view(P, x, jr)
Base.unsafe_getindex(P::OrthogonalPolynomial{T}, x::Number, jr::AbstractInfUnitRange{Int}) where T =
    BroadcastVector{T}(Base.unsafe_getindex, Ref(P), x, jr)

function getindex(P::BroadcastVector{<:Any,typeof(Base.unsafe_getindex), <:Tuple{Ref{<:OrthogonalPolynomial},Number,AbstractVector{Int}}}, kr::AbstractVector)
    Pr, x, jr = P.args
    Base.unsafe_getindex(Pr[], x, jr[kr])
end

###
# Clenshaw
###

function unsafe_getindex(f::Expansion{<:Any,<:OrthogonalPolynomial}, x::Number)
    P,c = arguments(f)
    _p0(P)*clenshaw(paddeddata(c), recurrencecoefficients(P)..., x)
end

Base.@propagate_inbounds function getindex(f::Expansion{<:Any,<:OrthogonalPolynomial}, x::Number)
    @inbounds checkbounds(f, x)
    unsafe_getindex(f, x)
end

getindex(f::Expansion{T,<:OrthogonalPolynomial}, x::AbstractVector{<:Number}) where T =
    copyto!(Vector{T}(undef, length(x)), view(f, x))

function copyto!(dest::AbstractVector{T}, v::SubArray{<:Any,1,<:Expansion{<:Any,<:OrthogonalPolynomial}, <:Tuple{AbstractVector{<:Number}}}) where T
    f = parent(v)
    (x,) = parentindices(v)
    P,c = arguments(f)
    clenshaw!(paddeddata(c), recurrencecoefficients(P)..., x, Fill(_p0(P), length(x)), dest)
end

###
# Operator clenshaw
###

Base.@propagate_inbounds function _clenshaw_next!(n, A::AbstractFill, ::Zeros, C::Ones, x::AbstractMatrix, c, bn1::AbstractMatrix{T}, bn2::AbstractMatrix{T}) where T
    muladd!(getindex_value(A), x, bn1, -one(T), bn2)
    view(bn2,band(0)) .+= c[n]
    bn2
end

Base.@propagate_inbounds function _clenshaw_next!(n, A::AbstractVector, ::Zeros, C::AbstractVector, x::AbstractMatrix, c, bn1::AbstractMatrix{T}, bn2::AbstractMatrix{T}) where T
    muladd!(A[n], x, bn1, -C[n+1], bn2)
    view(bn2,band(0)) .+= c[n]
    bn2
end

Base.@propagate_inbounds function _clenshaw_next!(n, A::AbstractVector, B::AbstractVector, C::AbstractVector, x::AbstractMatrix, c, bn1::AbstractMatrix{T}, bn2::AbstractMatrix{T}) where T
    bn2 .= B[n] .* bn1 .- C[n+1] .* bn2
    muladd!(A[n], x, bn1, one(T), bn2)
    view(bn2,band(0)) .+= c[n]
    bn2
end

# Operator * f Clenshaw
Base.@propagate_inbounds function _clenshaw_next!(n, A::AbstractFill, ::Zeros, C::Ones, X::AbstractMatrix, c, f::AbstractVector, bn1::AbstractVector{T}, bn2::AbstractVector{T}) where T
    muladd!(getindex_value(A), X, bn1, -one(T), bn2)
    bn2 .+= c[n] .* f
    bn2
end

Base.@propagate_inbounds function _clenshaw_next!(n, A, ::Zeros, C, X::AbstractMatrix, c, f::AbstractVector, bn1::AbstractVector{T}, bn2::AbstractVector{T}) where T
    muladd!(A[n], X, bn1, -C[n+1], bn2)
    bn2 .+= c[n] .* f
    bn2
end

Base.@propagate_inbounds function _clenshaw_next!(n, A, B, C, X::AbstractMatrix, c, f::AbstractVector, bn1::AbstractVector{T}, bn2::AbstractVector{T}) where T
    bn2 .= B[n] .* bn1 .- C[n+1] .* bn2 .+ c[n] .* f
    muladd!(A[n], X, bn1, one(T), bn2)
    bn2
end

# allow special casing first arg, for ChebyshevT in ClassicalOrthogonalPolynomials
Base.@propagate_inbounds function _clenshaw_first!(A, ::Zeros, C, X, c, bn1, bn2)
    muladd!(A[1], X, bn1, -C[2], bn2)
    view(bn2,band(0)) .+= c[1]
    bn2
end

Base.@propagate_inbounds function _clenshaw_first!(A, B, C, X, c, bn1, bn2)
    bn2 .= B[1] .* bn1 .- C[2] .* bn2
    muladd!(A[1], X, bn1, one(eltype(bn2)), bn2)
    view(bn2,band(0)) .+= c[1]
    bn2
end

Base.@propagate_inbounds function _clenshaw_first!(A, ::Zeros, C, X, c, f::AbstractVector, bn1, bn2)
    muladd!(A[1], X, bn1, -C[2], bn2)
    bn2 .+= c[1] .* f
    bn2
end

Base.@propagate_inbounds function _clenshaw_first!(A, B, C, X, c, f::AbstractVector, bn1, bn2)
    bn2 .= B[1] .* bn1 .- C[2] .* bn2 .+ c[1] .* f
    muladd!(A[1], X, bn1, one(eltype(bn2)), bn2)
    bn2
end

_clenshaw_op(::AbstractBandedLayout, Z, N) = BandedMatrix(Z, (N-1,N-1))

function clenshaw(c::AbstractVector, A::AbstractVector, B::AbstractVector, C::AbstractVector, X::AbstractMatrix)
    N = length(c)
    T = promote_type(eltype(c),eltype(A),eltype(B),eltype(C),eltype(X))
    @boundscheck check_clenshaw_recurrences(N, A, B, C)
    m = size(X,1)
    m == size(X,2) || throw(DimensionMismatch("X must be square"))
    N == 0 && return zero(T)
    bn2 = _clenshaw_op(MemoryLayout(X), Zeros{T}(m, m), N)
    bn1 = _clenshaw_op(MemoryLayout(X), c[N]*Eye{T}(m), N)
    _clenshaw_op!(c, A, B, C, X, bn1, bn2)
end

function clenshaw(c::AbstractVector, A::AbstractVector, B::AbstractVector, C::AbstractVector, X::AbstractMatrix, f::AbstractVector)
    N = length(c)
    T = promote_type(eltype(c),eltype(A),eltype(B),eltype(C),eltype(X))
    @boundscheck check_clenshaw_recurrences(N, A, B, C)
    m = size(X,1)
    m == size(X,2) || throw(DimensionMismatch("X must be square"))
    m == length(f) || throw(DimensionMismatch("Dimensions must match"))
    N == 0 && return zero(T)
    bn2 = zeros(T,m)
    bn1 = Vector{T}(undef,m)
    bn1 .= c[N] .* f
    _clenshaw_op!(c, A, B, C, X, f, bn1, bn2)
end

function _clenshaw_op!(c, A, B, C, X, bn1, bn2)
    N = length(c)
    N == 1 && return bn1
    @inbounds begin
        for n = N-1:-1:2
            bn1,bn2 = _clenshaw_next!(n, A, B, C, X, c, bn1, bn2),bn1
        end
        bn1 = _clenshaw_first!(A, B, C, X, c, bn1, bn2)
    end
    bn1
end

function _clenshaw_op!(c, A, B, C, X, f::AbstractVector, bn1, bn2)
    N = length(c)
    N == 1 && return bn1
    @inbounds begin
        for n = N-1:-1:2
            bn1,bn2 = _clenshaw_next!(n, A, B, C, X, c, f, bn1, bn2),bn1
        end
        bn1 = _clenshaw_first!(A, B, C, X, c, f, bn1, bn2)
    end
    bn1
end

"""
    Clenshaw(a, X)

represents the operator `a(X)` where a is a polynomial.
Here `a` is to stored as a quasi-vector.
"""
struct Clenshaw{T, Coefs<:AbstractVector, AA<:AbstractVector, BB<:AbstractVector, CC<:AbstractVector, Jac<:AbstractMatrix} <: AbstractBandedMatrix{T}
    c::Coefs
    A::AA
    B::BB
    C::CC
    X::Jac
    p0::T
end

Clenshaw(c::AbstractVector{T}, A::AbstractVector, B::AbstractVector, C::AbstractVector, X::AbstractMatrix{T}, p0::T) where T =
    Clenshaw{T,typeof(c),typeof(A),typeof(B),typeof(C),typeof(X)}(c, A, B, C, X, p0)

Clenshaw(c::Number, A, B, C, X, p) = Clenshaw([c], A, B, C, X, p)

function Clenshaw(a::AbstractQuasiVector, X::AbstractQuasiMatrix)
    P,c = arguments(a)
    Clenshaw(paddeddata(c), recurrencecoefficients(P)..., jacobimatrix(X), _p0(P))
end

copy(M::Clenshaw) = M
size(M::Clenshaw) = size(M.X)
axes(M::Clenshaw) = axes(M.X)
bandwidths(M::Clenshaw) = (length(M.c)-1,length(M.c)-1)

Base.array_summary(io::IO, C::Clenshaw{T}, inds::Tuple{Vararg{OneToInf{Int}}}) where T =
    print(io, Base.dims2string(length.(inds)), " Clenshaw{$T} with $(length(C.c)) degree polynomial")

struct ClenshawLayout <: AbstractLazyBandedLayout end
MemoryLayout(::Type{<:Clenshaw}) = ClenshawLayout()
sublayout(::ClenshawLayout, ::Type{<:NTuple{2,AbstractUnitRange{Int}}}) = ClenshawLayout()
sublayout(::ClenshawLayout, ::Type{<:Tuple{AbstractUnitRange{Int},Union{Slice,AbstractInfUnitRange{Int}}}}) = LazyBandedLayout()
sublayout(::ClenshawLayout, ::Type{<:Tuple{Union{Slice,AbstractInfUnitRange{Int}},AbstractUnitRange{Int}}}) = LazyBandedLayout()
sublayout(::ClenshawLayout, ::Type{<:Tuple{Union{Slice,AbstractInfUnitRange{Int}},Union{Slice,AbstractInfUnitRange{Int}}}}) = LazyBandedLayout()
sub_materialize(::ClenshawLayout, V) = BandedMatrix(V)

function _BandedMatrix(::ClenshawLayout, V::SubArray{<:Any,2})
    M = parent(V)
    kr,jr = parentindices(V)
    b = bandwidth(M,1)
    jkr=max(1,min(jr[1],kr[1])-b÷2):max(jr[end],kr[end])+b÷2
    # relationship between jkr and kr, jr
    kr2,jr2 = kr.-jkr[1].+1,jr.-jkr[1].+1
    lmul!(M.p0, clenshaw(M.c, M.A, M.B, M.C, M.X[jkr, jkr])[kr2,jr2])
end

function getindex(M::Clenshaw{T}, kr::AbstractUnitRange, j::Integer) where T
    b = bandwidth(M,1)
    jkr=max(1,min(j,kr[1])-b÷2):max(j,kr[end])+b÷2
    # relationship between jkr and kr, jr
    kr2,j2 = kr.-jkr[1].+1,j-jkr[1]+1
    f = [Zeros{T}(j2-1); one(T); Zeros{T}(length(jkr)-j2)]
    lmul!(M.p0, clenshaw(M.c, M.A, M.B, M.C, M.X[jkr, jkr], f)[kr2])
end

getindex(M::Clenshaw, k::Int, j::Int) = M[k:k,j][1]

transposelayout(M::ClenshawLayout) = LazyBandedMatrices.LazyBandedLayout()
# TODO: generalise for layout, use Base.PermutedDimsArray
Base.permutedims(M::Clenshaw{<:Number}) = transpose(M)

function materialize!(M::MatMulVecAdd{<:ClenshawLayout,<:PaddedLayout,<:PaddedLayout})
    α,A,x,β,y = M.α,M.A,M.B,M.β,M.C
    length(y) == size(A,1) || throw(DimensionMismatch("Dimensions must match"))
    length(x) == size(A,2) || throw(DimensionMismatch("Dimensions must match"))
    x̃ = paddeddata(x);
    m = length(x̃)
    b = bandwidth(A,1)
    jkr=1:m+b
    p = [x̃; zeros(eltype(x̃),length(jkr)-m)];
    Ax = lmul!(A.p0, clenshaw(A.c, A.A, A.B, A.C, A.X[jkr, jkr], p))
    _fill_lmul!(β,y)
    resizedata!(y, last(jkr))
    v = view(paddeddata(y),jkr)
    BLAS.axpy!(α, Ax, v)
    y
end

function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), a::Expansion{<:Any,<:OrthogonalPolynomial}, P::OrthogonalPolynomial)
    axes(a,1) == axes(P,1) || throw(DimensionMismatch())
    P * Clenshaw(a, P)
end

function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), a::Expansion{<:Any,<:WeightedOrthogonalPolynomial}, P::OrthogonalPolynomial)
    axes(a,1) == axes(P,1) || throw(DimensionMismatch())
    wQ,c = arguments(a)
    w,Q = arguments(wQ)
    (w .* P) * Clenshaw(Q * c, P)
end

function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), wv::BroadcastQuasiVector{<:Any,typeof(*),<:Tuple{Weight,AffineQuasiVector}}, P::OrthogonalPolynomial)
    w,v = arguments(wv)
    Q = OrthogonalPolynomial(w)
    a = (w .* Q) * (Q \ v)
    a .* P
end

##
# Banded dot is slow
###

LinearAlgebra.dot(x::AbstractVector, A::Clenshaw, y::AbstractVector) = dot(x, mul(A, y))


text:
struct CapacityExpansionData{T <: AbstractFloat} <: AbstractModelData
    hydrodata::HydroPlantCollection{T,2}
    resolution::Resolution
    investment_levels::Vector{T}
    expansion_cost::T
    discount_rate::T
    investment_period::Int

    function CapacityExpansionData(plantdata::HydroPlantCollection{T,2},
                                   investment_levels::Vector{T};
                                   hours_in_period::Integer = 24,
                                   expansion_cost::T = 0.79,
                                   discount_rate::T = 0.05,
                                   investment_period::Integer = 40) where T <: AbstractFloat
        return new{T}(plantdata, Resolution(hours_in_period), investment_levels, expansion_cost, discount_rate, investment_period)
    end
end

function CapacityExpansionData(plantfilename::String, investment_levels::AbstractVector; kw...)
    return CapacityExpansionData(HydroPlantCollection(plantfilename), investment_levels; kw...)
end

function CapacityExpansionData(plantfilename::String; kw...)
    investment_levels = collect(0.0:10:100.0)
    return CapacityExpansionData(HydroPlantCollection(plantfilename), investment_levels; kw...)
end

function equivalent_cost(data::CapacityExpansionData, horizon::Horizon)
    # C = data.investment_levels[level_index] * data.expansion_cost * 1e6
    # if iszero(C)
    #     return C
    # end
    r = data.discount_rate
    P = data.investment_period
    T = num_days(horizon)
    # Calculate equivalent rate
    rᴱ = (1 + r)^(T/365) - 1
    # Calculate equivalent cost
    Cᴱ = data.expansion_cost * rᴱ / (1 - (1 + rᴱ)^(-P * 365 / T))
    # Return equivalent cost
    return Cᴱ
end


text:
using DataFrames, DataFramesMeta, CSV
using PGFPlots
using Statistics
using DelimitedFiles

# Create a table of log-evidence x subject
# in order to perform a random effects analysis (Rigoux et al., 2014; Stephan et al., 2009)
# in MATLAB
# Like the one that now is in someresults/KlaasAnalysis/
function createklaastable_crossvalLL(;
        filepath =  joinpath(dirname(pathof(SinergiafMRI_datafit)),"mcmc_rl_fit/projects/fmri/someresults/crossvalidation/"),
        filestr = "",
        )

    savepath = makehomesavepath("createklaastable")
    mkdir(savepath)

    nSubj = 21
    prefix = split(filepath, "/")[end-1]
    runs = getsubfolders(filepath)
    algorithms = getsubfolders(joinpath(filepath, runs[1]))

    LLallalgos = []
    LLalgo = []
    LLallalgosmean = zeros(nSubj, length(algorithms))
    LLallalgosstd = zeros(nSubj, length(algorithms))
    LLallalgosstderror = zeros(nSubj, length(algorithms))
    for i in 1:length(algorithms)

        for irun in 1:length(runs)

            if irun == 1 # Initialize
                 LLalgo = zeros(length(runs), nSubj)
            end

            dfalgo = CSV.read(joinpath(filepath, runs[irun], algorithms[i],
                    "cross_validation/LL_out_subjectwise.csv"), header=false)
            LLalgo[irun, :] = [sum(dfalgo[:,j]) for j in 1:size(dfalgo,2)]

        end
        push!(LLallalgos, LLalgo)
        LLallalgosmean[:, i] = mean(LLalgo, dims=1)
        LLallalgosstd[:, i] = std(LLalgo, dims=1)
        LLallalgosstderror[:, i] = std(LLalgo, dims=1) ./sqrt(length(runs))
    end
    open(joinpath(savepath, "LLmean_algoxsubj.csv"), "a") do io
       writedlm(io, transpose(LLallalgosmean))
       end
    open(joinpath(savepath, "LLstd_algoxsubj.csv"), "a") do io
       writedlm(io, transpose(LLallalgosstd))
       end
    open(joinpath(savepath, "LLstderror_algoxsubj.csv"), "a") do io
        writedlm(io, transpose(LLallalgosstderror))
        end
    open(joinpath(savepath, "algolist.csv"), "a") do io
       writedlm(io, algorithms)
       end

    LLallalgosmean, LLallalgosstd, LLallalgosstderror# LLallalgos
end