text:
"""
As the cylinder flow blows up somehow,
here is the code which uses the split boundary under a lower Mach number.
The initial still gas is driven by the far-field flow.
"""

using KitBase, FluxReconstruction, LinearAlgebra, OrdinaryDiffEq, Plots
using KitBase.OffsetArrays
using KitBase.ProgressMeter: @showprogress
using Base.Threads: @threads

pyplot()
cd(@__DIR__)

begin
    set = Setup(
        case = "cylinder",
        space = "2d0f",
        flux = "hll",
        collision = "nothing",
        interpOrder = 3,
        limiter = "positivity",
        boundary = "euler",
        cfl = 0.1,
        maxTime = 1.0,
    )
    ps = begin
        ps0 = KitBase.CSpace2D(1.0, 6.0, 10, 0.0, π, 30, 0, 1)
        deg = set.interpOrder - 1
        FRPSpace2D(ps0, deg)
    end
    vs = nothing
    gas = Gas(Kn = 1e-6, Ma = 0.5, K = 1.0)
    ib = nothing

    ks = SolverSet(set, ps0, vs, gas, ib)
end

u0 = OffsetArray{Float64}(undef, 1:ps.nr, 0:ps.nθ+1, deg + 1, deg + 1, 4)
for i in axes(u0, 1), j in axes(u0, 2), k in axes(u0, 3), l in axes(u0, 4)
    sos = sound_speed(1.0, ks.gas.γ)
    prim = begin
        if i == ps.nr
            [1.0, ks.gas.Ma * sos, 0.0, 1.0]
        else
            [1.0, 0.0, 0.0, 1.0]
        end
    end
    u0[i, j, k, l, :] .= prim_conserve(prim, ks.gas.γ)
end

n1 = [[0.0, 0.0] for i = 1:ps.nr+1, j = 1:ps.nθ]
for i = 1:ps.nr+1, j = 1:ps.nθ
    angle = sum(ps.dθ[1, 1:j-1]) + 0.5 * ps.dθ[1, j]
    n1[i, j] .= [cos(angle), sin(angle)]
end

n2 = [[0.0, 0.0] for i = 1:ps.nr, j = 1:ps.nθ+1]
for i = 1:ps.nr, j = 1:ps.nθ+1
    angle = π / 2 + sum(ps.dθ[1, 1:j-1])
    n2[i, j] .= [cos(angle), sin(angle)]
end

function dudt!(du, u, p, t)
    du .= 0.0

    iJ, ll, lr, dhl, dhr, lpdm, γ = p

    nx = size(u, 1)
    ny = size(u, 2) - 2
    nsp = size(u, 3)

    f = OffsetArray{Float64}(undef, 1:nx, 0:ny+1, nsp, nsp, 4, 2)
    @inbounds @threads for i in axes(f, 1)
        for j in axes(f, 2), k = 1:nsp, l = 1:nsp
            fg, gg = euler_flux(u[i, j, k, l, :], γ)
            for m = 1:4
                f[i, j, k, l, m, :] .= iJ[i, j][k, l] * [fg[m], gg[m]]
            end
        end
    end

    u_face = OffsetArray{Float64}(undef, 1:nx, 0:ny+1, 4, nsp, 4)
    f_face = OffsetArray{Float64}(undef, 1:nx, 0:ny+1, 4, nsp, 4, 2)
    @inbounds @threads for i in axes(u_face, 1)
        for j in axes(u_face, 2), l = 1:nsp, m = 1:4
            u_face[i, j, 1, l, m] = dot(u[i, j, l, :, m], ll)
            u_face[i, j, 2, l, m] = dot(u[i, j, :, l, m], lr)
            u_face[i, j, 3, l, m] = dot(u[i, j, l, :, m], lr)
            u_face[i, j, 4, l, m] = dot(u[i, j, :, l, m], ll)

            for n = 1:2
                f_face[i, j, 1, l, m, n] = dot(f[i, j, l, :, m, n], ll)
                f_face[i, j, 2, l, m, n] = dot(f[i, j, :, l, m, n], lr)
                f_face[i, j, 3, l, m, n] = dot(f[i, j, l, :, m, n], lr)
                f_face[i, j, 4, l, m, n] = dot(f[i, j, :, l, m, n], ll)
            end
        end
    end

    fx_interaction = zeros(nx + 1, ny, nsp, 4)
    @inbounds @threads for i = 2:nx
        for j = 1:ny, k = 1:nsp
            fw = @view fx_interaction[i, j, k, :]
            uL = local_frame(u_face[i-1, j, 2, k, :], n1[i, j][1], n1[i, j][2])
            uR = local_frame(u_face[i, j, 4, k, :], n1[i, j][1], n1[i, j][2])
            flux_hll!(fw, uL, uR, γ, 1.0)
            fw .= global_frame(fw, n1[i, j][1], n1[i, j][2])
        end
    end

    @inbounds @threads for j = 1:ny
        for k = 1:nsp
            ul = local_frame(u_face[1, j, 4, k, :], n1[1, j][1], n1[1, j][2])
            fw = @view fx_interaction[1, j, k, :]
            flux_boundary_maxwell!(
                fw,
                [1., 0., 0., 1.],
                ul,
                1,
                γ,
                1.0,
                1.0,
                1,
            )
            fw .= global_frame(fw, n1[1, j][1], n1[1, j][2])
        end
    end

    fy_interaction = zeros(nx, ny + 1, nsp, 4)
    @inbounds @threads for i = 1:nx
        for j = 1:ny+1, k = 1:nsp
            fw = @view fy_interaction[i, j, k, :]
            uL = local_frame(u_face[i, j-1, 3, k, :], n2[i, j][1], n2[i, j][2])
            uR = local_frame(u_face[i, j, 1, k, :], n2[i, j][1], n2[i, j][2])
            flux_hll!(fw, uL, uR, γ, 1.0)
            fw .= global_frame(fw, n2[i, j][1], n2[i, j][2])
        end
    end

    rhs1 = zeros(nx, ny, nsp, nsp, 4)
    @inbounds @threads for i = 1:nx
        for j = 1:ny, k = 1:nsp, l = 1:nsp, m = 1:4
            rhs1[i, j, k, l, m] = dot(f[i, j, :, l, m, 1], lpdm[k, :])
        end
    end
    rhs2 = zeros(nx, ny, nsp, nsp, 4)
    @inbounds @threads for i = 1:nx
        for j = 1:ny, k = 1:nsp, l = 1:nsp, m = 1:4
            rhs2[i, j, k, l, m] = dot(f[i, j, k, :, m, 2], lpdm[l, :])
        end
    end

    @inbounds @threads for i = 2:nx-1
        for j = 1:ny, k = 1:nsp, l = 1:nsp, m = 1:4
            fxL = (inv(ps.Ji[i, j][4, l])*n1[i, j])[1] * fx_interaction[i, j, l, m]
            fxR = (inv(ps.Ji[i, j][2, l])*n1[i+1, j])[1] * fx_interaction[i+1, j, l, m]
            fyL = (inv(ps.Ji[i, j][1, k])*n2[i, j])[2] * fy_interaction[i, j, l, m]
            fyR = (inv(ps.Ji[i, j][3, k])*n2[i, j+1])[2] * fy_interaction[i, j+1, l, m]
            du[i, j, k, l, m] = -(
                rhs1[i, j, k, l, m] +
                rhs2[i, j, k, l, m] +
                (fxL - f_face[i, j, 4, l, m, 1]) * dhl[k] +
                (fxR - f_face[i, j, 2, l, m, 1]) * dhr[k] +
                (fyL - f_face[i, j, 1, k, m, 2]) * dhl[l] +
                (fyR - f_face[i, j, 3, k, m, 2]) * dhr[l]
            )
        end
    end

    return nothing
end

tspan = (0.0, 1.0)
p = (ps.iJ, ps.ll, ps.lr, ps.dhl, ps.dhr, ps.dl, ks.gas.γ)
prob = ODEProblem(dudt!, u0, tspan, p)

dt = 0.0005
nt = tspan[2] ÷ dt |> Int
itg = init(prob, Euler(), save_everystep = false, adaptive = false, dt = dt)

@showprogress for iter = 1:20#nt
    @inbounds for i = 1:ps.nr, k = 1:ps.deg+1, l = 1:ps.deg+1
        u1 = itg.u[i, 1, 4-k, 4-l, :]
        ug1 = [u1[1], u1[2], -u1[3], u1[4]]
        itg.u[i, 0, k, l, :] .= ug1

        u2 = itg.u[i, ps.nθ, 4-k, 4-l, :]
        ug2 = [u2[1], u2[2], -u2[3], u2[4]]
        itg.u[i, ps.nθ+1, k, l, :] .= ug2
    end
    @inbounds for j = 1:ps.nθ÷2, k = 1:ps.deg+1, l = 1:ps.deg+1
        itg.u[ps.nr, j, k, l, :] .= itg.u[ps.nr-1, j, k, l, :]
    end

    step!(itg)
end

sol = zeros(ps.nr, ps.nθ, 4)
for i = 1:ps.nr, j = 1:ps.nθ
    sol[i, j, :] .= conserve_prim(itg.u[i, j, 1, 1, :], ks.gas.γ)
    sol[i, j, 4] = 1 / sol[i, j, 4]
end

contourf(
    ps.x[1:ps.nr-1, 1:ps.nθ],
    ps.y[1:ps.nr-1, 1:ps.nθ],
    sol[1:ps.nr-1, :, 2],
    aspect_ratio = 1,
    legend = true,
)


text:
SRC_DIR = joinpath(Base.source_dir(), "..", "..", "src")
OUT_DIR = joinpath(Base.source_dir(), "..", "..", "output")
mkpath(OUT_DIR)
#NOTE additional packages that need installing:
#=
    CSV, DataFrames, RollingFunctions

=#

#include(joinpath(SRC_DIR, "BridgeSDEInference.jl"))
#using Main.BridgeSDEInference
include(joinpath(SRC_DIR, "BridgeSDEInference_for_tests.jl"))
#using BridgeSDEInference
#const BSI = BridgeSDEInference
using Distributions # to define priors
using Random        # to seed the random number generator
using DataFrames
using CSV
using LinearAlgebra

DIR = "auxiliary"
include(joinpath(SRC_DIR, DIR, "read_and_write_data.jl"))
include(joinpath(SRC_DIR, DIR, "transforms.jl"))
include(joinpath(SRC_DIR, DIR, "utility_functions.jl"))

# decide if first passage time observations or partially observed diffusion
fptObsFlag = false

# pick dataset
filename = "path_part_obs_conj.csv"

# fetch the data
(df, x0, obs, obs_time, fpt,
      fptOrPartObs) = readData(Val(fptObsFlag), joinpath(OUT_DIR, filename))
param = :complexConjug
θ_init = [10.0, -8.0, 15.0, 0.0, 3.0]
P˟ = FitzhughDiffusion(param, θ_init...)
P̃ = std_aux_laws(FitzhughDiffusionAux, param, θ_init, obs, obs_time, 1)
L = @SMatrix [1. 0.]
Σ = @SMatrix [10^(-10)]

model_setup = DiffusionSetup(P˟, P̃, fptOrPartObs)
set_observations!(model_setup, [L for _ in P̃], [Σ for _ in P̃], obs, obs_time, fpt)
set_imputation_grid!(model_setup, 1/1000)
set_x0_prior!(model_setup, KnownStartingPt(x0))
initialise!(eltype(x0), model_setup, Vern7(), false, NoChangePt(100))
set_auxiliary!(model_setup; skip_for_save=10^0, adaptive_prop=NoAdaptation())

mcmc_setup = MCMCSetup(
      Imputation(NoBlocking(), 0.975, Vern7()),
      ParamUpdate(MetropolisHastingsUpdt(), 5, θ_init,
                  UniformRandomWalk(0.5, true), ImproperPosPrior(),
                  UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 5))
                  ),
      ParamUpdate(ConjugateUpdt(), [1,2,3], θ_init, nothing,
                  MvNormal(fill(0.0, 3), diagm(0=>fill(1000.0, 3))),
                  UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [1,2,3]))
                  ))

schedule = MCMCSchedule(1*10^4, [[1,2,3]],
                        (save=3*10^2, verbose=10^2, warm_up=100,
                         readjust=(x->x%100==0), fuse=(x->false)))

Random.seed!(4)
out, elapsed = @timeit mcmc(mcmc_setup, schedule, model_setup)
#out = mcmc(mcmc_setup, schedule, model_setup)

include(joinpath(SRC_DIR, DIR, "plotting_fns.jl"))
plot_chains(out[2]; truth=[10.0, -8.0, 15.0, 0.0, 3.0])
plot_paths(out[1], out[2], schedule; transf=[(x,θ)->x, (x,θ)->conjugToRegular(x, θ[1], 0)],
           obs=(times=obs_time[2:end], vals=obs[2:end], indices=1))

ext: jl


text:
__precompile__(true)

module SupStable

  using Distributions
  using StatsBase
  using LambertW
  using FastGaussQuadrature

  importall Distributions

  export
    # distribution types
    PositiveStable,
    SupremumStable,
    Stable,

    # methods
    rand,
    maximum,
    minimum,
    insupport,
    pdf,
    cdf,
    mgf,
    cf,
    mellin,
    mean,
    params

  # Source Files

  include("positivestable.jl")
  include("supremumstable.jl")
  include("stable.jl")

end


text:
using Main.Test
using Turing: _leapfrog

D = 10

lp_grad_func(x) = nothing, x

step_size = 0.1

# Test 1 - single mutiple-step call v.s. plain Julia

@testset "1. single mutiple-step call v.s. plain Julia" begin

    for _ = 1:100

        theta = rand(D)
        p = rand(D)

        theta_turing, p_turing, _ = _leapfrog(theta, p, 1, step_size, lp_grad_func)

        p -= step_size .* lp_grad_func(theta)[2] / 2
        theta += step_size .* p
        p -= step_size .* lp_grad_func(theta)[2] / 2

        @test theta == theta_turing
        @test p == p_turing

    end

end

# Test 2 - multiple single-step call v.s. single mutiple-step call

@testset "2. multiple single-step call v.s. single mutiple-step call" begin

    for _ = 1:100

        theta_0 = rand(2)
        p_0 = rand(2)

        theta_1, p_1, _ = _leapfrog(theta_0, p_0, 1, step_size, lp_grad_func)

        theta_2, p_2, _ = _leapfrog(theta_1, p_1, 1, step_size, lp_grad_func)

        theta_2_one_call, p_2_one_call, _ = _leapfrog(theta_0, p_0, 2, step_size, lp_grad_func)

        @test theta_2 == theta_2_one_call
        @test p_2 == p_2_one_call

    end

end


text:
# Quick testing function:
using Revise
using RadiativeTransfer.Scattering
using ForwardDiff, DiffResults
using Distributions

x = [0.3, 6.82, 1.33,0.00001]
x2 = copy(x);
δ = 1e-8
x2[3] = x[3]+δ

result = DiffResults.JacobianResult(zeros(4568),x);

result = ForwardDiff.jacobian!(result, Scattering.f_test, x);

a1= Scattering.f_test(x)

@time a2= Scattering.f_test(x2)

plot(DiffResults.jacobian(result)[:,3], label="ForwardDiff")
plot!((a2 - a1)/δ, label="Finite Diff")

# Test CUDA Matrix stuff:
elty = Float32
m = 32
A = [rand(elty,m,m) for i in 1:10000]
# move to device
d_A = CuArray{elty, 2}[]
for i in 1:length(A)
    push!(d_A,CuArray(A[i]))
end

function invCUDA(d_A)
    pivot, info = CUBLAS.getrf_batched!(d_A, true)
    pivot, info, d_C = CUBLAS.getri_batched(d_A, pivot)
    return d_C
end


text:
using MAT, Meganet, BenchmarkTools, Compat, JLD, ProfileView

# Macro Benchmark on CIFAR10
n = 64
miniBatchSize = 32

path2data = Pkg.dir("Meganet")*"/data/CIFAR10/"
history = Pkg.dir("Meganet")*"/benchmarks/CIFAR10/cifar10_512_64.jld"

Y_train,C_train,Y_test,C_test = getCIFAR10(n, path2data)

nImg = [32; 32]
cin  = 3
nc   = [16;32;64;64]
nt   = 2*[1;1;1]
h    = [1.;1.;1.]

TYPE = Float32;

getConvKernel = (nImg,sK) -> getConvGEMMKernel(TYPE,nImg,sK);

# opening layer
K1 = getConvKernel(nImg,[3,3,cin,nc[1]]);

nL = getBatchNormLayer(TYPE,[prod(nImg);nc[1]],isTrainable=true);
blocks = [getSingleLayer(TYPE,K1,nL)]

for k=1:length(nt)
    # ResNN layers
    K2 = getConvKernel(nImg,[3,3,nc[k],nc[k]])
    nL = getBatchNormLayer(TYPE,[prod(nImg);nc[k]],isTrainable=true)
    L2 = getDoubleSymLayer(TYPE,K2,nL)
    RN  = getResNN(TYPE,L2,nt[k],h[k])

    if k<length(nt)
        RN.outTimes *=0
    end
    blocks = [blocks;RN]
    # change channels
    Kc = getConvKernel(nImg,[1,1,nc[k],nc[k+1]]);

    nL = getBatchNormLayer(TYPE,[prod(nImg);nc[k+1]],isTrainable=true)
    blocks = [blocks; getSingleLayer(TYPE,Kc,nL)]

    if k<length(nt)
        Kp = getAverageMatrix(TYPE,nImg,nc[k+1])
       	blocks = [blocks; getConnector(TYPE,Kp)]
        # nImg ./=2
       	nImg = div.(nImg,2)
    end
end

# Connector block
B      = kron(speye(TYPE,nc[end]),ones(TYPE, prod(nImg)))/prod(nImg);
blocks = [blocks; getConnector(TYPE,B')];
blocks[end].outTimes=1

net = getNN(blocks)
theta = initTheta(net);

display(net)

# regularizers
pRegTh = getTikhonovReg(TYPE;alpha=4e-4)
pRegW = getTikhonovReg(TYPE;alpha=4e-4)
pLoss = getSoftMaxLoss(TYPE);
objFun = dnnObjFctn(net,pLoss,pRegTh,pRegW)
opt = getSGDsolver(TYPE,learningRate=1e-2,maxEpochs=1,miniBatch=miniBatchSize,out=true)

W      = 0.1*vec(randn(TYPE,10,nFeatOut(net)+1));
W = min.(W,.2);
W = max.(W,-.2);
W = convert(Array{TYPE},W);

# Save benchmarks
solve(opt,objFun::dnnObjFctn,[vec(theta);vec(W)],Y_train,C_train,Y_test,C_test)

#@profiler solve(opt,objFun::dnnObjFctn,[vec(theta);vec(W)],Y_train,C_train,Y_test,C_test)
Profile.clear()
Profile.init(100000000, 0.001)
@profile solve(opt,objFun::dnnObjFctn,[vec(theta);vec(W)],Y_train,C_train,Y_test,C_test)
ProfileView.view()


text:
# This file is a part of SimilaritySearch.jl
struct KnnResultState
    shift::Int
end

export initialstate, KnnResultShift

"""
    KnnResultShift(ksearch::Integer)

Creates a priority queue with fixed capacity (`ksearch`) representing a knn result set.
It starts with zero items and grows with [`push!(res, id, dist)`](@ref) calls until `ksearch`
size is reached. After this only the smallest items based on distance are preserved.
"""
struct KnnResultShift
    id::Vector{Int32}
    dist::Vector{Float32}
    k::Int # number of neighbors
end

function KnnResultShift(k::Integer)
    @assert k > 0
    res = KnnResultShift(Vector{Int32}(undef, 0), Vector{Float32}(undef, 0), k)
    sizehint!(res.id, k)
    sizehint!(res.dist, k)
    res
end

function initialstate(::KnnResultShift)
    KnnResultState(0)
end

"""
    push!(res::KnnResultShift, item::Pair)
    push!(res::KnnResultShift, id::Integer, dist::Real)

Appends an item into the result set
"""
@inline function Base.push!(res::KnnResultShift, st::KnnResultState, id::Integer, dist::Real)
    if length(res, st) < maxlength(res)
        k = res.k
        if length(res.id) >= 2k-1
            compact!(res, st, 1)
            st = KnnResultState(0)
            @inbounds res.id[end], res.dist[end] = id, dist
        else
            push!(res.id, id)
            push!(res.dist, dist)
        end

        _shifted_fixorder!(res, st.shift)
        #_shifted_fixorder!(res.shift, res.id, res.dist)
        return st
    end

    dist >= last(res.dist) && return st

    @inbounds res.id[end], res.dist[end] = id, dist
    _shifted_fixorder!(res, st.shift)
    #_shifted_fixorder!(res.shift, res.id, res.dist)
    st
end

#@inline Base.push!(res::KnnResultShift, st::KnnResultState, id::Integer, dist::Real) = push!(res, st, convert(Int32, id), convert(Float32, dist))
@inline Base.push!(res::KnnResultShift, st::KnnResultState, p::Pair) = push!(res, st, p.first, p.second)

function compact!(res::KnnResultShift, st::KnnResultState, resize_extra)
    shift = st.shift
    if shift > 0
        n = length(res, st)
        j = shift
        @inbounds for i in 1:n
            j += 1
            res.id[i] = res.id[j]
            res.dist[i] = res.dist[j]
        end

        resize!(res.id, n+resize_extra)
        resize!(res.dist, n+resize_extra)
    end

    res
end

"""
    popfirst!(p::KnnResultShift, st::KnnResultState)

Removes and returns the nearest neeighboor pair from the pool, an O(length(p.pool)) operation
"""
@inline function Base.popfirst!(res::KnnResultShift, st::KnnResultState)
    p = argmin(res, st) => minimum(res, st)
    res.id[1] = 0  # mark as deleted
    p, KnnResultState(st.shift+1)
end

"""
    pop!(res::KnnResultShift, st::KnnResultState)

Removes and returns the last item in the pool, it is an O(1) operation
"""
@inline function Base.pop!(res::KnnResultShift, st::KnnResultState)
    pop!(res.id) => pop!(res.dist), st
end

"""
    maxlength(res::KnnResultShift)

The maximum allowed cardinality (the k of knn)
"""
@inline maxlength(res::KnnResultShift) = res.k
@inline Base.length(res::KnnResultShift, st::KnnResultState) = length(res.id) - st.shift

function Base.length(res::KnnResultShift)
    i = 1
    n = length(res.id)
    while i < n && res.id[i] == 0
        i += 1
    end

    n - i + 1
end

"""
    reuse!(res::KnnResultShift)
    reuse!(res::KnnResultShift, k::Integer)

Returns a result set and a new initial state; reuse the memory buffers
"""
@inline function reuse!(res::KnnResultShift, k::Integer=res.k)
    @assert k > 0
    empty!(res.id)
    empty!(res.dist)
    if k > res.k
        sizehint!(res.id, k)
        sizehint!(res.dist, k)
    end
    KnnResultShift(res.id, res.dist, k)
end

"""
    getindex(res::KnnResultShift, st::KnnResultState, i)

Access the i-th item in `res`
"""
@inline function getpair(res::KnnResultShift, st::KnnResultState, i)
    i += st.shift
    @inbounds res.id[i] => res.dist[i]
end

@inline getid(res::KnnResultShift, st::KnnResultState, i) = @inbounds res.id[i+st.shift]
@inline getdist(res::KnnResultShift, st::KnnResultState, i) = @inbounds res.dist[i+st.shift]

@inline Base.last(res::KnnResultShift, st::KnnResultState) = last(res.id) => last(res.dist)
@inline Base.first(res::KnnResultShift, st::KnnResultState) = res.id[st.shift+1] => res.dist[st.shift+1]
@inline Base.maximum(res::KnnResultShift, st::KnnResultState) = last(res.dist)
@inline Base.minimum(res::KnnResultShift, st::KnnResultState) = res.dist[1+st.shift]
@inline Base.argmax(res::KnnResultShift, st::KnnResultState) = last(res.id)
@inline Base.argmin(res::KnnResultShift, st::KnnResultState) = res.id[1+st.shift]

Base.maximum(res::KnnResultShift) = last(res.dist)
Base.argmax(res::KnnResultShift) = last(res.id)
Base.minimum(res::KnnResultShift) = res.dist[_find_start_position(res)]
Base.argmin(res::KnnResultShift) = res.id[_find_start_position(res)]

@inline idview(res::KnnResultShift, st::KnnResultState) = @view res.id[st.shift+1:end]
@inline distview(res::KnnResultShift, st::KnnResultState) = @view res.dist[st.shift+1:end]

@inline Base.eachindex(res::KnnResultShift, st::KnnResultState) = 1:length(res, st)
Base.eltype(res::KnnResultShift) = Pair{Int32,Float32}

##### iterator interface
### KnnResultShift
"""
    Base.iterate(res::KnnResultShift, state::Int=1)

Support for iteration
"""
@inline function _find_start_position(res::KnnResultShift)
    i = 1
    id = res.id
    n = length(id)
    @inbounds while i <= n && id[i] == 0
        i += 1
    end

    i
end

function Base.iterate(res::KnnResultShift, i::Int=-1)
    n = length(res.id)
    n == 0 && return nothing
    if i == -1
        i = _find_start_position(res)
    end

    i > n && return nothing
    @inbounds res.id[i] => res.dist[i], i+1
end

max_stars_repo_name: sadit/NNS

max_issues_repo_name: sadit/NNS

max_forks_repo_name: sadit/NNS


text:
# SPDX-License-Identifier: MIT

#######################################################################
# Model definition

"""
    struct PPCA{T,D,Q}
        hprior  # Prior over the latent embeddings
        trans   # Affine transform
        λ       # Precision parameter
    end

Standard PPCA model.
"""
struct PPCA{T,D,Q}
    hprior::Normal{T,Q}
    trans::AffineTransform{T,D,Q}
    λ::BayesParam{Gamma{T}}
end

function PPCA(T::Type{<:AbstractFloat}; datadim, latentdim,
                   pstrength = 1e-3, W_MAP = false)

    D, Q = datadim, latentdim
    hprior = Normal(zeros(T, Q), Symmetric(Matrix{T}(I, Q, Q)))

    trans = AffineTransform(T, outputdim = D, inputdim = Q,
                            pstrength = pstrength, W_MAP = W_MAP)
    λprior = Gamma{T}(pstrength, pstrength)
    λposterior = Gamma{T}(pstrength, pstrength)
    λ = BayesParam(λprior, λposterior)

    PPCA{T,datadim,latentdim}(hprior, trans, λ)
end

PPCA(;datadim, latentdim, pstrength = 1e-3,
          W_MAP = false) = PPCA(Float64;
                                         datadim = datadim,
                                         latentdim = latentdim,
                                         pstrength = pstrength,
                                         W_MAP = W_MAP)

#######################################################################
# Estimate the latent variables

function (m::PPCA{T,D,Q})(X::AbstractVector) where {T,D,Q}
    S₁, S₂ = hstats(m.trans, X)
    λ̄ = mean(m.λ.posterior)
    S₁, S₂ = λ̄*S₁, λ̄*S₂

    Λ₀ = inv(m.hprior.Σ)
    Λ₀μ₀ = Λ₀ * m.hprior.μ
    Σ = Symmetric(inv(Λ₀ + S₂))
    [Normal(Σ * (Λ₀μ₀ + mᵢ), Σ) for mᵢ in S₁]
end

#######################################################################
# Pretty print

function Base.show(io::IO, ::MIME"text/plain", model::PPCA)
    println(io, typeof(model), ":")
    println(io, "  hprior: $(typeof(model.hprior))")
    println(io, "  trans: $(typeof(model.trans))")
    println(io, "  λ: $(typeof(model.λ))")
end

#######################################################################
# Log-likelihood

function _llh_d(::PPCA, x, Tŵ, Tλ, Th)
    λ, lnλ = Tλ
    h, hhᵀ = Th

    # Extract the bias parameter
    wwᵀ = Tŵ[2][1:end-1, 1:end-1]
    w = Tŵ[1][1:end-1]
    μ = Tŵ[1][end]
    μ² = Tŵ[2][end, end]

    x̄ = dot(w, h) + μ
    lognorm = (-log(2π) + lnλ - λ*(x^2))/2
    K = λ*(x̄*x - dot(w, h)*μ - (dot(vec(hhᵀ), vec(wwᵀ)) + μ²)/2)

    lognorm + K
end

function _llh(m::PPCA, x, Tŵs, Tλ, Th)
    f = (a,b) -> begin
        xᵢ, Tŵᵢ = b
        a + _llh_d(m, xᵢ, Tŵᵢ, Tλ, Th)
    end
    foldl(f, zip(x, Tŵs), init = 0)
end

function loglikelihood(m::PPCA{T,D,Q}, X) where {T,D,Q}
    hposts = X |> m
    _llh.(
        [m],
        X,
        [[gradlognorm(w.posterior, vectorize = false) for w in m.trans.W ]],
        [gradlognorm(m.λ.posterior, vectorize = false)],
        [gradlognorm(p, vectorize = false) for p in hposts]
   ) - kldiv.(hposts, [m.hprior])
end

#######################################################################
# Update of the precision parameter λ

function _λstats_d(::PPCA, x::Real, Tŵ, Th)
    h, hhᵀ = Th

    # Extract the bias parameter
    w = Tŵ[1][1:end-1]
    wwᵀ = Tŵ[2][1:end-1, 1:end-1]
    μ = Tŵ[1][end]
    μ² = Tŵ[2][end, end]

    x̄ = dot(w, h) + μ
    Tλ₁ = -.5*x^2 + x*x̄ - dot(w, h)*μ - .5*(dot(vec(hhᵀ), vec(wwᵀ)) + μ²)
    Tλ₂ = 1/2
    vcat(Tλ₁, Tλ₂)
end

function _λstats(m::PPCA, x::AbstractVector, Tŵ, Th)
    sum(_λstats_d.([m], x, Tŵ, [Th]))
end

function λstats(m::PPCA, X, hposts)
    Tŵ = [gradlognorm(w.posterior, vectorize = false) for w in m.trans.W]
    Th = [gradlognorm(p, vectorize = false) for p in hposts]
    sum(_λstats.([m], X, [Tŵ], Th))
end

function update_λ!(m::PPCA, accstats)::Nothing
    η₀ = naturalparam(m.λ.prior)
    update!(m.λ.posterior, η₀ + accstats)
    nothing
end

#######################################################################
# Update of the bases W

function wstats(m::PPCA, X, hposts)
    λ, _ = gradlognorm(m.λ.posterior, vectorize = false)
    S₁, S₂ = wstats(m.trans, X, hposts)
    [λ*S₁, λ*S₂]
end

update_W!(m::PPCA, accstats) = update_W!(m.trans, accstats)

#######################################################################
# Training

"""
    fit!(model, dataloader, [, epochs = 1, callback = x -> x])

Fit a PPCA model to a data set by estimating the variational posteriors
over the parameters.
"""
function fit!(model::PPCA, dataloader; epochs = 1, callback = x -> x)

    @everywhere dataloader = $dataloader

    # NOTE: By 1 epoch we mean TWO passes over the data, one pass to
    # update the bases and the other to update the precision parameter

    for e in 1:epochs
        # Propagate the model to all the workers
        @everywhere model = $model

        ###############################################################
        # Step 1: update the posterior of the bases
        waccstats = @distributed (+) for X in dataloader
            # E-step: estimate the posterior of the embeddings
            hposts = X |> model

            # Accumulate statistics for the bases w
            wstats(model, X, hposts)
        end
        update_W!(model, waccstats)

        # Propagate the model to all the workers
        @everywhere model = $model

        ###############################################################
        # Step 2: update the posterior of the precision λ
        λaccstats = @distributed (+) for X in dataloader
            # E-step: estimate the posterior of the embeddings
            hposts = X |> model

            # Accumulate statistics for λ
            λstats(model, X, hposts)
        end

        # M-step 2: update the posterior of the precision parameter λ
        update_λ!(model, λaccstats)

        # Notify the caller
        callback(e)
    end
end

max_stars_repo_name: BUTSpeechFIT/PPCA

max_issues_repo_name: BUTSpeechFIT/PPCA

max_forks_repo_name: BUTSpeechFIT/PPCA


text:
export step!, dummystep

"""
    step!(model, agent_step! [, model_step!], n::Integer = 1)

Update agents `n` steps. Agents will be updated as specified by the `model.scheduler`.
If given the optional function `model_step!`, it is triggered _after_ every scheduled
agent has acted.

    step!(model, agent_step! [, model_step!], n::Function)

`n` can be also be a function.
Then `step!` runs the model until `n(model, s)` returns `true`, where `s` is the
current amount of steps taken (starting from 0).
"""
function step! end

"""
    dummystep(model)

Ignore the model dynamics. Use instead of `model_step!`.
"""
dummystep(model) = nothing
"""
    dummystep(agent, model)

Ignore the agent dynamics. Use instead of `agent_step!`.
"""
dummystep(agent, model) = nothing

until(ss, n::Int, model) = ss < n
until(ss, n, model) = !n(model, ss)

step!(model::ABM, agent_step!, n = 1) = step!(model, agent_step!, dummystep, n)

function step!(model::ABM, agent_step!, model_step!, n)
  s = 0
  while until(s, n, model)
    activation_order = model.scheduler(model)
    for index in activation_order
      haskey(model.agents, index) || continue
      agent_step!(model.agents[index], model)
    end
    model_step!(model)
    s += 1
  end
end

function step!(model, agent_step!, model_step!, n; kwargs...)
  @warn "`step!` with keyword arguments is deprecated. Use `run!` instead."
  run!(model, agent_step!, model_step!, n; kwargs...)
end


text:
# SECTION 3: other visualization in BENCHMARK & fertility experiments
# ----------------------------------------------
# 3.1 benchmark: correlations (during plotting years)
   println(
      "Cor between Aging Population Ratio & gap/pool expenditures: ",
      Statistics.cor( dict_Demog_base["AgingPopuRatio"][idx_plot], (DatPkg_base.Dt[:LI] ./ DatPkg_base.Dt[:AggPoolExp])[idx_plot] )
   )
   println(
      "Cor between Aging Population Ratio & gap/pool incomes: ",
      Statistics.cor( dict_Demog_base["AgingPopuRatio"][idx_plot], (DatPkg_base.Dt[:LI] ./ DatPkg_base.Dt[:AggPoolIn])[idx_plot] )
   )
   println(
      "Cor between Aging Population Ratio & gap/GDP: ",
      Statistics.cor( dict_Demog_base["AgingPopuRatio"][idx_plot], (DatPkg_base.Dt[:LI] ./ DatPkg_base.Dt[:Y])[idx_plot] )
   )
   println(
      "Cor between Aging Population Ratio & gap/tax revenues: ",
      Statistics.cor( dict_Demog_base["AgingPopuRatio"][idx_plot], (DatPkg_base.Dt[:LI] ./ ( DatPkg_base.Dt[:TRw] .+ DatPkg_base.Dt[:TRc] ) )[idx_plot] )
   )

#


text:
using Statistics
using NCDatasets

using ArgParse
using JSON

include("nanop.jl")
include("PCA.jl")

using .PCA

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--output-file"
            help = "Output file."
            arg_type = String
            required = true

        "--data-file"
            help = "PSLA data file."
            arg_type = String
            required = true

        "--domain-file"
            help = "Domain file."
            arg_type = String
            required = true

        "--sparsity"
            help = "Because solving PCA matrix is time consuming. This parameter let user to use coarse data to derive PCA. Sparseity `n` means to skip every `n` other grid point in lon / lat. So sparsity 0 (default) means do not skip any data point. Sparsity 1 means means data density will drop to 1/4 of the original (1/2 * 1/2)."
            arg_type = Int64
            default = 0

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

sparsity = parsed["sparsity"]
inc = sparsity + 1

Dataset(parsed["domain-file"], "r") do ds
    global lon = replace(ds["xc"][:, 1], missing=>NaN)
    global lat = replace(ds["yc"][1, :], missing=>NaN)

    global Nx = length(lon)
    global Ny = length(lat)

    global valid_grids = []
    global share_grids = []

    # detect the first lat
    global southest_j = 1
    for j = 1:Ny
        if 20.0 <= lat[j]
            global southest_j = j
            break
        end
    end

    for i = 1:inc:Nx, j = southest_j:inc:Ny
        push!(valid_grids, (i, j))
    end
end

Dataset(parsed["data-file"], "r") do ds

    global PSLA  = replace(ds["PSL_MA"][:, :, 1, :], missing=>NaN)
    global months  = ds.dim["time"]

    global EOF_input = zeros(Float64, length(valid_grids), months)
    for i = 1:length(valid_grids)
        EOF_input[i, :] = PSLA[valid_grids[i][1], valid_grids[i][2], :]
    end
end

modes = 2

println("Solving for PCA...")
eigen_vectors = PCA.findPCAs(EOF_input, num=modes)
println("done.")

PCAs = zeros(Float64, Nx, Ny, modes)
PCAs .= NaN

for i = 1:length(valid_grids)
    PCAs[valid_grids[i][1], valid_grids[i][2], :] = eigen_vectors[i, :]
end

# Map grid back because of sparsity
if sparsity > 0
    for m = 1:modes, i = 1:inc:Nx, j=southest_j:inc:Ny
        PCAs[i:min(Nx, i+sparsity), j:min(Ny, j+sparsity), m] .= PCAs[i, j, m]
    end
end

AO = zeros(Float64, months)

AO_mode = view(PCAs, :, :, 1)
for i = 1:length(AO)
    AO[i] = nansum(AO_mode .* view(PSLA, :, :, i))
end

Dataset(parsed["output-file"], "c") do ds

    defDim(ds, "Nx", Nx)
    defDim(ds, "Ny", Ny)
    defDim(ds, "time", Inf)
    defDim(ds, "modes", modes)

    for (varname, vardata, vardim, attrib) in [
        ("PCAs",  PCAs, ("Nx", "Ny", "modes",), Dict()),
        ("AO",   AO, ("time",), Dict()),
    ]
        println("Doing varname:", varname)
        var = defVar(ds, varname, Float64, vardim)
        var.attrib["_FillValue"] = 1e20

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
function test_func(x)
    max = NighttimeLights.max_missing(x)
    return x.+ max
end

@testset "cross_apply" begin
    for i in 1:10
        x = rand(10:30)
        y = rand(10:30)
        z = rand(10:30)
        rad = rand(20:100.0, x, y, z)
        @test size(cross_apply(test_func, rad)) == size(rad)
    end
end

@testset "long_apply" begin
    for i in 1:10
        x = rand(10:30)
        y = rand(10:30)
        z = rand(10:30)
        rad = rand(20:100.0, x, y, z)
        mask = rand(0:1, x, y)
        @test size(long_apply(test_func, rad, mask)) == size(rad)
    end
end


text:
module MeltingParam

# If you want to add a new method here, feel free to do so.
# Remember to also export the function name in GeoParams.jl (in addition to here)

using Parameters, LaTeXStrings, Unitful
using ..Units
using GeoParams:
    AbstractMaterialParam, PhaseDiagram_LookupTable, AbstractMaterialParamsStruct
import Base.show, GeoParams.param_info
using ..MaterialParameters: MaterialParamsInfo
using Setfield # allows modifying fields in immutable struct

abstract type AbstractMeltingParam{T} <: AbstractMaterialParam end

export  compute_meltfraction,
        compute_meltfraction!,    # calculation routines
        compute_dϕdT,             # derivative of melt fraction versus
        compute_dϕdT!,
        param_info,
        MeltingParam_Caricchi,
        MeltingParam_4thOrder,
        MeltingParam_5thOrder,
        MeltingParam_Quadratic,
        MeltingParam_Assimilation,
        SmoothMelting

include("../Utils.jl")
include("../Computations.jl")

# Caricchi  -------------------------------------------------------
"""
    MeltingParam_Caricchi()

Implements the T-dependent melting parameterisation used by Caricchi, Simpson et al. (as for example described in Simpson)
```math
    \\theta = {(a - (T + c)) \\over b}
```
```math
    \\phi_{melt} = {1.0 \\over (1.0 + e^\\theta)}
```

Note that T is in Kelvin. As default parameters we employ:
```math
b=23K,  a=800K , c=273.15K
```
Which gives a reasonable fit to experimental data of granodioritic composition (Piwinskii and Wyllie, 1968):

![MeltParam_Carrichi](./assets/img/MeltingParam_Caricchi.png)

References
====
- Simpson G. (2017) Practical finite element modelling in Earth Sciences Using MATLAB.
"""
@with_kw_noshow struct MeltingParam_Caricchi{T,U} <: AbstractMeltingParam{T}
    a::GeoUnit{T,U}     = 800.0K
    b::GeoUnit{T,U}     = 23.0K
    c::GeoUnit{T,U}     = 273.15K # shift from C to K
    apply_bounds::Bool  = true
end
MeltingParam_Caricchi(args...) = MeltingParam_Caricchi(convert.(GeoUnit, args)...)

function param_info(s::MeltingParam_Caricchi) # info about the struct
    return MaterialParamsInfo(;
        Equation=L"\phi = {1 \over {1 + \exp( {800-T[^oC] \over 23})}}"
    )
end

# Calculation routines
function (p::MeltingParam_Caricchi)(; T, kwargs...)
    if T isa Quantity
        @unpack_units a, b, c = p
    else
        @unpack_val a, b, c = p
    end
    θ = (a - (T - c)) / b
    return 1.0 / (1.0 + exp(θ))
end

function compute_meltfraction!(ϕ::AbstractArray, p::MeltingParam_Caricchi; T, kwargs...)
    for i in eachindex(T)
        @inbounds ϕ[i] = p(; T=T[i])
    end
    return nothing
end

function compute_dϕdT(p::MeltingParam_Caricchi; T, kwargs...)
    if T isa Quantity
        @unpack_units a, b, c = p
    else
        @unpack_val a, b, c = p
    end

    dϕdT = exp((a + c - T) / b) / (b * ((1.0 + exp((a + c - T) / b))^2))

    return dϕdT
end

function compute_dϕdT!(dϕdT::AbstractArray, p::MeltingParam_Caricchi; T, kwargs...)
    for i in eachindex(T)
        @inbounds dϕdT[i] = compute_dϕdT(p, T=T[i])
    end

    return nothing
end

# Print info
function show(io::IO, g::MeltingParam_Caricchi)
    return print(io, "Caricchi et al. melting parameterization")
end
#-------------------------------------------------------------------------

# MeltingParam_5thOrder  -------------------------------------------------------
"""
    MeltingParam_5thOrder(a,b,c,d,e,f,T_s,T_l)

Uses a 5th order polynomial to describe the melt fraction `phi` between solidus temperature `T_s` and liquidus temperature `T_l`
```math
    \\phi = a T^5 + b T^4 + c T^3 + d T^2 + e T + f  \\textrm{   for   } T_s ≤ T ≤ T_l
```
```math
    \\phi = 1  \\textrm{   if   } T>T_l
```
```math
    \\phi = 0  \\textrm{   if   } T<T_s
```
Temperature `T` is in Kelvin.

![MeltingParam_5thOrder](./assets/img/MeltingParam_5thorder.png)

The default values are for a composite liquid-line-of-descent:
- the upper part is for Andesite from: (Blatter, D. L. & Carmichael, I. S. (2001) Hydrous phase equilibria of a Mexican highsilica andesite: a candidate for a mantle origin? Geochim. Cosmochim. Acta 65, 4043–4065
- the lower part is extrapolated to the granitic minimum using the Marxer & Ulmer LLD for Andesite (Marxer, F. & Ulmer, P. (2019) Crystallisation and zircon saturation of calc-alkaline tonalite from the Adamello Batholith at upper crustal conditions: an experimental study. *Contributions Mineral. Petrol.* 174, 84)

"""
@with_kw_noshow struct MeltingParam_5thOrder{T,U,U1,U2,U3,U4,U5,U6} <:
                       AbstractMeltingParam{T}
    a::GeoUnit{T,U1} = 2.083291971482524e-12 / K^5
    b::GeoUnit{T,U2} = -1.239502833666574e-08 / K^4
    c::GeoUnit{T,U3} = 2.938887604687626e-05 / K^3
    d::GeoUnit{T,U4} = -0.034711533077108 / K^2
    e::GeoUnit{T,U5} = 20.425403874539178 / K
    f::GeoUnit{T,U6} = -4.790664658179178e+03 * NoUnits
    T_s::GeoUnit{T,U} = 963.15K
    T_l::GeoUnit{T,U} = 1388.2K
    apply_bounds::Bool  = true
end
MeltingParam_5thOrder(args...) = MeltingParam_5thOrder(convert.(GeoUnit, args)...)

function param_info(s::MeltingParam_5thOrder) # info about the struct
    return MaterialParamsInfo(; Equation=L"\phi = aT^5 + bT^4 + cT^3 + dT^2 + eT + f")
end

# Calculation routines
function (p::MeltingParam_5thOrder)(; T, kwargs...)
    if T isa Quantity
        @unpack_units a, b, c, d, e, f, T_s, T_l  = p
    else
        @unpack_val a, b, c, d, e, f, T_s, T_l = p
    end

    ϕ = a * T^5 + b * T^4 + c * T^3 + d * T^2 + e * T + f
    if p.apply_bounds
        if T < T_s
            ϕ = 0.0
        elseif T > T_l
            ϕ = 1.0
        end
    end

    return ϕ
end

function compute_meltfraction!(
    ϕ::AbstractArray, p::MeltingParam_5thOrder; T::AbstractArray, kwargs...
)
    for i in eachindex(T)
        @inbounds ϕ[i] = p(; T=T[i])
    end

    return nothing
end

function compute_dϕdT(p::MeltingParam_5thOrder; T::Real, kwargs...)
    if T isa Quantity
        @unpack_units a, b, c, d, e, T_s, T_l = p
    else
        @unpack_val a, b, c, d, e, T_s, T_l = p
    end

    dϕdT = 5 * a * T^4 + 4 * b * T^3 + 3 * c * T^2 + 2 * d * T + e
    if (T < T_s || T > T_l) && p.apply_bounds
        dϕdT = 0.0
    end

    return dϕdT
end

function compute_dϕdT!(
    dϕdT::AbstractArray, p::MeltingParam_5thOrder; T::AbstractArray, kwargs...
)

    for i in eachindex(T)
        @inbounds dϕdT[i] = compute_dϕdT(p, T=T[i])
    end

    return nothing
end

# Print info
function show(io::IO, g::MeltingParam_5thOrder)
    return print(
        io,
        "5th order polynomial melting curve: phi = $(NumValue(g.a)) T^5 + $(NumValue(g.b))T^4 + $(NumValue(g.c))T^3 + $(NumValue(g.d))T^2 + $(NumValue(g.e))T + $(NumValue(g.f))  $(Value(g.T_s)) ≤ T ≤ $(Value(g.T_l))",
    )
end
#-------------------------------------------------------------------------

# MeltingParam_4thOrder  -------------------------------------------------------
"""
    MeltingParam_4thOrder(b,c,d,e,f,T_s,T_l)

Uses a 4th order polynomial to describe the melt fraction `phi` between solidus temperature `T_s` and liquidus temperature `T_l`
```math
    \\phi = b T^4 + c T^3 + d T^2 + e T + f  \\textrm{   for   } T_s ≤ T ≤ T_l
```
```math
    \\phi = 1 \\textrm{   if   } T>T_l
```
```math
    \\phi = 0  \\textrm{   if   } T<T_s
```
Temperature `T` is in Kelvin.

![MeltingParam_4thOrder](./assets/img/MeltingParam_4thorder.png)

The default values are for Tonalite experiments from Marxer and Ulmer (2019):
- Marxer, F. & Ulmer, P. (2019) Crystallisation and zircon saturation of calc-alkaline tonalite from the Adamello Batholith at upper crustal conditions: an experimental study. *Contributions Mineral. Petrol.* 174, 84

"""
@with_kw_noshow struct MeltingParam_4thOrder{T,U,U2,U3,U4,U5,U6} <: AbstractMeltingParam{T}
    b::GeoUnit{T,U2} = -7.594512597174117e-10 / K^4
    c::GeoUnit{T,U3} = 3.469192091489447e-06 / K^3
    d::GeoUnit{T,U4} = -0.005923529809260 / K^2
    e::GeoUnit{T,U5} = 4.482855645604745 / K
    f::GeoUnit{T,U6} = -1.268730161921053e+03 * NoUnits
    T_s::GeoUnit{T,U} = 963.15K
    T_l::GeoUnit{T,U} = 1270.15K
    apply_bounds::Bool  = true
end
MeltingParam_4thOrder(args...) = MeltingParam_4thOrder(convert.(GeoUnit, args)...)

function param_info(s::MeltingParam_4thOrder) # info about the struct
    return MaterialParamsInfo(; Equation=L"\phi = bT^4 + cT^3 + dT^2 + eT + f")
end

# Calculation routines
function (p::MeltingParam_4thOrder)(; T::Real, kwargs...)
    if T isa Quantity
        @unpack_units b, c, d, e, f, T_s, T_l = p
    else
        @unpack_val  b, c, d, e, f, T_s, T_l = p
    end

    ϕ = b * T^4 + c * T^3 + d * T^2 + e * T + f
    if p.apply_bounds
        if T < T_s
            ϕ = 0.0
        elseif T > T_l
            ϕ = 1.0
        end
    end

    return ϕ
end

function compute_meltfraction!(
    ϕ::AbstractArray, p::MeltingParam_4thOrder; T::AbstractArray, kwargs...
)
    for i in eachindex(T)
        @inbounds ϕ[i] = p(; T=T[i])
    end

    return nothing
end

function compute_dϕdT(p::MeltingParam_4thOrder; T::Real, kwargs...)
    if T isa Quantity
        @unpack_units  b, c, d, e, T_s, T_l = p
    else
        @unpack_val b, c, d, e, T_s, T_l = p
    end

    dϕdT = 4 * b * T^3 + 3 * c * T^2 + 2 * d * T + e
    if (T < T_s || T > T_l) && p.apply_bounds
        dϕdT = 0.0
    end
    return dϕdT
end

function compute_dϕdT!(
    dϕdT::AbstractArray, p::MeltingParam_4thOrder; T::AbstractArray, kwargs...
)

    for i in eachindex(T)
        @inbounds dϕdT[i] = compute_dϕdT(p, T=T[i])
    end

    return nothing
end

# Print info
function show(io::IO, g::MeltingParam_4thOrder)
    return print(
        io,
        "4th order polynomial melting curve: phi = $(NumValue(g.b))T^4 + $(NumValue(g.c))T^3 + $(NumValue(g.d))T^2 + $(NumValue(g.e))T + $(NumValue(g.f))  $(Value(g.T_s)) ≤ T ≤ $(Value(g.T_l))",
    )
end
#-------------------------------------------------------------------------

# MeltingParam_Quadratic  -------------------------------------------------------
"""
    MeltingParam_Quadratic(T_s,T_l)

Quadratic melt fraction parameterisation where melt fraction ``\\phi`` depends only on solidus (``T_s``) and liquidus (``T_l``) temperature:
```math
    \\phi = 1.0 - \\left( {T_l - T} \\over {T_l - T_s} \\right)^2
```
```math
    \\phi = 1.0 \\textrm{ if } T>T_l
```
```math
    \\phi = 0.0 \\textrm{ if } T<T_s
```
Temperature `T` is in Kelvin.

![MeltingParam_Quadratic](./assets/img/MeltingParam_Quadratic.png)

This was used, among others, in Tierney et al. (2016) Geology

"""
@with_kw_noshow struct MeltingParam_Quadratic{T,U} <: AbstractMeltingParam{T}
    T_s::GeoUnit{T,U}   = 963.15K
    T_l::GeoUnit{T,U}   = 1273.15K
    apply_bounds::Bool  = true
end
MeltingParam_Quadratic(args...) = MeltingParam_Quadratic(convert.(GeoUnit, args)...)

function param_info(s::MeltingParam_Quadratic) # info about the struct
    return MaterialParamsInfo(; Equation=L"\phi = 1.0 - ((T_l - T)/(T_l - T_s))^2")
end

# Calculation routines
function (p::MeltingParam_Quadratic)(; T, kwargs...)
    if T isa Quantity
        @unpack_units T_s, T_l = p
    else
        @unpack_val T_s, T_l = p
    end

    ϕ = 1.0 - ((T_l - T) / (T_l - T_s))^2
    if p.apply_bounds
        if T > T_l
            ϕ = 1.0
        elseif T < T_s
            ϕ = 0.0
        end
    end
    return ϕ
end

function compute_meltfraction!(
    ϕ::AbstractArray, p::MeltingParam_Quadratic; T::AbstractArray, kwargs...
)

    for i in eachindex(T)
        @inbounds ϕ[i] = p(; T=T[i])
    end

    return nothing
end

function compute_dϕdT(p::MeltingParam_Quadratic; T, kwargs...)
    if T isa Quantity
        @unpack_units T_s, T_l = p
    else
        @unpack_val T_s, T_l = p
    end

    dϕdT = (2T_l - 2T) / ((T_l - T_s)^2)
    if (T > T_l || T < T_s) && p.apply_bounds
        dϕdT = 0.0
    end
    return dϕdT
end

function compute_dϕdT!(
    dϕdT::AbstractArray, p::MeltingParam_Quadratic; T::AbstractArray, kwargs...
)

    for i in eachindex(T)
        @inbounds dϕdT[i] = compute_dϕdT(p, T=T[i])
    end

    return nothing
end

# Print info
function show(io::IO, g::MeltingParam_Quadratic)
    return print(
        io,
        "Quadratic melting curve:  ϕ = 1.0 - ((Tₗ-T)/(Tₗ-Tₛ))² with Tₛ=$(Value(g.T_s)), Tₗ=$(Value(g.T_l)) ",
    )
end
#-------------------------------------------------------------------------

# MeltingParam_Assimilation  -------------------------------------------------------
"""
    MeltingParam_Assimilation(T_s,T_l,a)

Melt fraction parameterisation that takes the assimilation of crustal host rocks into account, as used by Tierney et al. (2016) based upon a parameterisation of Spera and Bohrson (2001)

Here, the fraction of molten and assimilated host rocks ``\\phi`` depends on the solidus (``T_s``) and liquidus (``T_l``) temperatures of the rocks, as well as on a parameter ``a=0.005``
```math
    X = \\left( {T - T_s} \\over {T_l - T_s} \\right)
```
```math
    \\phi = a \\cdot \\left( \\exp^{2ln(100)X} - 1.0 \\right) \\textrm{ if } X ≤ 0.5
```
```math
    \\phi = 1- a \\cdot \\exp^{2ln(100)(1-X)}  \\textrm{ if } X > 0.5
```
```math
    \\phi = 1.0 \\textrm{ if } T>T_l
```
```math
    \\phi = 0.0 \\textrm{ if } T<T_s
```
Temperature `T` is in Kelvin.

![MeltingParam_Assimilation](./assets/img/MeltingParam_Assimilation.png)

This was used, among others, in Tierney et al. (2016), who employed as default parameters:
```math
   T_s=973.15, T_l=1173.15, a=0.005
```

References
==========
- Spera, F.J., and Bohrson, W.A., 2001, Energy-Constrained Open-System Magmatic Processes I: General Model and Energy-Constrained Assimilation and Fractional Crystallization (EC- AFC) Formulation: Journal of Petrology, v. 42, p. 999–1018.
- Tierney, C.R., Schmitt, A.K., Lovera, O.M., de Silva, S.L., 2016. Voluminous plutonism during volcanic quiescence revealed by thermochemical modeling of zircon. Geology 44, 683–686. https://doi.org/10.1130/G37968.1

"""
@with_kw_noshow struct MeltingParam_Assimilation{T,U,U1} <: AbstractMeltingParam{T}
    T_s::GeoUnit{T,U}   =   973.15K
    T_l::GeoUnit{T,U}   =   1173.15K
    a::GeoUnit{T,U1}    =   0.005NoUnits
    apply_bounds::Bool  =   true
end
MeltingParam_Assimilation(args...) = MeltingParam_Assimilation(convert.(GeoUnit,args)...)

function param_info(s::MeltingParam_Assimilation) # info about the struct
    return MaterialParamsInfo(Equation =  L"\phi = f(T_S,T_l,a) taking crustal assimilation into account.")
end

# Calculation routines
function (p::MeltingParam_Assimilation)(; T::Real, kwargs...)
    if T isa Quantity
        @unpack_units T_s,T_l,a   = p
    else
        @unpack_val T_s,T_l,a   = p
    end

    X = (T - T_s)/(T_l - T_s)

    if X <= 0.5
        ϕ   =    a * (exp(2*log(100)*X) - 1.0 )
    else
        ϕ   =    1.0 - a * exp(2*log(100)*(1-X) )
    end
    if p.apply_bounds
        if T>T_l
            ϕ = 1.0
        elseif T<T_s
            ϕ = 0.0
        end
    end
    return ϕ
end

function compute_meltfraction!(
    ϕ::AbstractArray, p::MeltingParam_Assimilation; T::AbstractArray, kwargs...
)

    for i in eachindex(T)
        @inbounds ϕ[i] = p(; T=T[i])
    end

    return nothing
end

function compute_dϕdT(p::MeltingParam_Assimilation; T::Real, kwargs...)
    if T isa Quantity
        @unpack_units T_s,T_l,a   = p
    else
        @unpack_val T_s,T_l,a   = p
    end

    X      =   (T - T_s)/(T_l - T_s)
    dϕdT   =   (9.210340371976184*a*exp((9.210340371976184*T - 9.210340371976184*T_s) / (T_l - T_s))) / (T_l - T_s)
    if X>0.5
        dϕdT   = (9.210340371976184*a*exp(9.210340371976184 + (9.210340371976184*T_s - 9.210340371976184*T) / (T_l - T_s))) / (T_l - T_s)
    end

    if (T>T_l || T<T_s) && p.apply_bounds
        dϕdT = 0.0
    end
    return dϕdT
end

function compute_dϕdT!(
    dϕdT::AbstractArray, p::MeltingParam_Assimilation; T::AbstractArray, kwargs...
)

    for i in eachindex(T)
        @inbounds dϕdT[i] = compute_dϕdT(p, T=T[i])
    end

    return nothing
end

# Print info
function show(io::IO, g::MeltingParam_Assimilation)
    return print(
        io,
        "Quadratic melting assimilation parameterisation after Spera & Bohrson (2001)",
    )
end

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

# Smooth melting function ------------------------------------------------

"""
    SmoothMelting(; p=MeltingParam_4thOrder(), k_sol=0.2/K,  k_liq=0.2/K)

This smoothens the melting parameterisation ``p`` around the solidus ``T_{sol}`` and liquidus ``T_{liq}``
using a smoothened Heaviside step functions for the solidus:

```math
    H_{sol} =  {1.0 \\over { 1 + \\exp( -2 k_{sol} (T - T_{sol} - {2 \\over k_{sol}}) )  }}
```
and liquidus:
```math
    H_{liq} =  1.0 - {1.0 \\over { 1 + \\exp( -2 k_{liq} (T - T_{liq} + {2 \\over k_{liq}}) )  }}
```
The resulting melt fraction ``\\phi`` is computed from the original melt fraction ``\\phi_0`` (computed using one of the methods above) as:
```math
    \\phi =  \\phi_0 H_{sol} H_{liq} + 1.0 - H_{liq}
```
The width of the smoothening zones is controlled by ``k_{sol}, k_{liq}`` (larger values = sharper boundary).

This is important, as jumps in the derivative ``dϕ/dT`` can cause numerical instabilities in latent heat computations, which is prevented with this smoothening.

Example
====

Let's consider a 4th order parameterisation:
```julia
julia> using GeoParams, Plots
julia> p = MeltingParam_4thOrder();
julia> T= collect(650.0:1:1050.) .+ 273.15;
julia> T,phi,dϕdT =  PlotMeltFraction(p,T=T);
```

The same but with smoothening:
```julia
julia> p_s = SmoothMelting(p=MeltingParam_4thOrder(), k_liq=0.21/K);
4th order polynomial melting curve: phi = -7.594512597174117e-10T^4 + 3.469192091489447e-6T^3 + -0.00592352980926T^2 + 4.482855645604745T + -1268.730161921053  963.15 K ≤ T ≤ 1270.15 K with smooth Heaviside function smoothening using k_sol=0.1 K⁻¹·⁰, k_liq=0.11 K⁻¹·⁰
julia> T_s,phi_s,dϕdT_s =  PlotMeltFraction(p_s,T=T);
```

We can create plots of this with:
```julia
julia> plt1 = plot(T.-273.15, phi, ylabel="Melt Fraction ϕ", color=:red, label="original", xlabel="Temperature [C]")
julia> plt1 = plot(plt1, T.-273.15, phi_s,  color=:black, label="smoothened", legend=:bottomright)
julia> plt2 = plot(T.-273.15, dϕdT, ylabel="dϕ/dT", color=:red, label="original", xlabel="Temperature [C]")
julia> plt2 = plot(plt2, T.-273.15, dϕdT_s,  color=:black, label="smoothened", legend=:topright)
julia> plot!(plt1,plt2,   xlabel="Temperature [C]", layout=(2,1))
```
The derivative no longer has a jump now:

![MeltingParam_Smooth](./assets/img/MeltingParam_Smooth.png)

"""
struct SmoothMelting{P,T,U} <: AbstractMeltingParam{T}
    p::P
    k_sol::GeoUnit{T,U}
    k_liq::GeoUnit{T,U}
end

# Set default values:
function SmoothMelting(; p=MeltingParam_4thOrder(), k_sol=0.2/K,  k_liq=0.2/K)
    k_sol = convert(GeoUnit,k_sol)
    k_liq = convert(GeoUnit,k_liq)
    p = @set p.apply_bounds = false
    SmoothMelting(p, k_sol, k_liq)
end

SmoothMelting(p::AbstractMeltingParam) =  SmoothMelting(p=p)
#SmoothMelting(p::AbstractMeltingParam; k_liq) =  SmoothMelting(p=p, k_liq=convert.(GeoUnit,k_liq))

# Calculation routines
function (param::SmoothMelting)(; T, kwargs...)
    if T isa Quantity
        @unpack_units k_sol,k_liq  = param
    else
        @unpack_val k_sol,k_liq  = param
        k_sol = param.k_sol.val
    end

    ϕ = param.p(; T, kwargs...)     # Melt fraction computed in usual manner

    T_s = param.p.T_s
    H_s = 1.0/( 1.0 + exp(-2*k_sol*(T-T_s-(2/k_sol))));

    T_l = param.p.T_l
    H_l = 1.0 - 1.0/( 1.0 + exp(-2*k_liq*(T-T_l+(2/k_liq))));

    # Apply heaviside smoothening above liquidus & below solidus
    ϕ = ϕ*H_s*H_l + 1.0 - H_l

    return ϕ
end

function compute_meltfraction!(
    ϕ::AbstractArray, param::SmoothMelting;  T::AbstractArray, kwargs...
)
    for i in eachindex(T)
        @inbounds ϕ[i] = param(; T=T[i], kwargs=kwargs)
    end

    return nothing
end

function compute_dϕdT(param::SmoothMelting; T, kwargs...)

    if T isa Quantity
        @unpack_units k_sol,k_liq  = param
    else
        @unpack_val k_sol,k_liq  = param
    end

    # compute heaviside functions & derivatives of that vs. T
    T_s = param.p.T_s
    T_l = param.p.T_l
    H_s = 1.0/( 1.0 + exp(-2*k_sol*(T-T_s-(2/k_sol))));
    H_l = 1.0 - 1.0/( 1.0 + exp(-2*k_liq*(T-T_l+(2/k_liq))));

    dHs_dT = 2k_sol*( 1.0 / ((1.0 + exp(-2k_sol*(T + -2 / k_sol - T_s)))^2))*exp(-2k_sol*(T + -2 / k_sol - T_s))
    dHl_dT = 2k_liq*(-1.0 / ((1.0 + exp(-2k_liq*(T +  2 / k_liq - T_l)))^2))*exp(-2k_liq*(T +  2 / k_liq - T_l))

    # melt fraction & derivative
    dϕdT = compute_dϕdT(param.p, T=T)
    ϕ    = param.p(; T, kwargs...)

    # The derivative of the function
    # ϕ = ϕ(T)*H_s(T)*H_l(T)  + 1.0 - H_l
    # versus T is

    dϕdT_tot = dϕdT*H_s*H_l + ϕ*dHs_dT*H_l + ϕ*H_s*dHl_dT - dHl_dT

    return dϕdT_tot
end

function compute_dϕdT!(
    dϕdT::AbstractArray, p::SmoothMelting; T::AbstractArray, kwargs...
)

    for i in eachindex(T)
        @inbounds dϕdT[i] = compute_dϕdT(p, T=T[i])
    end

    return nothing
end

# Print info
function show(io::IO, g::SmoothMelting)
    param = show(io,g.p);
    return print(
        io,
        " with smooth Heaviside function smoothening using k_sol=$(Value(g.k_sol)), k_liq=$(Value(g.k_liq))",
    )
end
#-------------------------------------------------------------------------

"""
    compute_meltfraction(P,T, p::AbstractPhaseDiagramsStruct)

Computes melt fraction in case we use a phase diagram lookup table. The table should have the column `:meltFrac` specified.
"""
function compute_meltfraction(
    p::PhaseDiagram_LookupTable; P::_T, T::_T, kwargs...
) where {_T}
    return p.meltFrac.(T, P)
end

compute_meltfraction(p::PhaseDiagram_LookupTable, args) = compute_meltfraction(p; args...)
"""
    compute_meltfraction!(ϕ::AbstractArray{<:AbstractFloat}, P::AbstractArray{<:AbstractFloat},T:AbstractArray{<:AbstractFloat}, p::PhaseDiagram_LookupTable)

In-place computation of melt fraction in case we use a phase diagram lookup table. The table should have the column `:meltFrac` specified.
"""
function compute_meltfraction!(
    ϕ::AbstractArray{_T},
    p::PhaseDiagram_LookupTable;
    P::AbstractArray{_T},
    T::AbstractArray{_T},
    kwargs...,
) where {_T}
    ϕ[:] = p.meltFrac.(T, P)

    return nothing
end

compute_meltfraction!( ϕ::AbstractArray, p::PhaseDiagram_LookupTable, args) = compute_meltfraction!(p; args...)

"""
    compute_dϕdT(P,T, p::AbstractPhaseDiagramsStruct)

Computes derivative of melt fraction vs T in case we use a phase diagram lookup table. The table should have the column `:meltFrac` specified.
The derivative is computed by finite differencing.
"""
function compute_dϕdT(p::PhaseDiagram_LookupTable; P::_T, T::_T, kwargs...) where {_T}
    dT = (maximum(T) - minimum(T)) / 2.0 * 1e-6 + 1e-6   # 1e-6 of the average T value
    ϕ1 = p.meltFrac.(T .+ dT, P)
    ϕ0 = p.meltFrac.(T, P)
    dϕdT = (ϕ1 - ϕ0) / dT

    return dϕdT
end

compute_dϕdT(p::PhaseDiagram_LookupTable, args) = compute_dϕdT(p; args...)

"""
    compute_dϕdT!(dϕdT::AbstractArray{<:AbstractFloat}, P::AbstractArray{<:AbstractFloat},T:AbstractArray{<:AbstractFloat}, p::PhaseDiagram_LookupTable)

In-place computation of melt fraction in case we use a phase diagram lookup table. The table should have the column `:meltFrac` specified.
The derivative is computed by finite differencing.
"""
function compute_dϕdT!(
    dϕdT::AbstractArray{_T},
    p::PhaseDiagram_LookupTable;
    P::AbstractArray{_T},
    T::AbstractArray{_T},
    kwargs...,
) where {_T}
    dT = (maximum(T) - minimum(T)) / 2.0 * 1e-6 + 1e-6   # 1e-6 of the average T value
    ϕ1 = p.meltFrac.(T .+ dT, P)
    ϕ0 = p.meltFrac.(T, P)
    dϕdT = (ϕ1 - ϕ0) / dT

    return nothing
end

compute_dϕdT!( ϕ::AbstractArray, p::PhaseDiagram_LookupTable, args) = compute_dϕdT!(p; args...)

# fill methods programatically
for myType in (
    :MeltingParam_Caricchi,
    :MeltingParam_5thOrder,
    :MeltingParam_4thOrder,
    :MeltingParam_Quadratic,
    :MeltingParam_Assimilation,
    :SmoothMelting,
)
    @eval begin
        (p::$(myType))(args) = p(; args...)
        compute_meltfraction(p::$(myType), args) = p(; args...)
        function compute_meltfraction!(ϕ::AbstractArray, p::$(myType), args)
            return compute_meltfraction!(ϕ, p; args...)
        end
        compute_dϕdT(p::$(myType), args) = compute_dϕdT(p; args...)
        function compute_dϕdT!(dϕdT::AbstractArray, p::$(myType), args)
            return compute_dϕdT!(dϕdT, p; args...)
        end
    end
end

# Computational routines needed for computations with the MaterialParams structure
function compute_meltfraction(s::AbstractMaterialParamsStruct, args)
    if isempty(s.Melting) #in case there is a phase with no melting parametrization
        return zero(typeof(args).types[1])
    else
        return compute_meltfraction(s.Melting[1], args)
    end
end

function compute_dϕdT(s::AbstractMaterialParamsStruct, args)
    if isempty(s.Melting) #in case there is a phase with no melting parametrization
        return zero(typeof(args).types[1])
    else
        return compute_dϕdT(s.Melting[1], args)
    end
end

"""
    ϕ = compute_meltfraction(Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})

Computation of melt fraction ϕ for the whole domain and all phases, in case an array with phase properties `MatParam` is provided, along with `P` and `T` arrays.
"""
compute_meltfraction(args...) = compute_param(compute_meltfraction, args...)

"""
    compute_meltfraction(ϕ::AbstractArray{<:AbstractFloat}, Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})

In-place computation of melt fraction ϕ for the whole domain and all phases, in case an array with phase properties `MatParam` is provided, along with `P` and `T` arrays.
"""
compute_meltfraction!(args...) = compute_param!(compute_meltfraction, args...)

"""
    ϕ = compute_dϕdT(Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})

Computates the derivative of melt fraction ϕ versus temperature `T` for the whole domain and all phases, in case an array with phase properties `MatParam` is provided, along with `P` and `T` arrays.
This is employed in computing latent heat terms in an implicit manner, for example
"""
compute_dϕdT(args...) = compute_param(compute_dϕdT, args...)

"""
    compute_dϕdT!(ϕ::AbstractArray{<:AbstractFloat}, Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})

Computates the derivative of melt fraction `ϕ` versus temperature `T`, ``{\\partial \\phi} \\over {\\partial T}`` for the whole domain and all phases, in case an array with phase properties `MatParam` is provided, along with `P` and `T` arrays.
This is employed, for example, in computing latent heat terms in an implicit manner.
"""
compute_dϕdT!(args...) = compute_param!(compute_dϕdT, args...)

end


text:
println()
println("Testing MathProgNLSModel")

@printf("%-15s  %4s  %4s  %4s  %10s  %10s  %10s\n",
        "Problem", "nequ", "nvar", "ncon", "‖F(x₀)‖²", "‖JᵀF‖",
        "‖c(x₀)‖")
# Test that every problem can be instantiated.
for prob in [:mgh07,:mgh35,:hs30,:hs43]
  prob_fn = eval(prob)
  nls = prob_fn()
  N, n, m = nls.nls_meta.nequ, nls.meta.nvar, nls.meta.ncon
  x = nls.meta.x0
  Fx = residual(nls, x)
  Jx = jac_op_residual(nls, x)
  nFx = dot(Fx, Fx)
  JtF = norm(Jx' * Fx)
  ncx = m > 0 ? @sprintf("%10.4e", norm(cons(nls, x))) : "NA"
  @printf("%-15s  %4d  %4d  %4d  %10.4e  %10.4e  %10s\n",
          prob, N, n, m, nFx, JtF, ncx)
end
println()


text:
#################
# Shift & Scale #
#################
struct Shift{T, N} <: Bijector{N}
    a::T
end

Shift(a::T; dim::Val{D} = Val(0)) where {T<:Real, D} = Shift{T, D}(a)
Shift(a::A; dim::Val{D} = Val(N)) where {T, D, N, A<:AbstractArray{T, N}} = Shift{A, N}(a)

(b::Shift)(x) = b.a + x
(b::Shift{<:Real})(x::AbstractArray) = b.a .+ x
(b::Shift{<:AbstractVector})(x::AbstractMatrix) = b.a .+ x

inv(b::Shift{T, N}) where {T, N} = Shift{T, N}(-b.a)

# FIXME: implement custom adjoint to ensure we don't get tracking
logabsdetjac(b::Shift{T, N}, x) where {T, N} = _logabsdetjac_shift(b.a, x, Val(N))

_logabsdetjac_shift(a::Real, x::Real, ::Val{0}) = zero(eltype(x))
_logabsdetjac_shift(a::Real, x::AbstractVector{T}, ::Val{0}) where {T<:Real} = zeros(T, length(x))
_logabsdetjac_shift(a::T1, x::AbstractVector{T2}, ::Val{1}) where {T1<:Union{Real, AbstractVector}, T2<:Real} = zero(T2)
_logabsdetjac_shift(a::T1, x::AbstractMatrix{T2}, ::Val{1}) where {T1<:Union{Real, AbstractVector}, T2<:Real} = zeros(T2, size(x, 2))

function _logabsdetjac_shift(a::TrackedReal, x::Real, ::Val{0})
    return Tracker.param(_logabsdetjac_shift(data(a), data(x), Val(0)))
end
function _logabsdetjac_shift(a::TrackedReal, x::AbstractVector{T}, ::Val{0}) where {T<:Real}
    return Tracker.param(_logabsdetjac_shift(data(a), data(x), Val(0)))
end
function _logabsdetjac_shift(a::T1, x::AbstractVector{T2}, ::Val{1}) where {T1<:Union{TrackedReal, TrackedVector}, T2<:Real}
    return Tracker.param(_logabsdetjac_shift(data(a), data(x), Val(1)))
end
function _logabsdetjac_shift(a::T1, x::AbstractMatrix{T2}, ::Val{1}) where {T1<:Union{TrackedReal, TrackedVector}, T2<:Real}
    return Tracker.param(_logabsdetjac_shift(data(a), data(x), Val(1)))
end


text:
module TestProductRF

include("preamble.jl")

@testset "count_sum" begin
    @testset "$fold" for fold in [foldl, simple_reduce, foldxt_bs1, foldxt]
        foldl = nothing
        count_sum(xs, xf = Map(identity)) =
            fold(ProductRF(+, +), opcompose(xf, Map(x -> (1, x))), xs)

        @test count_sum([5, 2, 6, 8, 3]) == (5, 24)
        @test count_sum([5, 2, 6, 8, 3], Filter(isodd)) == (2, 8)
        @test_throws EmptyResultError count_sum(1:0)
        @test_throws EmptyResultError count_sum(1:10, Filter(_ -> false))
    end
end

@testset "partially started" begin
    @testset "$fold" for fold in [foldl, simple_reduce, foldxt_bs1, foldxt]
        foldl = nothing
        f(xs) =
            fold(ProductRF(min, reducingfunction(Filter(isodd), max)), Map(x -> (x, x)), xs)
        @test f([5, 2, 6, 8, 3]) == (2, 5)
        @test_throws EmptyResultError f(2:2:8)
    end
end

@testset "init" begin
    @testset "$fold" for fold in [foldl, simple_reduce, foldxt_bs1, foldxt]
        foldl = nothing
        @test fold(ProductRF(+, +), Map(x -> (x, x)), 1:9) == (45, 45)
        @test fold(ProductRF(+, +), Map(x -> (x, x)), 1:9, init = (0, 0)) == (45, 45)
    end
    @test foldl(ProductRF(+, +), Map(x -> (x, x)), 1:9; init = (1, 2)) == (46, 47)
end

@testset "standard fold" begin
    @test mapfoldl(x -> (x, x), ProductRF(min, max), [5, 2, 6, 8, 3]) == (2, 8)
end

@testset "error handling" begin
    rf = ProductRF(min, max)
    err = @test_error rf(0, (1, 2))
    msg = sprint(showerror, err)
    @test occursin("incompatible accumulator", msg)
end

end  # module


text:
using PyCall
@pyimport bs4

soup=bs4.BeautifulSoup(test_txt,"html.parser")
list=soup[:find](class="suggestions")[:find_all]("li");
for li in list
    word=li[:find](class="word")[:get_text]()
    def=li[:find](class="definition")[:get_text]()
    println("\033[1m $word \033[0m")
    println(def)
    print("\n")

end

max_stars_repo_name: shencebebetterme/my-sanity-arxiv

max_issues_repo_name: shencebebetterme/my-sanity-arxiv

max_forks_repo_name: shencebebetterme/my-sanity-arxiv


text:
# Make a graphWidth x graphWidth grid with random weights on edges from 1 : maxWeight
function makeGrid(graphWidth::Int64, maxWeight::Int64)
    vertexAmount = graphWidth ^ 2
    graph = SimpleWeightedGraphs.SimpleWeightedGraph(vertexAmount)
    for i::Int64 in 1 : vertexAmount
        if (i % graphWidth > 0)
            LightGraphs.add_edge!(graph, i, i + 1, rand(1:maxWeight))
        end
        if ((i - 1) / graphWidth < graphWidth - 1)
            LightGraphs.add_edge!(graph, i, i + graphWidth, rand(1:maxWeight))
        end
    end
    return graph
end;

# Create an A* graph square graph with size graphWidth x graphWidth
function makeGrid(graphWidth::Int64)
    vertexAmount = graphWidth ^ 2
    graph = SimpleGraph(vertexAmount, 0; seed = -1)
    for i::Int64 in 1 : vertexAmount
        if (i % graphWidth > 0)
            LightGraphs.add_edge!(graph, i, i + 1)
        end
        if ((i - 1) / graphWidth < graphWidth - 1)
            LightGraphs.add_edge!(graph, i, i + graphWidth)
        end
    end
    return graph
end;

# Remove a bunch of edges on g to make a maze
function makeMaze(g::AbstractGraph)
    numberOfVertices = nv(g)
    vertices = LightGraphs.vertices(g)

    # Vertex array is a list of unexplored things
    vertexArray = Array{Int64}(undef, numberOfVertices)
    for vertex in vertices
        vertexArray[vertex] = vertex;
    end

    # Create a set which will hold visited vertices
    closedVertices = []

    currentVertex = vertexArray[rand(1:length(vertexArray))]
    push!(closedVertices, currentVertex)
    filter!(x -> x != currentVertex, vertexArray)

    # Create an array of nodes which we can reach from current set
    reachable = []
    for reachableVertex in LightGraphs.outneighbors(g, currentVertex)
        push!(reachable, reachableVertex)
    end

    # Make a maze without any edges
    maze = SimpleGraph(numberOfVertices, 0; seed = -1)

    # Intersection is vertices which are both reachable and not visited
    intersection = []

    for i in 1 : (numberOfVertices - 1)
        # Find reachable unvisited nodes
        intersection = intersect(reachable, vertexArray)

        # Grab a random intersection vertex and close it / remove it from open
        randomIntersection = rand(1:length(intersection))
        currentVertex = intersection[randomIntersection]
        push!(closedVertices, currentVertex)
        filter!(x -> x != currentVertex, vertexArray)

        # Add the reachable things
        for reachableVertex in LightGraphs.outneighbors(g, currentVertex)
            push!(reachable, reachableVertex)
        end

        # Grab a random neighbour of our current one and add an edge to it
        neighbours = LightGraphs.outneighbors(g, currentVertex)
        while true
            random = rand(1:length(neighbours))
            if issubset(neighbours[random], closedVertices)
                add_edge!(maze, currentVertex, neighbours[random])
                break
            end
        end
    end
    return maze
end;

# Draw a graphWidth x graphWidth with connections defined by maze edges, saves as pictureName.png
function drawEdges(mazeEdges, graphWidth::Int64, pictureName::String)
    Drawing(10 * graphWidth, 10 * graphWidth, pictureName * ".png")
    background("black")
    sethue("purple")
    for j in 0 : graphWidth - 1
        for i in 0 : graphWidth - 1
            rect((10 * i) + 1, (10 * j) + 1, 8, 8, :fill)
        end
    end

    for edge in mazeEdges
        source = edge.src
        dest = edge.dst
        sourceX = ((source - 1) % graphWidth)
        sourceY = convert(Int64, div((source - 1), graphWidth))
        destX = ((dest - 1) % graphWidth)
        destY = convert(Int64, div((dest - 1), graphWidth))

        rect((10 * sourceX) + 1, (10 * sourceY) + 1, 8 + (10 * (destX - sourceX)), 8 + (10 * (destY - sourceY)), :fill)
    end

    Luxor.finish()
end;

# Draw a graphWidth x graphWidth with connections defined by maze edges and path overlayed, saves as pictureName.png
function drawEdges(mazeEdges, path, graphWidth::Int64, pictureName::String)
    Drawing(10 * graphWidth, 10 * graphWidth, pictureName * ".png")
    background("black")

    sethue("purple")
    for j in 0 : graphWidth - 1
        for i in 0 : graphWidth - 1
            rect((10 * i) + 1, (10 * j) + 1, 8, 8, :fill)
        end
    end

    for edge in mazeEdges
        source = edge.src
        dest = edge.dst
        sourceX = ((source - 1) % graphWidth)
        sourceY = convert(Int64, div((source - 1), graphWidth))
        destX = ((dest - 1) % graphWidth)
        destY = convert(Int64, div((dest - 1), graphWidth))

        rect((10 * sourceX) + 1, (10 * sourceY) + 1, 8 + (10 * (destX - sourceX)), 8 + (10 * (destY - sourceY)), :fill)
    end

    sethue("silver")
    for edge in path
        source = edge.src
        dest = edge.dst
        if source < dest
            sourceX = ((source - 1) % graphWidth)
            sourceY = convert(Int64, div((source - 1), graphWidth))
            destX = ((dest - 1) % graphWidth)
            destY = convert(Int64, div((dest - 1), graphWidth))
        else
            destX = ((source - 1) % graphWidth)
            destY = convert(Int64, div((source - 1), graphWidth))
            sourceX = ((dest - 1) % graphWidth)
            sourceY = convert(Int64, div((dest - 1), graphWidth))
        end

        rect((10 * sourceX) + 1, (10 * sourceY) + 1, 8 + (10 * (destX - sourceX)), 8 + (10 * (destY - sourceY)), :fill)
    end

    Luxor.finish()
end;

# Create a graph with graphLength number of vertices and connections number of edges from each vertex
function randomWeightedGraph(graphLength, connections)
    g = LightGraphs.SimpleGraph(connections)
    g = barabasi_albert!(g, graphLength, connections)

    sources = Array{Int64, 1}(undef, ne(g))
    dests = Array{Int64, 1}(undef, ne(g))
    weights = Array{Float64, 1}(undef, ne(g))

    unweightedEdges = LightGraphs.edges(g)

    i = 1
    for edge in unweightedEdges
        sources[i] = edge.src
        dests[i] = edge.dst
        weights[i] = rand(Float64) * 100
        i += 1
    end

    weightedG = SimpleWeightedGraphs.SimpleWeightedGraph(sources, dests, weights)
end

# Example code to draw a maze and its path
# graphWidth = 100;
# graph = makeGrid(graphWidth);
# maze = makeMaze(graph);
# path = a_star(maze, 1, 100);
# mazeEdges = LightGraphs.edges(maze);
# drawEdges(mazeEdges, graphWidth, "maze");
# drawEdges(mazeEdges, path, graphWidth, "path");

# graphWidth = 100
# g = makeGrid(graphWidth, 10)
# gEdges = LightGraphs.edges(g)
# primsEdges = LightGraphs.prim_mst(g);
# drawEdges(gEdges, graphWidth, "weightedGraph");
# drawEdges(primsEdges, graphWidth, "mst");


text:
import TPT
using Base.Test
using Plots; pyplot()

println("--- Unit NFETB ---")

#
# Liquid Fe
# The parameters are taken from: Bretonnet et al.: J. Phys.: Condens.
# Matter, 7 (1995), 3803-3815
# a is modified so that WCA can produce realistic g(r).
#
Ω = 89.35698
ρ = 1 / Ω
T = 1833.
zs = 1.4
rc = 1.540
a = 0.333
zd = 6.6
rd = 1.512

ahs = TPT.AHS(σ = 4.25, ρ = ρ, approx=:RFA)
lwca = TPT.LWCA(ahs, T, struct=:full)
wca = TPT.WCA(ahs, T)

nfe = TPT.NFE(ρ, zs, TPT.BretonnetSilbert(zs, rc, a))
tb = TPT.WHTB(zd, rd, version=:original)
nfetb = TPT.NFETB(nfe, tb)

sys1 = TPT.TPTSystem(ahs, nfetb, T = 1833., m = 1.018e+5)
sys2 = TPT.TPTSystem(wca, nfetb, m = 1.018e+5)
sys3 = TPT.TPTSystem(lwca, nfetb, m = 1.018e+5)

resid2 = sys2.ref.residue
resid3 = sys2.ref.residue

σ_estim = TPT.hsdiameter_estimate(nfetb, 1883.)[1,1]
σ1 = TPT.hsdiameter(sys1.ref)[1,1]
σ2 = TPT.hsdiameter(sys2.ref)[1,1]
σ3 = TPT.hsdiameter(sys3.ref)[1,1]

u_nfe = TPT.pairpotential(nfe)[1,1]
u_tb = TPT.pairpotential(tb)[1,1]
u_tot = TPT.pairpotential(nfetb)[1,1]
u′_tot = TPT.pairpotential_derivative(nfetb)[1,1]
rmin = TPT.pairpotential_minimizer(nfetb)[1,1]

plot([u_nfe, u_tb, u_tot, u′_tot], 2, 10, ylims = (-0.06, 0.1),
     labels = ["NFE" "TB" "Total" "derivative"])

g1 = TPT.paircorrelation(sys1)[1,1]
g2 = TPT.paircorrelation(sys2)[1,1]
g3 = TPT.paircorrelation(sys3)[1,1]
plot([g1, g2, g3], 3, 7, labels = ["AHS" "WCA" "LWCA"], ylims=:auto)

B2 = TPT.blipfunction(sys2.ref)[1,1]
B3 = TPT.blipfunction(sys3.ref)[1,1]
plot([r -> B2(r)*r^2, r -> B3(r)*r^2], 3.0, 5.5, ylims=:auto, labels=["WCA" "LWCA"])

S1 = TPT.structurefactor(sys1.ref)[1,1]
S2 = TPT.structurefactor(sys2.ref)[1,1]
S3 = TPT.structurefactor(sys3.ref)[1,1]
plot([S1, S2, S3], 0.1, 6, labels = ["AHS" "WCA" "LWCA"], ylims=:auto)

# Total free energy
F1 = TPT.helmholtz(sys1)
F2 = TPT.helmholtz(sys2)

@testset "Unit NFETB" begin
  @test isapprox(σ2, 4.24, atol=1e-2)
  @test isapprox(σ3, 4.26, atol=1e-2)
  @test isapprox(F1, -0.876, atol=1e-3)
  @test isapprox(F2, -0.884, atol=1e-3)
end


text:
export FixedParameterHomotopy, fix_parameters

"""
    FixedParameterHomotopy(H:AbstractHomotopy, parameters)

Construct a homotopy from the given [`AbstractHomotopy`](@ref) `H` with the given `parameters`
fixed.
"""
struct FixedParameterHomotopy{S<:AbstractHomotopy,T} <: AbstractHomotopy
    homotopy::S
    parameters::Vector{T}
end
FixedParameterHomotopy(H::Homotopy, p; compile::Union{Bool,Symbol} = COMPILE_DEFAULT[]) =
    FixedParameterHomotopy(fixed(H; compile = compile), p)
FixedParameterHomotopy(
    H::AbstractHomotopy,
    p;
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
) = FixedParameterHomotopy(H, p)
Base.size(H::FixedParameterHomotopy) = size(H.homotopy)

ModelKit.variables(H::FixedParameterHomotopy) = variables(H.homotopy)
ModelKit.parameters(H::FixedParameterHomotopy) = Variable[]
ModelKit.variable_groups(H::FixedParameterHomotopy) = variable_groups(H.homotopy)

(H::FixedParameterHomotopy)(x, t, p = nothing) = H.homotopy(x, t, H.parameters)

ModelKit.evaluate!(u, H::FixedParameterHomotopy, x, t) =
    evaluate!(u, H.homotopy, x, t, H.parameters)
ModelKit.evaluate_and_jacobian!(u, U, H::FixedParameterHomotopy, x, t) =
    evaluate_and_jacobian!(u, U, H.homotopy, x, t, H.parameters)
ModelKit.taylor!(u, v::Val, H::FixedParameterHomotopy, tx, t) =
    taylor!(u, v, H.homotopy, tx, t, H.parameters)

"""
    fix_parameters(H::Union{Homotopy,AbstractHomotopy}, p; compile::Union{Bool,Symbol} = $(COMPILE_DEFAULT[]))

Fix the parameters of the given homotopy `H`. Returns a [`FixedParameterHomotopy`](@ref).
"""
fix_parameters(
    H::Union{Homotopy,AbstractHomotopy},
    p;
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
) = FixedParameterHomotopy(H, p; compile = compile)


text:
#
# Defaults
#
CIRCLE_DEFAULT_PLOT_POINTS = 720

#
# Plotting Recipe – Poincaré Ball
#
@recipe function f(
    M::Hyperbolic{2},
    pts::AbstractVector{P};
    circle_points = CIRCLE_DEFAULT_PLOT_POINTS,
    geodesic_interpolation = -1,
    hyperbolic_border_color = RGBA(0.0, 0.0, 0.0, 1.0),
) where {P<:PoincareBallPoint}
    @series begin
        φr = range(0, stop = 2 * π, length = circle_points)
        x = [cos(φ) for φ in φr]
        y = [sin(φ) for φ in φr]
        seriestype := :path
        seriescolor := hyperbolic_border_color
        label := ""
        x, y
    end
    show_axis --> false
    framestyle -> :none
    axis --> false
    xlims --> (-1.01, 1.01)
    ylims --> (-1.01, 1.01)
    grid --> false
    aspect_ratio --> :equal
    tickfontcolor --> RGBA(1.0, 1.0, 1.0, 1.0)
    if geodesic_interpolation < 0
        seriestype --> :scatter
        return [p.value[1] for p in pts], [p.value[2] for p in pts]
    else
        lpts = empty(pts)
        for i in 1:(length(pts) - 1)
            # push interims points on geodesics between two points.
            push!(
                lpts,
                shortest_geodesic(
                    M,
                    pts[i],
                    pts[i + 1],
                    collect(range(0, 1, length = geodesic_interpolation + 2))[1:(end - 1)], # omit end point
                )...,
            )
        end
        push!(lpts, last(pts)) # add last end point
        # split into x, y and plot as curve
        seriestype --> :path
        return [p.value[1] for p in lpts], [p.value[2] for p in lpts]
    end
end
# Tangents as quiver
@recipe function f(
    ::Hyperbolic{2},
    pts::AbstractVector{P},
    vecs::AbstractVector{T},
    circle_points = CIRCLE_DEFAULT_PLOT_POINTS,
    hyperbolic_border_color = RGBA(0.0, 0.0, 0.0, 1.0),
) where {P<:PoincareBallPoint,T<:PoincareBallTVector}
    @series begin
        φr = range(0, stop = 2 * π, length = circle_points)
        x = [cos(φ) for φ in φr]
        y = [sin(φ) for φ in φr]
        seriestype := :path
        seriescolor := hyperbolic_border_color
        label := ""
        x, y
    end
    show_axis --> false
    framestyle -> :none
    axis --> false
    xlims --> (-1.01, 1.01)
    ylims --> (-1.01, 1.01)
    grid --> false
    aspect_ratio --> :equal
    tickfontcolor --> RGBA(1.0, 1.0, 1.0, 1.0)
    seriestype := :quiver
    quiver := ([v.value[1] for v in vecs], [v.value[2] for v in vecs])
    return [p.value[1] for p in pts], [p.value[2] for p in pts]
end
#
# Plotting Recipe – Poincaré Half plane
#
@recipe function f(
    M::Hyperbolic{2},
    pts::AbstractVector{T};
    geodesic_interpolation = -1,
) where {T<:PoincareHalfSpacePoint}
    aspect_ratio --> :equal
    framestyle --> :origin
    if geodesic_interpolation < 0
        seriestype --> :scatter
        return [p.value[1] for p in pts], [p.value[2] for p in pts]
    else
        lpts = empty(pts)
        for i in 1:(length(pts) - 1)
            # push interims points on geodesics between two points.
            push!(
                lpts,
                shortest_geodesic(
                    M,
                    pts[i],
                    pts[i + 1],
                    collect(range(0, 1, length = geodesic_interpolation + 2))[1:(end - 1)], # omit end point
                )...,
            )
        end
        push!(lpts, last(pts)) # add last end point
        # split into x, y, z and plot as curve
        seriestype --> :path
        return [p.value[1] for p in lpts], [p.value[2] for p in lpts]
    end
end
# Tangents as quiver
@recipe function f(
    ::Hyperbolic{2},
    pts::AbstractVector{P},
    vecs::AbstractVector{T},
) where {P<:PoincareHalfSpacePoint,T<:PoincareHalfSpaceTVector}
    aspect_ratio --> :equal
    framestyle --> :origin
    seriestype := :quiver
    quiver := ([v.value[1] for v in vecs], [v.value[2] for v in vecs])
    return [p.value[1] for p in pts], [p.value[2] for p in pts]
end


text:
using FFTW,Plots;N=128;P=64N;dt=1/4N;NT=1024;W=200;w=W/P*N;K=zeros(NT,4);
x=rand(P);v=collect(1:P.>P/2).*2 .-1;u()=(x.=mod.(x.+v/2*dt,1));n=zeros(N);
k=im.*2π*vcat(1,1:N/2,-N/2+1:-1);f(x)=Int(mod1(round(x*N),N));Σ=sum
@gif for t in 1:NT;
  u();E=real.(ifft((E=fft((n.*=0;for j∈x;n[f(j)]+=w;end;n))./k;E[1]*=0;E)));
  u();v+=E[f.(x)]*dt;scatter(x,v,ylims=(-3,3),legend=false);K[t,4]=sum(v)/P;
  K[t,1:2].=Σ.((E.^2N,v.^2*W/P))./W;K[t,3]=Σ(K[t,1:2]);plot!((1:t)/NT,K[1:t,:])
end every 8


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct EffectiveNetworkSecurityGroupAssociation <: SwaggerModel
    subnet::Any # spec type: Union{ Nothing, SubResource } # spec name: subnet
    networkInterface::Any # spec type: Union{ Nothing, SubResource } # spec name: networkInterface

    function EffectiveNetworkSecurityGroupAssociation(;subnet=nothing, networkInterface=nothing)
        o = new()
        validate_property(EffectiveNetworkSecurityGroupAssociation, Symbol("subnet"), subnet)
        setfield!(o, Symbol("subnet"), subnet)
        validate_property(EffectiveNetworkSecurityGroupAssociation, Symbol("networkInterface"), networkInterface)
        setfield!(o, Symbol("networkInterface"), networkInterface)
        o
    end
end # type EffectiveNetworkSecurityGroupAssociation

const _property_map_EffectiveNetworkSecurityGroupAssociation = Dict{Symbol,Symbol}(Symbol("subnet")=>Symbol("subnet"), Symbol("networkInterface")=>Symbol("networkInterface"))
const _property_types_EffectiveNetworkSecurityGroupAssociation = Dict{Symbol,String}(Symbol("subnet")=>"SubResource", Symbol("networkInterface")=>"SubResource")
Base.propertynames(::Type{ EffectiveNetworkSecurityGroupAssociation }) = collect(keys(_property_map_EffectiveNetworkSecurityGroupAssociation))
Swagger.property_type(::Type{ EffectiveNetworkSecurityGroupAssociation }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_EffectiveNetworkSecurityGroupAssociation[name]))}
Swagger.field_name(::Type{ EffectiveNetworkSecurityGroupAssociation }, property_name::Symbol) =  _property_map_EffectiveNetworkSecurityGroupAssociation[property_name]

function check_required(o::EffectiveNetworkSecurityGroupAssociation)
    true
end

function validate_property(::Type{ EffectiveNetworkSecurityGroupAssociation }, name::Symbol, val)
end


text:
using JuMP
using CPLEX
using Combinatorics

include("LSP_PLNE.jl")
include("InstanceLoader.jl")
include("ResolvePlne.jl")
include("Helper.jl")
include("VRP_PLNE.jl")

function createPDI_Bard_Nananukul_compacte(params, nodes, demands, costs)
    model = Model(CPLEX.Optimizer)
    n = params["n"] #nombre de revendeur (exclu le fournisseur donc)
    m = params["k"] #nb de véhicules
    l = params["l"] #nombre de pas de temps
    Q = params["Q"] #capacité d'un vahicule
    C = params["C"] #production max du fournisseur à un certain pas de temps
    u = params["u"] #coût de production d'une unité
    f = params["f"] #coût de setup

    M=Dict(t=>min(C,sum(sum( demands[i,j] for i in 1:n) for j in t:l)) for t in 1:l )
    M_tilde=Dict((i,t)=> min( nodes[i]["L"], Q, sum(demands[i,j] for j in t:l) ) for i in 1:n,t in 1:l)

    @variable(model, p[t=1:l]>=0) #(13)

    @variable(model, I[0:n, 0:l] >= 0) #(13)
    for i in 0:params["n"]
        @constraint(model, I[i, 0] == nodes[i]["L0"])
    end

    @variable(model, q[1:n, 1:l] >= 0) #(13)
    @variable(model, y[1:l], Bin) #(14)

    @variable(model, x[i=0:n, j=0:n,t=1:l], Bin) #(14)
    for i in 0:n, t in 1:l
        delete(model, x[i, i, t]) #on enlève les variables qui correspondent aux arêtes en trop (les (i, i))
    end

    @variable(model, z[i=1:n, t=1:l], Bin) #(15)
    @variable(model, 0 <= z0[ 0 , t=1:l ] <= m, Int) #(16) et (10) #ATTENTION IL Y A DEUX z À GERER: z0=z QUAND i=0 CAR JuMP NE NOUS LAISSE PAS INSTANCIER EN DEUX FOIS UNE VARIABLE, OR C'EST NECESSAIRE CAR LES z0 SONT DES ENTIERS POSITIFS ET LES z SONT DES NOMBRE BINAIRES
    # @variable(model, 0<=w[i=1:n,t=1:l]) #(12 partie 1)

    for t in 1:l

        @constraint(model, I[0, t-1] + p[t] == I[0, t] + sum(q[i, t] for i in 1:n)) #(2)
        @constraint(model, p[t] <= M[t] * y[t]) #(4)
        @constraint(model, I[0, t] <= nodes[0]["L"]) #(5)

        @constraint(model, sum( x[j,0,t] + x[0,j,t] for j in 1:n ) == 2*z0[0,t]) #(9 que pour i=0)

        for i in 1:n
			@constraint(model, I[i, t-1] + q[i, t] == demands[i, t] + I[i, t]) #(3)
			@constraint(model, I[i, t-1] + q[i, t] <= nodes[i]["L"]) #(6)
            @constraint(model,q[i,t]<=M_tilde[i,t]*z[i,t]) #(7) la qté à livrer ne peut être positive que si on visite le noeud à ce pas de temps et ne doit pas dépasser la capacité d'un seul véhicule (un noeud n'est servi que par un seul véhicule)

            nodesIndexWithoutI = filter(e -> e != i, 0:n)
            @constraint(model, sum( x[i,j,t] for j in nodesIndexWithoutI ) == z[i,t]) #(8)
            @constraint(model, sum( x[j,i,t] + x[i,j,t] for j in nodesIndexWithoutI ) == 2*z[i,t]) #(9 avec i in 1:n)

            # @constraint(model,w[i , t] <= Q*z[i,t]) # (12 partie 2)
		end

        # for ((i,j),edgeCost) in costs
        #     if(i!=0 && j!=0)
        #         @constraint(model,w[i,t]-w[j,t]>=q[i,t]-M_tilde[i,t]*(1-x[i,j,t])) #(11)
        #     end
        # end
    end

    @objective(model,Min,sum(u*p[t]+f*y[t]+sum(nodes[i]["h"]*I[i,t] for i in 0:n)+sum(edgeCost*x[i,j,t] for ((i,j),edgeCost) in costs) for t in 1:l)) #(1)

    return model

end

#=
Formulation non compacte, la contrainte (29) contient 2^n inégalités
=#
function createPDI_Boudia(params, nodes, demands, costs)
    model = Model(CPLEX.Optimizer)
    n = params["n"] #nombre de revendeur (exclu le fournisseur donc)
    m = params["k"] #nb de véhicules
    l=params["l"] #nombre de pas de temps
    Q = params["Q"] #capacité d'un vahicule
    C=params["C"] #production max du fournisseur à un certain pas de temps
    u=params["u"] #coût de production d'une unité
    f=params["f"] #coût de setup
    a=1

    M=Dict(t=>min(C,sum(sum( demands[i,j] for i in 1:n) for j in t:l)) for t in 1:l )
    M_tilde=Dict((i,t)=> min( nodes[i]["L"], Q, sum(demands[i,j] for j in t:l) ) for i in 1:n,t in 1:l)

    @variable(model,p[t=1:l]>=0) #(31)

    @variable(model, I[0:n, 0:l] >= 0) #(31)
    for i in 0:params["n"]
        @constraint(model, I[i, 0] == nodes[i]["L0"])
    end

    @variable(model, q[1:n, k=1:m, 1:l] >= 0) #(31)
    @variable(model, y[1:l], Bin) #(32)

    @variable(model, x[i=0:n, j=0:n, k=1:m, t=1:l], Bin) #(32)
    for i in 0:n, t in 1:l, k=1:m
        delete(model, x[i, i, k, t]) #on enlève les variables qui correspondent aux arêtes en trop (les (i, i))
    end

    @variable(model, z[i=0:n, k=1:m, t=1:l], Bin) #(32)

    for t in 1:l
        @constraint(model, I[0,t-1]+p[t] == sum(sum(q[i,k,t]+I[0,t] for k in 1:m) for i in 1:n)) #(21)
        @constraint(model,p[t]<=M[t]*y[t]) #(23)
        @constraint(model, I[0,t]<=nodes[0]["L"]) #(24)
        for i in 1:n
            @constraint(model,I[i,t-1]+sum(q[i,k,t] for k in 1:m)==demands[i,t]+I[i,t]) #(22)
            #@constraint(model, I[i,t-1]+sum(q[k,i,t] for k in 1:m)<=nodes[i]["L"]) #(25) #erreur sur le PDF ?
            @constraint(model, I[i,t-1]+sum(q[i,k,t] for k in 1:m)<=nodes[i]["L"]) #(25)
            @constraint(model,sum(z[i,k,t] for k in 1:m)<=1) #(27)
            for k in 1:m
                @constraint(model,q[i,k,t]<=M_tilde[i,t]*z[i,k,t]) #(26)
                nodesIndexWithoutI = filter(e -> e != i, 0:n)
                @constraint(model,sum(x[j,i,k,t] for j in nodesIndexWithoutI)+sum(x[i,j,k,t] for j in nodesIndexWithoutI)==2*z[i,k,t]) #(28 pour i>=1)
            end
        end

        for k in 1:m
            @constraint(model,sum(x[j,0,k,t] for j in 0:n)+sum(x[0,j,k,t] for j in 1:n)==2*z[0,k,t]) #(28 pour i=0)
            @constraint(model,sum(q[i,k,t] for i in 1:n)<=Q*z[0,k,t]) #(30)

            for S in powerset(1:n, 2, n)
                @constraint(model,sum(sum(x[i,j,k,t] for j in S) for i in S)<=length(S)-1) #(29)
            end
        end

    end
    @objective(model,Min,sum(u*p[t]+f*y[t]+sum(nodes[i]["h"]*I[i,t] for i in 0:n)+sum(edgeCost*sum(x[i,j,k,t] for k in 1:m) for ((i,j),edgeCost) in costs) for t in 1:l)) #(20)

    return model

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VpnDeviceScriptParameters <: SwaggerModel
    vendor::Any # spec type: Union{ Nothing, String } # spec name: vendor
    deviceFamily::Any # spec type: Union{ Nothing, String } # spec name: deviceFamily
    firmwareVersion::Any # spec type: Union{ Nothing, String } # spec name: firmwareVersion

    function VpnDeviceScriptParameters(;vendor=nothing, deviceFamily=nothing, firmwareVersion=nothing)
        o = new()
        validate_property(VpnDeviceScriptParameters, Symbol("vendor"), vendor)
        setfield!(o, Symbol("vendor"), vendor)
        validate_property(VpnDeviceScriptParameters, Symbol("deviceFamily"), deviceFamily)
        setfield!(o, Symbol("deviceFamily"), deviceFamily)
        validate_property(VpnDeviceScriptParameters, Symbol("firmwareVersion"), firmwareVersion)
        setfield!(o, Symbol("firmwareVersion"), firmwareVersion)
        o
    end
end # type VpnDeviceScriptParameters

const _property_map_VpnDeviceScriptParameters = Dict{Symbol,Symbol}(Symbol("vendor")=>Symbol("vendor"), Symbol("deviceFamily")=>Symbol("deviceFamily"), Symbol("firmwareVersion")=>Symbol("firmwareVersion"))
const _property_types_VpnDeviceScriptParameters = Dict{Symbol,String}(Symbol("vendor")=>"String", Symbol("deviceFamily")=>"String", Symbol("firmwareVersion")=>"String")
Base.propertynames(::Type{ VpnDeviceScriptParameters }) = collect(keys(_property_map_VpnDeviceScriptParameters))
Swagger.property_type(::Type{ VpnDeviceScriptParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VpnDeviceScriptParameters[name]))}
Swagger.field_name(::Type{ VpnDeviceScriptParameters }, property_name::Symbol) =  _property_map_VpnDeviceScriptParameters[property_name]

function check_required(o::VpnDeviceScriptParameters)
    true
end

function validate_property(::Type{ VpnDeviceScriptParameters }, name::Symbol, val)
end


text:
@testset "Wiggle" begin
  x = randn(8)
  y = acwt(x, wavelet(WT.db4))         # Normal ACDWT
  @test typeof(wiggle(y)) == Plots.Plot{Plots.GRBackend}
  @test typeof(wiggle!(y)) == Plots.Plot{Plots.GRBackend}
end

max_stars_repo_name: ShozenD/AutocorrelationShell

max_issues_repo_name: ShozenD/AutocorrelationShell

max_forks_repo_name: ShozenD/AutocorrelationShell


text:
function gss_vech(x)
# PURPOSE: creates a column vector by stacking columns of x
#          on and below the diagonal
#----------------------------------------------------------
# USAGE:  v = vech(x)
# where:  x = an input matrix
#---------------------------------------------------------
# RETURNS:
#         v = output vector containing stacked columns of x
#----------------------------------------------------------

# Written by Mike Cliff, UNC Finance  mcliff@unc.edu
# CREATED: 12/08/98
  if (nargin != 1)
    error('Wrong # of arguments to vech');
  end

  r, c = size(x);
  v = [];
  for i = 1:c
    v = [v ; x(i:r,i)];
  end
  return v
end