text:
ENV["GKSwstype"] = "100"

using Printf
using Statistics
using Random
using JLD2

using Oceananigans
using Oceananigans.Units
using Oceananigans: fields
using Oceananigans.TurbulenceClosures
using Oceananigans.TurbulenceClosures: FluxTapering

filename = "zonally_averaged_baroclinic_adjustment_withGM"

# Architecture
architecture = CPU()

# Domain
Ly = 1000kilometers  # north-south extent [m]
Lz = 1kilometers     # depth [m]

Ny = 128
Nz = 40

save_fields_interval = 0.5day
stop_time = 60days
Δt₀ = 5minutes

# We choose a regular grid though because of numerical issues that yet need to be resolved
grid = RectilinearGrid(architecture;
                       topology = (Flat, Bounded, Bounded),
                       size = (Ny, Nz),
                       y = (-Ly/2, Ly/2),
                       z = (-Lz, 0),
                       halo = (3, 3))

coriolis = BetaPlane(latitude = -45)

Δy, Δz = Ly/Ny, Lz/Nz

𝒜 = Δz/Δy   # Grid cell aspect ratio.

κh = 0.1    # [m² s⁻¹] horizontal diffusivity
νh = 0.1    # [m² s⁻¹] horizontal viscosity
κz = 𝒜 * κh # [m² s⁻¹] vertical diffusivity
νz = 𝒜 * νh # [m² s⁻¹] vertical viscosity

vertical_closure = VerticalScalarDiffusivity(ν = νz, κ = κz)

horizontal_closure = HorizontalScalarDiffusivity(ν = νh, κ = κh)

diffusive_closures = (vertical_closure, horizontal_closure)

convective_adjustment = ConvectiveAdjustmentVerticalDiffusivity(convective_κz = 1.0,
                                                                convective_νz = 0.0)

gerdes_koberle_willebrand_tapering = FluxTapering(1e-2)

gent_mcwilliams_diffusivity = IsopycnalSkewSymmetricDiffusivity(κ_skew = 1000,
                                                                κ_symmetric = 900,
                                                                slope_limiter = gerdes_koberle_willebrand_tapering)
#####
##### Model building
#####

@info "Building a model..."

closures = (diffusive_closures..., convective_adjustment, gent_mcwilliams_diffusivity)

model = HydrostaticFreeSurfaceModel(grid = grid,
                                    coriolis = coriolis,
                                    buoyancy = BuoyancyTracer(),
                                    closure = closures,
                                    tracers = (:b, :c),
                                    momentum_advection = WENO5(),
                                    tracer_advection = WENO5(),
                                    free_surface = ImplicitFreeSurface())

@info "Built $model."

#####
##### Initial conditions
#####

"""
Linear ramp from 0 to 1 between -Δy/2 and +Δy/2.

For example:

y < y₀           => ramp = 0
y₀ < y < y₀ + Δy => ramp = y / Δy
y > y₀ + Δy      => ramp = 1
"""
ramp(y, Δy) = min(max(0, y/Δy + 1/2), 1)

# Parameters
N² = 4e-6 # [s⁻²] buoyancy frequency / stratification
M² = 8e-8 # [s⁻²] horizontal buoyancy gradient

Δy = 100kilometers
Δz = 100

Δc = 2Δy
Δb = Δy * M²
ϵb = 1e-2 * Δb # noise amplitude

bᵢ(x, y, z) = N² * z + Δb * ramp(y, Δy)
cᵢ(x, y, z) = exp(-y^2 / 2Δc^2) * exp(-(z + Lz/4)^2 / 2Δz^2)

set!(model, b=bᵢ, c=cᵢ)

#####
##### Simulation building
#####

simulation = Simulation(model, Δt=Δt₀, stop_time=stop_time)

# add timestep wizard callback
wizard = TimeStepWizard(cfl=0.1, max_change=1.1, max_Δt=20minutes)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(20))

# add progress callback
wall_clock = [time_ns()]

function print_progress(sim)
    @printf("[%05.2f%%] i: %d, t: %s, wall time: %s, max(u): (%6.3e, %6.3e, %6.3e) m/s, next Δt: %s\n",
            100 * (sim.model.clock.time / sim.stop_time),
            sim.model.clock.iteration,
            prettytime(sim.model.clock.time),
            prettytime(1e-9 * (time_ns() - wall_clock[1])),
            maximum(abs, sim.model.velocities.u),
            maximum(abs, sim.model.velocities.v),
            maximum(abs, sim.model.velocities.w),
            prettytime(sim.Δt))

    wall_clock[1] = time_ns()

    return nothing
end

simulation.callbacks[:print_progress] = Callback(print_progress, IterationInterval(20))

#####
##### Output
#####

Redi_diffusivity = IsopycnalSkewSymmetricDiffusivity(κ_skew = (b=0, c=0),
                                                     κ_symmetric = (b=1, c=0),
                                                     slope_limiter = gerdes_koberle_willebrand_tapering)

dependencies = (Redi_diffusivity,
                model.tracers.b,
                Val(1),
                model.clock,
                model.diffusivity_fields,
                model.tracers,
                model.buoyancy,
                model.velocities)

using Oceananigans.TurbulenceClosures: ∇_dot_qᶜ

∇_q_op = KernelFunctionOperation{Center, Center, Center}(∇_dot_qᶜ,
                                                         grid,
                                                         computed_dependencies = dependencies)

# R(b) eg the Redi operator applied to buoyancy
Rb = Field(∇_q_op)

outputs = merge(fields(model), (; Rb))

simulation.output_writers[:fields] = JLD2OutputWriter(model, outputs,
                                                      schedule = TimeInterval(save_fields_interval),
                                                      field_slicer = nothing,
                                                      prefix = filename * "_fields",
                                                      force = true)

@info "Running the simulation..."

run!(simulation, pickup=false)

@info "Simulation completed in " * prettytime(simulation.run_wall_time)

#=

#####
##### Visualize
#####

using GLMakie

fig = Figure(resolution = (1400, 700))

filepath = filename * "_fields.jld2"

ut = FieldTimeSeries(filepath, "u")
bt = FieldTimeSeries(filepath, "b")
ct = FieldTimeSeries(filepath, "c")
rt = FieldTimeSeries(filepath, "Rb")

# Build coordinates, rescaling the vertical coordinate
x, y, z = nodes((Center, Center, Center), grid)

zscale = 1
z = z .* zscale

#####
##### Plot buoyancy...
#####

times = bt.times
Nt = length(times)

un(n) = interior(ut[n])[1, :, :]
bn(n) = interior(bt[n])[1, :, :]
cn(n) = interior(ct[n])[1, :, :]
rn(n) = interior(rt[n])[1, :, :]

@show min_c = 0
@show max_c = 1
@show max_u = maximum(abs, un(Nt))
min_u = - max_u

@show max_r = maximum(abs, rn(Nt))
@show min_r = - max_r

n = Node(1)
u = @lift un($n)
b = @lift bn($n)
c = @lift cn($n)
r = @lift rn($n)

ax = Axis(fig[1, 1], title="Zonal velocity")
hm = heatmap!(ax, y * 1e-3, z * 1e-3, u, colorrange=(min_u, max_u), colormap=:balance)
contour!(ax, y * 1e-3, z * 1e-3, b, levels = 25, color=:black, linewidth=2)
cb = Colorbar(fig[1, 2], hm)

ax = Axis(fig[2, 1], title="Tracer concentration")
hm = heatmap!(ax, y * 1e-3, z * 1e-3, c, colorrange=(0, 0.5), colormap=:speed)
contour!(ax, y * 1e-3, z * 1e-3, b, levels = 25, color=:black, linewidth=2)
cb = Colorbar(fig[2, 2], hm)

ax = Axis(fig[3, 1], title="R(b)")
hm = heatmap!(ax, y * 1e-3, z * 1e-3, r, colorrange=(min_r, max_r), colormap=:balance)
contour!(ax, y * 1e-3, z * 1e-3, b, levels = 25, color=:black, linewidth=2)
cb = Colorbar(fig[3, 2], hm)

title_str = @lift "Baroclinic adjustment with GM at t = " * prettytime(times[$n])
ax_t = fig[0, :] = Label(fig, title_str)

display(fig)

record(fig, filename * ".mp4", 1:Nt, framerate=8) do i
    @info "Plotting frame $i of $Nt"
    n[] = i
end

=#


text:
#!/usr/bin/env julia

using SplineHazard
using Test

t = Tuner(0)
@test SplineHazard.get_tune_param(t) == 2.38

@test (t.cur_blk == 1) & (t.cur_iter == 1)
update!(t, true)
@test (t.cur_blk == 1) & (t.cur_iter == 2)

t = Tuner(1000)
@test SplineHazard.get_tune_param(t) == 2.38

ext: jl


text:
"""
Déclare les nouveaux types utilisés dans ce projet.
"""

import Base.show

# objet Station modélisant une station de navette
struct Station
    id::Int
    name::String
    latitude::Float64
    longitude::Float64
end

Base.show(io::IO, x::Station) = print(io, x.name)

# objet représentant une requête utilisateur
struct Request
    t::Int
    departureStation::String
    arrivalStation::String
end

# objet représentant un élément du planning
struct Elt
    id::Int
    state::Etat
    idReq::Int
    isSource::Bool
    sourceRemaining::Bool
    limit::Float64 # si elt est un s alors limit dépend de w et r, si elt est un e alors limit dépend de s et epsilon et e (1+epsilon)d(s,e)
    name::String
    cMin::Float64
end

# objet représentant un noeud de l'arbre des plannings possibles (pour le branch and bound)
mutable struct Node
    empty::Bool # si le noeud est supprimé
    probed::Bool # si le noeud est sondé (que tous ses enfants sont sondés ou qu'il ne peut plus posséder d'enfants)
    passedElts::Vector{Elt} # éléments passés triés dans l'ordre de parcours
    passedDT::Vector{Float64} # dT passés correspondant aux temps des éléments passés
    futureElts::Vector{Elt} # éléments futures non triés
    sumFutureCMin::Float64 # somme des CMin des éléments futures
    e::Elt # élément du noeud
    h::Int # hauteur du noeud
    dT::Float64 # distance du noeud à la racine
    lowerBound::Float64 # plus petite distance trouvée en développant ce noeud
    children::Vector{Node} # noeuds enfants de ce noeud
end

# retourne un nouvel élément vide (élément supprimé par exemple)
function emptyElt()
    return Elt(0, Non, 0, true,false, 0., "", 0.)
end

# retourne un nouveau noeud vide (noeud supprimé car trop mauvais par exemple)
function emptyNode()
    passedElts = Vector{Elt}(undef,0)
    children= Vector{Node}(undef,0)
    return Node(true,true,passedElts,Vector{Float64}(undef,0),passedElts,0.,emptyElt(),0,0.,0.,children)
end

# retourne un nouveau noeud
function newNode(passedElts::Vector{Elt},passedDT::Vector{Float64},futureElts::Vector{Elt},sumFutureCMin::Float64,e::Elt,h::Int,dT::Float64,lowerBound::Float64)
    children = Vector{Node}(undef,0)
    return Node(false,false,passedElts,passedDT,futureElts,sumFutureCMin,e,h,dT,lowerBound,children)
end


text:
export
    LaneFollowingDriver,
    StaticLaneFollowingDriver,
    PrincetonDriver,
    IntelligentDriverModel,
    ProportionalSpeedTracker,
    track_longitudinal!

include("lane_following_drivers.jl")
include("static_lane_following_drivers.jl")
include("princeton_drivers.jl")
include("speed_trackers.jl")
include("intelligent_driver_models.jl")


text:
"""

Module `MultivariateOnlineStatistics` deals with statistics of multi-variate
data.  Computation are performed *on-line* that is in one pass, as the data is
obtained.

"""
module MultivariateOnlineStatistics

export
    IndependentStatistics,
    mean,
    nobs,
    order,
    std,
    var

using Base: @propagate_inbounds, OneTo

using Statistics
import Statistics: mean, std, var

using StatsBase
import StatsBase: nobs

"""
    IndependentStatistics{L,T,N}(dims)

yields an object to store the `L` first statistical moments of independent
`N`-dimensional data of dimensions `dims`.  Type parameter `T` is the
floating-point type for computed statistics.  The number `N` of dimensions may
be omitted as it is equal to `length(dims)`.

For fine tuning the type of storage used by the object, the arrays `s1`, `s2`,
..., and `sL` storing the statistics may be provided as an `L`-tuple:

    IndependentStatistics((s1, s2, ..., sL))

where `s1`, `s2`, ..., and `sL` are arrays having the same floating-point
element type and the same indices.  The storage arrays will be zero-filled.

If storage arrays already contain statistical moments collected from a number
of independent samples, then call the constructor with the number `n` of
samples specified as the last argument:

    IndependentStatistics((s1, s2, ..., sL), n)

In that case, input arrays `s1`, `s2`, ..., and `sL` must have been correctly
intialized, typically as follows (∀i):

    s1[i] = (x_1[i] + ... + x_n[i])/n
    s2[i] = (x_1[i] - s1[i])^2 + ... + (x_n[i] - s1[i])^2
    ...
    sL[i] = (x_1[i] - s1[i])^L + ... + (x_n[i] - s1[i])^L

with `x_j` the `j`-th data sample and where index `i` may be multi-dimensional.
That is to say that `s1` is the element-wise empirical mean of the samples
while `s2`, ..., and `sL` are the sum over the samples of the element-wise
sample difference with their element-wise empirical mean raised to the
corresponding power.

Assuming `A` is an instance of `IndependentStatistics`, then collecting
statistics from more samples is done by:

    push!(A, x...) -> A
    merge!(A, itr) -> A

where each `x...` is a single data sample, an (abstract) array of suitable
dimension, while `itr` is an iterable object which yields independent data
samples.  The recurrence formula of Welford (1962) is used to avoid loss of
precision due to rounding errors.

!!! note
    It is assumed that data samples are mutually independent.

If `B` is another instance of `IndependentStatistics`, the statistics collected
by `B` can be merged into `A` by:

    merge!(A, B) -> A

In this case, the recurrence formula of Chan, Golub and LeVeque (1979) is used
to avoid loss of precision due to rounding errors.

To retrieve statistics, a number of methods from the `Statistics` and
`StatsBase` package are re-exported:

    nobs(A)                # the number of independent samples
    mean(A)                # the element-wise sample mean
    var(A; corrected=true) # the element-wise sample variance
    std(A; corrected=true) # the element-wise sample standard deviation

If keyword `corrected` is true (the default) then an unbiased estimator
is returned; otherwise, the maximum-likelihood estimator is returned.

The following basic methods are also applicable to an instance of
`IndependentStatistics`:

    ndims(A)   # the number of dimensions of a data sample
    size(A)    # the dimensions of a data sample
    size(A, k) # the k-th dimension of a data sample
    axes(A)    # the axes of a data sample
    axes(A, k) # the k-th axis of a data sample
    eltype(A)  # the floating-point type of the collected statistics
    order(A)   # the maximum order of statistical moments

"""
mutable struct IndependentStatistics{L,T<:AbstractFloat,N,
                                     A<:AbstractArray{T,N}}
    s::NTuple{L,A} # collected statistical moments
    n::Int         # number of samples
    function IndependentStatistics{L,T,N,A}(
        s::NTuple{L,A},
        n::Integer) where {T<:AbstractFloat,N,L,
                           A<:AbstractArray{T,N}}
        have_same_axes(s...) ||
            dimension_mismatch("storage arrays have different indices")
        n ≥ 0 || throw(ArgumentError("bad number of samples"))
        return new{L,T,N,A}(s, n)
    end
end

# When number of samples is not provided, assume 0 samples and zero-fill after
# building with the inner constructor (to check correctness of arguments).
function IndependentStatistics{L,T,N,A}(
    s::NTuple{L,A}) where {T<:AbstractFloat,N,L,A<:AbstractArray{T,N}}
    empty!(IndependentStatistics{L,T,N,A}(s, 0))
end

# For basic constructors, just provide the missing type parameters and call one
# of the fully qualified constructors.
for f in (:(IndependentStatistics{L,T,N}),
          :(IndependentStatistics{L,T}),
          :(IndependentStatistics{L}),
          :(IndependentStatistics))
    @eval begin
        function $f(s::NTuple{L,A}) where {T<:AbstractFloat,N,L,
                                           A<:AbstractArray{T,N}}
            IndependentStatistics{L,T,N,A}(s)
        end
        function $f(s::NTuple{L,A}, n::Integer) where {T<:AbstractFloat,N,L,
                                                       A<:AbstractArray{T,N}}
            IndependentStatistics{L,T,N,A}(s, n)
        end
    end
end

IndependentStatistics{L,T}(dims::Integer...) where {L,T} =
    IndependentStatistics{L,T}(dims)
IndependentStatistics{L,T}(dims::NTuple{N,Integer}) where {L,T,N} =
    IndependentStatistics{L,T}(Dims(dims))
IndependentStatistics{L,T}(dims::Dims{N}) where {L,T,N} =
    IndependentStatistics(ntuple(x -> zeros(T, dims), L), 0)

IndependentStatistics{L,T,N}(dims::Integer...) where {L,T,N} =
    IndependentStatistics{L,T,N}(dims)
IndependentStatistics{L,T,N}(dims::NTuple{N,Integer}) where {L,T,N} =
    IndependentStatistics{L,T}(Dims(dims))
IndependentStatistics{L,T,N}(dims::NTuple{Np,Integer}) where {L,T,N,Np} =
    throw(ArgumentError("number of dimensions is not equal to type parameter"))

"""
    MultivariateOnlineStatistics.storage(A)

yields a tuple of the arrays storing the statistics collected by `A`.

---
    MultivariateOnlineStatistics.storage(A, k)

yields the `k`-th array storing the statistics collected by `A`.

"""
storage(A::IndependentStatistics) = A.s
storage(A::IndependentStatistics, k) = A.s[k]

"""
    order(A)

yields the highest order of the statistical moments collected by `A`.

"""
order(A::IndependentStatistics) = order(typeof(A))
order(::Type{<:IndependentStatistics{L,T,N}}) where {L,T,N} = L

# Extend simple Base methods.
Base.eltype(A::IndependentStatistics) = eltype(typeof(A))
Base.ndims(A::IndependentStatistics) = ndims(typeof(A))
Base.eltype(::Type{<:IndependentStatistics{L,T}}) where {L,T} = T
Base.ndims(::Type{<:IndependentStatistics{L,T,N}}) where {L,T,N} = N
Base.size(A::IndependentStatistics) = size(storage(A, 1))
Base.size(A::IndependentStatistics, k) = size(storage(A, 1), k)
Base.axes(A::IndependentStatistics) = axes(storage(A, 1))
Base.axes(A::IndependentStatistics, k) = axes(storage(A, 1), k)

"""
     MultivariateOnlineStatistics.checked_axes(A) -> I

yields the indices of the data integrated by `A`.  Compared to `axes(A)`, this
method also checks the number of samples and the indices of all the storage
arrays used by `A`.

"""
@inline function checked_axes(A::IndependentStatistics{L}) where {L}
    A.n ≥ 0 || throw_bad_number_of_samples()
    I = axes(storage(A, 1))
    L < 2 || all_match(I, axes, storage(A, 2:L)...) ||
        throw_bad_storage_indices()
    return I
end

function Base.empty!(A::IndependentStatistics{L,T}) where {L,T}
    for s in storage(A)
        fill!(s, zero(eltype(s)))
    end
    A.n = 0
    return A
end

# Extend StatsBase methods.
nobs(A::IndependentStatistics) = A.n

# Extend Statistics methods.
mean(A::IndependentStatistics) = storage(A, 1)
var(A::IndependentStatistics{L,T,N}; kwds...) where {L,T,N} =
    var!(Array{T,N}(undef, size(A)), A; kwds...)
std(A::IndependentStatistics{L,T,N}; kwds...) where {L,T,N} =
    std!(Array{T,N}(undef, size(A)), A; kwds...)

# Idem but with indices.  Index argument must be a bit specific to avoid
# calling statistics method with an iterator.  NOTE: `getindex` and
# `checkbounds` methods may be called with a linear index (an `Integer`), a
# `CartesianIndex`, or a tuple of Cartesian indices.
@inline @propagate_inbounds function mean(A::IndependentStatistics,
                                          i::Union{Integer,CartesianIndex})
    return getindex(mean(A), i)
end

@inline @propagate_inbounds function mean(A::IndependentStatistics,
                                          i::Tuple{Vararg{Integer}})
    return getindex(mean(A), i...)
end

for (f, g) in ((:var, :identity),
               (:std, :sqrt))
    @eval begin
        @inline @propagate_inbounds function $f(A::IndependentStatistics,
                                                i::Union{Integer,
                                                         CartesianIndex};
                                                kwds...)
            return _mapvar($g, A, i; kwds...)
        end
        @inline @propagate_inbounds function $f(A::IndependentStatistics,
                                                i::Tuple{Vararg{Integer}};
                                                kwds...)
            return _mapvar($g, A, i...; kwds...)
        end
    end
end

function var!(dst::AbstractArray{<:AbstractFloat,N},
              A::IndependentStatistics{<:Any,<:AbstractFloat,N};
              kwds...) where {N}
    return _mapvar!(identity, dst, A; kwds...)
end

function std!(dst::AbstractArray{<:AbstractFloat,N},
              A::IndependentStatistics{<:Any,<:AbstractFloat,N};
              kwds...) where {N}
   return _mapvar!(sqrt, dst, A; kwds...)
end

# Map the variance.
function _mapvar!(f, dst::AbstractArray{<:AbstractFloat,N},
                  A::IndependentStatistics{L,T,N};
                  corrected::Bool=true) where {L,T,N}
    L ≥ 2 || error("2nd order statistical moment is not available")
    s2 = storage(A, 2)
    axes(dst) == axes(s2) ||
        dimension_mismatch("destination has incompatible indices")
    n = nobs(A)
    n ≥ 0 || throw_bad_number_of_samples()
    if corrected
        n > 1 || error("not enough samples")
        n -= 1
    end
    if n > 1
        a = one(T)/n
        @inbounds @simd for i in eachindex(dst, s2)
            dst[i] = f(a*s2[i])
        end
    elseif n == 1
        @inbounds @simd for i in eachindex(dst, s2)
            dst[i] = f(s2[i])
        end
    else
        fill!(dst, f(zero(eltype(dst))))
    end
    return dst
end

# This version does not check the number of samples for efficiency reasons.
@inline function _mapvar(f, A::IndependentStatistics{L},
                         i1, inds...; corrected::Bool=true) where {L}
    L ≥ 2 || error("2nd order statistical moment is not available")
    S2 = storage(A, 2)
    @boundscheck checkbounds(S2, i1, inds...)
    @inbounds s2 = getindex(S2, i1, inds...)
    n = nobs(A)
    return f(s2/ifelse(corrected, n - 1, n))
end

function Base.push!(A::IndependentStatistics{L,T,N},
                    x::AbstractArray{<:Real,N}) where {L,T<:AbstractFloat,N}
    axes(x) == checked_axes(A) ||
        dimension_mismatch("data sample has incompatible indices")
    return unsafe_push!(A, x)
end

# Push without checking arguments.
function unsafe_push!(A::IndependentStatistics{L,T,N},
                      x::AbstractArray{<:Real,N}) where {L,T<:AbstractFloat,N}
    n = A.n
    if n == 0
        # No data have yet been collected by A.
        s1 = storage(A, 1)
        @inbounds @simd for i in eachindex(s1, x)
            s1[i] = x[i]
        end
    elseif L == 1
        s1, = storage(A)
        w1 = T(1)/(n + 1)
        @inbounds @simd for i in eachindex(s1, x)
            # Apply the recurrence given in Welford (1962).
            s1[i] += w1*(x[i] - s1[i])
        end
    elseif L == 2
        s1, s2 = storage(A)
        w1 = T(1)/(n + 1)
        wn = T(n)/(n + 1)
        @inbounds @simd for i in eachindex(s1, s2, x)
            # Apply the recurrence given in Welford (1962).  Note that this
            # guarantees that s2[i] ≥ O always hold.
            u = x[i] - s1[i]
            s1[i] += w1*u
            s2[i] += wn*u^2
        end
    else
        error("statistical moments of order higher than 2 not yet implemented")
    end
    A.n += 1
    return A
end

function Base.copyto!(dst::IndependentStatistics{L,<:AbstractFloat,N},
                      src::IndependentStatistics{L,<:AbstractFloat,N}) where {L,N}
    copy!(dst, src)
end

function Base.copy!(dst::IndependentStatistics{L,<:AbstractFloat,N},
                    src::IndependentStatistics{L,<:AbstractFloat,N}) where {L,N}
    checked_axes(dst) == checked_axes(src) ||
        dimension_mismatch("statistics have incompatible indices")
    return unsafe_copy!(dst, src)
end

# Copy without checking arguments.
function unsafe_copy!(dst::IndependentStatistics{L,<:AbstractFloat,N},
                      src::IndependentStatistics{L,<:AbstractFloat,N}) where {L,N}
    for k in 1:L
        copy!(storage(dst, k), storage(src, k))
    end
    dst.n = src.n
    return dst
end

function Base.merge!(A::IndependentStatistics, itr)
    I = checked_axes(A)
    for x in itr
        axes(x) == I ||
            dimension_mismatch("data sample has incompatible indices")
        unsafe_push!(A, x)
    end
    return A
end

function Base.merge!(A::IndependentStatistics{L,<:AbstractFloat,N},
                     B::IndependentStatistics{L,<:AbstractFloat,N}) where {L,N}
    checked_axes(A) == checked_axes(B) ||
        dimension_mismatch("statistics have incompatible indices")
    if A.n == 0
        # No data have been collected in A, just copy B if any data there.
        if B.n > 0
            unsafe_copy!(A, B)
        end
    elseif B.n == 1
        # A single sample has been collected by B, it is sufficient to push it
        # into A.
        unsafe_push!(A, storage(B, 1))
    else
        # Some samples have been integrated in A and B.
        _merge!(A, B)
    end
    return A
end

function _merge!(A::IndependentStatistics{1,T,N},
                 B::IndependentStatistics{1,<:AbstractFloat,N}) where {T,N}
    # Apply the formula given in Chan, Golub and LeVeque (1979) at the
    # numerical precision of the destination to combine the collected
    # statistics.
    As1, = storage(A)
    Bs1, = storage(B)
    n = A.n + B.n
    wa = T(A.n)/n
    wb = T(B.n)/n
    wab = wa*B.n
    @inbounds @simd for i in eachindex(As1, As2, Bs1, Bs2)
        As1[i] = wa*As1[i] + wb*T(Bs1[i])
    end
    A.n = n
    nothing
end

function _merge!(A::IndependentStatistics{2,T,N},
                 B::IndependentStatistics{2,<:AbstractFloat,N}) where {T,N}
    # Apply the formulae given in Chan, Golub and LeVeque (1979) at the
    # numerical precision of the destination to combine the collected
    # statistics.
    As1, As2 = storage(A)
    Bs1, Bs2 = storage(B)
    n = A.n + B.n
    wa = T(A.n)/n
    wb = T(B.n)/n
    wab = wa*B.n
    @inbounds @simd for i in eachindex(As1, As2, Bs1, Bs2)
        a, b = As1[i], T(Bs1[i])
        As1[i] = wa*a + wb*b
        As2[i] += T(Bs2[i]) + wab*(a - b)^2
    end
    A.n = n
    nothing
end

function _merge!(A::IndependentStatistics{L,T,N},
                 B::IndependentStatistics{L,<:AbstractFloat,N}) where {L,T,N}
    error("statistical moments of order higher than 2 not yet implemented")
end

# This method borrowed from ArrayTools.
"""
    MultivariateOnlineStatistics.all_match(val, f, args...) -> bool

yields as soon as possible (short-circuit) whether `f(arg) == val` for each
argument `arg` in `args...`.  The returned value is `true` if there are no
arguments after `f`.

"""
all_match(val, f::Function) = true
all_match(val, f::Function, A) = f(A) == val
@inline all_match(val, f::Function, A, B...) =
    all_match(val, f, A) && all_match(val, f::Function, B...)

"""
    MultivariateOnlineStatistics.same_axes(A...) -> inds

yields the axes of arrays `A...` throwing an error if these arrays have
not the same axes.

"""
same_axes(A::AbstractArray) = axes(A)
@inline same_axes(A::AbstractArray, B::AbstractArray...) = begin
    I = axes(A)
    all_match(I, axes, B...) ||
        dimension_mismatch("arguments have different indices")
    return I
end

"""
    MultivariateOnlineStatistics.have_same_axes(A...) -> bool

yields whether arrays `A...` have the same axes.

"""
have_same_axes(A::AbstractArray) = true
@inline have_same_axes(A::AbstractArray, B::AbstractArray...) =
    all_match(axes(A), axes, B...)

"""
    MultivariateOnlineStatistics.to_int(x)

yields integer `x` converted to an `Int`.

"""
to_int(x::Int) = x
to_int(x::Integer) = Int(x)

@noinline dimension_mismatch(args...) = dimension_mismatch(string(args...))
@noinline dimension_mismatch(msg::AbstractString) =
    throw(DimensionMismatch(msg))

@noinline assertion_error(args...) = assertion_error(string(args...))
@noinline assertion_error(msg::AbstractString) = throw(AssertionError(msg))

# Throw an AssertionError when contents of a structure is inconsistent.
@noinline throw_bad_number_of_samples() =
    assertion_error("bad number of samples")
@noinline throw_bad_storage_indices() =
    assertion_error("storage arrays have different indices")

end # module

max_issues_repo_name: emmt/MultivariateOnlineStatistics

max_forks_repo_name: emmt/MultivariateOnlineStatistics


text:
module AnalyticCartesianField

    import ..ElectromagneticFields
    import ..ElectromagneticFields: AnalyticEquilibrium, AnalyticPerturbation

    export CartesianField, CartesianEquilibrium, CartesianPerturbation

    abstract type CartesianEquilibrium <: AnalyticEquilibrium end
    abstract type CartesianPerturbation <: AnalyticPerturbation end

    CartesianField = Union{CartesianEquilibrium, CartesianPerturbation}

    X(x::AbstractVector, ::CartesianField) = x[1]
    Y(x::AbstractVector, ::CartesianField) = x[2]
    Z(x::AbstractVector, ::CartesianField) = x[3]

    ElectromagneticFields.J(x::AbstractVector, ::CartesianField) = one(eltype(x))

    ElectromagneticFields.x¹(ξ::AbstractVector, ::CartesianField) = ξ[1]
    ElectromagneticFields.x²(ξ::AbstractVector, ::CartesianField) = ξ[2]
    ElectromagneticFields.x³(ξ::AbstractVector, ::CartesianField) = ξ[3]

    ElectromagneticFields.ξ¹(x::AbstractVector, ::CartesianField) = x[1]
    ElectromagneticFields.ξ²(x::AbstractVector, ::CartesianField) = x[2]
    ElectromagneticFields.ξ³(x::AbstractVector, ::CartesianField) = x[3]

    ElectromagneticFields.g₁₁(x::AbstractVector, ::CartesianField) = one(eltype(x))
    ElectromagneticFields.g₂₂(x::AbstractVector, ::CartesianField) = one(eltype(x))
    ElectromagneticFields.g₃₃(x::AbstractVector, ::CartesianField) = one(eltype(x))

end


text:
"""
Abstract type for all components used to compose a [`DynamicInjection`](@ref) device
"""
abstract type DynamicComponent <: DeviceParameter end

"""
Abstract type for all dynamic injection types
"""
abstract type DynamicInjection <: Device end

"""
Return all the dynamic components of a [`DynamicInjection`](@ref) device
"""
function get_dynamic_components(device::T) where {T <: DynamicInjection}
    return (
        getfield(device, x) for
        (x, y) in zip(fieldnames(T), fieldtypes(T)) if y <: DynamicComponent
    )
end

supports_services(::DynamicInjection) = false
get_states(::DynamicInjection) = Vector{Symbol}()
"""
    Default implementation of get_state_types for dynamic components. Assumes all states are
    Differential
"""
function get_states_types(d::DynamicComponent)
    return fill(StateTypes.Differential, get_n_states(d))
end


text:
# vim: set ts=8 sw=4 sts=4 et:

module Unicode2LaTeX

export unicode2latex

# Mapping generated from the LaTeX distribution.
include("u2latex.jl")

unicode2latex(c::Char) = get(u2latex, c, string(c))
function unicode2latex(s::String)
    o = IOBuffer()
    for (i, c) in enumerate(s)
        write(o, unicode2latex(c))
    end
    String(take!(o))
end

function generate(;nowarn=false)
    nowarn || @warn "Do not call unless you are the package maintainer"
    texp       = "base/utf8ienc.dtx"
    (ver, tex) = mktempdir() do d
        cd(d) do
            run(`git clone --quiet https://github.com/latex3/latex2e.git`)
            cd("latex2e")
            ver = chomp(read(`git log --pretty='%h' -n 1 --no-merges $texp`,
                String))
            tex = read(texp, String)
            (ver, tex)
        end
    end
    mktemp() do p, f
        print(f, """
            # This file is programmatically generated, do not edit it!
            #   Generated from `$texp` from the LaTeX 2e source tree.
            #   Available under the LaTeX Project Public License (LPPL).
            #   File commit hash: $ver
            """)
        println(f, "const u2latex = Dict(")
        d = Dict{Char,String}()
        i = 1
        while true
            m = match(
                r"""
                \\DeclareUnicodeCharacter{
                        ([0-9A-Fa-f]+)
                    }{
                        ((?:[^{}]|(?R))+)
                    }
                """x,
                tex, i)
            m !== nothing || break
            i = last(m.offsets)
            c = Char(parse(Int, m[1], base=16))
            # Remove redefinition protection.
            s = replace(m[2], "@tabacckludge" => "")
            haskey(d, c) == ((c => s) in d) || error("Multiple mappings: $c")
            d[c] = s
            println(f, "    $(:($c => $s)),")
        end
        println(f, ")")
        close(f)
        cp(p, "$(dirname(@__FILE__))/u2latex.jl", force=true)
        length(d)
    end
end

end # module


text:
using IndexedGraphs
using Test

include("indexeddigraph.jl")
include("indexedbidigraph.jl")
include("indexedgraph.jl")
include("factorgraph.jl")
include("dijkstra.jl")


text:
#------------------------------------------------------------------------------

# rational number reconstruction implementation borrowed from CLUE
# and modified a bit to suit the 'Modern Computer Algebra' definitions
# returns a rational r // h of QQ field in a canonical form such that
#   r // h ≡ a (mod m)
#
# let n = max( λ(a), λ(m) ) , where λ(x) is a number of bits for x
# O(n^2)
"""
    rational_reconstruction(a, m)

Rational number reconstruction implementation borrowed from CLUE
and modified a bit to suit the 'Modern Computer Algebra' definitions.
Returns a rational of QQ field in a canonical form that
is congruent a modulo m

a, m are integers
"""
function rational_reconstruction(a::I, m::I) where {I<:Union{Int, BigInt}}
    a = mod(a, m)
    if a == 0 || m == 0
        return Nemo.QQ(0, 1)
    end
    if m < 0
        m = -m
    end
    if a < 0
        a = m - a
    end
    if a == 1
        return Nemo.QQ(1, 1)
    end
    bnd = sqrt(float(m) / 2)

    U = I[1, 0, m]
    V = I[0, 1, a]
    while abs(V[3]) >= bnd
        q = div(U[3], V[3])
        T = U .- q .* V
        U = V
        V = T
    end

    t = abs(V[2])
    r = V[3] * sign(V[2])
    # changed from `<= bnd` to `<= m / bnd`
    if t <= m / bnd && gcd(r, t) == 1
        return Nemo.QQ(r, t)
    end

    throw(DomainError(
        :($a//$m), "rational reconstruction of $a (mod $m) does not exist"
    ))
end

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

function modular_reduction(x::FracElem, field)
    n, d = field(numerator(x)), field(denominator(x))
    if iszero(d)
        throw(DomainError(
            :($x), "modular reduction of $x (to $field) does not exist"
        ))
    end
    n // d
end

function modular_reduction(x::gfp_fmpz_elem, field)
    return field(convert(BigInt, x))
end

function modular_reduction(x::gfp_elem, field)
    return field(x)
end

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

# Returns the minimal polynomial of the sequence `h` modulo `m`,
# the sequence `h` is formed as the given polynomial coefficients
# If n is max(deg(h), deg(m))
# O(n^2)
function minimal_polynomial(h::PolyElem, m::PolyElem)
    # The function implements Algorithm 12.9 from
    #   "Modern Computer Algebra", second edition,
    #   Joachim von zur Gathen, Jürgen Gerhard

    n = div(degree(m), 2)
    FF = parent(h)

    U = (FF(1), FF(0), m)
    V = (FF(0), FF(1), h)
    while degree(V[3]) >= n
        q = div(U[3], V[3])
        T = U .- q .* V
        U = V
        V = T
    end

    s, t = V[3], V[2]

    if iszero(coeff(t, 0))
        throw(SingularException(0))
    end

    t *= inv(coeff(t, 0))
    d = max(degree(s) + 1, degree(t))

    if gcd(s, t) == 1 && degree(t) <= n
        return reverse(t, d + 1)
    end

    throw(DomainError(
        h, "minimal polynomial of $h (mod $m) does not exist"
    ))
end

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

# Returns such vector, that it is orthogonal to the given vector set
# The given `vectors` set must be linearly independent
# and for any vect from `vectors` the following must hold
#   length(vectors) = dim(vect) - 1
#
# vectors :: dict of elems (row_idx, row)
#
# Note that vectors may be modified, CAREFUL ABOUT CONCURRENCY!
# It is guaranteed that after the function yields,
# `vectors` will be in the original state
#
# NOT GUARANTEED FOR RETURNED VECTOR TO BE ORTHOGONAL
function find_orthogonal!(vectors::AbstractDict)
    first_vector = first(values(vectors))
    field = base_ring(first_vector)
    n = dim(first_vector)

    i = 1
    found = false
    nnz_rows = Array(1 : n)
    f = zero(first_vector)

    # some complexity estimations should be done
    while ! found
        # taking the i-th unit vector
        u = unit_sparsik(n, i, field)

        # the only change of `vectors`
        vectors[n] = u

        b = unit_sparsik(n, n, field)

        B = from_rows(n, n, field, nnz_rows, vectors)

        try
            # if the selected `u` is not independent with `vectors`
            # then the wiedemann must fail and we should choose a new vector
            # Gleb: we could use probabilistic here
            # Alex: we could, so, the answer could be "almost orthogonal"
            # (only orthogonal to some of the vectors)
            f = square_nonsingular_randomized_wiedemann(B, b)
            found = true
        catch e
            if isa(e, SingularException)
                @info e
            else
                rethrow()
            end
        end

        i += 1
        i % 10 == 0 && error("something is wrong")
    end

    # undo the change
    delete!(vectors, n)

    return f
end

# A slower analogue of the upper function
#
# vectors :: array of sparsiks
function find_orthogonal!(vectors::AbstractVector)
    return find_orthogonal!(
        Dict(i => vect for (i, vect) in enumerate(vectors))
    )
end

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

#=
function construct_jacobians(system)
    nnz = [
        (i, var_index(v), derivative(f, v))
        for (i, f) in enumerate(system)
            for v in vars(f)
    ]
    m, n = length(system), length(degrees(first(system)))
    ring = parent(first(system))
    return from_COO(m, n, nnz, ring)
end
=#

# for the given system of polynomials in variables xi
# consturucts a set of matrices Ai over number field
# such that the Jacobian J of the provided system can represented as the sum
# J = Aᵢxⁱ
"""
    construct_jacobians(system)

For the given system of polynomials in variables xi
consturucts a set of matrices Aᵢ over number field
such that the Jacobian J of the provided system can represented as the sum
J = Aᵢxⁱ
"""
function construct_jacobians(system::AbstractArray{T}) where {T<:MPolyElem}
    domain = base_ring(first(system))
    poly_ring = parent(first(system))
    gen2idx = Dict(x=>i for (i, x) in enumerate(gens(poly_ring)))

    jacobians = Dict()

    for (p_idx, poly) in enumerate(system)
        for var in vars(poly)
            v_idx = gen2idx[var]
            # term is of form α*monomial
            for term in terms(derivative(poly, var))
                monom = monomial(term, 1)
                cf = coeff(term, 1)

                # how can we do it in a normal way..
                !haskey(jacobians, monom) && (jacobians[monom] = Dict())
                idx = (v_idx, p_idx)
                if !haskey(jacobians[monom], idx)
                    jacobians[monom][idx] = zero(domain)
                end
                jacobians[monom][idx] += cf
            end
        end
    end

    m, n = length(system), length(gens(poly_ring))
    factors = [
        from_COO(m, n, jac, domain)
        for jac in values(jacobians)
    ]

    @info "constructed a set of $(length(factors)) matrices $m×$n from the system Jacobian"

    return factors
end

function construct_jacobians(system)
    ring = parent(first(keys(system)))
    construct_jacobians([system[x] for x in gens(ring)])
end

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

# vectors - an array of 1D iterables
# domain - an MPoly domain
#
# converts the given vectors into the polynomial form
# with respect to the generators of the given domain
function polynormalize(vectors, domain) where {T}
    polynômes = []
    vars = Nemo.gens(domain)

    for v in vectors
        push!(
            polynômes,
            sum(map(prod, zip(vars, to_dense(v))))
        )
    end

    polynômes
end

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

# for As[i] in array As compute dok_sparsik Tr: Tr_ij = trace(A_i, A_j)
function gram_matrix(As::Array{AbstractSparseObject{FlintRationalField},1})
    F = base_ring(As[1])
    n = length(As)

    traces = Dict{Tuple{Int, Int}, elem_type(F)}(
        (i, j) => tr(As[i] * As[j])
        for i in 1 : n
            for j in i : n
    )

    nnz_coords = [
        (i, j, traces[min(i, j), max(i, j)])
        for i in 1 : n
            for j in 1 : n
                if ! iszero(traces[min(i, j), max(i, j)])
    ]

    return from_COO(n, n, nnz_coords, F)
end


text:
# A simple model
struct MySimpleModel <: MCMC.Model
  K::Int
end
MCMC.make_init_state(m::MySimpleModel) = (a=0, b=0.0, c=zeros(m.K))
function MCMC.update(::MySimpleModel, s::T) where T
  newstate = (a=s.a + 1,
              b=s.b + randn(),
              c=s.c .+ 1)
  return setproperties!!(s, newstate)
end


text:
using Base.Test

try
    if Pkg.installed("AIMASamples") == nothing
        Pkg.add("AIMASamples")
    end
catch
    Pkg.clone("https://github.com/sambitdash/AIMASamples.jl")
end

include(joinpath(Pkg.dir("AIMASamples"), "test", "runtests.jl"))


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

using Markdown
using InteractiveUtils

# ╔═╡ 59750850-b0ea-11eb-2086-87707f76e327
md"""
# [Day 24: It Hangs in the Balance](https://adventofcode.com/2015/day/24)

It's Christmas Eve, and Santa is loading up the sleigh for this year's deliveries. However, there's one small problem: he can't get the sleigh to balance. If it isn't balanced, he can't defy physics, and nobody gets presents this year.

No pressure.

Santa has provided you a list of the weights of every package he needs to fit on the sleigh. The packages need to be split into three groups of exactly the same weight, and every package has to fit. The first group goes in the passenger compartment of the sleigh, and the second and third go in containers on either side. Only when all three groups weigh exactly the same amount will the sleigh be able to fly. Defying physics has rules, you know!

Of course, that's not the only problem. The first group - the one going in the passenger compartment - needs as few packages as possible so that Santa has some legroom left over. It doesn't matter how many packages are in either of the other two groups, so long as all of the groups weigh the same.

Furthermore, Santa tells you, if there are multiple ways to arrange the packages such that the fewest possible are in the first group, you need to choose the way where the first group has the smallest quantum entanglement to reduce the chance of any "complications". The quantum entanglement of a group of packages is the [product](https://en.wikipedia.org/wiki/Product_%28mathematics%29) of their weights, that is, the value you get when you multiply their weights together. Only consider quantum entanglement if the first group has the fewest possible number of packages in it and all groups weigh the same amount.

For example, suppose you have ten packages with weights `1` through `5` and `7` through `11`. For this situation, some of the unique first groups, their quantum entanglements, and a way to divide the remaining packages are as follows:

```
Group 1;             Group 2; Group 3
11 9       (QE= 99); 10 8 2;  7 5 4 3 1
10 9 1     (QE= 90); 11 7 2;  8 5 4 3
10 8 2     (QE=160); 11 9;    7 5 4 3 1
10 7 3     (QE=210); 11 9;    8 5 4 2 1
10 5 4 1   (QE=200); 11 9;    8 7 3 2
10 5 3 2   (QE=300); 11 9;    8 7 4 1
10 4 3 2 1 (QE=240); 11 9;    8 7 5
9 8 3      (QE=216); 11 7 2;  10 5 4 1
9 7 4      (QE=252); 11 8 1;  10 5 3 2
9 5 4 2    (QE=360); 11 8 1;  10 7 3
8 7 5      (QE=280); 11 9;    10 4 3 2 1
8 5 4 3    (QE=480); 11 9;    10 7 2 1
7 5 4 3 1  (QE=420); 11 9;    10 8 2

```

Of these, although `10 9 1` has the smallest quantum entanglement (`90`), the configuration with only two packages, `11 9`, in the passenger compartment gives Santa the most legroom and wins. In this situation, the quantum entanglement for the ideal configuration is therefore `99`. Had there been two configurations with only two packages in the first group, the one with the smaller quantum entanglement would be chosen.

**What is the quantum entanglement of the first group of packages in the ideal configuration?**
"""

# ╔═╡ 45a666a1-ffaa-4517-b79b-cc6b182505b1
puzzle_input = "1
2
3
7
11
13
17
19
23
31
37
41
43
47
53
59
61
67
71
73
79
83
89
97
101
103
107
109
113"

# ╔═╡ Cell order:
# ╟─59750850-b0ea-11eb-2086-87707f76e327
# ╟─45a666a1-ffaa-4517-b79b-cc6b182505b1


text:
import LightXML: free, parse_string
import SeisIO: ChanSpec,
  KW,
  TimeSpec,
  check_for_gap!,
  dtconst,
  endtime,
  fastread,
  fwrite_note!,
  mk_t,
  mk_xml!,
  parsetimewin,
  read_station_xml!,
  split_id,
  sxml_mergehdr!,
  sμ,
  t_win,
  trunc_x!,
  x_inds,
  μs
import SeisIO.Quake: event_xml!,
  fwrite_note_quake!,
  new_qml!,
  write_qml!


text:
module JungleHelperPlankJumpthru

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Jump Through (Plank) (Jungle Helper)" => Ahorn.EntityPlacement(
        Maple.JumpThru,
        "rectangle",
        Dict{String, Any}(
          "texture" => "JungleHelper/Plank"
        )
    )
)

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
module exactQmodule

using Base, SparseArrays, Arpack, LinearAlgebra
using Combinatorics, KrylovKit, JLD

export AbstractQBasis

#from (1)
export bitSwap, bitCount, findFirstBit, getBitPos, FermionicBitSwap

#from (2)
export QState, shiftStateBasis, getCijMatrix, projectOntoStates, entanglementHalfCut

#from (3)
export QBosonicBasis, addParticleBosonic, constructBosonicState, bosonicHashfunc

#from 4
export fermionicCr, QFermionicBasis, fermionicCrMatrix, addParticleFermionic, constructFermionicState, fermionicHashfunc

#from 5
export loadFermionicCr, saveFermionicCr

abstract type AbstractQBasis end

include("./bitOperations.jl")
include("./Qstate.jl")
include("./bosonic.jl")
include("./fermionic.jl")
include("./preBuild.jl")

end #modlue

max_stars_repo_name: Gawatz/exactQmodule

max_issues_repo_name: Gawatz/exactQmodule

max_forks_repo_name: Gawatz/exactQmodule


text:
import Revise
import GLRenderer
import PoseComposition
import Rotations
import FileIO

R = Rotations
P = PoseComposition
GL = GLRenderer

obj_path = joinpath(@__DIR__, "035_power_drill/textured_simple.obj")
texture_path = joinpath(@__DIR__, "035_power_drill/texture_map.png")

camera_intrinsics = GL.CameraIntrinsics(
    640, 480,
    1000.0, 1000.0,
    320.0, 240.0,
    0.01, 5.0
)
camera_intrinsics = GL.scale_down_camera(camera_intrinsics, 4)

renderer = GL.setup_renderer(camera_intrinsics, GL.DepthMode())

mesh_data = GL.get_mesh_data_from_obj_file(obj_path)
GL.load_object!(renderer, mesh_data)

depth_image = GL.gl_render(renderer, [1], [P.Pose([0.0, 0.0, 1.0], R.RotXYZ(0.1, 0.4, 0.9))], P.IDENTITY_POSE)

times = [
    let
        r = rand() * 2*π
        t = @elapsed GL.gl_render(renderer, [1], [P.Pose([0.0, 0.0, 1.0], R.RotXYZ(0.1, 0.4, r))], P.IDENTITY_POSE)
        t
    end
    for _ in 1:1000
]
avg_time = sum(times)/length(times)
@show avg_time

FileIO.save("imgs/speed_test_img.png", GL.view_depth_image(depth_image))


text:
using YAML
using Logging
using ArgParse

include("src/network.jl")
include("src/plot.jl")

argtable = ArgParseSettings()
@add_arg_table argtable begin
    "--seed"
    arg_type = Int
    default = 0
    "--log"
    arg_type = String
    default = "timings.log"
    "--plot"
    arg_type = String
    default = "timings.pdf"
    "--cfg"
    arg_type = String
    default = "cfg/default.yaml"
    "--vinput"
    arg_type = Int64
    default = 15
    "--connectivity"
    arg_type = Float64
    default = 0.1
end
args = parse_args(argtable)

cfg = YAML.load_file(args["cfg"])
cfg["vinput"] = args["vinput"]

Logging.configure(filename=args["log"], level=INFO)

nn = 100
nin = 10
nout = 10
tsim = 1000

srand(args["seed"])
n = Network(nn, nin, nout, args["connectivity"], cfg)

inputs = rand(nin, tsim) .< (0.5 / nin)
reward = 0

for t in 1:tsim
    input!(n, inputs[:, t])
    s = spike!(n)
    learn!(n, s)
    if inputs[1, t] && inputs[4, t] && inputs[7, t]
        reward = t + rand(1:10)
    end
    if t == reward
        n.da[1] += 1.0
    end
    if ((t > 0) && (mod(t, 10) == 0))
        Logging.info(@sprintf("E: %d %0.4f %0.4f %0.4f %0.4f %0.4f %0.4f", t,
                              mean(n.ge), mean(n.gi), mean(n.neurons[:, 1]),
                              mean(n.neurons[:, 2]), mean(n.neurons[:, 3]),
                              mean(n.weights)))
    end
end

plot_spikes(n, args["plot"])


text:
"""
Heat transfer coefficients.

* [`therm_cond_dryair`](@ref)
* [`heat_cap_dryair`](@ref)
* [`therm_cond_vapor`](@ref)
* [`heat_cap_vapor`](@ref)
* [`therm_cond_moistair`](@ref)
* [`heat_cap_moistair`](@ref)
* [`heat_cap_mass_moistair`](@ref)
* [`therm_diff_moistair`](@ref)
* [`prandtl`](@ref)

# References

* [T08] Tsilingiris, P. T. (2008). Thermophysical and transport properties of
  humid air at temperature range between 0 and 100°C. *Energy Conversion and
  Management*, 49(5), 1098--1110.
  <https://doi.org/10.1016/j.enconman.2007.09.015>
"""
module Heat

include("../docstring_style.jl")

using Canopy.Constants: T_0, M_w, M_d
using Canopy.Water: vapor_mole_frac
using Canopy.Air: air_density
using Canopy.Transfer.Momentum: dyn_visc_moistair, viscosity_interaction_terms

export therm_cond_dryair,
    heat_cap_dryair,
    therm_cond_vapor,
    heat_cap_vapor,
    therm_cond_moistair,
    heat_cap_moistair,
    heat_cap_mass_moistair,
    therm_diff_moistair,
    prandtl

# == dry air ==

"""
Calculate the thermal conductivity of dry air [W m^-1 K^-1] from temperature
[K]. Applicable between 250 and 1050 K.
"""
therm_cond_dryair(temp) = @evalpoly(
    temp, -2.276501e-3, 1.2598485e-4, -1.4815235e-7,
    1.73550646e-10, -1.066657e-13, 2.47663035e-17)

"""
Calculate the isobaric heat capacity of dry air [J mol^-1 K^-1] from
temperature [K]. Applicable between 250 and 1050 K.
"""
heat_cap_dryair(temp) = @evalpoly(
    temp, 1.03409, -0.284887e-3, 0.7816818e-6,
    -0.4970786e-9, 0.1077024e-12) * 1e3 * M_d

# == water vapor ==

"""
Calculate the themal conductivity of water vapor [W m^-1 K^-1] from temperature
[K]. Applicable between 273 and 393 K.
"""
therm_cond_vapor(temp) = @evalpoly(
    temp - T_0,  # <- the empirical eqn takes Celsius values
    1.761758242e-2, 5.558941059e-5, 1.663336663e-7)

"""
Calculate the isobaric heat capacity of water vapor [J mol^-1 K^-1] from
temperature [K]. Applicable between 273 and 393 K.
"""
heat_cap_vapor(temp) = @evalpoly(
    temp - T_0,  # <- the empirical eqn takes Celsius values
    1.86910989, -2.578421578e-4, 1.941058941e-5) * 1e3 * M_w

# == moist air ==

"""
Calculate the thermal conductivity of moist air [W m^-1 K^-1] from temperature
[K], pressure [Pa], and relative humidity [0--1].
"""
function therm_cond_moistair(temp, pressure, RH)
    x_w = vapor_mole_frac(temp, pressure, RH)
    kappa_d = therm_cond_dryair(temp)
    kappa_w = therm_cond_vapor(temp)
    phi_dw, phi_wd, _, _ = viscosity_interaction_terms(temp)
    return (1. - x_w) * kappa_d / ((1. - x_w) + x_w * phi_dw) +
        x_w * kappa_w / (x_w + (1. - x_w) * phi_wd)
end

"""
Calculate the isobaric heat capacity of moist air [J mol^-1 K^-1] from
temperature [K], pressure [Pa], and relative humidity [0--1].
"""
function heat_cap_moistair(temp, pressure, RH)
    cp_d = heat_cap_dryair(temp)
    cp_w = heat_cap_vapor(temp)
    x_w = vapor_mole_frac(temp, pressure, RH)
    return (1. - x_w) * cp_d + x_w * cp_w
end

"""
Calculate the isobaric heat capacity (per unit mass) of moist air [J kg^-1
K^-1] from temperature [K], pressure [Pa], and relative humidity [0--1].
"""
function heat_cap_mass_moistair(temp, pressure, RH)
    cp_d = heat_cap_dryair(temp)
    cp_w = heat_cap_vapor(temp)
    x_w = vapor_mole_frac(temp, pressure, RH)
    return ((1. - x_w) * cp_d + x_w * cp_w) / ((1. - x_w) * M_d + x_w * M_w)
end

"""
Calculate the thermal diffusivity of moist air [m^2 s^-1] from temperature [K],
pressure [Pa], and relative humidity [0--1].
"""
function therm_diff_moistair(temp, pressure, RH)
    rho_a = air_density(temp, pressure, RH)
    # isobaric specific heat capacity (per unit mass) of moist air
    cpm_a = heat_cap_mass_moistair(temp, pressure, RH)  # [J kg^-1 K^-1]
    return therm_cond_moistair(temp, pressure, RH) / (rho_a * cpm_a)
end

"""
Calculate the Prandtl number of moist air from temperature [K], pressure [Pa],
and relative humidity [0--1].
"""
prandtl(temp, pressure, RH) = dyn_visc_moistair(temp, pressure, RH) *
    heat_cap_mass_moistair(temp, pressure, RH) /
    therm_cond_moistair(temp, pressure, RH)

end  # module


text:
export PerspectiveMap3D

function PerspectiveMap3D(kx::Real, ky::Real; rx = 0.00001, ry = 0.0001)
    pmap(p) = (p[1] + kx * p[3] + rx * p[2], p[2] + ky * p[3] + ry * p[1])
    return pmap
end

# rx ry make every lines not horizontal or vertical
function PerspectiveMap3D(array::Array, kx::Real, ky::Real; rx = 0.00001, ry = 0.0001)
    tmap = PerspectiveMap3D(kx, ky; rx = rx, ry = ry)
    N = size(array)[1]
    new_locations = zeros(Float64, N, 2)
    for i = 1:N
        new_locations[i, :] .= tmap(array[i, :])
    end
    return stretch2D!(new_locations)
end

PerspectiveMap3D(lattice::Lattice, kx::Real, ky::Real) =
    typeof(lattice)(lattice.graph, PerspectiveMap3D(lattice.locations, kx, ky))

function stretch2D!(array::Array; xytol = 0.1)
    xmax = maximum(array[:, 1]) + xytol
    xmin = minimum(array[:, 1]) - xytol
    ymax = maximum(array[:, 2]) + xytol
    ymin = minimum(array[:, 2]) - xytol
    xspan, yspan = xmax - xmin, ymax - ymin
    array[:, 1] .= (array[:, 1] .- xmin) / xspan
    array[:, 2] .= (array[:, 2] .- ymin) / yspan
    return array
end

max_stars_repo_name: qiyang-ustc/FTensorKit

max_issues_repo_name: qiyang-ustc/FTensorKit

max_forks_repo_name: qiyang-ustc/FTensorKit


text:
using CompileBot

botconfig = BotConfig(
  "AcuteML";
  yml_path = "SnoopCompile.yml"        # parse `os` and `version` from `SnoopCompile.yml`
)

println("Benchmarking the inference time of `using AcuteML`")
snoop_bench(
  botconfig,
  :(using AcuteML),
)

println("Benchmarking the inference time of `using AcuteML` and runtests")
snoop_bench(
  botconfig,
  "$(@__DIR__)/example_script.jl",
)

max_stars_repo_name: aminya/AML

max_forks_repo_name: aminya/AML


text:
#=
Copyright 2022 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

"""
    setInletBC(t :: Float64, dt :: Float64, v :: Vessel,  h  :: Heart)

Flow or pressure value is set from the inlet file stored in the `Heart` structure.
"""
function setInletBC(t :: Float64, dt :: Float64, v :: Vessel)
    h = v.heart

  	if h.inlet_type == "Q"
		v.Q[1] = inputFromData(t, h)
  	else
		v.P[1] = inputFromData(t, h)
  	end

	inletCompatibility(dt, v, h)
end

"""
    inputFromData(t :: Float64, h :: Heart)

Inlet flow waveform is stored in `h` data structure inside `input_data` matrix.
`input_data` has two columns and as many rows as time steps.
"""
function inputFromData(t :: Float64, h :: Heart)
	idt = h.input_data[:,1]
	idq = h.input_data[:,2]

	t_hat = div(t,h.cardiac_T)

	t -= t_hat*h.cardiac_T

	idx = 0
	@inbounds for i in 1:length(idt)
		if ((t >= idt[i]) && (t <= idt[i+1]))
			idx=i
			break
		end
	end

	@inbounds qu = idq[idx] + (t - idt[idx]) * (idq[idx+1] - idq[idx]) /
       (idt[idx+1] - idt[idx])

	return qu
end

"""
    inletCompatibility(dt :: Float64, v :: Vessel, h :: Heart)

Compute compatibility relations at inlet node.
"""
function inletCompatibility(dt :: Float64, v :: Vessel, h :: Heart)
	W11, W21 = riemannInvariants(1, v)
	W12, W22 = riemannInvariants(2, v)

	@fastmath @inbounds W11 += (W12 - W11)*(v.c[1] - v.u[1])*dt*v.invDx
	@fastmath @inbounds W21 = 2.0*v.Q[1]/v.A[1] - W11

	v.u[1], v.c[1] = inverseRiemannInvariants(W11, W21)

	if h.inlet_type == "Q"
		@fastmath @inbounds v.A[1] = v.Q[1]/v.u[1]
		v.P[1] = pressure(v.A[1], v.A0[1], v.beta[1], v.Pext)
	else
		v.A[1] = areaFromPressure(v.P[1], v.A0[1], v.beta[1], v.Pext)
		@fastmath @inbounds v.Q[1] = v.u[1]*v.A[1]
	end

end

"""
    riemannInvariants(i :: Int, v :: Vessel)

Calculate Riemann invariants at the node `i` from `u` and `c`.
"""
function riemannInvariants(i :: Int, v :: Vessel)
  @fastmath @inbounds W1 = v.u[i] - 4.0*v.c[i]
  @fastmath @inbounds W2 = v.u[i] + 4.0*v.c[i]

  return W1, W2
end

"""
    inverseRiemannInvariants(W1 :: Float64, W2 :: Float64)

Calculate `u` and `c` given `W1` and `W2`
"""
function inverseRiemannInvariants(W1 :: Float64, W2 :: Float64)
  @fastmath u = 0.5*(W1 + W2)
  @fastmath c = (W2 - W1)*0.125

  return u, c
end

"""
    areaFromPressure(P :: Float64, A0 :: Float64, beta :: Float64, Pext :: Float64)

Inverse constitutive equation. This is used only when a pressure inlet-time-function is
imposed (not recommended).
"""
function areaFromPressure(P :: Float64, A0 :: Float64, beta :: Float64, Pext :: Float64)
   return A0 * ((P-Pext)/beta + 1.0)*((P-Pext)/beta + 1.0)
end

"""
    setOutletBC(dt :: Float64, v :: Vessel)

Outlet boundary condition is applied to the last node of the involved vessel. Two
boundary conditions are available. Either a reflection coefficient is provided or a
three-element windkessel is coupled.
"""
function setOutletBC(dt :: Float64, v :: Vessel)
    if v.outlet == "reflection"
        v.P[end] = 2.0*v.P[end-1] - v.P[end-2]
		outletCompatibility(dt, v)
    elseif v.outlet == "wk3"
        threeElementWindkessel(dt, v)
	end
end

"""
    function outletCompatibility(dt :: Float64, v :: Vessel)

Outlet compatibility relations compute all the quantities not directly assigned by
the outlet boundary condition.
"""
function outletCompatibility(dt :: Float64, v :: Vessel)
	W1M1, W2M1 = riemannInvariants(v.M-1, v)
	W1M, W2M   = riemannInvariants(v.M, v)

	W2M += (W2M1 - W2M)*(v.u[end] + v.c[end])*dt/v.dx
	W1M = v.W1M0 - v.Rt * (W2M - v.W2M0)

	v.u[end], v.c[end] = inverseRiemannInvariants(W1M, W2M)
	v.Q[end] = v.A[end]*v.u[end]
end

"""
    function threeElementWindkessel(dt :: Float64, v :: Vessel)

The three element windkessel simulates the perfusion of downstream vessels.
This 0D model is coupled by `wk3` function to 1D model terminal branches via the
solution of a Riemann problem at the 0D/1D interface.
"""
function threeElementWindkessel(dt :: Float64, v :: Vessel)
	Pout = 0.0

	Al = v.A[end]
	ul = v.u[end]

	v.Pc += dt/v.Cc * (Al*ul - (v.Pc - Pout)/v.R2)

	As = Al

	ssAl = sqrt(sqrt(Al))
	sgamma =  2*sqrt(6*v.gamma[end])
	sA0 = sqrt(v.A0[end])
	bA0 = v.beta[end]/sA0

	fun(As) = As*v.R1*(ul + sgamma * (ssAl - sqrt(sqrt(As)))) -
		(v.Pext + bA0*(sqrt(As) - sA0)) + v.Pc

	dfun(As) = v.R1*(ul + sgamma * (ssAl - 1.25*sqrt(sqrt(As)))) - bA0*0.5/sqrt(As)

    try
        As = newtonSolver(fun, dfun, As)
    catch e
        vlab = v.label
        println("\nNewton solver doesn't converge at $vlab outlet!")
        throw(e)
    end

	us = (pressure(As, v.A0[end], v.beta[end], v.Pext) - Pout)/(As*v.R1)

	v.A[end] = As
	v.u[end] = us
end

"""
    function newtonSolver(f, df, x0 :: Float64)

Solve windkessel equation by means of Newton method.
"""
function newtonSolver(f, df, x0 :: Float64)
	xn = x0 - f(x0)/df(x0)
	if abs(xn-x0)<= 1e-5
		return xn
	else
		newtonSolver(f, df, xn)
	end
end

"""
    updateGhostCells(v :: Vessel)

Ghost cells are updated by assuming both ends as transmissive boundaries.
"""
function updateGhostCells(v :: Vessel)
	v.U00A = v.A[1]
	v.U00Q = v.Q[1]
	v.U01A = v.A[2]
	v.U01Q = v.Q[2]

	v.UM1A = v.A[v.M]
	v.UM1Q = v.Q[v.M]
	v.UM2A = v.A[v.M-1]
	v.UM2Q = v.Q[v.M-1]
end

function updateGhostCells(vessels :: Array{Vessel,1})
	for vessel in vessels
		updateGhostCells(vessel)
	end
end


text:
using Pkg
pkg"activate ."
using Revise
using Layered
using Colors

function example()
    c, tl = canvas(5, 5, color=Gray(0.5))

    l = rectlayer!(tl, c.rect, :w, :norm, margin=10)

    gaussian(x, μ, σ²) = 1 / √(2 * π * σ²) * exp(-(x - μ)^2 / (2σ²))
    scaled_gaussian(x) = (gaussian(x, 0, 1) - gaussian(1, 0, 1)) / (gaussian(0, 0, 1) - gaussian(1, 0, 1))

    nfreqs=6
    basefreqs = range(2, 10, length=nfreqs)
    ndegs = 6
    basedegs = range(0, 180, length=ndegs+1)[1:end-1]

    freqs = [f for f in basefreqs, d in basedegs]
    degs = [d for f in basefreqs, d in basedegs]

    rlayers = [layer!(l, Transform()) for f in freqs]

    rs = rect!.((i, j, freq, d) -> begin
        p = P((i - 0.5) / nfreqs, (j - 0.5) / ndegs)
        r = Rect(p, 1/nfreqs, 1/ndegs, deg(0))

        gfrom = r.center + P(deg(d)) * r.width/2
        gto = r.center - P(deg(d)) * r.width/2
        grad = Gradient(gfrom, gto, Gray.(sin.(range(0, freq*pi, length=50)) .* 0.5 .+ 0.5)...)
        r, Fill(grad)
    end, rlayers, 1:nfreqs, reshape(collect(1:ndegs), 1, :), freqs, degs)

    circs = circle!.((r, freq, d) -> begin
        radius = r.width * 0.4
        c = Circle(r.center, radius)

        alphas = scaled_gaussian.(range(1, 0, length=50))

        grad = RadialGradient(Circle(r.center, 0), Circle(r.center, radius), GrayA.(0.5, alphas)...)

        c, Fill(grad)

    end, rlayers, rs, freqs, degs) .+ Stroke(nothing)

    clipcircs = circle!.((c) -> begin
        scale(c, 0.99)
    end, rlayers, circs) .+ Visible(false)

    rs .+ Clip.(clipcircs)

    txt!.((r, f, d) -> begin
        Txt(center(topline(r)), "$f, $(d)°", 6, :c, :t)
    end, tl, rs, freqs, degs)

    # rlayers .+ Operator(:mult)

    c

    # draw_svg(c, "./examples/02_gabor_patches/gabors.svg")

end; example()


text:
@compat abstract type ABCtype end
@compat abstract type Particle end

#type Prior2
#  distribution
#  Prior2(distributionarray::Array{Distributions.Distribution{Distributions.Univariate,Distributions.Continuous},1}) = new(tuple(distributionarray...))
#end

"""
    Prior(distributions)

    Create Prior type for ABC algorithm specifying priors for each parameters. This is an array of Distribution types from Distribution.jl, each element corresponding to a parameter.
"""
type Prior
  distribution
  Prior(distributionarray) = new(tuple(distributionarray...))
end

type ParticleRejection <: Particle
  params::Array{Float64, 1}
  distance::Float64
  other::Any
end

type ParticleRejectionModel <: Particle
  params::Array{Float64, 1}
  model::Int64
  distance::Float64
  other::Any
end

type ParticleSMC <: Particle
  params::Array{Float64, 1}
  weight::Float64
  scales::Array{Float64, 1}
  distance::Float64
  other::Any
end

type ParticleSMCModel <: Particle
  params::Array{Float64, 1}
  weight::Float64
  scales::Array{Float64, 1}
  model::Int64
  distance::Float64
  other::Any
end

"""
    ABCRejection(sim_func::Function, nparams::Int64, ϵ::Float64, prior::Prior; <keyword arguments>)

Create an ABCRejection type which will simulate data with sim_func. nparams is the number of parameters inputted into sim_func, ϵ is the target tolerance and prior sets the priors for the parameters. sim_func needs to take in 3 values, the parameters (in an array), constants (array) and target data in that order and needs to return 2 values, the first being the distance between the target data and simulated data and the second can be anything but is useful if for example you want to record some additional information about the simulations.
...
## Arguments
- `maxiterations = 10^5`: Maximum number of samples before the ABC algorithm terminates.
- `constants = []`: Any constants needed to simulate from sim_func
- `nparticles = 100`: Number of particles (ie samples) of ABC algorithm
...
"""
type ABCRejection <: ABCtype

  simfunc::Function
  nparams::Int64
  ϵ::Float64
  nparticles::Int64
  constants::Array{Any,1}
  maxiterations::Int64
  prior::Prior

  ABCRejection(sim_func::Function, nparams::Int64, ϵ::Float64, prior::Prior;
    maxiterations = 10000,
    constants = [1.0],
    nparticles = 100,
    ) =
  new(sim_func, nparams, ϵ, nparticles, constants, maxiterations, prior)

end

"""
    ABCRejection(sim_func::Function, nparams::Int64, ϵT::Float64, prior::Prior; <keyword arguments>)

Create an ABCSMC type which will simulate data with sim_func. nparams is the number of parameters inputted into sim_func, ϵT is the target tolerance and prior sets the priors for the parameters. sim_func needs to take in 3 values, the parameters (in an array), constants (array) and target data in that order and needs to return 2 values, the first being the distance between the target data and simulated data and the second can be anything but is useful if for example you want to record some additional information about the simulations.
...
## Arguments
- `maxiterations = 10^5`: Maximum number of samples before the ABC algorithm terminates.
- `constants = []`: Any constants needed to simulate from sim_func
- `nparticles = 100`: Number of particles (ie samples) of ABC algorithm
- `α = 0.3`: The αth quantile of population i is chosen as the ϵ for population i + 1
- `ϵ1 = 10^5`: Starting ϵ for first ABC SMC populations
- `convergence = 0.05`: ABC SMC stops when ϵ in population i + 1 is within 0.05 of populations i
- `scalefactor = 2`: : Parameter for perturbation kernel for parameter values. Larger values means space will be explored more slowly but fewer particles will be perturbed outside prior range.
...
"""
type ABCSMC <: ABCtype

  simfunc::Function
  nparams::Int64
  ϵ1::Float64
  ϵT::Float64
  nparticles::Int64
  constants::Array{Any,1}
  maxiterations::Int64
  prior::Prior
  α::Float64
  convergence::Float64
  scalefactor::Int64

  ABCSMC(sim_func::Function, nparams::Int64, ϵT::Float64, prior::Prior;
    maxiterations = 10^5,
    constants = [1.0],
    nparticles = 100,
    α = 0.3,
    ϵ1 = 10000.0,
    convergence = 0.05,
    scalefactor = 2,
    ) =
  new(sim_func, nparams, ϵ1, ϵT, nparticles, constants, maxiterations, prior, α, convergence, scalefactor)

end

"""
    ABCRejectionModel(sim_func::Array{Function, 1}, nparams::Array{Int64, 1}, ϵ::Float64, prior::Array{Prior, 1}; <keyword arguments>)

Create an ABCRejectionModel type which will create a type to run ABC with model selection. Each model is specified with a function, first input is an array of functions. nparams and priors are arrays for the number of parameters and priors for each model. each sim_func needs to take in 3 values, the parameters (in an array), constants (array) and target data in that order and needs to return 2 values, the first being the distance between the target data and simulated data and the second can be anything but is useful if for example you want to record some additional information about the simulations.
...
## Arguments
- `maxiterations = 10^5`: Maximum number of samples before the ABC algorithm terminates.
- `constants = [[]]`: Any constants needed to simulate from sim_func, needs to be an array of arrays, each one corresponding to a model function.
- `nparticles = 100`: Number of particles (ie samples) of ABC algorithm
...
"""
type ABCRejectionModel <: ABCtype

  Models::Array{ABCRejection, 1}
  nmodels::Int64

  ABCRejectionModel(sim_func::Array{Function, 1}, nparams::Array{Int64, 1}, ϵ::Float64, prior::Array{Prior, 1};
    constants = repeat([[1.0]], outer = length(sim_func)),
    maxiterations = 10000,
    nparticles = 100,
    ) =
  new([ABCRejection(sim_func[i], nparams[i], ϵ, prior[i],  maxiterations = maxiterations, constants = constants[i], nparticles = nparticles) for i in 1:length(sim_func)], length(sim_func))

end

"""
    ABCSMCModel(sim_func::Array{Function, 1}, nparams::Array{Int64, 1}, ϵT::Float64, prior::Array{Prior, 1}; <keyword arguments>)

Create an ABCSMCModel type which will create a type to run the ABC SMC with model selection algorithm. Each model is specified with a function, first input is an array of functions. nparams and priors are arrays for the number of parameters and priors for each model, ϵT is the target tolerance. Each sim_func needs to take in 3 values, the parameters (in an array), constants (array) and target data in that order and needs to return 2 values, the first being the distance between the target data and simulated data and the second can be anything but is useful if for example you want to record some additional information about the simulations.
...
## Arguments
- `maxiterations = 10^5`: Maximum number of samples before the ABC algorithm terminates.
- `constants = []`: Any constants needed to simulate from sim_func
- `nparticles = 100`: Number of particles (ie samples) of ABC algorithm
- `α = 0.3`: The αth quantile of population i is chosen as the ϵ for population i + 1
- `ϵ1 = 10^5`: Starting ϵ for first ABC SMC populations
- `convergence = 0.05`: ABC SMC stops when ϵ in population i + 1 is within 0.05 of populations i
- `scalefactor = 2`: Parameter for perturbation kernel for parameter values. Larger values means space will be explored more slowly but fewer particles will be perturbed outside prior range.
- `modelkern = 0.7`: Probability model stays the same in model perturbation kernel, ie 70% of the time the model perturbation kernel will leave the model the same.
...
"""
type ABCSMCModel <: ABCtype

  Models::Array{ABCSMC, 1}
  nmodels::Int64
  modelkern::Float64
  nparticles::Int64
  α::Float64
  ϵT::Float64
  maxiterations::Int64
  convergence::Float64
  scalefactor::Int64
  other::Any

  function ABCSMCModel(sim_func::Array{Function, 1}, nparams::Array{Int64, 1}, ϵT::Float64, prior::Array{Prior, 1};
    constants = repeat([[1.0]], outer = length(sim_func)),
    maxiterations = 10^5,
    nparticles = 100,
    α = 0.3,
    ϵ1 = 10000.0,
    modelkern = 0.7,
    convergence = 0.05,
    scalefactor = 2,
    other = []
    )
    smcarray = [ABCSMC(sim_func[i], nparams[i], ϵT, prior[i],  maxiterations = maxiterations, constants = constants[i], nparticles = nparticles, α = α, ϵ1 = ϵ1, convergence = convergence) for i in 1:length(sim_func)]
    nmodels = length(sim_func)
    new(smcarray, nmodels, modelkern, nparticles, α, ϵT, maxiterations, convergence, scalefactor, other)
  end

end

type ABCrejectionresults

  parameters::Array{Float64,2}
  accratio::Float64
  numsims::Int64
  dist::Array{Float64,1}
  particles::Array{ParticleRejection, 1}
  setup::ABCRejection

   function ABCrejectionresults(particles, its, ABCsetup, dist)
      parameters = hcat(map(x -> x.params, particles)...)'
      accratio = ABCsetup.nparticles/its
      new(parameters, accratio, its, dist, particles, ABCsetup)
   end
end

type ABCrejectionmodelresults

   parameters
   accratio::Float64
   numsims::Int64
   dist::Array{Float64,1}
   particles::Array{ParticleRejectionModel, 1}
   modelfreq::Array{Float64, 1}
   setup::ABCRejectionModel

   function ABCrejectionmodelresults(particles, its, ABCsetup, dist)
     parameters = []
     modelfreq = []
     for i in 1:ABCsetup.nmodels
        push!(modelfreq, sum(map(x -> x.model, particles) .== i))
        models = map(x -> x.model, particles)
        parameters = push!(parameters,
        hcat(map(x -> x.params, particles[map(x -> x.model, particles) .== i])...)')
     end
     accratio = ABCsetup.Models[1].nparticles/its
     modelfreq = modelfreq ./ sum(modelfreq)

     new(parameters, accratio, its, dist, particles, modelfreq, ABCsetup)
   end
end

type ABCSMCresults

  parameters::Array{Float64,2}
  weights::Array{Float64, 1}
  accratio::Float64
  numsims::Array{Int64,1}
  ϵ::Array{Float64,1}
  particles::Array{ParticleSMC, 1}
  setup::ABCSMC

   function ABCSMCresults(particles, numsims, ABCsetup, epsvec)
      parameters = hcat(map(x -> x.params, particles)...)'
      weights = map(x -> x.weight, particles)
      accratio = ABCsetup.nparticles/sum(numsims)
      new(parameters, weights, accratio, numsims, epsvec, particles, ABCsetup)
   end
end

type ABCSMCmodelresults

   parameters
   weights
   accratio::Float64
   numsims::Array{Int64, 1}
   ϵ::Array{Float64,1}
   particles::Array{ParticleSMCModel, 1}
   modelfreq::Array{Float64, 1}
   modelprob::Array{Float64, 1}
   setup::ABCSMCModel

   function ABCSMCmodelresults(particles, its, ABCsetup, dist)

     parameters = []
     weights = []
     modelfreq = []
     modelweights = map(x -> x.weight, particles)
     models = map(x -> x.model, particles)
     modelprob = []
     for i in 1:ABCsetup.nmodels
        push!(modelfreq, sum(map(x -> x.model, particles) .== i))
        push!(modelprob, sum(modelweights[models .== i]))
        parameters = push!(parameters, hcat(map(x -> x.params, particles[map(x -> x.model, particles) .== i])...)')
        weights = push!(weights, map(x -> x.weight, particles[map(x -> x.model, particles) .== i]))
     end
     accratio = ABCsetup.nparticles ./ sum(its)
     modelfreq = modelfreq ./ sum(modelfreq)

     new(parameters, weights, accratio, its, dist, particles, modelfreq, modelprob, ABCsetup)
   end
end

type SimData
  params::Array{Float64, 1}
  dist::Float64
end


text:
import ACE.Transforms: GetVal

@doc raw"""
`Trig1pBasis`

One-particle basis component of the form
```math
\phi_{l}(\theta) = \exp(i l \theta)
```
Fundamental building block of ACE basis sets of the form
```math
\begin{aligned}
   \phi_{nl}({\bm r}) &= R_{n}(r_{ij}) \exp(i l \theta) \\
   \phi_{nlm}({\bm r}) &= R_{n}(r_{ij}) \exp(i l \theta) T_{m}(z) \\
```
However, for various reasons the implementation takes as input not an angle
``\theta`` but a 2-dimensional vector ``{\bm r}`` taken to be synonymous with
```math
   {\bm r} = r e^{i \theta}.
```
"""
function Trig1pBasis(L::Integer; varsym = :rr, lsym = :l, label = "E$lsym")
   spec = [ (l = l,) for l = -L:L ]
   degrees = [ abs(b.l) for b in spec ]
   return B1pComponent(Trig(L), GetVal{varsym}(), spec, degrees, label)
end

# -------- prototype implementation
#          a more sophisticated implementation could go into the polynomials

using ForwardDiff

struct Trig
   L::Int
end

function cossin(rr)
   r = norm(rr)
   return SVector(rr[1]/r, rr[2]/r)
end

function cossin_ed(rr)
   r = norm(rr)
   c = rr[1]/r
   s = rr[2]/r
   return SVector(c, s), SVector(-s/r, c/r)
end

_theta(rr) = atan(rr[2], rr[1])

function _theta_ed(rr)
   rsq = rr[1]^2+rr[2]^2
   return _theta(rr), SVector(-rr[2]/rsq, rr[1]/rsq)
end

valtype(E::Trig) = ComplexF64
valtype(E::Trig, T) = ComplexF64

evaluate(E::Trig, θ::Real) = [ exp(im * θ * l) for l = -E.L:E.L ]

evaluate(E::Trig, rr::AbstractVector) = evaluate(E, _theta(rr))

ACE.evaluate_d(E::Trig, θ::Real) = [ im * l * exp(im * θ * l)  for l = -E.L:E.L ]

function ACE.evaluate_d(E::Trig, rr)
   θ, dθ = _theta_ed(rr)
   return [ (im * l * exp(im * θ * l)) * dθ  for l = -E.L:E.L ]
end

evaluate_ed(E::Trig, x) = evaluate(E, x), evaluate_d(E, x)


text:
using DataStructures: nlargest

"""
    GlobalPool(aggr, dim...)

Global pooling layer.

It pools all features with `aggr` operation.

# Arguments

- `aggr`: An aggregate function applied to pool all features.
"""
struct GlobalPool{A}
    aggr
    cluster::A
    function GlobalPool(aggr, dim...)
        cluster = ones(Int64, dim)
        new{typeof(cluster)}(aggr, cluster)
    end
end

(g::GlobalPool)(X::AbstractArray) = NNlib.scatter(g.aggr, X, g.cluster)

"""
    LocalPool(aggr, cluster)

Local pooling layer.

It pools features with `aggr` operation accroding to `cluster`. It is implemented with `scatter` operation.

# Arguments

- `aggr`: An aggregate function applied to pool all features.
- `cluster`: An index structure which indicates what features to aggregate with.
"""
struct LocalPool{A<:AbstractArray}
    aggr
    cluster::A
end

(l::LocalPool)(X::AbstractArray) = NNlib.scatter(l.aggr, X, l.cluster)

"""
    TopKPool(adj, k, in_channel)

Top-k pooling layer.

# Arguments

- `adj`: Adjacency matrix  of a graph.
- `k`: Top-k nodes are selected to pool together.
- `in_channel`: The dimension of input channel.
"""
struct TopKPool{T,S}
    A::AbstractMatrix{T}
    k::Int
    p::AbstractVector{S}
    Ã::AbstractMatrix{T}
end

function TopKPool(adj::AbstractMatrix, k::Int, in_channel::Int; init=glorot_uniform)
    TopKPool(adj, k, init(in_channel), similar(adj, k, k))
end

function (t::TopKPool)(X::AbstractArray)
    y = t.p' * X / norm(t.p)
    idx = topk_index(y, t.k)
    t.Ã .= view(t.A, idx, idx)
    X_ = view(X, :, idx) .* σ.(view(y, idx)')
    return X_
end

function topk_index(y::AbstractVector, k::Int)
    v = nlargest(k, y)
    return collect(1:length(y))[y .>= v[end]]
end

topk_index(y::Adjoint, k::Int) = topk_index(y', k)


text:
using transD_GP.GP, Random, PyPlot, Statistics, LinearAlgebra, PositiveFactorizations
##
Random.seed!(3)
f = 1.0/0.5
t = 0:(10*f)^-1:2/f
ntrain = 1
ntries = 10
λ = [0.1, 0.1]
δ = 0.01
δtry = 5δ
p = 2
demean = true
findmin = true
switchflag = findmin ? -1. : 1.
knownvalue = NaN
acqfun = GP.EI()

y = [100000*switchflag*(2x1+ sin.(2*pi*f*x1) + sinc.(2*f*(x1.-middle(x1))))*
                       (2x2 + sin.(2*pi*f*x2) + sinc.(2*f*(x2.-middle(x2)))) for x1 in t, x2 in t]
ynoisy = y + δ*randn(size(y))

X1, X2 = [x1 for x1 in t, x2 in t], [x2 for x1 in t, x2 in t]
T = [X1[:]';X2[:]']
tperm = randperm(length(y))
trainidx = tperm[1:ntrain]
ytrain = ynoisy[trainidx]
ttrain = T[:,trainidx]

vmin, vmax = extrema(ynoisy)
ytest, σ2, σ_prior  = GP.GPfit(GP.SqEuclidean(), ytrain', ttrain, T, λ.^2, δtry, p=p, demean=demean)
AF = GP.getAF(acqfun, ytrain, vec(ytest), diag(σ2), findmin=findmin, knownvalue=knownvalue)
f, ax = plt.subplots(3, ntries+1, figsize=(22,4), sharex=true, sharey=true)
ax[1].imshow(ynoisy, extent=[t[1], t[end], t[end], t[1]])
ax[1].plot(ttrain[2,:][:],ttrain[1,:][:], ".m",markersize=5)
ax[2].imshow(reshape(ytest, size(y)), extent=[t[1], t[end], t[end], t[1]], vmin=vmin, vmax=vmax)
ax[3].imshow(reshape(AF, size(y)), extent=[t[1], t[end], t[end], t[1]])
ax[1].set_title("0")
##
for i = 1:ntries
    global AF
    nextpos = argmax(AF)
    push!(ytrain, ynoisy[nextpos])
    ttrain = hcat(ttrain, T[:,nextpos])
    ytest, σ2, σ_prior  = GP.GPfit(GP.SqEuclidean(), ytrain', ttrain, T, λ.^2, δtry, p=p, demean=demean)
    AF = GP.getAF(acqfun, ytrain, vec(ytest), diag(σ2), findmin=findmin, knownvalue=knownvalue)
    ax[3i+1].imshow(y, extent=[t[1], t[end], t[end], t[1]])
    ax[3i+1].plot(ttrain[2,:][:],ttrain[1,:][:], ".m",markersize=5)
    ax[3i+2].imshow(reshape(ytest, size(y)), extent=[t[1], t[end], t[end], t[1]], vmin=vmin, vmax=vmax)
    ax[3i+3].imshow(reshape(AF, size(y)), extent=[t[1], t[end], t[end], t[1]])
    ax[3i+1].set_title("$i")
end
plt.tight_layout()


text:
abstract type AbstractTerm{𝒯} end

struct PressureDivergence{𝒯} <: AbstractTerm{𝒯} end

@inline calc_component!(flux, ::Nothing, _...) = nothing
@inline calc_component!(flux, ::AbstractTerm, _...) = nothing

@inline function calc_component!(flux, ::PressureDivergence, state, aux, physics)
    eos = physics.eos

    p = calc_pressure(eos, state, aux)

    flux.ρu += p * I

    nothing
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


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

using Markdown
using InteractiveUtils

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

# ╔═╡ e8a5afbe-dc1e-11ea-023e-31002d1c7688
using Pkg, DrWatson, PlutoUI

# ╔═╡ 397f31b2-eba9-11ea-3fd5-1704e3a50a96
md"## Clip-00-01-03t.jl"

# ╔═╡ 1b540d64-dc1f-11ea-321a-f1e0d2bab3e8
md"### snippet 0.1"

# ╔═╡ 624b8daa-dc1f-11ea-114f-19896cd76aa5
"All models are wrong, but some are useful."

# ╔═╡ 8fa79b3e-dc1f-11ea-1437-d7bfe6cdea51
md"### snippet 0.2"

# ╔═╡ c4321a66-de4f-11ea-3449-774edd51eb10
@bind N Slider(1:5, default=3)

# ╔═╡ 35330a6a-dc20-11ea-2929-9b904165e207
md"##### Variable x initially is a StepRange, not a vector. The log.(x) notation `broadcast` the log function to all steprange elements in x and returns a vector."

# ╔═╡ df73cd1c-dc1f-11ea-0e14-c7e381f49c4a
begin
	x = 1:N
	x = x*10
	x = log.(x)
	x = sum(x)
	x = exp(x)
	x = x*10
	x = log(x)
	x = sum(x)
	x = exp(x)
end

# ╔═╡ 975e15f4-122c-11eb-2a4d-ed836c7bade0
md"##### Notice that in Pluto notebooks varables can only be defined once. This is needed in Pluto to make the `reactivity` work (e.g. when updating above slider setting). Bracketing above sequence of assignments with a `begin` and `end` is allowed."

# ╔═╡ b1d6b4ac-dc20-11ea-3237-a515964f9176
md"### snippet 0.3"

# ╔═╡ 88e3bb64-dc1e-11ea-1136-bd822afb72e4
[log(0.01^200) 200 * log(0.01)]

# ╔═╡ c4403280-dc20-11ea-0550-5f2f606a9fea
md"## End of clip-00-01-03t.jl"

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
DrWatson = "634d3b9d-ee7a-5ddf-bec9-22491ea816e1"
Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
DrWatson = "~2.4.1"
PlutoUI = "~0.7.9"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[DrWatson]]
deps = ["Dates", "FileIO", "LibGit2", "MacroTools", "Pkg", "Random", "Requires", "UnPack"]
git-tree-sha1 = "023a12e0004225f5b6860f94f3c0df5e485f54b4"
uuid = "634d3b9d-ee7a-5ddf-bec9-22491ea816e1"
version = "2.4.1"

[[FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "3c041d2ac0a52a12a27af2782b34900d9c3ee68c"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.11.1"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "0fb723cd8c45858c22169b2e42269e53271a6df7"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.7"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "438d35d2d95ae2c5e8780b330592b6de8494e779"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.3"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[UnPack]]
git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
version = "1.0.2"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─397f31b2-eba9-11ea-3fd5-1704e3a50a96
# ╠═e8a5afbe-dc1e-11ea-023e-31002d1c7688
# ╟─1b540d64-dc1f-11ea-321a-f1e0d2bab3e8
# ╠═624b8daa-dc1f-11ea-114f-19896cd76aa5
# ╟─8fa79b3e-dc1f-11ea-1437-d7bfe6cdea51
# ╠═c4321a66-de4f-11ea-3449-774edd51eb10
# ╟─35330a6a-dc20-11ea-2929-9b904165e207
# ╠═df73cd1c-dc1f-11ea-0e14-c7e381f49c4a
# ╟─975e15f4-122c-11eb-2a4d-ed836c7bade0
# ╟─b1d6b4ac-dc20-11ea-3237-a515964f9176
# ╠═88e3bb64-dc1e-11ea-1136-bd822afb72e4
# ╟─c4403280-dc20-11ea-0550-5f2f606a9fea
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = test_page_model()
    include("../src/components/CO2forcing.jl")

    add_comp!(m, co2forcing)

    set_param!(m, :co2forcing, :c_CO2concentration, readpagedata(m, "test/validationdata/$valdir/c_co2concentration.csv"))

    ##running Model
    run(m)

    forcing = m[:co2forcing,:f_CO2forcing]
    forcing_compare = readpagedata(m, "test/validationdata/$valdir/f_co2forcing.csv")

    @test forcing ≈ forcing_compare rtol = 1e-3
end


text:
module Socp
using LinearAlgebra
using SparseArrays
using StaticArrays
using Unrolled
using SuiteSparse

abstract type Cone{D} end
struct POC{D} <: Cone{D} # positive orthant cone
	offs::Int
	POC(offs,dim) = new{dim}(offs)
end
struct SOC{D} <: Cone{D} # second order cone
	offs::Int
	SOC(offs,dim) = new{dim}(offs)
end
@generated conedim(::POC{d}) where d = d
@generated conedim(::SOC{d}) where d = d

struct Problem{C <: Tuple{Vararg{Cone}}, n, m, k, sing}
	#minimize c' x
	c::Vector{Float64} # dim n

	# s.t. A x = b
	A::SparseMatrixCSC{Float64, Int} # dim m x n
	b::Vector{Float64} # dim m

	# s.t. G x + s = h
	G::SparseMatrixCSC{Float64, Int} # dim k x n
	h::Vector{Float64} # dim k

	# wrt cones
	cones::C

	#dimensions
	n::Int
	m::Int
	k::Int

	Problem(c, A, b, G, h, cones::C) where C <: Tuple{Vararg{Cone}} = begin
	    n = length(c)
	    m = size(A)[1]
	    @assert length(b) == m
	    @assert size(A)[2] == n
	    @assert size(G)[2] == n
	    k = size(G)[1]
	    @assert length(h) == k

	    sing = false
	    try
	    	cholesky(G' * G)
		catch e
			if isa(e, PosDefException) || isa(e, SingularException) || isa(e, ArgumentError)
				sing = true
			end
		end
	    #todo: check that each variable appears in a type of cone only once
	    return new{C, n, m, k, sing}(c, A, b, G, h, cones, n, m, k)
	end
end

struct State
	x::Vector{Float64}
	y::Vector{Float64}
	z::Vector{Float64}
	s::Vector{Float64}
	State(prob::Problem, x, y, z, s) = begin
		#todo: initalization proc
		@assert length(x) == prob.n
		@assert length(y) == prob.m
		@assert length(z) == prob.k
		@assert length(s) == prob.k
		return new(x,y,z,s)
	end
end

abstract type AbstractScaling end
abstract type KKTSolver{T} end

include("cholutils.jl")
include("sqrscalings.jl")
include("scalings.jl")
include("vectors.jl")
include("mats.jl")
include("spsolver.jl")
include("densesolver.jl")
include("solver.jl")
include("moi.jl")

end # module


text:
export Lens, set, get, modify
export @lens
export set, get, modify
using ConstructionBase
export setproperties
export constructorof

import Base: get
using Base: setindex, getproperty

"""
    Lens

A `Lens` allows to access or replace deeply nested parts of complicated objects.

# Example
```jldoctest
julia> using Setfield

julia> struct T;a;b; end

julia> obj = T("AA", "BB")
T("AA", "BB")

julia> lens = @lens _.a
(@lens _.a)

julia> get(obj, lens)
"AA"

julia> set(obj, lens, 2)
T(2, "BB")

julia> obj
T("AA", "BB")

julia> modify(lowercase, obj, lens)
T("aa", "BB")
```

# Interface
Concrete subtypes of `Lens` have to implement
* `set(obj, lens, val)`
* `get(obj, lens)`

These must be pure functions, that satisfy the three lens laws:

```jldoctest; output = false, setup = :(using Setfield; obj = (a="A", b="B"); lens = @lens _.a; val = 2; val1 = 10; val2 = 20)
@assert get(set(obj, lens, val), lens) == val
        # You get what you set.
@assert set(obj, lens, get(obj, lens)) == obj
        # Setting what was already there changes nothing.
@assert set(set(obj, lens, val1), lens, val2) == set(obj, lens, val2)
        # The last set wins.

# output

```

See also [`@lens`](@ref), [`set`](@ref), [`get`](@ref), [`modify`](@ref).
"""
abstract type Lens end

"""
    modify(f, obj, l::Lens)

Replace a deeply nested part `x` of `obj` by `f(x)`. See also [`Lens`](@ref).
"""
function modify end

"""
    get(obj, l::Lens)

Access a deeply nested part of `obj`. See also [`Lens`](@ref).
"""
function get end

"""
    set(obj, l::Lens, val)

Replace a deeply nested part of `obj` by `val`. See also [`Lens`](@ref).
"""
function set end

@inline function modify(f, obj, l::Lens)
    old_val = get(obj, l)
    new_val = f(old_val)
    set(obj, l, new_val)
end

struct IdentityLens <: Lens end
get(obj, ::IdentityLens) = obj
set(obj, ::IdentityLens, val) = val

struct PropertyLens{fieldname} <: Lens end

function get(obj, l::PropertyLens{field}) where {field}
    getproperty(obj, field)
end

@generated function set(obj, l::PropertyLens{field}, val) where {field}
    Expr(:block,
         Expr(:meta, :inline),
        :(setproperties(obj, ($field=val,)))
       )
end

struct ComposedLens{LO, LI} <: Lens
    outer::LO
    inner::LI
end

compose() = IdentityLens()
compose(l::Lens) = l
compose(::IdentityLens, ::IdentityLens) = IdentityLens()
compose(::IdentityLens, l::Lens) = l
compose(l::Lens, ::IdentityLens) = l
compose(outer::Lens, inner::Lens) = ComposedLens(outer, inner)
function compose(l1::Lens, ls::Lens...)
    # We can build _.a.b.c as (_.a.b).c or _.a.(b.c)
    # The compiler prefers (_.a.b).c
    compose(l1, compose(ls...))
end

"""
    lens₁ ∘ lens₂
    compose([lens₁, [lens₂, [lens₃, ...]]])

Compose lenses `lens₁`, `lens₂`, ..., `lensₙ` to access nested objects.

# Example
```jldoctest
julia> using Setfield

julia> obj = (a = (b = (c = 1,),),);

julia> la = @lens _.a
       lb = @lens _.b
       lc = @lens _.c
       lens = la ∘ lb ∘ lc
(@lens _.a.b.c)

julia> get(obj, lens)
1
```
"""
Base.:∘(l1::Lens, l2::Lens) = compose(l1, l2)

function get(obj, l::ComposedLens)
    inner_obj = get(obj, l.outer)
    get(inner_obj, l.inner)
end

function set(obj,l::ComposedLens, val)
    inner_obj = get(obj, l.outer)
    inner_val = set(inner_obj, l.inner, val)
    set(obj, l.outer, inner_val)
end

struct IndexLens{I <: Tuple} <: Lens
    indices::I
end

Base.@propagate_inbounds function get(obj, l::IndexLens)
    getindex(obj, l.indices...)
end
Base.@propagate_inbounds function set(obj, l::IndexLens, val)
    setindex(obj, val, l.indices...)
end

"""
    ConstIndexLens{I}

Lens with index stored in type parameter.  This is useful for type-stable
[`get`](@ref) and [`set`](@ref) operations on tuples and named tuples.

This lens can be constructed by, e.g., `@lens _[\$1]`.  Complex expression
must be wrapped with `\$(...)` like `@lens _[\$(length(xs))]`.

# Examples
```jldoctest
julia> using Setfield

julia> get((1, 2.0), @lens _[\$1])
1

julia> Base.promote_op(get, typeof.(((1, 2.0), @lens _[\$1]))...)
Int64

julia> Base.promote_op(get, typeof.(((1, 2.0), @lens _[1]))...) !== Int
true
```
"""
struct ConstIndexLens{I} <: Lens end

Base.@propagate_inbounds get(obj, ::ConstIndexLens{I}) where I = obj[I...]

Base.@propagate_inbounds set(obj, ::ConstIndexLens{I}, val) where I =
    setindex(obj, val, I...)

@generated function set(obj::Union{Tuple, NamedTuple},
                        ::ConstIndexLens{I},
                        val) where I
    if length(I) == 1
        n, = I
        args = map(1:length(obj.types)) do i
            i == n ? :val : :(obj[$i])
        end
        quote
            $(Expr(:meta, :inline))
            ($(args...),)
        end
    else
        quote
            throw(ArgumentError($(string(
                "A `Tuple` and `NamedTuple` can only be indexed with one ",
                "integer.  Given: $I"))))
        end
    end
end

struct DynamicIndexLens{F} <: Lens
    f::F
end

Base.@propagate_inbounds get(obj, I::DynamicIndexLens) = obj[I.f(obj)...]

Base.@propagate_inbounds set(obj, I::DynamicIndexLens, val) =
    setindex(obj, val, I.f(obj)...)

"""
    FunctionLens(f)
    @lens f(_)

Lens with [`get`](@ref) method definition that simply calls `f`.
[`set`](@ref) method for each function `f` must be implemented manually.
Use `methods(set, (Any, Setfield.FunctionLens, Any))` to get a list of
supported functions.

Note that `FunctionLens` flips the order of composition; i.e.,
`(@lens f(_)) ∘ (@lens g(_)) == @lens g(f(_))`.

# Example
```jldoctest
julia> using Setfield

julia> obj = ((1, 2), (3, 4));

julia> lens = (@lens first(_)) ∘ (@lens last(_))
(@lens last(first(_)))

julia> get(obj, lens)
2

julia> set(obj, lens, '2')
((1, '2'), (3, 4))
```

# Implementation

To use `myfunction` as a lens, define a `set` method with the following
signature:

```julia
Setfield.set(obj, ::typeof(@lens myfunction(_)), val) = ...
```

`typeof` is used above instead of `FunctionLens` because how actual
type of `@lens myfunction(_)` is implemented is not the part of stable
API.
"""
struct FunctionLens{f} <: Lens end
FunctionLens(f) = FunctionLens{f}()

get(obj, ::FunctionLens{f}) where f = f(obj)

Base.@deprecate constructor_of(T) constructorof(T)
Base.@deprecate get(lens::Lens, obj)       get(obj, lens)
Base.@deprecate set(lens::Lens, obj, val)  set(obj, lens, val)
Base.@deprecate modify(f, lens::Lens, obj) modify(f, obj, lens)


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

# This file is primarily to be included from runtest.jl. We tried to cover various
# features of MXNet.jl in this example in order to detect regression errors.

module MNISTTest

using MXNet
using Test

include("mnist-data.jl")

function get_mnist_mlp()
  @mx.chain mx.Variable(:data)                   =>
    mx.FullyConnected(name=:fc1, num_hidden=128) =>
    mx.Activation(name=:relu1, act_type=:relu)   =>
    mx.FullyConnected(name=:fc2, num_hidden=64)  =>
    mx.Activation(name=:relu2, act_type=:relu)   =>
    mx.FullyConnected(name=:fc3, num_hidden=10)  =>
    mx.SoftmaxOutput(name=:softmax)
end

get_mnist_data(batch_size = 100) = get_mnist_providers(batch_size)

function mnist_fit_and_predict(optimizer, initializer, n_epoch)
  mlp = get_mnist_mlp()
  train_provider, eval_provider = get_mnist_data()

  # setup model
  model = mx.FeedForward(mlp, context = mx.cpu())

  # fit parameters
  cp_prefix = "mnist-test-cp"
  mx.fit(model, optimizer, train_provider, eval_data=eval_provider, n_epoch=n_epoch,
         initializer=initializer, callbacks=[mx.speedometer(), mx.do_checkpoint(cp_prefix, save_epoch_0=true)])

  # make sure the checkpoints are saved
  @test isfile("$cp_prefix-symbol.json")
  for i_epoch = 0:n_epoch
    @test isfile(mx.format("{1}-{2:04d}.params", cp_prefix, i_epoch))
  end
  mlp_load = mx.load("$cp_prefix-symbol.json", mx.SymbolicNode)
  @test mx.to_json(mlp_load) == mx.to_json(mlp)
  mlp_load = mx.from_json(read("$cp_prefix-symbol.json", String), mx.SymbolicNode)
  @test mx.to_json(mlp_load) == mx.to_json(mlp)

  #--------------------------------------------------------------------------------
  # the predict API
  probs = mx.predict(model, eval_provider)

  # collect all labels from eval data
  labels = Array[]
  for batch in eval_provider
    push!(labels, copy(mx.get(eval_provider, batch, :softmax_label)))
  end
  labels = cat(labels..., dims = 1)

  # Now we use compute the accuracy
  correct = 0
  for i = 1:length(labels)
    # labels are 0...9
    if argmax(probs[:,i]) == labels[i]+1
      correct += 1
    end
  end
  accuracy = 100correct/length(labels)
  println(mx.format("Accuracy on eval set: {1:.2f}%", accuracy))

  # try to call visualization
  dot_code = mx.to_graphviz(mlp)

  return accuracy
end

function test_mnist_mlp()
  @info("MNIST::SGD")
  @test mnist_fit_and_predict(mx.SGD(η=.2), mx.UniformInitializer(.01), 2) > 90

  @info("MNIST::SGD::η scheduler")
  @test mnist_fit_and_predict(mx.SGD(η_sched=mx.LearningRate.Inv(.25)),
                              mx.UniformInitializer(.01), 2) > 90

  @info("MNIST::SGD::momentum μ")
  @test mnist_fit_and_predict(mx.SGD(η=.1, μ=.9), mx.UniformInitializer(.01), 2) > 90

  @info("MNIST::ADAM")
  @test mnist_fit_and_predict(mx.ADAM(), mx.NormalInitializer(), 2) > 90

  @info("MNIST::AdaGrad")
  @test mnist_fit_and_predict(mx.AdaGrad(), mx.NormalInitializer(), 2) > 90

  @info("MNIST::AdaDelta")
  @test mnist_fit_and_predict(mx.AdaDelta(), mx.NormalInitializer(), 2) > 90

  @info("MNIST::AdaMax")
  @test mnist_fit_and_predict(mx.AdaMax(), mx.NormalInitializer(), 2) > 90

  @info("MNIST::RMSProp")
  @test mnist_fit_and_predict(mx.RMSProp(), mx.NormalInitializer(), 2) > 90

  @info("MNIST::Nadam")
  @test mnist_fit_and_predict(mx.Nadam(), mx.NormalInitializer(), 2) > 90
end

test_mnist_mlp()

end # module MNISTTest