text:
using OrdinaryDiffEq, StochasticDiffEq, DelayDiffEq, Flux, DiffEqFlux,
      Zygote, Test, DiffEqSensitivity

mp = Float32[0.1,0.1]
x = Float32[2.; 0.]
xs = Float32.(hcat([0.; 0.], [1.; 0.], [2.; 0.]))
tspan = (0.0f0,1.0f0)
dudt = Chain(Dense(2,50,tanh),Dense(50,2))
Slowdudt = SlowChain(SlowDense(2,50,tanh),SlowDense(50,2))
staticdudt = SlowChain(StaticDense(2,50,tanh),StaticDense(50,2))

NeuralODE(dudt,tspan,Tsit5(),save_everystep=false,save_start=false)(x)
NeuralODE(dudt,tspan,Tsit5(),saveat=0.1)(x)
NeuralODE(dudt,tspan,Tsit5(),saveat=0.1,sensealg=TrackerAdjoint())(x)

NeuralODE(dudt,tspan,Tsit5(),save_everystep=false,save_start=false)(xs)
NeuralODE(dudt,tspan,Tsit5(),saveat=0.1)(xs)
NeuralODE(dudt,tspan,Tsit5(),saveat=0.1,sensealg=TrackerAdjoint())(xs)

@info "Test some gradients"

node = NeuralODE(dudt,tspan,Tsit5(),save_everystep=false,save_start=false)
grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
@test ! iszero(grads[x])
@test ! iszero(grads[node.p])

grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
@test ! iszero(grads[xs])
@test ! iszero(grads[node.p])

node = NeuralODE(dudt,tspan,Tsit5(),save_everystep=false,save_start=false,sensealg=TrackerAdjoint())
grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
@test ! iszero(grads[x])
@test ! iszero(grads[node.p])

grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
@test ! iszero(grads[xs])
@test ! iszero(grads[node.p])

node = NeuralODE(dudt,tspan,Tsit5(),save_everystep=false,save_start=false,sensealg=BacksolveAdjoint())
grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
@test ! iszero(grads[x])
@test ! iszero(grads[node.p])

grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
@test ! iszero(grads[xs])
@test ! iszero(grads[node.p])

node = NeuralODE(dudt,tspan,Tsit5(),save_everystep=false,save_start=false,sensealg=BacksolveAdjoint(autojacvec=ZygoteVJP()))
grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
@test ! iszero(grads[x])
@test ! iszero(grads[node.p])

grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
@test ! iszero(grads[xs])
@test ! iszero(grads[node.p])

## Slow

@info "Test some Slow layers"

node = NeuralODE(Slowdudt,tspan,Tsit5(),save_everystep=false,save_start=false)
grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
@test ! iszero(grads[x])
@test ! iszero(grads[node.p])

grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
@test ! iszero(grads[xs])
@test ! iszero(grads[node.p])

node = NeuralODE(Slowdudt,tspan,Tsit5(),save_everystep=false,save_start=false,sensealg=TrackerAdjoint())
grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
@test ! iszero(grads[x])
@test ! iszero(grads[node.p])

grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
@test ! iszero(grads[xs])
@test ! iszero(grads[node.p])

goodgrad = grads[node.p]
p = node.p

node = NeuralODE(Slowdudt,tspan,Tsit5(),save_everystep=false,save_start=false,sensealg=BacksolveAdjoint(),p=p)
grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
@test ! iszero(grads[x])
@test ! iszero(grads[node.p])

grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
@test ! iszero(grads[xs])
@test ! iszero(grads[node.p])
goodgrad2 = grads[node.p]
@test goodgrad ≈ goodgrad2 # Make sure adjoint overloads are correct

node = NeuralODE(staticdudt,tspan,Tsit5(),save_everystep=false,save_start=false,sensealg=BacksolveAdjoint(),p=p)
grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
@test ! iszero(grads[x])
@test ! iszero(grads[node.p])

@test_throws ErrorException grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))

@info "Test some adjoints"

# Adjoint
@testset "adjoint mode" begin
    node = NeuralODE(dudt,tspan,Tsit5(),save_everystep=false,save_start=false)
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
    @test ! iszero(grads[xs])
    @test ! iszero(grads[node.p])

    node = NeuralODE(dudt,tspan,Tsit5(),saveat=0.0:0.1:1.0)
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
    @test ! iszero(grads[xs])
    @test ! iszero(grads[node.p])

    node = NeuralODE(dudt,tspan,Tsit5(),saveat=0.1)
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
    @test ! iszero(grads[xs])
    @test ! iszero(grads[node.p])

    node = NeuralODE(Slowdudt,tspan,Tsit5(),save_everystep=false,save_start=false)
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    @test_broken grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node)) isa Tuple
    #@test ! iszero(grads[xs])
    #@test ! iszero(grads[node.p])

    node = NeuralODE(Slowdudt,tspan,Tsit5(),saveat=0.0:0.1:1.0)
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    @test_broken grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node)) isa Tuple
    #@test ! iszero(grads[xs])
    #@test ! iszero(grads[node.p])

    node = NeuralODE(Slowdudt,tspan,Tsit5(),saveat=0.1)
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    @test_broken grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node)) isa Tuple
    #@test ! iszero(grads[xs])
    #@test ! iszero(grads[node.p])
end

@info "Test Tracker"

# RD
@testset "Tracker mode" begin
    node = NeuralODE(dudt,tspan,Tsit5(),save_everystep=false,save_start=false,sensealg=TrackerAdjoint())
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
    @test ! iszero(grads[xs])
    @test ! iszero(grads[node.p])

    NeuralODE(dudt,tspan,Tsit5(),saveat=0.0:0.1:1.0,sensealg=TrackerAdjoint())
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
    @test ! iszero(grads[xs])
    @test ! iszero(grads[node.p])

    node = NeuralODE(dudt,tspan,Tsit5(),saveat=0.1,sensealg=TrackerAdjoint())
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
    @test ! iszero(grads[xs])
    @test ! iszero(grads[node.p])

    node = NeuralODE(Slowdudt,tspan,Tsit5(),save_everystep=false,save_start=false,sensealg=TrackerAdjoint())
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    @test_broken grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node)) isa Tuple
    #@test ! iszero(grads[xs])
    #@test ! iszero(grads[node.p])

    NeuralODE(Slowdudt,tspan,Tsit5(),saveat=0.0:0.1:1.0,sensealg=TrackerAdjoint())
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    @test_broken grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node)) isa Tuple
    #@test ! iszero(grads[xs])
    #@test ! iszero(grads[node.p])

    node = NeuralODE(Slowdudt,tspan,Tsit5(),saveat=0.1,sensealg=TrackerAdjoint())
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    @test_broken grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node)) isa Tuple
    #@test ! iszero(grads[xs])
    #@test ! iszero(grads[node.p])
end

@info "Test non-ODEs"

dudt2 = Chain(Dense(2,50,tanh),Dense(50,2))
Slowdudt2 = SlowChain(SlowDense(2,50,tanh),SlowDense(50,2))
NeuralDSDE(dudt,dudt2,(0.0f0,.1f0),SOSRI(),saveat=0.1)(x)
sode = NeuralDSDE(dudt,dudt2,(0.0f0,.1f0),SOSRI(),saveat=0.0:0.01:0.1)

grads = Zygote.gradient(()->sum(sode(x)),Flux.params(x,sode))
@test ! iszero(grads[x])
@test ! iszero(grads[sode.p])
@test ! iszero(grads[sode.p][end])

grads = Zygote.gradient(()->sum(sode(xs)),Flux.params(xs,sode))
@test ! iszero(grads[xs])
@test ! iszero(grads[sode.p])
@test ! iszero(grads[sode.p][end])

sode = NeuralDSDE(Slowdudt,Slowdudt2,(0.0f0,.1f0),SOSRI(),saveat=0.0:0.01:0.1)

grads = Zygote.gradient(()->sum(sode(x)),Flux.params(x,sode))
@test ! iszero(grads[x])
@test ! iszero(grads[sode.p])
@test ! iszero(grads[sode.p][end])

@test_broken grads = Zygote.gradient(()->sum(sode(xs)),Flux.params(xs,sode)) isa Tuple
#@test ! iszero(grads[xs])
#@test ! iszero(grads[sode.p])
#@test ! iszero(grads[sode.p][end])

dudt22 = Chain(Dense(2,50,tanh),Dense(50,4),x->reshape(x,2,2))
Slowdudt22 = SlowChain(SlowDense(2,50,tanh),SlowDense(50,4),(x,p)->reshape(x,2,2))
NeuralSDE(dudt,dudt22,(0.0f0,.1f0),2,LambaEM(),saveat=0.01)(x)

sode = NeuralSDE(dudt,dudt22,(0.0f0,0.1f0),2,LambaEM(),saveat=0.0:0.01:0.1)

grads = Zygote.gradient(()->sum(sode(x)),Flux.params(x,sode))
@test ! iszero(grads[x])
@test ! iszero(grads[sode.p])
@test ! iszero(grads[sode.p][end])

@test_broken grads = Zygote.gradient(()->sum(sode(xs)),Flux.params(xs,sode))
#@test ! iszero(grads[xs])
#@test ! iszero(grads[sode.p])
#@test ! iszero(grads[sode.p][end])

sode = NeuralSDE(Slowdudt,Slowdudt22,(0.0f0,0.1f0),2,LambaEM(),saveat=0.0:0.01:0.1)

grads = Zygote.gradient(()->sum(sode(x)),Flux.params(x,sode))
@test ! iszero(grads[x])
@test ! iszero(grads[sode.p])
@test ! iszero(grads[sode.p][end])

@test_broken grads = Zygote.gradient(()->sum(sode(xs)),Flux.params(xs,sode))
#@test ! iszero(grads[xs])
#@test ! iszero(grads[sode.p])
#@test ! iszero(grads[sode.p][end])

ddudt = Chain(Dense(6,50,tanh),Dense(50,2))
NeuralCDDE(ddudt,(0.0f0,2.0f0),(p,t)->zero(x),(1f-1,2f-1),MethodOfSteps(Tsit5()),saveat=0.1)(x)
dode = NeuralCDDE(ddudt,(0.0f0,2.0f0),(p,t)->zero(x),(1f-1,2f-1),MethodOfSteps(Tsit5()),saveat=0.0:0.1:2.0)

grads = Zygote.gradient(()->sum(dode(x)),Flux.params(x,dode))
@test ! iszero(grads[x])
@test ! iszero(grads[dode.p])

@test_broken grads = Zygote.gradient(()->sum(dode(xs)),Flux.params(xs,dode)) isa Tuple
#@test ! iszero(grads[xs])
#@test ! iszero(grads[dode.p])

Slowddudt = SlowChain(SlowDense(6,50,tanh),SlowDense(50,2))
NeuralCDDE(Slowddudt,(0.0f0,2.0f0),(p,t)->zero(x),(1f-1,2f-1),MethodOfSteps(Tsit5()),saveat=0.1)(x)
dode = NeuralCDDE(Slowddudt,(0.0f0,2.0f0),(p,t)->zero(x),(1f-1,2f-1),MethodOfSteps(Tsit5()),saveat=0.0:0.1:2.0)

grads = Zygote.gradient(()->sum(dode(x)),Flux.params(x,dode))
@test ! iszero(grads[x])
@test ! iszero(grads[dode.p])

@test_broken grads = Zygote.gradient(()->sum(dode(xs)),Flux.params(xs,dode)) isa Tuple
#@test ! iszero(grads[xs])
#@test ! iszero(grads[dode.p])


text:
N=10000000
cnt=0
for i = 1:N
    x=rand()
    y=rand()
    if x^2+y^2<1.0
        cnt+=1
    end
end

@printf "Estimate of PI for %d trials is %8.5f\n" N 4.0*(cnt / N)

ext: jl

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
### QTP Example ###

using ModelPredictiveControl
using DifferentialEquations
using Plots
using LinearAlgebra

### QTP with DifferentialEquation.jl
a1      = 1.31e-4; #m²
a2      = 1.51e-4; #m²
a3      = 9.27e-5; #m²
a4      = 8.82e-5; #m²
S       = 0.06; #m
g       = 9.81; #m²/s
gamma_a = 0.3;
gamma_b = 0.4;

qa = 0; #m3/h
qb = 0; #m3/h

function QTP(dh, h, p, t)
    dh[1] = -a1 / S * sqrt(2 * g * h[1]) + a3 / S * sqrt(2 * g * h[3]) + gamma_a / S * qa / 3600;
    dh[2] = -a2 / S * sqrt(2 * g * h[2]) + a4 / S * sqrt(2 * g * h[4]) + gamma_b / S * qb / 3600;
    dh[3] = -a3 / S * sqrt(2 * g * h[3]) + (1 - gamma_b) * qb / (3600 * S);
    dh[4] = -a4 / S * sqrt(2 * g * h[4]) + (1 - gamma_a) * qa / (3600 * S);
 end

tspan = (0.0,60.0);
tsteps = range(tspan[1], tspan[2], step = 0.1);
u0 = [0.68; 0.68; 0.68; 0.68];
prob = ODEProblem(QTP,u0,tspan);
sol = solve(prob, BS3(), saveat = tsteps);

plot(sol)

### Design the MPC controller to regulate the water level

xl = [0.65;0.66;0.65;0.66]
ul = [1.63 ; 2.00]

xcons = [0.2 1.36;
         0.2 1.36;
         0.2 1.30;
         0.2 1.30] #m

ucons = [0  3.26;
         0  4] #m3/h

Ts = 5.0;#s
method = "l"; #linear MPC
Q =  Matrix(I, 4, 4)
R = 0.01 *  Matrix(I, 2, 2)
N = 15; #horizon length

 f((x1, x2, x3, x4, u1, u2)) =
 [ -a1 / S * sqrt(2 * g * x1) + a3 / S * sqrt(2 * g * x3) + gamma_a / S * u1 / 3600,
   -a2 / S * sqrt(2 * g * x2) + a4 / S * sqrt(2 * g * x4) + gamma_b / S * u2 / 3600,
   -a3 / S * sqrt(2 * g * x3) + (1 - gamma_b) * u2 / (3600 * S),
   -a4 / S * sqrt(2 * g * x4) + (1 - gamma_a) * u1 / (3600 * S),
 ]

RegulatorMPC = mpc_problem(f, xl, ul, xcons, ucons, Ts, method, Q, R, N)

state_reference = [0.3;0.3;0.301;0.306]
input_reference = [1.112;1.351]

    RegulatorMPC.x_t = xl
    RegulatorMPC.x_ref = state_reference .* ones(4, N+1)
    RegulatorMPC.u_ref = input_reference .* ones(2,N)
    println("toto")
    u, x, utilde, xtilde, objc = solve_mpc(RegulatorMPC)

obj = []
function ctrl_cb(int)
    x_meas = int.u
    RegulatorMPC.x_t = convert(Vector{Float64}, x_meas)
    RegulatorMPC.x_ref = state_reference .* ones(4, N+1)
    RegulatorMPC.u_ref = input_reference .* ones(2,N)
    u, x, utilde, xtilde, objc = solve_mpc(RegulatorMPC)
    global qa = u[1,1]
    global qb = u[2,1]
    println("qa : $qa")
    println("qb : $qb")
    push!(obj, objc)
end

pcb  = PeriodicCallback(ctrl_cb, 5.0)

tspan = (0.0,1500.0);
tsteps2 = range(tspan[1], tspan[2], step = 0.01);

prob2 = ODEProblem(QTP,u0, tspan, callback = pcb);

### Run the QTP problem with the regulator
sol2 = DifferentialEquations.solve(prob2, Tsit5());

plot(sol2)
plot(obj)


text:
"""
    LocalApproximation{N, VN<:AbstractVector{N}}

Type that represents a local approximation in 2D.

### Fields

- `p1`        -- first inner point
- `d1`        -- first direction
- `p2`        -- second inner point
- `d2`        -- second direction
- `q`         -- intersection of the lines l1 ⟂ d1 at p1 and l2 ⟂ d2 at p2
- `refinable` -- states if this approximation is refinable
- `err`       -- error upper bound

### Notes

The criteria for being refinable are determined in the method `new_approx`.
"""
struct LocalApproximation{N, VN<:AbstractVector{N}}
    p1::VN
    d1::VN
    p2::VN
    d2::VN
    q::VN
    refinable::Bool
    err::N
end

"""
    constraint(approx::LocalApproximation)

Convert a local approximation to a linear constraint.

### Input

- `approx` -- local approximation

### Output

A linear constraint.
"""
function constraint(approx::LocalApproximation)
    return LinearConstraint(approx.d1, dot(approx.d1, approx.p1))
end

"""
    PolygonalOverapproximation{N, SN<:LazySet{N}, VN<:AbstractVector{N}}

Type that represents the polygonal approximation of a convex set.

### Fields

- `S`            -- convex set
- `approx_stack` -- stack of local approximations that still need to be examined
- `constraints`  -- vector of linear constraints that are already finalized
                    (i.e., they satisfy the given error bound)
"""
struct PolygonalOverapproximation{N, SN<:LazySet{N}, VN<:AbstractVector{N}}
    S::SN
    approx_stack::Vector{LocalApproximation{N, VN}}
    constraints::Vector{LinearConstraint{N, VN}}
end

function PolygonalOverapproximation(S::SN) where {N, SN<:LazySet{N}}
    empty_local_approx = Vector{LocalApproximation{N, Vector{N}}}()
    empty_constraints = Vector{LinearConstraint{N,Vector{N}}}()
    return PolygonalOverapproximation(S, empty_local_approx, empty_constraints)
end

"""
    new_approx(S::LazySet, p1::VN, d1::VN,
               p2::VN, d2::VN) where {N<:AbstractFloat, VN<:AbstractVector{N}}

Create a `LocalApproximation` instance for the given excerpt of a polygonal
approximation.

### Input

- `S`  -- convex set
- `p1` -- first inner point
- `d1` -- first direction
- `p2` -- second inner point
- `d2` -- second direction

### Output

A local approximation of `S` in the given directions.
"""
function new_approx(S::LazySet, p1::VN, d1::VN,
                    p2::VN, d2::VN) where {N<:AbstractFloat, VN<:AbstractVector{N}}
    if norm(p1-p2, 2) <= _rtol(N)
        # this approximation cannot be refined and we set q = p1 by convention
        q = p1
        refinable = false
        err = zero(N)
    else
        ndir = normalize([p2[2]-p1[2], p1[1]-p2[1]])
        q = element(intersection(Line2D(d1, dot(d1, p1)), Line2D(d2, dot(d2, p2))))
        err = min(norm(q - σ(ndir, S)), dot(ndir, q - p1))
        refinable = (err > _rtol(N)) && (norm(p1-q, 2) > _rtol(N)) &&
                    (norm(q-p2, 2) > _rtol(N))
    end
    return LocalApproximation(p1, d1, p2, d2, q, refinable, err)
end

"""
    addapproximation!(Ω::PolygonalOverapproximation, p1::VN, d1::VN,
                      p2::VN, d2::VN) where {N, VN<:AbstractVector{N}}

### Input

- `Ω`  -- polygonal overapproximation of a convex set
- `p1` -- first inner point
- `d1` -- first direction
- `p2` -- second inner point
- `d2` -- second direction

### Output

The list of local approximations in `Ω` of the set `Ω.S` is updated in-place and
the new approximation is returned by this function.
"""
function addapproximation!(Ω::PolygonalOverapproximation, p1::VN, d1::VN,
                           p2::VN, d2::VN) where {N, VN<:AbstractVector{N}}
    approx = new_approx(Ω.S, p1, d1, p2, d2)
    push!(Ω.approx_stack, approx)
    return approx
end

"""
    refine(approx::LocalApproximation, S::LazySet)

Refine a given local approximation of the polygonal approximation of a convex
set by splitting along the normal direction of the approximation.

### Input

- `approx` -- local approximation to be refined
- `S`      -- 2D convex set

### Output

The tuple consisting of the refined right and left local approximations.
"""
function refine(approx::LocalApproximation, S::LazySet)
    @assert approx.refinable
    ndir = normalize([approx.p2[2]-approx.p1[2], approx.p1[1]-approx.p2[1]])
    s = σ(ndir, S)
    ap1 = new_approx(S, approx.p1, approx.d1, s, ndir)
    ap2 = new_approx(S, s, ndir, approx.p2, approx.d2)
    return (ap1, ap2)
end

"""
    tohrep(Ω::PolygonalOverapproximation)

Convert a polygonal overapproximation into a concrete polygon.

### Input

- `Ω` -- polygonal overapproximation of a convex set

### Output

A polygon in constraint representation.

### Algorithm

Internally we keep the constraints sorted.
Hence we do not need to use `addconstraint!` when creating the `HPolygon`.
"""
function tohrep(Ω::PolygonalOverapproximation)
    # already finalized
    if isempty(Ω.approx_stack)
        return HPolygon(Ω.constraints, sort_constraints=false)
    end
    # some constraints not finalized yet
    constraints = copy(Ω.constraints)
    for approx in Ω.approx_stack
        push!(constraints, constraint(approx))
    end
    return HPolygon(constraints, sort_constraints=false)
end

"""
    _approximate(S::LazySet{N}, ε::Real) where {N<:AbstractFloat}

Return an ε-close approximation of the given 2D convex set (in terms of
Hausdorff distance) as an inner and an outer approximation composed by sorted
local `Approximation2D`.

### Input

- `S` -- 2D convex set
- `ε` -- error bound

### Output

An ε-close approximation of the given 2D convex set.
"""
function _approximate(S::LazySet{N}, ε::Real) where {N<:AbstractFloat}
    # initialize box directions
    pe = σ(DIR_EAST(N), S)
    pn = σ(DIR_NORTH(N), S)
    pw = σ(DIR_WEST(N), S)
    ps = σ(DIR_SOUTH(N), S)
    east = dir_east(N, pe)
    north = dir_north(N, pe)
    west = dir_west(N, pe)
    south = dir_south(N, pe)

    Ω = PolygonalOverapproximation(S)

    # add constraints in reverse (i.e., clockwise) order to the stack
    addapproximation!(Ω, ps, south, pe, east)
    addapproximation!(Ω, pw, west, ps, south)
    addapproximation!(Ω, pn, north, pw, west)
    addapproximation!(Ω, pe, east, pn, north)

    approx_stack = Ω.approx_stack
    while !isempty(approx_stack)
        approx = pop!(approx_stack)

        if !approx.refinable || approx.err <= ε
            # if the approximation is not refinable => continue
            push!(Ω.constraints, constraint(approx))
            continue
        end

        # refine
        (la1, la2) = refine(approx, Ω.S)

        if isempty(approx_stack)
            redundant = false
        else
            # check if the next local approximation became redundant
            next = approx_stack[end]
            redundant = (norm(la2.p1-next.p1) <= _rtol(N)) &&
                (norm(la2.q-next.q) <= _rtol(N))
        end
        if redundant
            # replace redundant old constraint
            approx_stack[end] = la2
        else
            # add new constraint
            push!(approx_stack, la2)
        end
        push!(approx_stack, la1)
    end
    return Ω
end


text:
module ReferenceFrameRotations

import Base: +, -, *, /, \, ≈, ==
import Base: conj, convert, copy, display, eltype, firstindex, getindex, imag
import Base: inv, iterate, lastindex, ndims, one, real, setindex!, show
import Base: zero, zeros
import Base: Broadcast.broadcastable
import LinearAlgebra: norm

using Crayons
using LinearAlgebra
using Printf
using StaticArrays

# Re-export `I` from LinearAlgebra.
export I

################################################################################
#                                    Types
################################################################################

include("types.jl")

################################################################################
#                                  Constants
################################################################################

# Pre-defined crayons.
const _reset_crayon = Crayon(reset = true)
const _crayon_bold  = crayon"bold"
const _crayon_g     = crayon"bold green"
const _crayon_u     = crayon"bold blue"
const _crayon_y     = crayon"bold yellow"

# Escape sequences related to the crayons.
const _b = _crayon_bold
const _d = _reset_crayon
const _g = _crayon_g
const _y = _crayon_y
const _u = _crayon_u

################################################################################
#                                   Includes
################################################################################

include("angle.jl")
include("angleaxis.jl")
include("compose_rotations.jl")
include("dcm.jl")
include("inv_rotations.jl")
include("quaternion.jl")

include("deprecations.jl")

include("./conversions/angle_to_angle.jl")
include("./conversions/angle_to_angleaxis.jl")
include("./conversions/angle_to_dcm.jl")
include("./conversions/angle_to_quat.jl")
include("./conversions/angle_to_rot.jl")
include("./conversions/angleaxis_to_angle.jl")
include("./conversions/angleaxis_to_dcm.jl")
include("./conversions/angleaxis_to_quat.jl")
include("./conversions/dcm_to_angle.jl")
include("./conversions/dcm_to_angleaxis.jl")
include("./conversions/dcm_to_quat.jl")
include("./conversions/quat_to_angle.jl")
include("./conversions/quat_to_angleaxis.jl")
include("./conversions/quat_to_dcm.jl")
include("./conversions/smallangle_to_dcm.jl")
include("./conversions/smallangle_to_quat.jl")
include("./conversions/smallangle_to_rot.jl")

end # module


text:
# Autogenerated wrapper script for Gloo_jll for aarch64-apple-darwin
export libgloo

JLLWrappers.@generate_wrapper_header("Gloo")
JLLWrappers.@declare_library_product(libgloo, "@rpath/libgloo.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libgloo,
        "lib/libgloo.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function is_reproducing(genotype, cutoff_reproductiblity, certainty)
    max_samples = log(1 - certainty) / log(1 - cutoff_reproductiblity)

    for _ in 1:max_samples
        snapshot_to_test = sample_snapshot_id() |> get_snapshot

        if is_reproducing(genotype, snapshot_to_test)
            return true
        end
    end

    return false
end

function is_reproducing(genotype, snapshot)
    organism_to_replace = get_organisms(snapshot) |> rand
    id_to_replace = get_id(snapshot, organism_to_replace)

    fitness = get_fitness(snapshot, id_to_replace, genotype)

    return 0 < fitness
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
@testset "model" begin

    @test_nowarn begin
        m = @agent(
            1,

            l::Local,
            g::Global,
            ga::GlobalArray,
            m::Medium,

            Equation = d_v = 34*dt ,
            UpdateLocal = l += 1,
            UpdateGlobal = begin
                g += Normal(1,2) + Normal(1,2) + Uniform(1,2)
                ga[3,2] += 3
            end,
            UpdateInteraction = l.i += 1,
            UpdateLocalInteraction = l.i += 1,

            Boundary = BoundaryFlat(1)
        )
        compile(m, platform="cpu",debug=false)
    end

end


text:
function gaussiansample(μ, σ2)
	ϵ = randn!(similar(μ,size(μ)))
	μ .+ sqrt.(σ2) .* ϵ
end

"""
		Implementation of VAE_ard with Gaussian prior and posterior.

		VAE_ard(q,g,β,s)

		q --- encoder
		g --- decoder

		β --- (variance of the posterior distribution) or the strength on KL-divergence ``D_{KL}(p(z)\\|q(z|x))``
		s = :unit  p(x|z) has a distribution ``N(\\mu(z),I)``
		s = :sigmadiag p(x|z) has a distribution ``N(\\mu(z),\\sigma^2(z))``

		Let's assume that latent dimension should be `n`, then encoder should have output dimension `2n`.
		Similarly, if `s = Val{:unit}`, the g should gave output dimension `d`, as it simply codes the mean
		```julia
			m = VAE_ard(layerbuilder(idim,hidden,2*zdim,3,"relu","linear","Dense"),
	      layerbuilder(zdim,hidden,idim,3,"relu","linear","Dense"),1.0,:unit)
		```
		but for inferring variance of the normal distribution `s = Val{sigmadiag}`, output dimension of `g` should be `2d`.

		```julia
			m = VAE_ard(layerbuilder(idim,hidden,2*zdim,3,"relu","linear","Dense"),
	      layerbuilder(zdim,hidden,2*idim,3,"relu","linear","Dense"),1.0,:sigma)
		```
"""
struct VAE_ard{T<:AbstractFloat,V<:Val}
	q  # encoder (inference modul)
	g  # decoder (generator)
    β::AbstractArray{T, 1} 	#penalization = log σ^2
    γ
	variant::V
end

VAE_ard(q,g,β,γ,s::Symbol = :unit) = VAE_ard(q,g,β,γ,Val(s))

function VAE_ard(inputDim, hiddenDim, latentDim, numLayers, nonlinearity, layerType, β, V = :unit, T = Float32)
	encoder = Adapt.adapt(T, FluxExtensions.layerbuilder(inputDim, hiddenDim, latentDim * 2, numLayers + 1, nonlinearity, "linear", layerType))
    decoder = Adapt.adapt(T, FluxExtensions.layerbuilder(latentDim, hiddenDim, inputDim, numLayers + 1, nonlinearity, "linear", layerType))
	return VAE_ard(encoder, decoder, param([T(β)]), param(ones(T, latentDim)), V)
end

Flux.@treelike(VAE_ard)

"""

	infer(m::VAE_ard{T,V},x)

	infer latent variables and sample output x

"""
# function infer(m::VAE_ard{T,V},x) where {T,V<:Val{:sigmadiag}}
# 	μz, σ2z = hsplitsoftp(m.q(x))
# 	z = gaussiansample(μz,σ2z)
# 	μx, σ2x = hsplitsoftp(m.g(z))
# 	μz, σ2z, μx, σ2x
# end

function infer(m::VAE_ard{T,V},x) where {T,V<:Val{:scalarsigma}}
	μz, σ2z = hsplitsoftp(m.q(x))
	z = gaussiansample(μz,σ2z)
	μx, σ2x = hsplit1softp(m.g(z))
	μz, σ2z, μx, σ2x
end

function infer(m::VAE_ard{T,V},x) where {T,V<:Val{:unit}}
	μz, σ2z = hsplitsoftp(m.q(x))
	z = gaussiansample(μz,σ2z)
	μx  = m.g(z)
	μz, σ2z, μx
end

"""

	loss(m::VAE_ard{T,V},x)

	loss of the Variational autoencoder ``\\mathbb{E}_{q(z|x)log(p(x|z) - KL(p(z)|q(z|x))``
	with the KL-divergence calculated analytically, since the p(z) and q(z|x) are both Gaussian distributions
"""
# function loss(m::VAE_ard{T,V},x) where {T,V<:Val{:sigmadiag}}
# 	μz, σ2z, μx, σ2x = infer(m,x)
# 	-mean(log_normal(x,μx,σ2x)) + m.β * mean(kldiv(μz,σ2z))
# end

function elbo_loss(m::VAE_ard{T,V},x) where {T,V<:Val{:scalarsigma}}
	μz, σ2z, μx, σ2x = infer(m,x)
	-mean(log_normal(x,μx,collect(σ2x'))) + mean(kldiv(μz,σ2z))
end

function elbo_loss(m::VAE_ard{T,V},x) where {T,V<:Val{:unit}}
	μz, σ2z, μx = infer(m,x)
	-mean(log_normal(x,μx)) + exp(m.β[1]) * mean(kldiv(μz,σ2z)) + size(x, 1) * size(x, 1) / 4f0 * m.β[1] * exp(m.β[1])
end

function elbo_loss_ard(m::VAE_ard{T,V}, x) where {T,V<:Val{:unit}}
    μz, σ2z = zparams(m, x)
    z = gaussiansample(μz, σ2z)
    μx = m.g(z .* m.γ)
	-mean(log_normal(x, μx)) + exp(m.β[1]) * mean(kldiv(μz,σ2z)) + size(x, 1) * size(x, 1) / 4f0 * m.β[1] * exp(m.β[1])
end

function x_from_z(m::VAE_ard{T,V}, z) where {T,V<:Val{:unit}}
	μx = m.g(z .* m.γ)
end

function decomposed_elbo_loss(m::VAE_ard{T,V},x) where {T,V<:Val{:unit}}
	μz, σ2z, μx = infer(m,x)
	-mean(log_normal(x,μx)), mean(kldiv(μz,σ2z))
end

function wass_dist(m::VAE_ard, x, d)
	μz, σ2z = zparams(m, x)
	z = gaussiansample(μz,σ2z)
	prior = randn(size(z))
	Ω = d(z, prior)
	return Ω, z
end

function wloss(m::VAE_ard{T,V}, x, d) where {T,V<:Val{:scalarsigma}}
	Ω, z = wass_dist(m, x, d)
	μx, σ2x = hsplit1softp(m.g(z))
	-mean(log_normal(x, μx, collect(σ2x'))) + Ω
end

function wloss(m::VAE_ard{T,V}, x, d) where {T,V<:Val{:unit}}
	Ω, z = wass_dist(m, x, d)
	μx = m.g(z)
	-mean(log_normal(x, μx)) + m.β * Ω
end

function rloss(m::VAE_ard{T,V}, x) where {T,V<:Val{:scalarsigma}}
	μz, σ2z = zparams(m, x)
	z = gaussiansample(μz,σ2z)
	μx, σ2x = hsplit1softp(m.g(z))
	-mean(log_normal(x, μx, collect(σ2x')))
end

function rloss(m::VAE_ard{T,V}, x) where {T,V<:Val{:unit}}
	μz, σ2z = zparams(m, x)
	z = gaussiansample(μz,σ2z)
	μx = m.g(z)
	-mean(log_normal(x, μx))
end

function samplez(m::VAE_ard{T,V}, x) where {T,V}
	z = gaussiansample(zparams(m, x)...)
end

function zparams(m::VAE_ard{T,V}, x) where {T,V}
	μz, σ2z = hsplitsoftp(m.q(x))
	return μz, σ2z
end

function printing_loss(m::VAE_ard{T,V},x) where {T,V<:Val{:scalarsigma}}
	μz, σ2z, μx, σ2x = infer(m,x)
	println("loglkl: $(-mean(log_normal(x,μx,collect(σ2x')))) | KL: $(mean(kldiv(μz,σ2z)))")
	-mean(log_normal(x,μx,collect(σ2x'))) + mean(kldiv(μz,σ2z))
end

function printing_loss(m::VAE_ard{T,V},x) where {T,V<:Val{:unit}}
	μz, σ2z, μx = infer(m,x)
	println("MSE:$(-mean(log_normal(x,μx))) Flux.mse: $(Flux.mse(x, μx)) KL: $(mean(kldiv(μz,σ2z)))")
	-mean(log_normal(x,μx)) + exp(m.β[1]) * mean(kldiv(μz,σ2z)) + size(x, 1) * size(x, 1) / 4f0 * m.β[1] * exp(m.β[1])
end

function printing_wloss(m::VAE_ard{T,V}, x, d) where {T,V<:Val{:scalarsigma}}
	Ω, z = wass_dist(m, x, d)
	μx, σ2x = hsplit1softp(m.g(z))
	lklh = -mean(log_normal(x, μx, collect(σ2x')))
	println("loglklh: $lklh | wass dist: $Ω")
	lklh + Ω
end

function printing_wloss(m::VAE_ard{T,V}, x, d) where {T,V<:Val{:unit}}
	Ω, z = wass_dist(m, x, d)
	lklh = -mean(log_normal(x, m.g(z)))
	println("loglklh: $lklh | wass dist: $Ω")
	lklh + m.β * Ω
end

log_pz(m::VAE_ard, x) = log_normal(hsplitsoftp(m.q(x))[1])

function log_pxexpectedz(m::VAE_ard{T,V},x, σ::AbstractFloat = 1.0) where {T, V<:Val{:unit}}
	μz, σ2z = hsplitsoftp(m.q(x))
	log_normal(x,m.g(μz .* m.γ),σ)
end

function log_pxexpectedz(m::VAE_ard{T,V},x, z::AbstractArray) where {T, V<:Val{:unit}}
	log_normal(x, m.g(z .* m.γ))
end

function log_pxexpectedz(m::VAE_ard{T,V}, x) where {T, V<:Val{:scalarsigma}}
	μz, σ2z = hsplitsoftp(m.q(x))
	μx, σ2x = hsplit1softp(m.g(μz .* m.γ))
	log_normal(x,μx,collect(σ2x'))
end

function log_pxexpectedz(m::VAE_ard{T,V}, x, z) where {T, V<:Val{:scalarsigma}}
	μx, σ2x = hsplit1softp(m.g(z .* m.γ))
	log_normal(x,μx,collect(σ2x'))
end

jacobian_encoder(m::VAE_ard, x) = Flux.Tracker.jacobian(a -> vec(hsplitsoftp(m.q(a))[1]), x)
jacobian_decoder(m::VAE_ard{T,V}, z) where {T, V<:Val{:scalarsigma}} = Flux.Tracker.jacobian(a -> vec(hsplit1softp(m.g(a))[1]), z)
jacobian_decoder(m::VAE_ard{T,V}, z) where {T, V<:Val{:unit}} = Flux.Tracker.jacobian(a -> m.g(a), z)

function log_det_jacobian_encoder(m::VAE_ard, x)
	if size(x, 2) > 1
		xs = [x[:, i] for i in 1:size(x, 2)]
		return map(x -> log_det_jacobian_encoder_singleinstance(m, x), xs)
	else
		return log_det_jacobian_encoder_singleinstance(m, x)
	end
end

function log_det_jacobian_encoder_singleinstance(m::VAE_ard, x)
	@assert size(x, 2) == 1
	s = svd(jacobian_encoder(m, x).data)
	d = reduce(+, log.(abs.(s.S))) * 2
end

function log_det_jacobian_decoder(m::VAE_ard, z)
	z = (z .* m.γ).data
	if size(z, 2) > 1
		zs = [z[:, i] for i in 1:size(z, 2)]
		return map(z -> log_det_jacobian_decoder_singleinstance(m, z), zs)
	else
		return log_det_jacobian_decoder_singleinstance(m, z)
	end
end

function log_det_jacobian_decoder_singleinstance(m::VAE_ard, z)
	@assert size(z, 2) == 1
	s = svd(jacobian_decoder(m, z).data)
	d = reduce(+, log.(abs.(s.S))) * 2
end

jacodeco(m::VAE_ard, x::Vector, z::Vector) = sum(log_normal((z .* m.γ).data)) + sum(log_pxexpectedz(m, x, z)) - det(m.g, (z .* m.γ).data)

function correctedjacodeco(model::VAE_ard, x::Vector, z::Vector)
	z = (z .* model.γ).data
	Σ = Flux.data(Flux.Tracker.jacobian(model.g, z))
	S = svd(Σ)
	J = inv(S)
	logd = 2*sum(log.(S.S .+ 1f-6))
	if (det(I + J * transpose(J)) < 0)
		# println("determinant = $(det(I + J * transpose(J))) was negative, returning jacodeco = -Inf")
		# println(J)
		# println(J * transpose(J))
		return -Inf32
	end
	logpz = log_normal(z, zeros(Float32, size(z)...), I + J * J')[1] # TODO I changed the transpose(J) to J' - might be an issue
	logpx = FluxExtensions.log_normal(x, x_from_z(model, z).data, 1, length(x) - length(z))[1]
	logpz + logpx - logd
end

max_stars_repo_name: honzabim/FewShotAnomalyDetection

max_issues_repo_name: honzabim/FewShotAnomalyDetection

max_forks_repo_name: honzabim/FewShotAnomalyDetection


text:
module MiniPascal
include("CodeGen.jl")

import Base: *,+,-,÷,<,!,&,==,showerror

DEBUG = false

end # module

max_stars_repo_name: juhamyllari/minipascal-julia

max_issues_repo_name: juhamyllari/minipascal-julia

max_forks_repo_name: juhamyllari/minipascal-julia


text:
module LagrangeTests
using ForwardDiff
using Morbit
using Test
#%%
mop = MixedMOP(2)
add_objective!(mop, x -> sum(x.^4), LagrangeConfig(;optimized_sampling=false))
x0 = rand(2)
smop, iter_data, data_base, sc, ac = Morbit.initialize_data(mop, x0, Float32[] )

mod = sc.surrogates[1].model
objf = Morbit.list_of_objectives(smop)[1]
meta = sc.surrogates[1].meta

# `Morbit.eval_objf` evaluates at a site from the original domain
_eval_objf_original = x -> Morbit.eval_objf(objf, x)
_eval_objf_scaled = x -> Morbit.eval_objf( objf, Morbit.unscale(x, smop) )

# `eval_models`, `get_gradient` etc., on the other hand,
# take input from [0,1]^n:
_eval_model = x -> Morbit.eval_models( mod, x )

x0_s = Morbit.scale( x0, smop )

#%%
# Test some configurations
# Only degrees 1 and 2 are allowed
@test_throws Any Morbit.LagrangeConfig(;degree=0)
@test_throws Any Morbit.LagrangeConfig(;degree=-1)
@test_throws Any Morbit.LagrangeConfig(;degree=3)

# Λ must be >1
@test_throws Any Morbit.LagrangeConfig(;LAMBDA=-1)
@test_throws Any Morbit.LagrangeConfig(;LAMBDA=1.0)

# Algorithm must be derivative-free
@test_throws Any Morbit.LagrangeConfig(; algo1_solver = :LD_MMA)
@test_throws Any Morbit.LagrangeConfig(; algo2_solver = :LD_MMA)
@test_throws Any Morbit.LagrangeConfig(; algo2_solver = :L)

#%%
# preparation

const scalar_test_funcs = [
	x -> sum(x),
	x -> sum( x.^2 ),
	x -> sum( x.^(2 : length(x) + 1) ),
	x -> exp( sum( x.^2 ) / 100 ),
	x -> 1 / (10 + abs(sum(x))) * sin( sum(x) )
]

const lagrange_configs = [
	LagrangeConfig(; degree = 1 ),
	LagrangeConfig(; degree = 1, LAMBDA = 10 ),
	LagrangeConfig(; degree = 2 ),
	LagrangeConfig(; degree = 2, LAMBDA = 10 ),
	LagrangeConfig(; optimized_sampling = false, degree = 1 ),
	LagrangeConfig(; optimized_sampling = false, degree = 2 ),
]

#%%
for cfg in lagrange_configs
	for n_vars = 1 : 3
		for n_out = 1 : 2
			for f1 in scalar_test_funcs
				for unconstrained in [true, false]

					lb = fill(-5.0, n_vars )
					ub = fill(6.0, n_vars )

					if unconstrained
						mop = Morbit.MixedMOP(n_vars)
					else
						mop = Morbit.MixedMOP(lb,ub)
					end

					if n_out == 1
						Morbit.add_objective!( mop, f1, cfg )
					else
						f_remaining = rand( scalar_test_funcs, n_out - 1 )
						test_func = x -> [ f1(x); [φ(x) for φ in f_remaining] ]
						Morbit.add_vector_objective!( mop, test_func, cfg; n_out )
					end

					x0 = Morbit._rand_box_point(lb, ub)

					smop, iter_data, data_base, sc, ac = Morbit.initialize_data(mop, x0, Float32[] )

					mod = sc.surrogates[1].model
					objf = Morbit.list_of_objectives(smop)[1]
					meta = sc.surrogates[1].meta

					# `Morbit.eval_objf` evaluates at a site from the original domain
					_eval_objf_original = x -> Morbit.eval_objf(objf, x)
					_eval_objf_scaled = x -> Morbit.eval_objf( objf, Morbit.unscale(x, smop) )

					# `eval_models`, `get_gradient` etc., on the other hand,
					# take input from [0,1]^n:
					_eval_model = x -> Morbit.eval_models( mod, x )

					x0_s = Morbit.scale( x0, smop )

					#%%
					# check if scaling went alright:
					@test x0_s ≈ Morbit.get_x( iter_data )
					@test _eval_objf_original( x0 ) ≈ _eval_objf_scaled( x0_s )

					# did we succeed in making an orthogonal basis?
					X = Morbit.get_site.( data_base, meta.interpolation_indices )
					Y = Morbit.get_value.( data_base, meta.interpolation_indices )

					num_polys = length(X)

					for (i,poly) in enumerate(meta.lagrange_basis)
						res_vec = zeros(num_polys)
						res_vec[i] = 1
						@test poly.(X) ≈ res_vec
					end

					for (i,poly) in enumerate(mod.basis)
						res_vec = zeros(num_polys)
						res_vec[i] = 1
						@test poly.(X) ≈ res_vec
					end

					# if the database is transformed, then the sites are from the
					# original domain => we must scale them for the model.
					# else they are in [0,1] and can be passed to the model.
					_eval_func = Morbit.is_transformed( data_base ) ? x -> _eval_model( Morbit.scale(x, smop) ) : _eval_model

					@test all([ isapprox(_eval_func( X[i] ), Y[i][ meta.out_indices ]; atol = 1e-2, rtol = 1e-3) for i = eachindex(X) ])

					#%%
					# does evaluation work alright?
					if 1 in meta.interpolation_indices
						@test _eval_model( x0_s ) ≈ _eval_objf_original( x0 )
					end
					@test Morbit.get_gradient( mod, x0_s, 1) ≈ ForwardDiff.gradient( x -> _eval_model(x)[1], x0_s )
					@test Morbit.get_jacobian( mod, x0_s ) ≈ ForwardDiff.jacobian( _eval_model, x0_s )
				end#for
			end#for
		end#for
	end#for
end#for

end#module

using .LagrangeTests


text:
using ForwardDiff: ForwardDiff

"""
    DarlingHollowayElbow()

Adiabatic elbow potential from Darling and Holloway: [Faraday Discuss., 1993, 96, 43-54](https://doi.org/10.1039/FD9939600043)
"""
Parameters.@with_kw struct DarlingHollowayElbow <: AdiabaticModel

    # Elbow parameters
    d::Float64 = austrip(4.76u"eV")
    α::Float64 = 1.028
    p::Float64 = 2.0
    zoff::Float64 = 0.3

    # Dissociation barrier parameters
    V₀::Float64 = austrip(0.7u"eV")
    βx::Float64 = 0.05
    xb::Float64 = 1.0
    βz::Float64 = 0.0
    zb::Float64 = 0.0

    # Physisorption parameters
    V₁::Float64 = austrip(0.85u"eV")
    βr::Float64 = 1.8
    zr::Float64 = 0.25
    Cvw::Float64 = austrip(5.306u"eV")
    kc::Float64 = 0.5
    zvw::Float64 = -1.9

    m₁::Float64 = 1.3
    m₂::Float64 = 0.0
    m₃::Float64 = -0.1

    w₁::Float64 = 1.5
    w₂::Float64 = 1.5
    w₃::Float64 = 2.0
end

NQCModels.ndofs(::DarlingHollowayElbow) = 1

function NQCModels.potential(model::DarlingHollowayElbow, R::AbstractVector)

    Parameters.@unpack d, α, p, zoff, V₀, βx, xb, βz, zb, V₁, βr, zr, Cvw, zvw, kc, m₁, m₂, m₃, w₁, w₂, w₃ = model

    y(x,z) = (1 + 1/x^4)*(1 + 1/z^4)
    f(x,z) = (y(x,z) - 1)^(-1/4)
    s(a,b) = (1 + tanh(a*b)) / 2
    function xₑ(x,z)
        xp = x + p
        zp = z + p
        fp = f(xp,zp)
        return fp + f(xp-fp, zp-fp) - p
    end

    zf(x,z) = z - zoff * s(x - m₁, w₁)

    Velbow(x,z) = d*(1 - exp(-α*xₑ(x,z)))^2
    Vbarrier(x,z) = V₀*exp(-βx*(x-xb)^2 - βz*(zf(x,z)-zb)^2) * s(x-zf(x,z)-m₂,w₂)

    Vᵣ(z) = V₁*exp(-βr*(z-zr))
    fc(x) = 1 - (2x*(1+x)+1)*exp(-2x)
    Vvw(z) = - Cvw/(z - zvw)^3 * fc(kc*(z-zvw))
    Vphys(x,z) = (Vᵣ(zf(x,z)) + Vvw(zf(x,z))) * s(zf(x,z)-x-m₂, w₂) * s(zf(x,z)-m₃, w₃)

    total_V(x,z) = Velbow(x,z) + Vbarrier(x,z) + Vphys(x,z)

    x = R[1]
    Z = R[2]

    return total_V(x, Z)
end

function NQCModels.derivative!(model::DarlingHollowayElbow, D, R::AbstractVector)
    f(x) = NQCModels.potential(model, x)
    copyto!(D, ForwardDiff.gradient(f, R))
    return D
end


text:
using Knet, Base.Test, BenchmarkTools

@show a = reshape(KnetArray([1.0:6.0...]), (2,3))
@show a[2]
@show a[2,3]
@show (a[2]=7;a)
@show (a[2,3]=8;a)
@show a[:]
@show (a[:]=9;a)
@show (a[:]=[11.0:16.0...];a)
@show (a[:]=[21:26...];a)
@show (b=KnetArray(rand(6));a[:]=b;a)
@show (b=KnetArray(rand(0:9,6));a[:]=b;a)
@show a[2:4]
@show (a.ptr, a[2:4].ptr)
@show (a[2:4]=10;a)
@show (a[1:2]=a[5:6];a)
@show (a[1:2]=ones(2);a)
@show (a[5:6]=2ones(Int,2);a)
@show a[:,:]
@show a[1,:]
@show a[:,1:2]
@show a[1,1:2]
@show a[1:2,1:2]
@show (a[:,:]=21;a)
@show (a[1,:]=22;a)
@show (a[:,1:2]=23;a)
@show (a[1,1:2]=24;a)
@show (a[1:2,1:2]=25;a)
@show (a[:,:]=31ones(6);a)
@show (a[1,:]=32ones(3);a)
@show (a[:,1:2]=33ones(4);a)
@show (a[1,1:2]=34ones(2);a)
@show (a[1:2,1:2]=35ones(4);a)
@show (a[:,:]=41ones(Int,6);a)
@show (a[1,:]=42ones(Int,3);a)
@show (a[:,1:2]=43ones(Int,4);a)
@show (a[1,1:2]=44ones(Int,2);a)
@show (a[1:2,1:2]=45ones(Int,4);a)

a = rand(3,3)
v = rand(3)
t = v'
ka = KnetArray(a)
kv = KnetArray(v)
kt = KnetArray(t)
@test hcat(a,a) == Array(hcat(ka,ka))
@test hcat(a,v) == Array(hcat(ka,kv))
@test hcat(v,a) == Array(hcat(kv,ka))
@test hcat(v,v) == Array(hcat(kv,kv))
@test vcat(a,a) == Array(vcat(ka,ka))
@test vcat(a,t) == Array(vcat(ka,kt))
@test vcat(t,a) == Array(vcat(kt,ka))
@test vcat(v,v) == Array(vcat(kv,kv))

a = rand(1000,1000)
v = rand(1000)
t = v'
ka = KnetArray(a)
kv = KnetArray(v)
kt = KnetArray(t)
@show @benchmark hcat(a,a)
@show @benchmark hcat(ka,ka)
@show @benchmark KnetArray(hcat(Array(ka),Array(ka)))
@show @benchmark hcat(a,v)
@show @benchmark hcat(ka,kv)
@show @benchmark KnetArray(hcat(Array(ka),Array(kv)))
@show @benchmark hcat(v,a)
@show @benchmark hcat(kv,ka)
@show @benchmark KnetArray(hcat(Array(kv),Array(ka)))
@show @benchmark hcat(v,v)
@show @benchmark hcat(kv,kv)
@show @benchmark KnetArray(hcat(Array(kv),Array(kv)))
@show @benchmark vcat(a,a)
@show @benchmark vcat(ka,ka)
@show @benchmark KnetArray(vcat(Array(ka),Array(ka)))
@show @benchmark vcat(a,t)
@show @benchmark vcat(ka,kt)
@show @benchmark KnetArray(vcat(Array(ka),Array(kt)))
@show @benchmark vcat(t,a)
@show @benchmark vcat(kt,ka)
@show @benchmark KnetArray(vcat(Array(kt),Array(ka)))
@show @benchmark vcat(v,v)
@show @benchmark vcat(kv,kv)
@show @benchmark KnetArray(vcat(Array(kv),Array(kv)))


text:
module NewtonsMethod

using LinearAlgebra, ForwardDiff

function newtonroot(f, f′; x₀, tolerance=1E-7, maxiter=1000)
    x_old = x₀
    normdiff = Inf
    iter = 1
    while normdiff > tolerance && iter <= maxiter
        x_new = x_old - (f(x_old) / f′(x_old))
        normdiff = norm(x_new - x_old)
        x_old = x_new
        iter = iter + 1
    end
    if iter > maxiter && normdiff > tolerance
        return NaN
    else
        return x_old
    end
end

function newtonroot(f; x₀, tolerance=1E-7, maxiter=1000)
    x_old = x₀
    f′ = x -> ForwardDiff.derivative(f, x)
    normdiff = Inf
    iter = 1
    while normdiff > tolerance && iter <= maxiter
        x_new = x_old - (f(x_old) / f′(x_old))
        normdiff = norm(x_new - x_old)
        x_old = x_new
        iter = iter + 1
    end
    if iter > maxiter && normdiff > tolerance
        return NaN
    else
        return x_old
    end
end

export newtonroot

end


text:
module PhysicalParameterizations

using OceanTurb
using OceanParameterizations.DataWrangling
using Oceananigans.Grids: Center, Face

export closure_kpp_full_evolution,
       closure_tke_full_evolution

include("k_profile_parameterization.jl")
include("turbulent_kinetic_energy_closure.jl")

end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.3.3
#   kernelspec:
#     display_name: Julia 1.3.0
#     language: julia
#     name: julia-1.3
# ---

# https://xavierbourretsicotte.github.io/loess.html

using LinearAlgebra, Statistics, Plots

# +
"""
    lowess_1d(x, y, tau = .005) -> yest

Locally weighted regression: fits a nonparametric regression curve
to a scatterplot.
The arrays x and y contain an equal number of elements; each pair
(x[i], y[i]) defines a data point in the scatterplot. The function
returns the estimated (smooth) values of y.
The kernel function is the bell shaped function with parameter tau.
Larger tau will result in a smoother curve.
"""
function lowess_1d(x, y; tau = 0.005)

    m = length(x)
    yest = zeros(Float64, m)

    #Initializing all weights from the bell shape kernel function
    w = [exp.(- (x .- x[i]).^2 ./(2*tau)) for i in 1:m]
    weights = zeros(Float64, m)

    #Looping through all x-points
    for i in 1:m
        weights .= w[i]
        b = [sum(weights .* y), sum(weights .* y .* x)]
        A = SymTridiagonal(
                [sum(weights), sum(weights .* x .* x)],
                [sum(weights .* x)])
        theta = A \ b
        yest[i] = theta[1] + theta[2] * x[i]
    end

    yest

end
# -

#Initializing noisy non linear data
x = range(0,stop=1,length=100) |> collect
noise = 0.2 * randn(100)
y = sin.(x * 3π/2 )
y_noise = y .+ noise
ypred = lowess_1d(x, y_noise; tau = 1e-2)
plot(x, ypred; label="lowess")
scatter!(x, y_noise)


text:
#formats the ccm table in a reasonable way
function prepccm(;
  ccmname::String = PARAM[:ccmname],
  ccmpath::String = PARAM[:ccmpath],
  ccmcolumns::Vector{Symbol} = PARAM[:ccmcolumns],
  lastdatestring::String = "20991231")

  ccm::DataFrame = CSV.File("$ccmpath\\$ccmname.csv") |> DataFrame

  #names ot lower case
  rename!(ccm, (s::String->lowercase(s)).(names(ccm)))
  rename!(ccm, :lpermno=>:permno)

  #fix the end date
  ccm.linkenddt .= (s::MString->
    coalesce(s,"")=="E" ? lastdatestring : s).(ccm.linkenddt)

  #below didn't work for some reason
  #filter!(r::DataFrameRow->((!ismissing(r.gvkey)) && (!ismissing(r.lpermno))), ccm)
  ccm = ccm[completecases(ccm, [:gvkey, :permno]),:] #only want valid links
  ccm.linkprim = (s::String->length(s)==1 ? s[1] : error("invalid linkprim in ccm")).(ccm.linkprim)
  ccm = ccm[(ccm.linkprim .== 'P') .| (ccm.linkprim .== 'C'), :] #avoid secondary issues

  #parse the dates
  ccm.ccmeffdate = (s->parseccm(Date, s)).(ccm.linkdt)
  ccm.ccmenddate = (s->parseccm(Date, s)).(ccm.linkenddt)

  ccm = select!(ccm, ccmcolumns)

  return ccm
end

#WARNING WARNING WARNING- a big difference between here and the leverage changes
#code is the use of inclusive eff and end dates (I think that the begin date was exclusive
#in the leverage changes code)
function mergecompccm(comp::DataFrame;
    maxadateinterval::DatePeriod =  PARAM[:compmaxadateinterval],
    defaultadateinterval::DatePeriod =  PARAM[:compdefaultadateinterval],
    )

  #first
  ccm::DataFrame = prepccm()
  #testoutput && CSV.write("output\\ccm.csv", ccm)

  (isempty(comp[nonunique(comp, [:gvkey, :fdate]), :])) || error("non-unique
  gvkey-fdate records detected, dedup prior to merging ccm")
  @assert comp[nonunique(comp, [:gvkey, :adate]),:] |> isempty

  #sort and preliminary contract
  @assert issorted(comp, [:gvkey, :fdate])
  @assert issorted(comp, [:gvkey, :adate])

  #lead the adate to form the announcement interval, and fix intervals where the end is missing
  #or the interval is too long
  comp.beginadate = comp.adate |> deepcopy
  leadwithin2!(comp, [:beginadate], :gvkey, date=:fdate)
  comp.endadate = incrementdt.(comp.Nbeginadate, Day(-1)) #use inclusive, non-overlapping ranges
  #rename!(comp, :Nbeginadate => :endadate)
  badenddates = view(comp, (comp.endadate .=== missing) .|
    (comp.beginadate .+ maxadateinterval .< comp.endadate), [:beginadate, :endadate])
  badenddates.endadate .=  badenddates.beginadate .+ defaultadateinterval

  #integrity checks
  @assert comp[nonunique(comp, [:gvkey, :endadate]),:] |> isempty
  @assert (comp.endadate .≥ comp.beginadate) |> all
  @assert (comp.endadate .≤ comp.beginadate .+ maxadateinterval) |> all
  @assert issorted(comp, [:gvkey, :endadate])

  comp = innerjoin(comp, ccm, on=:gvkey)
  sort!(comp, [:gvkey, :endadate])
  comp = filter!(comp) do r #drop rows with no interseciton between effective link ranges and adate interals
    (r.ccmeffdate ≤ r.beginadate) &&
    ((r.ccmenddate === missing) || (r.ccmenddate ≥ r.endadate))
  end
  #want the intersection of the link range and the announcement range
  comp.effdate = max.(comp.beginadate, comp.ccmeffdate)
  comp.enddate = min.(comp.endadate, comp.ccmenddate)
  comp.effrange = ((effdate,enddate)->effdate:Day(1):enddate).(comp.effdate, comp.enddate)
  @assert (comp.effdate .≤ comp.enddate) |> all

  #date range checks
  #use the intersection of the link range and the announcement range
  adateranges = ((beginadate, endadate)->
    beginadate:Day(1):endadate).(comp.beginadate,comp.endadate)
  ccmdateranges = ((ccmeffdate, ccmenddate)->
    ccmeffdate:Day(1):ccmenddate).(comp.ccmeffdate,comp.ccmenddate)
  effrangecheck = ((arange, ccmrange)->intersect(arange, ccmrange)).(adateranges, ccmdateranges)
  @assert (comp.effrange .== effrangecheck) |> all

  #integrity checks on the dates
  @assert view(comp, (comp.effdate .=== missing) .| (comp.enddate .=== missing), :) |> isempty
  @assert (comp.enddate .≥ comp.effdate) |> all
  @assert (comp.enddate .≤ comp.effdate .+ maxadateinterval) |> all

  #validates the exclusivity of dates within a grouping scheme
  function validatedateswithingrouping(Fgrp::Symbol)
    @assert issorted(comp, [Fgrp, :effdate])
    @assert issorted(comp, [Fgrp, :enddate])

    #check that the intervals are exclusive
    leadwithin2!(comp, [:effdate], Fgrp, date=:fdate)
    lagwithin2!(comp, [:enddate], Fgrp, date=:fdate)

    #check that the start of the next period is always greater than the end of the preceding
    completeNeffdate = view(comp, completecases(comp, [:Neffdate]), :)
    @assert (completeNeffdate.Neffdate .> completeNeffdate.enddate) |> all

    #same as before but in reverse- check end of preceding is less than the start of the current
    completeLenddate = view(comp, completecases(comp, [:Lenddate]), :)
    @assert (completeLenddate.Lenddate .< completeLenddate.effdate) |> all

    select!(comp, Not([:Neffdate, :Lenddate]))
  end

  #validate with the current gvkey groupings
  validatedateswithingrouping(:gvkey)

  #repeat the validaiton excercise, but this time with :permno
  sort!(comp, [:permno, :fdate])
  validatedateswithingrouping(:permno)

  return comp
end

function parseccm(::Type{T}, s::String;) where T

  #if it doesn't parse to the right type, set to missing
  v::Union{T,Missing} = something(tryparse(T, s), missing)

  return v
end

#the date case
parseccm(::Type{Date}, s::String;
    ccmdateformat::DateFormat = dateformat"yyyymmdd") = Dates.Date(s, ccmdateformat)

#helper methods and special cases
parseccm(::Type{<:Any}, v::Missing) = missing
parseccm(::Type{Date}, i::Int) = parseccm(Date, "$i")
parseccm(::Type{Symbol}, s::String) = Symbol(s)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# reminder, to type H*, do H\^+
cd("/Users/Frank/.julia/dev/Bridge")
outdir="/Users/Frank/Sync/DOCUMENTS/onderzoek/code/diffbridges/out_nclar/"

using Bridge, StaticArrays, Distributions
using Test, Statistics, Random, LinearAlgebra
using Bridge.Models
using DelimitedFiles
using DataFrames
using CSV

T = 1.0
dt = 1/5000
τ(T) = (x) ->  x * (2-x/T)
tt = τ(T).(0.:dt:T)

sk = 0 # skipped in evaluating loglikelihood

νHparam = false
generate_data = true

# settings in case of νH - parametrisation
ϵ = 10^(-3)
Σdiagel = 10^(-10)

# settings sampler
iterations = 5*10^4
skip_it = 1000
subsamples = 0:skip_it:iterations

ρ = obs_scheme=="full" ? 0.85 : 0.95

########
#########
n = 3 # number of particles

# q are (p)ositions, p are (m)oments ;-)
qs = 1:2:2n
ps = 2:2:2n

q(i) = 2i - 1
p(i) = 2i

# Gaussian kernel
kernel(x, P) = 1/(2*π*P.a)^(length(x)/2)*exp(-norm(x)^2/(2*P.a))

const Point = SArray{Tuple{2},Float64,1,2}       # point in 2
const Unc = SArray{Tuple{2,2},Float64,2,4}     # Matrix presenting uncertainty
const NP = SVector{16,Float64}

symmetrize!(H) = H .= (H+H')/2

fll(B::Matrix{Unc}) = [B[(i+1)÷2,(j+1)÷2][mod1(i,2),mod1(j,2)] for i in 1:4n, j in 1:4n]
fll(v::Vector{Point}, n = n) = [v[(i+1)÷2][mod1(i,2)] for i in 1:4n]
sps(v::Vector) = [Point(v[i], v[i+1]) for i in 1:2:4n-1]
sps(B::Matrix) = [Unc(B[i, j], B[i+1, j], B[i, j+1], B[i+1,j+1] ) for i in 1:2:4n-1, j in 1:2:4n-1]

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

m, d = size(L)
Σ = SMatrix{m,m}(Σdiagel*I)

struct Landmarks <: ContinuousTimeProcess{Point}
    a::Float64 # kernel parameter
    σ::Float64 # noise level
    λ::Float64 # mean reversion
    n::Int
end

# specify auxiliary process
struct LandmarksTilde <: ContinuousTimeProcess{Point}
    a::Float64 # kernel parameter
    σ::Float64 # noise level
    λ::Float64 # mean reversion
    qT::Vector{Point}
    n::Int
end

aa = 10.0 ; λ = 0.5 ; si = 2.0

P = Landmarks(aa, si, λ, n)

function Bridge.b!(t, x, out, P::Landmarks)
    @. out = zero(Point)
    for i in 1:P.n
        for j in 1:P.n
            out[q(i)] += 0.5*x[p(j)]*kernel(x[q(i)] - x[q(j)], P)
            # heath bath
            out[p(i)] += -P.λ*0.5*x[p(j)]*kernel(x[q(i)] - x[q(j)], P) +
                1/(2*P.a) * dot(x[p(i)], x[p(j)]) * (x[q(i)]-x[q(j)])*kernel(x[q(i)] - x[q(j)], P)
        end
    end
    fll(out)
end
Bridge.b(t, x, P::Union{Landmarks, LandmarksTilde}) = Bridge.b!(t, x, zeros(Point, 2*P.n), P)

Bridge.b(1.0,rand(n*2),P)

Bridge.σ(t, x, P::NclarDiffusion) = ℝ{3}(0.0, 0.0, P.σ)
Bridge.constdiff(::Landmarks) = true

function Bridge._b!((i,s), x, out, Po::LandmarksBridge)
    Bridge.b!(Po.tt[i], x, out, Bridge.P(Po)) # original drift
    out2 = zero(out)
    Bridge.ri!(i, x, out2, Po)
    for i in 1:npoints(P) #acts only on p coordinates
        out[p(i)] += Bridge.P(Po).σ^2 * out2[p(i)]
    end
    out
end

Bridge.bitilde!(i, x, out, Po::LandmarksBridge) = Bridge.b!(Po.tt[i], x, out, Bridge.Pt(Po))

function Bridge.β(t, P::LandmarksTilde)
    zeros(Point, 2npoints(P))
end
function Bridge.B(t, P::LandmarksTilde)
    out = zeros(Unc, 2npoints(P), 2npoints(P))
    for i in 1:npoints(P)
        for j in 1:npoints(P)
            out[q(i),p(j)] += 0.5*one(Unc)*kernel(P.qT[i] - P.qT[j], P)
            out[p(i),p(j)] += -0.5*one(Unc)*P.λ*kernel(P.qT[i] - P.qT[j], P)
        end
    end
    out
end

Bridge.σ(t, x, P::NclarDiffusionAux) = ℝ{3}(0.0,0.0, P.σ)
Bridge.constdiff(::LandmarksTilde) = true
Bridge.b(t, x, P::NclarDiffusionAux) = Bridge.B(t,P) * x + Bridge.β(t,P)
Bridge.a(t, P::NclarDiffusionAux) = Bridge.σ(t,0,P) * Bridge.σ(t, 0, P)'

q0 = x0[qs] .+ sqrt(Σ)*randn(Point, n)
qT = X.yy[qs, end] .+ sqrt(Σ)*randn(Point, n)
pT = X.yy[2:2:end, end]

Pt = LandmarksTilde(aa, si, λ, qT, n)

# Solve Backward Recursion
Po = νHparam ? Bridge.PartialBridgeνH(tt, P, Pt, L, ℝ{m}(v),ϵ, Σ) : Bridge.PartialBridge(tt, P, Pt, L, ℝ{m}(v), Σ)

####################### MH algorithm ###################
W = sample(tt, Wiener())
X = solve(Euler(), x0, W, P)
Xo = copy(X)
solve!(Euler(),Xo, x0, W, Po)

solve!(Euler(),X, x0, W, Po)
ll = llikelihood(Bridge.LeftRule(), X, Po,skip=sk)

# further initialisation
Wo = copy(W)
W2 = copy(W)
XX = Any[]
if 0 in subsamples
    push!(XX, copy(X))
end

acc = 0

for iter in 1:iterations
    # Proposal
    sample!(W2, Wiener())
    #ρ = rand(Uniform(0.95,1.0))
    Wo.yy .= ρ*W.yy + sqrt(1-ρ^2)*W2.yy
    solve!(Euler(),Xo, x0, Wo, Po)

    llo = llikelihood(Bridge.LeftRule(), Xo, Po,skip=sk)
    print("ll $ll $llo, diff_ll: ",round(llo-ll,3))

    if log(rand()) <= llo - ll
        X.yy .= Xo.yy
        W.yy .= Wo.yy
        ll = llo
        print("✓")
        acc +=1
    end
    println()
    if iter in subsamples
        push!(XX, copy(X))
    end
end

@info "Done."*"\x7"^6

# write mcmc iterates to csv file

fn = outdir*"iterates-"*obs_scheme*".csv"
f = open(fn,"w")
head = "iteration, time, component, value \n"
write(f, head)
iterates = [Any[s, tt[j], d, XX[i].yy[j][d]] for d in 1:3, j in 1:length(X), (i,s) in enumerate(subsamples) ][:]
writecsv(f,iterates)
close(f)

ave_acc_perc = 100*round(acc/iterations,2)

# write info to txt file
fn = outdir*"info-"*obs_scheme*".txt"
f = open(fn,"w")
write(f, "Choice of observation schemes: ",obs_scheme,"\n")
write(f, "Easy conditioning (means going up to 1 for the rough component instead of 2): ",string(easy_conditioning),"\n")
write(f, "Number of iterations: ",string(iterations),"\n")
write(f, "Skip every ",string(skip_it)," iterations, when saving to csv","\n\n")
write(f, "Starting point: ",string(x0),"\n")
write(f, "End time T: ", string(T),"\n")
write(f, "Endpoint v: ",string(v),"\n")
write(f, "Noise Sigma: ",string(Σ),"\n")
write(f, "L: ",string(L),"\n\n")
write(f,"Mesh width: ",string(dt),"\n")
write(f, "rho (Crank-Nicholsen parameter: ",string(ρ),"\n")
write(f, "Average acceptance percentage: ",string(ave_acc_perc),"\n\n")
write(f, "Backward type parametrisation in terms of nu and H? ",string(νHparam),"\n")
close(f)

println("Average acceptance percentage: ",ave_acc_perc,"\n")
println(obs_scheme)
println("Parametrisation of nu and H? ", νHparam)


text:
# auxiliary/error.jl
#
# An exception type for reporting errors in the index specificatino

immutable IndexError{S<:AbstractString} <: Exception
    msg::S
end

checkindices(A, IA) =
    length(IA) == numind(A) || throw(IndexError("invalid number of indices: $IA"))


text:
src = "../src"
if !in(src,LOAD_PATH)
    push!(LOAD_PATH, src)
end
using HTSLIB
using Base.Test

#=
include("bam.jl")
include("user/bam.jl")
include("user/header.jl")
include("user/sam.jl")
include("user/query.jl")
=#

bios = HTSLIB.bam_open("data/100_sort.bam","rb")

while !HTSLIB.eof(bios)
    line = HTSLIB.readline(bios)
    println(line)
end

HTSLIB.close(bios)


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

using Markdown
using InteractiveUtils

# ╔═╡ be75fc5e-3247-11eb-1c9f-6f30ea72331c
begin
	import Pkg; Pkg.add("ForwardDiff")
	Pkg.add("Plots")
	using Plots
end

# ╔═╡ 5f085b68-323f-11eb-3ee9-13b634011eac
md"""
# Julia Essentials
### Solutions to the Exercise
"""

# ╔═╡ e6599ddc-323f-11eb-3f27-a7f05eb26ae5
md"""
##### `Exercise 1 - 4 are quite naive and will only help you to get familiarity with the platform, you are advised to permorn on your own` """

# ╔═╡ e57c7f4c-323f-11eb-26dd-5b88be11bed5
html"""

"""

# ╔═╡ c994bf8e-3245-11eb-259f-a791ae0734dc
md"""
##### Exercise 5
The Julia libraries include functions for interpolation and approximation.

Nevertheless, let’s write our own function approximation routine as an exercise.

In particular, write a function linapprox that takes as arguments

A function f mapping some interval [a,b] into ℝ.
two scalars a and b providing the limits of this interval.
An integer n determining the number of grid points.
A number x satisfying a ≤ x ≤ b.
and returns the piecewise linear interpolation of f at x, based on n evenly spaced grid points a = point[1] < point[2] < ... < point[n] = b.

Aim for clarity, not efficiency.

Hint: use the function range to linearly space numbers."""

# ╔═╡ 20cd089a-3246-11eb-368f-fdbc5129f16c
function linapprox(f,a,b,n,x)
	Δ = (b-a)/(n-1)
	k = (x-a)÷Δ
	u,v = k*Δ +a,(k+1)*Δ +a
	δ= (f(v)-f(u))/(Δ)
	f(u) + (x-u)*δ
end

# ╔═╡ 4a812cde-3249-11eb-3d46-c1937d75ef1e
html"""

"""

# ╔═╡ 42fd5d3c-3248-11eb-0dac-43b98063620e
md"""
##### Exercise 6
`should be easy, try to do it yourself`
"""

# ╔═╡ 58b6b208-3249-11eb-110e-1f2d5f82756c
html"""

"""

# ╔═╡ 46a8a40c-3249-11eb-345a-8143e3ac48cc
md"""
##### Exercise 7
Redo Exercise 5 except

Pass in a range instead of the a, b, and n. Test with a range such as nodes = -1.0:0.5:1.0.
Instead of the while used in the solution to Exercise 5, find a better way to efficiently bracket the x in the nodes.
Hints:

Rather than the signature as function linapprox(f, a, b, n, x), it should be called as function linapprox(f, nodes, x).
step(nodes), length(nodes), nodes[1], and nodes[end] may be useful.
Type ?÷ into jupyter to explore quotients from Euclidean division for more efficient bracketing."""

# ╔═╡ 78f845ea-3249-11eb-3ac1-fb51649acd04
function linapprox(f,nodes,x)
	a,b = nodes[1],nodes[end]
	n = length(nodes)
	linapprox(f,a,b,n,x)
end

# ╔═╡ 9dbef2fe-3247-11eb-2827-1b0878e2cd0c
let
	f_ex5(x) = x^2
	g_ex5(x) = linapprox(f_ex5, -1, 1, 3, x)
	x_grid = range(-1.0, 1.0, length = 100)
	y_vals = f_ex5.(x_grid)
	y = g_ex5.(x_grid)
	plot(x_grid, y_vals, label = "true")
	plot!(x_grid, y, label = "approximation")
end

# ╔═╡ 9e334f98-324a-11eb-024f-93e55e491e58
let
	f_ex5(x) = x^2
	g_ex5(x) = linapprox(f_ex5, -1:5, x)
	x_grid = range(-1.0, 5.0, length = 100)
	y_vals = f_ex5.(x_grid)
	y = g_ex5.(x_grid)
	plot(x_grid, y_vals, label = "true")
	plot!(x_grid, y, label = "approximation")
end

# ╔═╡ Cell order:
# ╟─5f085b68-323f-11eb-3ee9-13b634011eac
# ╟─e6599ddc-323f-11eb-3f27-a7f05eb26ae5
# ╟─e57c7f4c-323f-11eb-26dd-5b88be11bed5
# ╟─c994bf8e-3245-11eb-259f-a791ae0734dc
# ╠═20cd089a-3246-11eb-368f-fdbc5129f16c
# ╟─9dbef2fe-3247-11eb-2827-1b0878e2cd0c
# ╟─4a812cde-3249-11eb-3d46-c1937d75ef1e
# ╟─42fd5d3c-3248-11eb-0dac-43b98063620e
# ╟─58b6b208-3249-11eb-110e-1f2d5f82756c
# ╟─46a8a40c-3249-11eb-345a-8143e3ac48cc
# ╠═78f845ea-3249-11eb-3ac1-fb51649acd04
# ╟─9e334f98-324a-11eb-024f-93e55e491e58
# ╟─be75fc5e-3247-11eb-1c9f-6f30ea72331c

max_stars_repo_path: getting started with julia/Julia Esswntials.jl

max_stars_repo_name: nehamsoni/Quantitative-Economics

max_issues_repo_name: nehamsoni/Quantitative-Economics

max_forks_repo_name: nehamsoni/Quantitative-Economics


text:
"""
    SrcDir(; file="$(contractuser(default_file("src", "module.jl")))")

Creates a module entrypoint.

## Keyword Arguments
- `file::AbstractString`: Template file for `src/<module>.jl`.
"""
@with_kw_noshow mutable struct SrcDir <: BasicPlugin
    file::String = default_file("src", "module.jl")
    destination::String = ""
end

Base.:(==)(a::SrcDir, b::SrcDir) = a.file == b.file

# Don't display the destination field.
function Base.show(io::IO, ::MIME"text/plain", p::SrcDir)
    indent = get(io, :indent, 0)
    print(io, repeat(' ', indent), "SrcDir:")
    print(io, "\n", repeat(' ', indent + 2), "file: ", show_field(p.file))
end

source(p::SrcDir) = p.file
destination(p::SrcDir) = p.destination
view(::SrcDir, ::Template, pkg::AbstractString) = Dict("PKG" => pkg)

# Update the destination now that we know the package name.
# Kind of hacky, but oh well.
function prehook(p::SrcDir, t::Template, pkg_dir::AbstractString)
    p.destination = joinpath("src", basename(pkg_dir) * ".jl")
end


text:
#Function for mimicking Matlab's combvec. Not named combvec because I have seen at least one Julia package with a function called combvec.
#25/07/2017

function combvechack(args...)
  numargs = length(args)
  Mvec = zeros(numargs)
  Nvec = zeros(numargs)
  kk = 1
  for arg in args
    if ndims(arg)==1
      Nvec[kk] = 1
      Mvec[kk] = length(arg)
    else
      sA = size(arg)
      Mvec[kk] = sA[1]
      Nvec[kk] = sA[2]
    end
    kk = kk + 1
  end
  Nprod = prod(Nvec)
  outmatrix = zeros(sum(Mvec),Nprod)
  Mstart = 1
  kk = 1
  for arg in args
    Mend = Int(sum(Mvec[1:kk]))
    val1 = rectpulse(arg',Int(Nprod/prod(Nvec[kk:end])))'
    outmatrix[Mstart:Mend,:] = repmat(val1,1,Int(Nprod/prod(Nvec[1:kk])))
    Mstart = Mend + 1
    kk = kk + 1
  end
  return outmatrix
end

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
function AbstractPlotting.backend_display(::GLBackend, scene::Scene)
    screen = global_gl_screen(size(scene), AbstractPlotting.use_display[])
    display_loading_image(screen)
    AbstractPlotting.backend_display(screen, scene)
    return screen
end

function AbstractPlotting.backend_display(screen::Screen, scene::Scene)
    empty!(screen)
    # So, the GLFW window events are not guarantee to fire
    # when we close a window, so we ensure this here!
    window_open = events(scene).window_open
    on(screen.window_open) do open
        window_open[] = open
    end
    register_callbacks(scene, screen)
    pollevents(screen)
    insertplots!(screen, scene)
    pollevents(screen)
    return
end

"""
    scene2image(scene::Scene)

Buffers the `scene` in an image buffer.
"""
function scene2image(scene::Scene)
    old = WINDOW_CONFIG.pause_rendering[]
    try
        WINDOW_CONFIG.pause_rendering[] = true
        screen = global_gl_screen(size(scene), false)
        AbstractPlotting.backend_display(screen, scene)
        return AbstractPlotting.colorbuffer(screen), screen
    finally
        WINDOW_CONFIG.pause_rendering[] = old
    end
end

raw_io(io::IO) = io
raw_io(io::IOContext) = raw_io(io.io)

function AbstractPlotting.backend_show(::GLBackend, io::IO, m::MIME"image/png", scene::Scene)
    img, screen = scene2image(scene)
    FileIO.save(FileIO.Stream(FileIO.format"PNG", raw_io(io)), img)
    return screen
end

function AbstractPlotting.backend_show(::GLBackend, io::IO, m::MIME"image/jpeg", scene::Scene)
    img, screen = scene2image(scene)
    FileIO.save(FileIO.Stream(FileIO.format"JPEG", raw_io(io)), img)
    return screen
end


text:
module NNlib
using Pkg
using Requires
using NNPACK_jll

# Include APIs
include("dim_helpers.jl")

if isdefined(NNPACK_jll, :libnnpack)
  include("nnpack/NNPACK.jl")
else
  @warn "NNPACK not available for your platform: " *
        "$( Pkg.BinaryPlatforms.platform_name(Pkg.BinaryPlatforms.platform_key_abi()))" *
        "($( Pkg.BinaryPlatforms.triplet(Pkg.BinaryPlatforms.platform_key_abi())))
        You will be able to use only the default Julia NNlib backend"
  is_nnpack_available() = false
end

include("activation.jl")
include("softmax.jl")
include("batched/batchedmul.jl")
include("gemm.jl")
include("conv.jl")
include("pooling.jl")

## Include implementations
include("impl/padding_edges.jl")

# Direct implementations of convolutional and depthwise-convolutional algorithms
include("impl/conv_direct.jl")
include("impl/depthwiseconv_direct.jl")
# im2col implementations of convolutional and depthwise-convolutional algorithms
include("impl/conv_im2col.jl")
include("impl/depthwiseconv_im2col.jl")

# Direct implementations of pooling
include("impl/pooling_direct.jl")

end # module NNlib


text:
# Point2D mutable struct

struct Point2 <: IncrementalInference.InferenceVariable
  dims::Int
  labels::Vector{String}
  Point2() = new(2, String["";])
end

mutable struct PriorPoint2D <: IncrementalInference.FunctorSingleton
  mv::MvNormal
  W::Array{Float64,1}
  PriorPoint2D() = new()
  PriorPoint2D(mu, cov, W) = new(MvNormal(mu, cov), W)
end
function getSample(p2::PriorPoint2D, N::Int=1)
  return (rand(p2.mv, N),)
end

mutable struct Point2DPoint2DRange <: IncrementalInference.FunctorPairwiseMinimize #Pairwise
    Zij::Vector{Float64} # bearing and range hypotheses as columns
    Cov::Float64
    W::Vector{Float64}
    Point2DPoint2DRange() = new()
    Point2DPoint2DRange(x...) = new(x[1],x[2],x[3])
end
function (pp2r::Point2DPoint2DRange)(
      res::Array{Float64},
      idx::Int,
      meas::Tuple, # Array{Float64,2},
      xi::Array{Float64,2},
      lm::Array{Float64,2} )
  #
  # TODO -- still need to add multi-hypotheses support here
  # this is the noisy range
  z = pp2r.Zij[1]+meas[1][1,idx]
  # theta = meas[2]
  XX = lm[1,idx] - (z*cos(meas[2][idx]) + xi[1,idx])
  YY = lm[2,idx] - (z*sin(meas[2][idx]) + xi[2,idx])
  res[1] = XX^2 + YY^2
  # nothing
end
function getSample(pp2::Point2DPoint2DRange, N::Int=1)
  return (pp2.Cov*randn(1,N),  2*pi*rand(N))
end

mutable struct Point2DPoint2D <: BetweenPoses
    Zij::Distribution
    Point2DPoint2D() = new()
    Point2DPoint2D(x) = new(x)
end
function (pp2r::Point2DPoint2D)(
      res::Array{Float64},
      idx::Int,
      meas::Tuple, # Array{Float64,2},
      xi::Array{Float64,2},
      xj::Array{Float64,2} )
  #
  # TODO -- still need to add multi-hypotheses support here
  res[1]  = meas[1][1,idx] - (xj[1,idx] - xi[1,idx])
  res[2]  = meas[1][2,idx] - (xj[2,idx] - xi[2,idx])
  nothing
end
function getSample(pp2::Point2DPoint2D, N::Int=1)
  return (rand(pp2.Zij,N),  )
end

mutable struct PriorPoint2DensityNH <: IncrementalInference.FunctorSingletonNH
  belief::BallTreeDensity
  nullhypothesis::Distributions.Categorical
  PriorPoint2DensityNH() = new()
  PriorPoint2DensityNH(belief, p) = new(belief, Distributions.Categorical(p))
end
function getSample(p2::PriorPoint2DensityNH, N::Int=1)
  return (rand(p2.belief, N), )
end
mutable struct PackedPriorPoint2DensityNH <: IncrementalInference.PackedInferenceType
    rpts::Vector{Float64} # 0rotations, 1translation in each column
    rbw::Vector{Float64}
    dims::Int
    nh::Vector{Float64}
    PackedPriorPoint2DensityNH() = new()
    PackedPriorPoint2DensityNH(x1,x2,x3, x4) = new(x1, x2, x3, x4)
end
function convert(::Type{PriorPoint2DensityNH}, d::PackedPriorPoint2DensityNH)
  return PriorPoint2DensityNH(
            kde!(EasyMessage( reshapeVec2Mat(d.rpts, d.dims), d.rbw)),
            Distributions.Categorical(d.nh)  )
end
function convert(::Type{PackedPriorPoint2DensityNH}, d::PriorPoint2DensityNH)
  return PackedPriorPoint2DensityNH( getPoints(d.belief)[:], getBW(d.belief)[:,1], Ndim(d.belief), d.nullhypothesis.p )
end

# Old evalPotential functions
function evalPotential(prior::PriorPoint2D, Xi::Array{Graphs.ExVertex,1}; N::Int=100)#, from::Int)
    return rand(prior.mv, N)
end

# Solve for Xid, given values from vertices [Xi] and measurement rho
function evalPotential(rho::Point2DPoint2DRange, Xi::Array{Graphs.ExVertex,1}, Xid::Int)
  fromX, ret = nothing, nothing
  if Xi[1].index == Xid
    fromX = getVal( Xi[2] )
    ret = deepcopy(getVal( Xi[1] )) # carry pose yaw row over if required
  elseif Xi[2].index == Xid
    fromX = getVal( Xi[1] )
    ret = deepcopy(getVal( Xi[2] )) # carry pose yaw row over if required
  end
  r,c = size(fromX)
  theta = 2*pi*rand(c)
  noisy = rho.Cov*randn(c) + rho.Zij[1]

  for i in 1:c
    ret[1,i] = noisy[i]*cos(theta[i]) + fromX[1,i]
    ret[2,i] = noisy[i]*sin(theta[i]) + fromX[2,i]
  end
  return ret
end

# ---------------------------------------------------------

mutable struct PackedPriorPoint2D  <: IncrementalInference.PackedInferenceType
    mu::Array{Float64,1}
    vecCov::Array{Float64,1}
    dimc::Int
    W::Array{Float64,1}
    PackedPriorPoint2D() = new()
    PackedPriorPoint2D(x...) = new(x[1], x[2], x[3], x[4])
end

passTypeThrough(d::FunctionNodeData{Point2DPoint2DRange}) = d

function convert(::Type{PriorPoint2D}, d::PackedPriorPoint2D)
  Cov = reshapeVec2Mat(d.vecCov, d.dimc)
  return PriorPoint2D(d.mu, Cov, d.W)
end
function convert(::Type{PackedPriorPoint2D}, d::PriorPoint2D)
  v2 = d.mv.Σ.mat[:];
  return PackedPriorPoint2D(d.mv.μ, v2, size(d.mv.Σ.mat,1), d.W)
end

# no longer needed -- using multiple dispatch
# function convert(::Type{PackedFunctionNodeData{PackedPriorPoint2D}}, d::FunctionNodeData{PriorPoint2D})
#   return PackedFunctionNodeData{PackedPriorPoint2D}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
#           string(d.frommodule), convert(PackedPriorPoint2D, d.fnc))
# end
# function convert(::Type{FunctionNodeData{PriorPoint2D}}, d::PackedFunctionNodeData{PackedPriorPoint2D})
#   return FunctionNodeData{PriorPoint2D}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
#           Symbol(d.frommodule), convert(PriorPoint2D, d.fnc))
# end
# function FNDencode(d::FunctionNodeData{PriorPoint2D})
#   return convert(PackedFunctionNodeData{PackedPriorPoint2D}, d)
# end
# function FNDdecode(d::PackedFunctionNodeData{PackedPriorPoint2D})
#   return convert(FunctionNodeData{PriorPoint2D}, d)
# end

mutable struct PackedPoint2DPoint2DRange  <: IncrementalInference.PackedInferenceType
    Zij::Vector{Float64} # bearing and range hypotheses as columns
    Cov::Float64
    W::Vector{Float64}
    PackedPoint2DPoint2DRange() = new()
    PackedPoint2DPoint2DRange(x...) = new(x[1],x[2],x[3])
    PackedPoint2DPoint2DRange(x::Point2DPoint2DRange) = new(x.Zij,x.Cov,x.W)
end
function convert(::Type{PackedPoint2DPoint2DRange}, d::Point2DPoint2DRange)
  return PackedPoint2DPoint2DRange(d)
end
function convert(::Type{Point2DPoint2DRange}, d::PackedPoint2DPoint2DRange)
  return Point2DPoint2DRange(d.Zij, d.Cov, d.W)
end

mutable struct PackedPoint2DPoint2D <: IncrementalInference.PackedInferenceType
    mu::Vector{Float64}
    sigma::Vector{Float64}
    sdim::Int
    PackedPoint2DPoint2D() = new()
    PackedPoint2DPoint2D(x, y, d) = new(x,y,d)
end
function convert(::Type{Point2DPoint2D}, d::PackedPoint2DPoint2D)
  return Point2DPoint2D( MvNormal(d.mu, reshapeVec2Mat(d.sigma, d.sdim)) )
end
function convert(::Type{PackedPoint2DPoint2D}, d::Point2DPoint2D)
  return PackedPoint2DPoint2D( d.Zij.μ, d.Zij.Σ.mat[:], size(d.Zij.Σ.mat,1) )
end


text:
struct LTInitialBelief
    robot_init::Union{Coord, Nothing}
    floor::Floor
end

sampletype(::Type{LTInitialBelief}) = LTState
function iterator(b::LTInitialBelief)
    states = LTState[]
    for (x,y) in product(1:b.floor.n_cols, 1:b.floor.n_rows)
        if isnull(b.robot_init)
            for (rx, ry) in product(1:b.floor.n_cols, 1:b.floor.n_rows)
                push!(states, LTState(Coord(rx,ry), Coord(x,y), false))
            end
        else
            push!(states, LTState(get(b.robot_init), Coord(x,y), false))
        end
    end
    return states
end

function Random.rand(rng::AbstractRNG, b::LTInitialBelief)
    opp = Coord(rand(rng, 1:b.floor.n_cols), rand(rng, 1:b.floor.n_rows))
    if b.robot_init === nothing
        rob = Coord(rand(rng, 1:b.floor.n_cols), rand(rng, 1:b.floor.n_rows))
    else
        rob = b.robot_init
    end
    return LTState(rob, opp, false)
end

function Distributions.pdf(b::LTInitialBelief, s::LTState)
    if s.terminal
        return 0.0
    else
        if b.robot_init === nothing
            return (1/n_pos(b.floor))^2
        else
            if s.robot == b.robot_init
                return 1/n_pos(b.floor)
            else
                return 0.0
            end
        end
    end
end

POMDPs.initialstate(p::LaserTagPOMDP) = LTInitialBelief(p.robot_init, p.floor)


text:
function save(f::Stream{format"STL_ASCII"}, mesh::AbstractMesh)
    io = stream(f)
    points = decompose(Point3f, mesh)
    faces = decompose(GLTriangleFace, mesh)
    normals = decompose_normals(mesh)

    n_points = length(points)
    n_faces = length(faces)

    # write the header
    write(io,"solid vcg\n")
    # write the data
    for i = 1:n_faces
        f = faces[i]
        n = normals[f][1] # TODO: properly compute normal(f)
        v1, v2, v3 = points[f]
        @printf io "  facet normal %e %e %e\n" n[1] n[2] n[3]
        write(io,"    outer loop\n")

        @printf io "      vertex  %e %e %e\n" v1[1] v1[2] v1[3]
        @printf io "      vertex  %e %e %e\n" v2[1] v2[2] v2[3]
        @printf io "      vertex  %e %e %e\n" v3[1] v3[2] v3[3]

        write(io,"    endloop\n")
        write(io,"  endfacet\n")
    end
    write(io,"endsolid vcg\n")
end

show(io::IO, ::MIME"model/stl+ascii", mesh::AbstractMesh) = save(io, mesh)

function save(f::Stream{format"STL_BINARY"}, mesh::AbstractMesh)
    io = stream(f)
    points = decompose(Point3f, mesh)
    faces = decompose(GLTriangleFace, mesh)
    normals = decompose_normals(mesh)
    n_faces = length(faces)
    # Implementation made according to https://en.wikipedia.org/wiki/STL_%28file_format%29#Binary_STL
    for i in 1:80 # write empty header
        write(io, 0x00)
    end

    write(io, UInt32(n_faces)) # write triangle count
    for i = 1:n_faces
        f = faces[i]
        n = normals[f][1] # TODO: properly compute normal(f)
        triangle = points[f]
        foreach(j-> write(io, n[j]), 1:3)
        for point in triangle
            foreach(p-> write(io, p), point)
        end
        write(io, 0x0000) # write 16bit empty bit
    end
end

function load(fs::Stream{format"STL_BINARY"}; facetype=GLTriangleFace,
              pointtype=Point3f, normaltype=Vec3f)
    #Binary STL
    #https://en.wikipedia.org/wiki/STL_%28file_format%29#Binary_STL
    io = stream(fs)
    read(io, 80) # throw out header
    triangle_count = read(io, UInt32)

    faces = Array{facetype}(undef, triangle_count)
    vertices = Array{pointtype}(undef, triangle_count * 3)
    normals = Array{normaltype}(undef, triangle_count * 3)

    i = 0
    while !eof(io)
        faces[i+1] = GLTriangleFace(i * 3 + 1, i * 3 + 2, i * 3 + 3)
        normal = (read(io, Float32), read(io, Float32), read(io, Float32))

        normals[i*3+1] = normaltype(normal...)
        normals[i*3+2] = normals[i*3+1] # hurts, but we need per vertex normals
        normals[i*3+3] = normals[i*3+1]

        vertices[i*3+1] = pointtype(read(io, Float32), read(io, Float32), read(io, Float32))
        vertices[i*3+2] = pointtype(read(io, Float32), read(io, Float32), read(io, Float32))
        vertices[i*3+3] = pointtype(read(io, Float32), read(io, Float32), read(io, Float32))

        skip(io, 2) # throwout 16bit attribute
        i += 1
    end

    return Mesh(meta(vertices; normals=normals), faces)
end

function load(fs::Stream{format"STL_ASCII"}; facetype=GLTriangleFace,
              pointtype=Point3f, normaltype=Vec3f, topology=false)
    #ASCII STL
    #https://en.wikipedia.org/wiki/STL_%28file_format%29#ASCII_STL
    io = stream(fs)

    points = pointtype[]
    faces = facetype[]
    normals = normaltype[]

    vert_count = 0
    vert_idx = [0, 0, 0]

    while !eof(io)
        line = split(lowercase(readline(io)))
        if !isempty(line) && line[1] == "facet"
            normal = normaltype(parse.(eltype(normaltype), line[3:5]))
            readline(io) # Throw away outerloop
            for i in 1:3
                vertex = pointtype(parse.(eltype(pointtype),
                                   split(readline(io))[2:4]))
                if topology
                    idx = findfirst(vertices(mesh), vertex)
                end
                if topology && idx != 0
                    vert_idx[i] = idx
                else
                    push!(points, vertex)
                    push!(normals, normal)
                    vert_count += 1
                    vert_idx[i] = vert_count
                end
            end
            readline(io) # throwout endloop
            readline(io) # throwout endfacet
            push!(faces, TriangleFace{Int}(vert_idx...))
        end
    end
    return Mesh(meta(points; normals=normals), faces)
end


text:
"""
a quantum link model of qed
"""
function nonsym_qed_qlm_ham(;m=0.8,a=1,e=1.0,l=1)
    (Xl,Yl,Zl,Il) = nonsym_spintensors(l);
    Pl = Xl + 1im*Yl; Ml = Xl - 1im*Yl;
    (Xf,Yf,Zf,If) = nonsym_spintensors(1//2);
    Pf = Xf + 1im*Yf; Mf = Xf - 1im*Yf;

    data = Array{Union{Missing,typeof(Pf)},3}(missing,2,4,4)

    data[1,1,1] = Il; data[1,end,end] = Il;
    data[2,1,1] = If; data[2,end,end] = If;

    data[1,1,end] = a*e*e*0.5*Zl*Zl
    data[2,1,end] = m*Zf

    pref = -1im/(a*2*sqrt(l*(l+1)))

    data[2,1,2] = pref*Pf
    data[1,2,2] = Pl
    data[2,2,end] = Pf

    data[2,1,3] = -pref*Mf
    data[1,3,3] = Ml
    data[2,3,end] = Mf

    return MPOHamiltonian(data)
end

"""
qed_qlm has a local gauge transform G_j, here we construct G_j^2
useful to find the groundstate of qed_qlm in the G = 0 sector
"""
function qed_qlm_G2(;l=1)
    (Xl,Yl,Zl,Il) = nonsym_spintensors(l);
    Pl = Xl + 1im*Yl; Ml = Xl - 1im*Yl;

    (Xf,Yf,Zf,If) = nonsym_spintensors(1//2);
    Pf = Xf + 1im*Yf; Mf = Xf - 1im*Yf;

    @tensor todecomp[-1 -2 -3;-4 -5 -6]:= Zl[-1 -4]*If[-2 -5]*Il[-3 -6]
    @tensor todecomp[-1 -2 -3;-4 -5 -6]+= Il[-1 -4]*If[-2 -5]*Zl[-3 -6]
    @tensor todecomp[-1 -2 -3;-4 -5 -6]-= Il[-1 -4]*(Zf+If*0.5)[-2 -5]*Il[-3 -6]
    todecomp = todecomp*todecomp

    string = MPSKit.decompose_localmpo(add_util_leg(todecomp))

    nIl = permute(add_util_leg(Il),(1,2),(4,3))
    nIf = permute(add_util_leg(If),(1,2),(4,3))

    data = Array{Union{Missing,typeof(string[1])},3}(missing,2,3,3)

    data[1,1,1] = nIl; data[1,end,end] = nIl;
    data[2,1,1] = nIf; data[2,end,end] = nIf;

    data[1,1,2] = string[1];
    data[2,2,2] = string[2];
    data[1,2,end] = string[3];

    return MPOHamiltonian(data)
end


text:
include("smear_FD.jl")
include("calc_entropy.jl")

function test_main(kT::Float64)
    Nstates = 5
    ev = Array{Float64}(Nstates)
    ev = [-2.4, -1.0, -0.5, -0.2, -0.19]
    efermi = 0.5*(ev[3] + ev[4])
    println("\nkT = ", kT)
    for spinpol in [true,false]
        println("\nspinpol = ", spinpol)
        Focc = smear_FD(ev, efermi, kT, is_spinpol=spinpol)
        for ist = 1:Nstates
            @printf("%18.10f %18.10f\n", ev[ist], Focc[ist])
        end
        @printf("sum(Focc) = %18.10f\n", sum(Focc))
        @printf("Entropy (-TS) = %18.10f\n", calc_entropy(Focc, kT, is_spinpol=spinpol))
    end
end

test_main(0.001)
test_main(0.01)
test_main(0.1)


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using LinearAlgebra
using ValueShapes, Distributions, ArraysOfArrays, ForwardDiff

@testset "test_distribution_transform" begin
    function test_back_and_forth(trg_d::Distribution, src_d::BAT.StdMvDist)
        @testset "transform $(typeof(trg_d).name) <-> $(typeof(src_d).name)" begin
            src_v = rand(src_d)
            prev_ladj = 7.9

            @test @inferred(BAT.apply_dist_trafo(trg_d, src_d, src_v, prev_ladj)) isa NamedTuple{(:v,:ladj)}
            trg_v, trg_ladj = BAT.apply_dist_trafo(trg_d, src_d, src_v, prev_ladj)

            @test BAT.apply_dist_trafo(src_d, trg_d, trg_v, trg_ladj) isa NamedTuple{(:v,:ladj)}
            src_v_reco, prev_ladj_reco = BAT.apply_dist_trafo(src_d, trg_d, trg_v, trg_ladj)

            @test src_v ≈ src_v_reco
            @test prev_ladj ≈ prev_ladj_reco
            @test trg_ladj ≈ logabsdet(ForwardDiff.jacobian(x -> unshaped(BAT.apply_dist_trafo(trg_d, src_d, x, prev_ladj).v), src_v))[1] + prev_ladj
        end
    end

    function test_dist_trafo_moments(trg_d::Distribution, src_d::BAT.StdMvDist)
        @testset "check moments of trafo $(typeof(trg_d).name) <- $(typeof(src_d).name)" begin
            let trg_d = trg_d, src_d = src_d
                X = rand(src_d, 10^5)
                trgxs = (x -> BAT.apply_dist_trafo(trg_d, src_d, x, 0.0).v).(nestedview(X))
                unshaped_trgxs = map(unshaped, trgxs)
                @test isapprox(mean(unshaped_trgxs), mean(unshaped(trg_d)), rtol = 0.1)
                @test isapprox(cov(unshaped_trgxs), cov(unshaped(trg_d)), rtol = 0.1)
            end
        end
    end

    standnorm1 = BAT.StandardMvNormal(1)
    standnorm2 = BAT.StandardMvNormal(2)

    mvnorm = MvNormal([0.3, -2.9], [1.7 0.5; 0.5 2.3])
    beta = Beta(3,1)
    gamma = Gamma(0.1,0.7)
    dirich = Dirichlet([0.1,4])

    test_back_and_forth(beta, standnorm1)
    test_back_and_forth(gamma, standnorm1)
    test_back_and_forth(mvnorm, standnorm2)

    test_dist_trafo_moments(beta, standnorm1)
    test_dist_trafo_moments(gamma, standnorm1)

    test_dist_trafo_moments(mvnorm, standnorm2)
    test_dist_trafo_moments(dirich, standnorm1)

    let
        primary_dist = NamedTupleDist(x = Normal(2), c = 5)
        f = x -> NamedTupleDist(y = Normal(x.x, 3), z = MvNormal([1.3 0.5; 0.5 2.2]))
        trg_d = @inferred(HierarchicalDistribution(f, primary_dist))
        src_d = BAT.StandardMvNormal(totalndof(varshape(trg_d)))
        test_back_and_forth(trg_d, src_d)
        test_dist_trafo_moments(trg_d, src_d)
    end

    #=
    using Cuba
    function integrate_over_unit(density::AbstractDensity)
        vs = varshape(density)
        f_cuba(source_x, y) = y[1] = exp(logvalof(density)(vs(source_x)))
        Cuba.vegas(f_cuba, 1, 1).integral[1]
    end
    =#
end

@testset "bat_transform_defaults" begin
    mvn = @inferred(product_distribution([Normal(-1), Normal(), Normal(1)]))
    uniform_prior = @inferred(product_distribution([Uniform(-3, 1), Uniform(-2, 2), Uniform(-1, 3)]))

    posterior_uniform_prior = @inferred(PosteriorDensity(mvn, uniform_prior))
    posterior_gaussian_prior = @inferred(PosteriorDensity(mvn, mvn))

    @test @inferred(bat_transform(PriorToGaussian(), posterior_uniform_prior)).result.prior.dist == @inferred(BAT.StandardMvNormal(3))
    @test @inferred(bat_transform(PriorToUniform(), posterior_gaussian_prior)).result.prior.dist == @inferred(BAT.StandardMvUniform(3))
    @test @inferred(bat_transform(NoDensityTransform(), posterior_uniform_prior)).result.prior.dist == uniform_prior
    pd = @inferred(product_distribution([Uniform() for i in 1:3]))
    density = @inferred(BAT.DistributionDensity(pd))
    @test @inferred(bat_transform(NoDensityTransform(), density)).result.dist == density.dist

    # ToDo: Improve comparison for bounds so `.dist` is not required here:
    @inferred(bat_transform(PriorToUniform(), convert(AbstractDensity, BAT.StandardUvUniform()))).result.dist == convert(AbstractDensity, BAT.StandardUvUniform()).dist
    @inferred(bat_transform(PriorToUniform(), convert(AbstractDensity, BAT.StandardMvUniform(4)))).result.dist == convert(AbstractDensity, BAT.StandardMvUniform(4)).dist
    @inferred(bat_transform(PriorToGaussian(), convert(AbstractDensity, BAT.StandardUvNormal()))).result.dist == convert(AbstractDensity, BAT.StandardUvNormal()).dist
    @inferred(bat_transform(PriorToGaussian(), convert(AbstractDensity, BAT.StandardMvNormal(4)))).result.dist == convert(AbstractDensity, BAT.StandardMvNormal(4)).dist
end


text:
"""
    compat(version_spec::String="auto", path_to_project::String=pwd(); package::String="", overwrite::Bool=false)

create and update compat in `Project.toml` via command line.

# Intro

By default the auto update the compat to latest non-breaking version by semantic version definition
(latest minor version for `0.x.y`, latest major version for other case). One can select only updating
specific packages by using the `package` option.

# Arguments

- `version_spec`: version_spec to be compatible with, can be `patch`, `minor`, `major` or `auto`.
- `path_to_project`: path to the project, default is the current working directory.

# Options

- `-p,--package=<name>`: package selector, use this option to update specified package compat only.

# Flags

- `--overwrite`: force update all compat by overwrite old compat.
"""
@cast function compat(version_spec::String="auto", path_to_project::String=pwd(); package::String="", overwrite::Bool=false)
    toml = Base.current_project(path_to_project)
    toml === nothing && cmd_error("cannot find (Julia)Project.toml in $path_to_project")
    d = TOML.parsefile(toml)
    compat = get!(Dict{String, String}, d, "compat")
    env = Pkg.Types.EnvCache(toml)
    for (name, uuid) in env.project.deps
        isempty(package) || package == name || continue
        env.manifest[uuid].version isa VersionNumber || continue # skip stdlib
        update_compat!(compat, name, env.manifest[uuid].version, version_spec, overwrite)
    end
    update_compat!(compat, "julia", VERSION, version_spec, overwrite)
    write_project_toml(toml, d)
    return
end

function update_compat!(compat::Dict{String, String}, name::String, version::String, version_spec::String, overwrite::Bool)
    compat_spec = compat_version(version, version_spec)
    compat[name] = if overwrite
        compat_spec
    else
        append_compat_spec(compat, name, compat_spec)
    end
    return compat
end

function append_compat_spec(compat::Dict{String, String}, name::String, spec::String)
    haskey(compat, name) || return spec
    # TODO: use more compact format, e.g x.y.z-a.b.c
    return string(compat[name], ", ", spec)
end

function compat_version(version::VersionNumber, version_spec::String)
    return if version_spec == "patch"
        string(version.major, ".", version.minor, ".", version.patch)
    elseif version_spec == "minor"
        string(version.major, ".", version.minor)
    elseif version_spec == "major"
        string(version.major)
    elseif version_spec == "auto"
        if version.major == 0
            string(version.major, ".", version.minor)
        else
            string(version.major)
        end
    else
        cmd_error("invalid version spec: $version_spec")
    end
end


text:
__precompile__()

module EmACQR

include("ACQR.jl")

export
	byrow, dare, dlyap,
	ACQR_kfilter, ACQR_kfilter_s,
	ACQR_ksmoother, ACQR_ksmoother_s,
	ACQR_em, ACQR_em_s,
	ACQR_em_starting,
	ACQRS_sub,
	simula01, ACQR_simula,
	ACQR_em_run, ACQR_em_s_run

end # module


text:
using Test

function AND(x1, x2)
    x = [x1, x2]
    w = [0.5, 0.5]
    b = -0.7
    tmp = sum(w.*x) + b
    return tmp > 0
end

function NAND(x1, x2)
    x = [x1, x2]
    w = [-0.5, -0.5]
    b = 0.7
    tmp = sum(w.*x) + b
    return tmp > 0
end

function OR(x1, x2)
    x = [x1, x2]
    w = [0.5, 0.5]
    b = -0.2
    tmp = sum(w.*x) + b
    return tmp > 0
end

function XOR(x1, x2)
    s1 = NAND(x1, x2)
    s2 = OR(x1, x2)
    y = AND(s1, s2)
    return y
end

function main()
    @testset "Test" begin
        @testset "AND" begin
            @test AND(0, 0) == 0
            @test AND(0, 1) == 0
            @test AND(1, 0) == 0
            @test AND(1, 1) == 1
        end
        @testset "NAND" begin
            @test NAND(0, 0) == 1
            @test NAND(0, 1) == 1
            @test NAND(1, 0) == 1
            @test NAND(1, 1) == 0
        end
        @testset "OR" begin
            @test OR(0, 0) == 0
            @test OR(0, 1) == 1
            @test OR(1, 0) == 1
            @test OR(1, 1) == 1
        end
        @testset "XOR" begin
            @test XOR(0, 0) == 0
            @test XOR(0, 1) == 1
            @test XOR(1, 0) == 1
            @test XOR(1, 1) == 0
        end
    end
end

main()


text:
# Code derived from metabolicEP (https://github.com/anna-pa-m/Metabolic-EP)

# TODO remove Y unused parameter
# TODO remove alpha Inf code
function prepareinput(S, lb, ub, alpha, verbose, solution, expval)

    M,N = size(S)
    verbose && M >= N && @warn("M = $M ≥ N = $N")
    all(lb .<= ub) || error("lower bound fluxes > upper bound fluxes. Consider swapping lower and upper bounds")

    verbose && println(stderr, "Analyzing a $M × $N stoichiometric matrix.")

    updatefunction = alpha == Inf ? eponesweepT0! : eponesweep!

    scalefact = max(maximum(abs.(lb)), maximum(abs.(ub)))
    epfields = isnothing(solution) ? epfields = EPFields(N, expval, eltype(S)) :
        epfields = deepcopy(solution.sol) # preserve the original solution!

    return updatefunction, scalefact, epfields
end


text:
using Random
using StatsBase
using Distributions
using StatsPlots
using DataFrames
using ShiftedArrays

Random.seed!(123)

arrival_rates = (1, 1)
# job_size_dist = DiscreteNonParametric([1], [1])
job_size = 1
private_service_rate(job_size) = 1 / (1.5 * job_size)
public_service_rate(job_size) = 1 / (0.7 * job_size)
p1_prob_public = 0.3
p2_prob_public = 0.3

p1_arrival_dist = Exponential(1 / arrival_rates[1])
p2_arrival_dist = Exponential(1 / arrival_rates[2])

private_service_dist = Exponential(1 / private_service_rate(job_size))
public_service_dist = Exponential(1 / public_service_rate(job_size))

p1_njobs = 5
p2_njobs = 5
total_njobs = p1_njobs + p2_njobs
p1_interarrival_times = rand(p1_arrival_dist, p1_njobs)
p2_interarrival_times = rand(p2_arrival_dist, p2_njobs)

p1_job_destinations = rand(Bernoulli(p1_prob_public), p1_njobs)
p2_job_destinations = rand(Bernoulli(p2_prob_public), p2_njobs)

p1_job_order = sort(DataFrame(dests = p1_job_destinations, order = 1:p1_njobs), :dests)
p2_job_order = sort(DataFrame(dests = p2_job_destinations, order = 1:p2_njobs), :dests)

function event_df(arrival_times, service_distribution, job_destinations, player, public, job_order)
    job_destinations = public ? job_destinations : map(!, job_destinations)
    filtered_job_order = filter(r -> r.dests == public, job_order)
    DataFrame(
        arrival_time = arrival_times[job_destinations],
        service_time = rand(service_distribution,
        sum(job_destinations)),
        player = player,
        public = public,
        order = filtered_job_order.order)
end

p1_public_arrivals_df = event_df(p1_interarrival_times, public_service_dist, p1_job_destinations, 1, true, p1_job_order)
p1_private_arrivals_df = event_df(p1_interarrival_times, private_service_dist, p1_job_destinations, 1, false, p1_job_order)
p2_public_arrivals_df = event_df(p2_interarrival_times, public_service_dist, p2_job_destinations, 2, true, p2_job_order)
p2_private_arrivals_df = event_df(p2_interarrival_times, private_service_dist, p2_job_destinations, 2, false, p2_job_order)

function calculate_events(public_arrivals, private_arrivals)
    function service_finishes(arrivals_df)
        sorted = sort(arrivals_df, :order)
        sorted[!, :cum_arrival_time] = cumsum(sorted.arrival_time)
        sorted[!, :service_finishes] = sorted[!, :cum_arrival_time] .+ sorted[!, :service_time]
        return sorted
    end

    arrivals_df = service_finishes(vcat(public_arrivals, private_arrivals))

    function f1(df, column, event_result, length_df = size(df)[1])
        times_df = select(df, [:player, :public, column])
        rename!(times_df, column => :time)
        times_df[!, :event] = repeat([event_result], length_df)
        return times_df
    end

    arrivaltimes_df = f1(arrivals_df, :cum_arrival_time, 1)
    servicetimes_df = f1(arrivals_df, :service_finishes, -1)

    return vcat(arrivaltimes_df, servicetimes_df)
end

p1_eventtimes_df = calculate_events(p1_public_arrivals_df, p1_private_arrivals_df)
p2_eventtimes_df = calculate_events(p2_public_arrivals_df, p2_private_arrivals_df)
eventtimes_df = sort(vcat(p1_eventtimes_df, p2_eventtimes_df), :time)
eventtimes_df[!, :public] = ifelse.(eventtimes_df[!, :public] .== true, "public", "private")
ped = select(deepcopy(eventtimes_df), [:time, :event])
ped[:, :event] = cumsum(ped[:, :event])
sort!(push!(ped, (0, 0)), :time)
ped[:, :endtime] = lag(ped[:, :time], -1)
ped[:, :nextevent] = lag(ped[:, :event], -1)

p = plot()
@inbounds for i = 1:(size(ped, 1) - 1)
    # plot!(p, [ped[i, :time]], [ped[i, :event]], st = :scatter, colour = :black, label = false)
    plot!(p, [ped[i, :time], ped[i, :endtime]], [ped[i, :event], ped[i, :event]], colour = :black, label = false)
    plot!(p, [ped[i, :endtime], ped[i, :endtime]], [ped[i, :event], ped[i, :nextevent]], colour = :black, label = false)
end
# plot!(p, [ped[size(ped, 1), :time]], [ped[size(ped, 1), :event]], st = :scatter, colour = :black, label = false)

ped
p

# public_njobs = sum(p1_job_destinations + p2_job_destinations)
# private_njobs = total_njobs - public_njobs
# public_service_times = rand(public_service_dist, public_njobs)
# private_service_times = rand(private_service_dist, private_njobs)

# public_service_times
# p1_interarrival_times[p1_job_destinations]
# p2_interarrival_times[p2_job_destinations]

# p1_arrivals_df = DataFrame(times = p1_interarrival_times, cumtimes = cumsum(p1_interarrival_times), source = 1, public_dest = p1_job_destinations)
# p2_arrivals_df = DataFrame(times = p2_interarrival_times, cumtimes = cumsum(p2_interarrival_times), source = 2, public_dest = p2_job_destinations)
# arrivals_df = vcat(p1_arrivals_df, p2_arrivals_df)

# public_services_df = DataFrame(service_time = public_service_times, public = true)
# private_services_df = DataFrame(service_time = private_service_times, public = false)
# services_df = vcat(public_services_df, private_services_df)
# hcat(arrivals_df, services_df)

# @df arrivals_df scatter(
#     :times,
#     1:total_njobs,
#     color = :source
# )

# plot(sort(vcat(cumsum(p1_interarrival_times), cumsum(p2_interarrival_times))), 1:20, seriestype = :scatter)


text:
module Day06

import ..data_dir # from parent module
input = read(joinpath(data_dir, "day06"), String)

export part1, part2

"""
Converts a string like "0,1,1,1,2,2" into a vector
[1
 3
 2]
Note the change in 0-based indexing to 1-based indexing.
"""
function collect_days_to_spawn(input::String)
    days_to_spawn = [parse(Int8, string(day)) for day in split(input, ",")]
    # the below line converts the 0-based indexing of the input to 1-based indexing
    [count(x -> x == (day - 1), days_to_spawn) for day = 1:size(age_matrix)[1]]
end

const age_matrix = [
    0 1 0 0 0 0 0 0 0
    0 0 1 0 0 0 0 0 0
    0 0 0 1 0 0 0 0 0
    0 0 0 0 1 0 0 0 0
    0 0 0 0 0 1 0 0 0
    0 0 0 0 0 0 1 0 0
    1 0 0 0 0 0 0 1 0
    0 0 0 0 0 0 0 0 1
    1 0 0 0 0 0 0 0 0]

function age(days_to_spawn::Vector{Int64}, days::Int64 = 1)
    age_matrix^days * days_to_spawn
end

function count_after_aging(days_to_spawn::Vector{Int64}, days::Int64)
    aged = age(days_to_spawn, days)
    sum(aged)
end

function part1(input = input)
    days_to_spawn = collect_days_to_spawn(input)
    count_after_aging(days_to_spawn, 80)
end

function part2(input = input)
    days_to_spawn = collect_days_to_spawn(input)
    count_after_aging(days_to_spawn, 256)
end

end # module


text:
using UNet
using Test

@info "Testing UNet architecture runs..."
@testset "UNet.jl" begin
    X = rand(256, 256, 1, 1)
    y = unet(X)
    @test size(y) == (68, 68, 1, 1)
end