text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from src/cholmod.jl, to deprecate
@eval SuiteSparseNext.CHOLMOD begin
    Base.Ac_ldiv_B(A::RealHermSymComplexHermF64SSL, B::StridedVecOrMat) = \(adjoint(A), B)
    Base.Ac_ldiv_B(L::Factor, B::Dense) = \(adjoint(L), B)
    Base.Ac_ldiv_B(L::Factor, B::VecOrMat) = \(adjoint(L), B)
    Base.Ac_ldiv_B(L::Factor, B::Sparse) = \(adjoint(L), B)
    Base.Ac_ldiv_B(L::Factor, B::SparseVecOrMat) = \(adjoint(L), B)
    Base.Ac_ldiv_B(L::FactorComponent, B) = \(adjoint(L), B)
    Base.Ac_ldiv_B(L::FactorComponent, B::RowVector) = \(adjoint(L), B)
    Base.Ac_mul_B(A::Sparse, B::Dense) = *(adjoint(A), B)
    Base.Ac_mul_B(A::Sparse, B::VecOrMat) =  *(adjoint(A), B)
    Base.Ac_mul_B(A::Sparse, B::Sparse) = *(adjoint(A), B)
    Base.A_mul_Bc(A::Sparse{Tv}, B::Sparse{Tv}) where {Tv<:VRealTypes} = *(A, adjoint(B))
end

# A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from src/umfpack.jl, to deprecate
@eval SuiteSparseNext.UMFPACK begin
    using LinearAlgebra: Adjoint, Transpose
    LinearAlgebra.A_ldiv_B!(X::StridedVecOrMat{T}, lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} =
        LinearAlgebra.ldiv!(X, lu, B)
    LinearAlgebra.At_ldiv_B!(X::StridedVecOrMat{T}, lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} =
        LinearAlgebra.ldiv!(X, transpose(lu), B)
    LinearAlgebra.Ac_ldiv_B!(X::StridedVecOrMat{T}, lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} =
        LinearAlgebra.ldiv!(X, adjoint(lu), B)
    LinearAlgebra.A_ldiv_B!(X::StridedVecOrMat{Tb}, lu::UmfpackLU{Float64}, B::StridedVecOrMat{Tb}) where {Tb<:Complex} =
        LinearAlgebra.ldiv!(X, lu, B)
    LinearAlgebra.At_ldiv_B!(X::StridedVecOrMat{Tb}, lu::UmfpackLU{Float64}, B::StridedVecOrMat{Tb}) where {Tb<:Complex} =
        LinearAlgebra.ldiv!(X, transpose(lu), B)
    LinearAlgebra.Ac_ldiv_B!(X::StridedVecOrMat{Tb}, lu::UmfpackLU{Float64}, B::StridedVecOrMat{Tb}) where {Tb<:Complex} =
        LinearAlgebra.ldiv!(X, adjoint(lu), B)
    LinearAlgebra.A_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = LinearAlgebra.ldiv!(B, lu, copy(B))
    LinearAlgebra.At_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = LinearAlgebra.ldiv!(B, transpose(lu), copy(B))
    LinearAlgebra.Ac_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = LinearAlgebra.ldiv!(B, adjoint(lu), copy(B))
    LinearAlgebra.A_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = LinearAlgebra.ldiv!(B, lu, copy(B))
    LinearAlgebra.At_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = LinearAlgebra.ldiv!(B, transpose(lu), copy(B))
    LinearAlgebra.Ac_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = LinearAlgebra.ldiv!(B, adjoint(lu), copy(B))
end

# A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from src/spqr.jl, to deprecate
@eval SuiteSparseNext.SPQR begin
    using LinearAlgebra: Adjoint, Transpose
    LinearAlgebra.A_mul_Bc!(A::StridedMatrix, Q::QRSparseQ) = LinearAlgebra.mul!(A, adjoint(Q))
    LinearAlgebra.Ac_mul_B!(Q::QRSparseQ, A::StridedVecOrMat) = LinearAlgebra.mul!(adjoint(Q), A)
    LinearAlgebra.A_mul_B!(A::StridedMatrix, Q::QRSparseQ) = LinearAlgebra.mul!(A, Q)
    LinearAlgebra.A_mul_B!(Q::QRSparseQ, A::StridedVecOrMat) = LinearAlgebra.mul!(Q, A)
end

# deprecate lufact to lu
@eval SuiteSparseNext.UMFPACK begin
    @deprecate(lufact(A::SparseMatrixCSC), lu(A))
    @deprecate(lufact(S::SparseMatrixCSC{<:UMFVTypes,<:UMFITypes}), lu(S))
    @deprecate(lufact(A::SparseMatrixCSC{<:Union{Float16,Float32},Ti}) where {Ti<:UMFITypes}, lu(A))
    @deprecate(lufact(A::SparseMatrixCSC{<:Union{ComplexF16,ComplexF32},Ti}) where {Ti<:UMFITypes}, lu(A))
    @deprecate(lufact(A::Union{SparseMatrixCSC{T},SparseMatrixCSC{Complex{T}}}) where {T<:AbstractFloat}, lu(A))
end

# deprecate qrfact to qr
@eval SuiteSparseNext.SPQR begin
    import LinearAlgebra: qrfact
    @deprecate(qrfact(A::SparseMatrixCSC{Tv}; tol = _default_tol(A)) where {Tv<:Union{ComplexF64,Float64}}, qr(A; tol=tol))
    @deprecate(qrfact(A::SparseMatrixCSC; tol = _default_tol(A)), qr(A; tol=tol))
end

# deprecate ldltfact to ldlt
@eval SuiteSparseNext.CHOLMOD begin
    import LinearAlgebra: ldltfact
    @deprecate(ldltfact(A::Sparse; shift::Real=0.0, perm::AbstractVector{SuiteSparse_long}=SuiteSparse_long[]), ldlt(A; shift=shift, perm=perm))
    @deprecate(ldltfact(A::Union{SparseMatrixCSC{T},SparseMatrixCSC{Complex{T}},
                        Symmetric{T,SparseMatrixCSC{T,SuiteSparse_long}},
                        Hermitian{Complex{T},SparseMatrixCSC{Complex{T},SuiteSparse_long}},
                        Hermitian{T,SparseMatrixCSC{T,SuiteSparse_long}}};
                        kws...) where {T<:Real},
                ldlt(A; kws...))
end

# deprecate ldltfact! to ldlt!
@eval SuiteSparseNext.CHOLMOD begin
    import LinearAlgebra: ldltfact!
    @deprecate(ldltfact!(F::Factor{Tv}, A::Sparse{Tv}; shift::Real=0.0) where Tv, ldlt!(F, A; shift=shift))
    @deprecate(ldltfact!(F::Factor, A::Union{SparseMatrixCSC{T},
                         SparseMatrixCSC{Complex{T}},
                         Symmetric{T,SparseMatrixCSC{T,SuiteSparse_long}},
                         Hermitian{Complex{T},SparseMatrixCSC{Complex{T},SuiteSparse_long}},
                         Hermitian{T,SparseMatrixCSC{T,SuiteSparse_long}}};
                         shift = 0.0) where {T<:Real},
               ldlt!(F, A; shift=shift))
end

# deprecate cholfact to cholesky
@eval SuiteSparseNext.CHOLMOD begin
    import LinearAlgebra: cholfact
    @deprecate(cholfact(A::Sparse; shift::Real=0.0, perm::AbstractVector{SuiteSparse_long}=SuiteSparse_long[]), cholesky(A; shift=shift, perm=perm))
    @deprecate(cholfact(A::Union{SparseMatrixCSC{T}, SparseMatrixCSC{Complex{T}},
                        Symmetric{T,SparseMatrixCSC{T,SuiteSparse_long}},
                        Hermitian{Complex{T},SparseMatrixCSC{Complex{T},SuiteSparse_long}},
                        Hermitian{T,SparseMatrixCSC{T,SuiteSparse_long}}};
                        kws...) where {T<:Real},
               cholesky(A; kws...))
end

# deprecate cholfact! to cholesky!
@eval SuiteSparseNext.CHOLMOD begin
    import LinearAlgebra: cholfact!
    @deprecate(cholfact!(F::Factor{Tv}, A::Sparse{Tv}; shift::Real=0.0) where Tv, cholesky!(F, A; shift=shift))
    @deprecate(cholfact!(F::Factor, A::Union{SparseMatrixCSC{T},
                        SparseMatrixCSC{Complex{T}},
                        Symmetric{T,SparseMatrixCSC{T,SuiteSparse_long}},
                        Hermitian{Complex{T},SparseMatrixCSC{Complex{T},SuiteSparse_long}},
                        Hermitian{T,SparseMatrixCSC{T,SuiteSparse_long}}};
                        shift = 0.0) where {T<:Real},
                cholesyk!(F, A; shift=shift))
end


text:
export UCBExplorer

using Random
using Flux

Base.@kwdef mutable struct UCBExplorer{R<:AbstractRNG} <: AbstractExplorer
    c::Float64
    actioncounts::Vector{Float64}
    step::Int
    rng::R
    is_training::Bool = true
end

Flux.testmode!(p::UCBExplorer, mode = true) = p.is_training = !mode

"""
    UCBExplorer(na; c=2.0, ϵ=1e-10, step=1, seed=nothing)

# Arguments
- `na` is the number of actions used to create a internal counter.
- `t` is used to store current time step.
- `c` is used to control the degree of exploration.
- `seed`, set the seed of inner RNG.
- `is_training=true`, in training mode, time step and counter will not be updated.
"""
UCBExplorer(na; c = 2.0, ϵ = 1e-10, step = 1, seed = nothing, is_training = true) =
    UCBExplorer(c, fill(ϵ, na), 1, MersenneTwister(seed), is_training)

@doc raw"""
    (ucb::UCBExplorer)(values::AbstractArray)
Unlike [`EpsilonGreedyExplorer`](@ref), uncertaintyies are considered in UCB.

!!! note
    If multiple values with the same maximum value are found.
    Then a random one will be returned!

```math
A_t = \underset{a}{\arg \max} \left[ Q_t(a) + c \sqrt{\frac{\ln t}{N_t(a)}} \right]
```

See more details at Section (2.7) on Page 35 of the book *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*
""" function (p::UCBExplorer)(values::AbstractArray)
    v, inds = find_all_max(@. values + p.c * sqrt(log(p.step + 1) / p.actioncounts))
    action = sample(p.rng, inds)
    if p.is_training
        p.actioncounts[action] += 1
        p.step += 1
    end
    action
end

ext: jl


text:
"""
```
Lifting_Aux::Dict{String, Dict{String, String}}
```
Is a dictionary of dictionaries containing non-essential information about exercises and progressions. Its default keys are,
```
julia> keys(Lifting_Aux)
Base.KeySet for a Dict{String,Dict{String,String}} with 5 entries. Keys:
  "modality"
  "setType"
  "muscles"
  "equipmentSize"
  "equipment"
```
where each entry is a dictionary with its own keys and values.
"""
Lifting_Aux = Dict{String,Dict{String,String}}()

"""
```
Lifting_Exercise_Names::Dict{String, Dict{String, String}}
```
Is a dictionary of dictionaries containing lifting exercises. Its default keys are,
```
julia> import: Lifting_Exercise_Names
julia> keys(Lifting_Exercise_Names)
Base.KeySet for a Dict{String,Dict{String,String}} with 10 entries. Keys:
  "grip"
  "bench"
  "squat"
  "press"
  "isometric"
  "core"
  "pull"
  "row"
  "deadlift"
  "lowBack"
```
where each entry is a dictionary with its own keys and values.
"""
Lifting_Exercise_Names = Dict{String,Dict{String,String}}()

modality = Dict(
    "default" => "Default",
    "belt" => "Belt",
    "tempo" => "Tempo",
    "paused" => "Paused",
    "speed" => "Speed",
    "deficit" => "Deficit",
    "chains" => "Chains",
    "bands" => "Bands",
    "reverse bands" => "Reverse Bands",
    "box" => "Box",
    "block" => "Block",
    "pin" => "Pin",
    "iso" => "Isometric",
    "standing" => "Standing",
    "seated" => "Seated",
    "unilateral" => "Unilateral",
    "reverse grip" => "Reverse Grip",
    "kneeling" => "Kneeling",
    "side" => "Side",
    "cheat" => "Cheat",
    "dragon" => "Dragon",
    "walking" => "Walking",
    "backstep" => "Backstep",
    "straddle" => "Straddle",
    "full" => "Full",
    "one leg out" => "One Leg Out",
    "tuck" => "Tuck",
    "adv tuck" => "Advanced Tuck",
)
setType = Dict(
    "default" => "Default",
    "amrap" => "AMRAP",
    "emom" => "EMOM",
    "amrap emom" => "AMRAP EMOM",
    "test" => "Test",
    "optional" => "Optional",
    "opt amrap" => "Optional AMRAP",
    "opt amrap emom" => "Opt AMRAP EMOM",
)
equipmentSize = Dict(
    "na" => "NA",
    "tiny" => "Tiny",
    "s" => "Small",
    "m" => "Medium",
    "m" => "Large",
    "xl" => "Extra Large",
)
equipment = Dict(
    "bb" => "Barbell",
    "db" => "Dumbell",
    "ssb" => "Safety Squat Bar",
    "trap bar" => "Trap Bar",
    "fat grips" => "Fat Grips",
    "axle" => "Axle",
    "dl bar" => "Deadlift Bar",
    "squat bar" => "Squat Bar",
    "stiff bar" => "Stiff Bar",
    "oly bar" => "Olympic Bar",
    "swiss bar" => "Swiss Bar",
    "duffalo bar" => "Duffalo Bar",
    "sling shot" => "SlingShot",
    "cable" => "Cable Machine",
    "other" => "Other Machine",
    "t bar" => "T-bar",
    "ez bar" => "EZ Bar",
    "none" => "None",
    "pullup bar" => "Pullup Bar",
    "chain belt" => "Chain Belt",
    "plate" => "Plate",
    "parallette" => "Parallette",
    "ab wheel" => "Ab Wheel",
)
muscles = Dict(
    "pecs" => "Pecs",
    "triceps" => "Triceps",
    "front delts" => "Front Delts",
    "side delts" => "Side Delts",
    "rear delts" => "Rear Delts",
    "rhomboids" => "Rhomboids",
    "traps" => "Traps",
    "lats" => "Lats",
    "biceps" => "Biceps",
    "wrist flexors" => "Wrist Flexors",
    "wrist extensors" => "Wrist Extensors",
    "finger flexors" => "Finger Flexors",
    "finger extensors" => "Finger Extensors",
    "hand" => "Hand Muscles",
    "core" => "Core",
    "rectus abdominis" => "Rectus Abdominis",
    "transverse abdominis" => "Transverse Abdominis",
    "obliques" => "Obliques",
    "quads" => "Quads",
    "glute major" => "Glute Major",
    "glute minor" => "Glute Minor",
    "glutes" => "Glutes",
    "hamstrings" => "Hamstrings",
    "calves" => "Calves",
    "spinal erectors" => "Spinal Erectors",
    "brachialis" => "Brachialis",
    "rot cuff" => "Rotator Cuff",
)
push!(Lifting_Aux, "modality" => modality)
push!(Lifting_Aux, "setType" => setType)
push!(Lifting_Aux, "equipmentSize" => equipmentSize)
push!(Lifting_Aux, "equipment" => equipment)
push!(Lifting_Aux, "muscles" => muscles)

grip = Dict(
    # Wrist
    "wrist curl" => "Wrist Curl",
    "wrist extension" => "Wrist Extension",
    "ulnar dev" => "Resisted Ulnar Deviation",
    "radial dev" => "Resisted Radial Deviation",
    "supination" => "Resisted Supination",
    "pronation" => "Resisted Pronation",
    # Support hold
    "dead hang" => "Dead Hang",
    "active hang" => "Active Hang",
    # Pinch
    "plate pinch" => "Plate Pinch",
    "plate flip" => "Plate Flip",
    # Dynamic
    "finger curl" => "Finger Curl",
    "crusher" => "Crusher",
    # Open Hand
    "blob" => "Blob",
    "hub" => "Hub",
    # Rolling
    "rolling thunder" => "Rolling Thunder",
)
squat = Dict(
    # Main variations
    "high" => "High Bar Squat",
    "low" => "Low Bar Squat",
    "front" => "Front Squat",
    # Novelty variations
    "zercher" => "Zercher Squat",
    "hatfield" => "Hatfield Squat",
    "goblet" => "Goblet Squat",
    # Unilateral
    "bulgarian" => "Bulgarian Split Squat",
    "lunge" => "Lunge",
    "pistol" => "Pistol Squat",
    "shrimp" => "Shrimp Squat",
    # Enforced bodyweight.
    "sissy" => "Sissy Squat",
    # Calf raises.
    "calf raise" => "Calf Raise",
)
bench = Dict(
    # Main variations.
    "flat" => "Flat Bench",
    "decline" => "Decline Bench",
    "incline" => "Incline Bench",
    "db incline" => "DB Incline Bench",
    "cg" => "Close Grip Bench",
    "feet up" => "Feet Up Bench",
    "spotto" => "Spoto Bench",
    # Novelty variations.
    "jm" => "JM Press",
    # Bodyweight
    "dips" => "Dips",
    "pushup" => "Pushup",
    "diamond" => "Diamond Pushup",
    "archer" => "Archer Pushup",
)
deadlift = Dict(
    # Main variations.
    "dl" => "Deadlift",
    "sumo" => "Sumo Deadlift",
    "rdl" => "Romanian Deadlift",
    "hip thrust" => "Hip Thrust",
    "single leg hip thrust" => "Single Leg Hip Thrust",
    # Novelty.
    "rack pull" => "Rack Pull",
    "jefferson" => "Jefferson",
    "zercher" => "Zercher",
)
press = Dict(
    "ohp" => "Overhead Press",
    "db ohp" => "Dumbell Overhead Press",
    "push" => "Push Press",
    "z" => "Z Press",
    "hspu" => "Handstand Push Up",
    "side raise" => "Side Raise",
)
row = Dict(
    "strict" => "Strict Row",
    "pendlay" => "Pendlay Row",
    "kroc" => "Kroc Row",
    "other" => "Other Row",
)
pull = Dict(
    "face" => "Face Pull",
    "pullup" => "Pullup",
    "archer" => "Archer Pullup",
    "mu" => "Muscle Up",
    "chin up" => "Chin Up",
    "hammer" => "Hammer Curl",
    "curl" => "Curl",
    "zottman" => "Zottman Curl",
)
lowBack = Dict(
    "nordic" => "Nordic Curl",
    "ghr" => "Glute-Ham Raise",
    "hyper extension" => "Hyper Extension",
)
core = Dict(
    "dragon flag" => "Dragon Flag",
    "hanging leg raise" => "Hanging Leg Raise",
    "toes to bar" => "Toes to Bar",
    "windshield wipers" => "Windshield Wipers",
    "ab roll" => "Ab Rollout",
    "parallette crunch" => "Parallette Crunch",
    "l sit" => "L-Sit",
    "v sit" => "V-Sit",
    "landmine" => "Land Mine Anti-Rotation",
    "plank" => "Plank",
    "leg raise" => "Leg Raise",
    "other" => "Other",
)
isometric = Dict(
    "back" => "Back Lever",
    "front" => "Front Lever",
    "planche" => "Planche",
    "other" => "Other",
)
push!(Lifting_Exercise_Names, "grip" => grip)
push!(Lifting_Exercise_Names, "squat" => squat)
push!(Lifting_Exercise_Names, "bench" => bench)
push!(Lifting_Exercise_Names, "deadlift" => deadlift)
push!(Lifting_Exercise_Names, "press" => press)
push!(Lifting_Exercise_Names, "row" => row)
push!(Lifting_Exercise_Names, "pull" => pull)
push!(Lifting_Exercise_Names, "lowBack" => lowBack)
push!(Lifting_Exercise_Names, "core" => core)
push!(Lifting_Exercise_Names, "isometric" => isometric)

max_stars_repo_name: dcelisgarza/Lifting


text:
# packing and unpacking of Graph related types

using KernelDensityEstimate
using IncrementalInference

using Base.Test

fg = emptyFactorGraph()

N=100

doors = zeros(1,4)
doors[1,1:4] = [-100.0;0.0;100.0;300.0]
pd = kde!(doors,[3.0])
pd = resample(pd,N);
bws = getBW(pd)[:,1]
bws2 = Array{Float64}(length(bws),1)
bws2[:,1] = bws[:]
doors2 = getPoints(pd);
v1 = addNode!(fg,:x1,doors,N=N)
f1  = addFactor!(fg,[v1],Obsv2( doors2, bws2, [1.0]))

tem = 2.0*randn(1,N)+getVal(v1)+50.0
v2 = addNode!(fg,:x2, tem, N=N)
f2 = addFactor!(fg, [:x1, :x2], Odo(50.0*ones(1,1),2.0*ones(1,1),[1.0]))

println("Testing conversion to packed function node data structure and back")

topack = getData(f1) #fg.f[4].attributes["data"]
dd = convert(PackedFunctionNodeData{PackedObsv2},topack)
upd = convert(FunctionNodeData{GenericWrapParam{Obsv2}}, dd)

@test compare(topack, upd)

topack = getData(f2) #fg.f[4].attributes["data"]
dd = convert(PackedFunctionNodeData{PackedOdo},topack)
upd = convert(FunctionNodeData{GenericWrapParam{Odo}}, dd)

@test compare(topack, upd)

packedv4data = FNDencode(IncrementalInference.PackedFunctionNodeData{PackedOdo}, getData(f2))
upv4data = FNDdecode(IncrementalInference.FunctionNodeData{GenericWrapParam{Odo}}, packedv4data)

@test compare(getData(f2), upv4data)

# data structure conversion tests for protobuffing
println("Testing conversion to packed variable node data structure and back")
dat = getData(fg,1)
pd = convert(PackedVariableNodeData, dat) #fg.v[1].attributes["data"]
# pd = convert(PackedVariableNodeData,fg.v[1].attributes["data"])
unpckd = convert(VariableNodeData, pd)

@test IncrementalInference.compare(getData(fg,1), unpckd)
@test getData(fg,1) == unpckd

packedv1data = VNDencoder(IncrementalInference.PackedVariableNodeData, getData(v1))
upv1data = VNDdecoder(IncrementalInference.VariableNodeData, packedv1data)

@test compare(getData(v1), upv1data)

println("Conversions and comparisons agree")

#


text:
module baffled_piston_examples
using FinEtools
using LinearAlgebra
using Gaston

function baffled_piston_H8_transient()
    rho = 1.21*phun("kg/m^3");# mass density
c  = 343.0*phun("m/s");# sound speed
bulk =  c^2*rho;
frequency=2000; # Hz
omega=2*pi*frequency;
P_piston = 1.0e-3*phun("MPa"); # amplitude of the piston pressure
R = 50.0*phun("mm");# radius of the piston
Ro = 400.0*phun("mm"); # radius of the enclosure
nref = 2;#number of refinements of the sphere around the piston
nlayers = 48;                     # number of layers of elements surrounding the piston
tolerance = R/(2^nref)/10
dt=1.0/frequency/20;
tfinal=50*dt;
nsteps=round(tfinal/dt)+1;

println("""

Baffled piston in a half-sphere domain with hard boundary condition.

Hexahedral mesh. Algorithm version.
""")

t0  =  time()

# Hexahedral mesh
fens,fes  =  H8sphere(R,nref);
bfes  =  meshboundary(fes)
# File  =   "baffledabc_boundary.vtk"
# vtkexportmesh(File, bfes.conn, fens.xyz, FinEtools.MeshExportModule.Q4)
#  @async run(`"paraview.exe" $File`)

l = selectelem(fens,bfes,facing = true,direction = [1.0 1.0  1.0])
ex(xyz, layer) = (R+layer/nlayers*(Ro-R))*xyz/norm(xyz)
fens1,fes1  =  H8extrudeQ4(fens, subset(bfes,l), nlayers, ex);
fens,newfes1,fes2 =  mergemeshes(fens1, fes1, fens, fes, tolerance)
fes = cat(newfes1,fes2)

# Piston surface mesh
bfes  =  meshboundary(fes)
l1 = selectelem(fens, bfes, facing = true, direction = [-1.0 0.0 0.0])
l2 = selectelem(fens, bfes, distance = R, from = [0.0 0.0 0.0], inflate = tolerance)
piston_fes = subset(bfes,intersect(l1,l2));
# File  =   "piston_fes.vtk"
# vtkexportmesh(File, piston_fes.conn, fens.xyz, FinEtools.MeshExportModule.Q4)
#  @async run(`"paraview.exe" $File`)

# Outer spherical boundary
louter = selectelem(fens, bfes, facing = true, direction = [1.0 1.0  1.0], dotmin= 0.001)
outer_fes = subset(bfes,louter);
# File  =   "outer_fes.vtk"
# vtkexportmesh(File, outer_fes.conn, fens.xyz, FinEtools.MeshExportModule.Q4)
#  @async run(`"paraview.exe" $File`)

println("Pre-processing time elapsed  =  ",time() - t0,"s")

t1  =  time()

material = MatAcoustFluid(bulk, rho)
# Region of the fluid
femm = FEMMAcoust(IntegData(fes, GaussRule(3, 2)), material)

geom = NodalField(fens.xyz);
P = NodalField(zeros(nnodes(geom),1));
piston_fenids = connectednodes(piston_fes);
setebc!(P, piston_fenids, true, 1, 0.0);
applyebc!(P);
numberdofs!(P);

S  =  acousticstiffness(femm, geom, P);
C  =  acousticmass(femm, geom, P);
D= (2.0/dt)*S +(dt/2.0)*C;

P0 = deepcopy(P)
Pdd0 = deepcopy(P)
P1 = deepcopy(P)
Pdd1 = deepcopy(P)
TMPF = deepcopy(P)
P = nothing # we don't need this field anymore
vP0 = zeros(P0.nfreedofs)
vP1 = deepcopy(vP0)
vPd0 = deepcopy(vP0)
vPd1 = deepcopy(vP0)

t = 0.0
P0.fixed_values[piston_fenids,1] = P_piston*sin(omega*t)
Pdd0.fixed_values[piston_fenids,1] = P_piston*(-omega^2)*sin(omega*t)
vP0 = gathersysvec!(P0, vP0)

nh = selectnode(fens, nearestto = [R+Ro/2, 0.0, 0.0] )
Pnh = [P1.values[nh, 1][1]]
ts = [t]
step =0;
while t <=tfinal
    step = step  +1;
    t=t+dt;
    P1.fixed_values[piston_fenids,1] = P_piston*sin(omega*t)
    Pdd1.fixed_values[piston_fenids,1] = P_piston*(-omega^2)*sin(omega*t)
    TMPF.fixed_values = P0.fixed_values + P1.fixed_values
    F = nzebcloadsacousticmass(femm, geom, TMPF);
    TMPF.fixed_values = Pdd0.fixed_values + Pdd1.fixed_values
    F = F + nzebcloadsacousticstiffness(femm, geom, TMPF);
    # println("$(norm(F))")
    vPd1 = D\((2/dt)*(S*vPd0) - C*(2*vP0+(dt/2)*vPd0) + F);
    vP1 = vP0 + (dt/2)*(vPd0+vPd1);
    scattersysvec!(P1, vP1); # store current pressure
    # Swap variables for the next step
    copyto!(vP0, vP1)
    copyto!(vPd0, vPd1)
    copyto!(P0, P1)
    copyto!(Pdd0, Pdd1)
    # Graphics output
    File  =   "baffled_piston-$(step).vtk"
    vtkexportmesh(File, fes.conn, fens.xyz, FinEtools.MeshExportModule.H8;
        scalars = [(  "P", P1.values)])
    #  @async run(`"paraview.exe" $File`)
    push!(Pnh, P1.values[nh, 1][1])
    push!(ts, t)
    println("step = $( step )")
end

set(axis="normal", plotstyle="linespoints", linewidth=4, pointsize = 2, color = "black", xlabel = "Time", ylabel = "Pnh", grid="on", title = "")
f = figure()
plot(vec(ts), vec(Pnh), legend = "", marker = "edmd")
figure(f)

# using Plots
# plotly()
# plot(vec(ts), vec(Pnh))
# gui()
println("$(Pnh[end])")

# # Surface for the ABC
# abc1  =  FDataDict("femm"=>FEMMAcoustSurf(IntegData(outer_fes, GaussRule(2, 2)),
#           material))
#
# # Surface of the piston
# flux1  =  FDataDict("femm"=>FEMMAcoustSurf(IntegData(piston_fes, GaussRule(2, 2)),
#           material),  "normal_flux"=> -rho*a_piston+0.0im);
#
# # Make model data
# modeldata =  FDataDict("fens"=>  fens,
#                  "omega"=>omega,
#                  "regions"=>[region1],
#                  "flux_bcs"=>[flux1], "ABCs"=>[abc1])
#
# # Call the solver
# modeldata = FinEtools.AlgoAcoustModule.steadystate(modeldata)
#
# println("Computing time elapsed  =  ",time() - t1,"s")
# println("Total time elapsed  =  ",time() - t0,"s")
#
# geom = modeldata["geom"]
# P = modeldata["P"]
#
# File  =   "baffledabc.vtk"
# vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.H8;
# scalars = [("absP", abs.(P.values))])
# @async run(`"paraview.exe" $File`)
#
# # using Winston
# # pl  =  FramedPlot(title = "Matrix",xlabel = "x",ylabel = "Re P, Im P")
# # setattr(pl.frame, draw_grid = true)
# # add(pl, Curve([1:length(C[:])],vec(C[:]), color = "blue"))
#
# # # pl = plot(geom.values[nLx,1][ix],scalars[nLx][ix])
# # # xlabel("x")
# # # ylabel("Pressure")
# # display(pl)
#
# true
end # baffled_piston_H8_transient

function allrun()
    println("#####################################################")
    println("# baffled_piston_H8_transient ")
    baffled_piston_H8_transient()
    return true
end # function allrun

end # module baffled_piston_examples


text:
#=
Make a table of latex substitution shortcuts for characters in ./LaText Substitutions.plist.
This can be dragged into the System Preferences > Keyboard > Text pane on Mac OS.
Uses julia latex symbols definitions in dictionary latex_symbols.jl

Simon de Szoeke
adapted from  UnicodeCompletion/latex_symbols.jl
=#

include(joinpath(JULIA_HOME, "..", "share", "julia", "base", "latex_symbols.jl"));

φ = open("LaText Substitutions.plist", "w") # that's \varphi

# plist header copied from dragging Text Shortcuts out of System Preferences > Keyboard > Text pane
println(φ, """<?xml version="1.0" encoding="UTF-8"?>""")
println(φ, """<!DOCTYPE plist PUBLIC "-//Apple//DTD PLIST 1.0//EN" "http://www.apple.com/DTDs/PropertyList-1.0.dtd">""")
println(φ, """<plist version="1.0">""")
println(φ, """<array>""")
# write each plist dictionary entry
for (ω, (α, β)) in enumerate(latex_symbols)
    println(φ, "   <dict>")
    println(φ, "      <key>phrase</key>")
    println(φ, "      <string>$β</string>")
    println(φ, "      <key>shortcut</key>")
    println(φ, "      <string>$α</string>")
    println(φ, "   </dict>")
end
# plist footer
println(φ, """</array>""")
println(φ, """</plist>""")

close(φ)

max_stars_repo_name: deszoeke/LaTeXSusbstitutionsPlist

max_issues_repo_name: deszoeke/LaTeXSusbstitutionsPlist

max_forks_repo_name: deszoeke/LaTeXSusbstitutionsPlist


text:
include("1.jl")
function angleofsquare4(r::Robot)
    angleofsquare(r::Robot)
    for side in (Nord,West,Sud,Ost)
        movehole(r,side)
        putmarker!(r)
    end

end

function angleofsquare(r::Robot)
    while ((!isborder(r,Sud)) || (!isborder(r,Ost)))
        if (!isborder(r,Sud))
            move!(r,Sud)
        end
        if (!isborder(r,Ost))
            move!(r,Ost)
        end
    end
end
function movehole(r::Robot,side::HorizonSide)
    while (!isborder(r,side))
        move!(r,side)
        end
    end

max_stars_repo_name: bubasooo/proggggggggg

max_issues_repo_name: bubasooo/proggggggggg

max_forks_repo_name: bubasooo/proggggggggg


text:
__precompile__()

module Blink

using Reexport
using Distributed: Future
using Sockets
using Logging
using Base64: stringmime
using RelocatableFolders

include("rpc/rpc.jl")
include("content/content.jl")

include("AtomShell/AtomShell.jl")
export AtomShell
@reexport using .AtomShell
import .AtomShell: resolve_blink_asset

end # module


text:
using SimilaritySearch, Random
using Test

#
# This file contains a set of tests for SearchGraph over databases of vectors (of Float32)
#

@testset "vector indexing with SearchGraph" begin
    # NOTE: The following algorithms are complex enough to say we are testing it doesn't have syntax errors, a more grained test functions are requiered
    Random.seed!(0)
    ksearch = 10
    n, m, dim = 100_000, 100, 8

    db = MatrixDatabase(rand(Float32, dim, n))
    queries = MatrixDatabase(rand(Float32, dim, m))

    dist = SqL2Distance()
    seq = ExhaustiveSearch(dist, db)
    goldtime = @elapsed goldI, goldD = searchbatch(seq, queries, ksearch)

    for bsize in [2, 12]
        search_algo = BeamSearch(; bsize)
        @info "=================== $search_algo"
        graph = SearchGraph(; db=DynamicMatrixDatabase(Float32, dim), dist, search_algo=search_algo)
        neighborhood = Neighborhood(reduce=IdentityNeighborhood())
        append!(graph, db; neighborhood, parallel_block=8, callbacks=SearchGraphCallbacks(hyperparameters=nothing))
        searchtime = @elapsed I, D= searchbatch(graph, queries, ksearch)
        #@info sort!(length.(graph.links), rev=true)
        @show goldD[:, 1]
        @show D[:, 1]
        @show goldI[:, 1]
        @show I[:, 1]
        recall = macrorecall(goldI, I)
        @info "testing search_algo: $(string(graph.search_algo)), time: $(searchtime)"
        @test recall >= 0.6
        @info "queries per second: $(m/searchtime), recall: $recall"
        @info "===="
    end

    @info "--- Optimizing parameters ParetoRadius ---"
    graph = SearchGraph(; dist, search_algo=BeamSearch(bsize=2), verbose=false)
    neighborhood = Neighborhood(reduce=DistalSatNeighborhood())
    append!(graph, db; neighborhood, callbacks=SearchGraphCallbacks(ParetoRadius()))
    @info "---- starting ParetoRadius optimization ---"
    optimize!(graph, ParetoRadius())
    searchtime = @elapsed I, D = searchbatch(graph, queries, ksearch)
    recall = macrorecall(goldI, I)
    @info "ParetoRadius:> queries per second: ", m/searchtime, ", recall:", recall
    @info graph.search_algo
    @test recall >= 0.3  # we don't expect high quality results on ParetoRadius

    @info "---- starting ParetoRecall optimization ---"
    optimize!(graph, ParetoRecall())
    searchtime = @elapsed I, D = searchbatch(graph, queries, ksearch)
    recall = macrorecall(goldI, I)
    @info "ParetoRecall:> queries per second: ", m/searchtime, ", recall:", recall
    @info graph.search_algo
    @test recall >= 0.6

    @info "========================= REBUILD TEST ======================"
    graph = SearchGraph(; db, dist, verbose=false)
    index!(graph; callbacks=SearchGraphCallbacks(MinRecall(0.9)))
    optimize!(graph, MinRecall(0.9); queries)
    searchtime = @elapsed I, D = searchbatch(graph, queries, ksearch)
    recall = macrorecall(goldI, I)
    @info "testing without additional optimizations: queries per second:", m/searchtime, ", recall: ", recall
    @info graph.search_algo
    @test recall >= 0.6

    @info "========================= rebuild process =========================="
    graph = rebuild(graph)
    optimize!(graph, MinRecall(0.9); queries)  # using the actual dataset makes prone to overfitting hyperparameters (more noticeable in rebuilt indexes)
    searchtime_ = @elapsed I, D = searchbatch(graph, queries, ksearch)
    recall_ = macrorecall(goldI, I)

    @info "-- old vs rebuild> searchtime: $searchtime vs $searchtime_; recall: $recall vs $recall_"

    @info "#############=========== Default parameters (useful as fast benchmark) ==========###########"
    dim = 4
    db = MatrixDatabase(randn(Float32, dim, n))
    queries = MatrixDatabase(randn(Float32, dim, m))
    seq = ExhaustiveSearch(; dist, db)
    goldI, goldD = searchbatch(seq, queries, ksearch)
    graph = SearchGraph(; db, dist, verbose=false)
    index!(graph)
    searchtime = @elapsed I, D = searchbatch(graph, queries, ksearch)
    recall = macrorecall(goldI, I)
    @info "testing without additional optimizations> queries per second:", m/searchtime, ", recall: ", recall
    @info graph.search_algo
    @test recall >= 0.7
end

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/SimilaritySearch

max_forks_repo_name: sadit/SimilaritySearch


text:
struct PDBijector <: Bijector{2} end

# This function has custom adjoints defined for Tracker, Zygote and ReverseDiff.
# I couldn't find a mutation-free implementation that maintains TrackedArrays in Tracker
# and ReverseDiff, hence the need for custom adjoints.
function replace_diag(f, X)
    g(i, j) = ifelse(i == j, f(X[i, i]), X[i, j])
    return g.(1:size(X, 1), (1:size(X, 2))')
end
function (b::PDBijector)(X::AbstractMatrix{<:Real})
    Y = cholesky(X; check = true).L
    return replace_diag(log, Y)
end
function (ib::Inverse{<:PDBijector})(Y::AbstractMatrix{<:Real})
    X = replace_diag(exp, Y)
    return LowerTriangular(X) * LowerTriangular(X)'
end


text:
"""
    div_rem(x, y)

Computes `quotient = fl(x/y)` and `remainder = x - quotient*y` such that `(x-remainder) / y = quotient`.
"""
@inline function div_rem(x::T, y::T) where {T}
   quotient = x / y
   remainder = fma(-y, quotient, x)
   return quotient, remainder
end

"""
    sqrt_rem(x)

Computes `root = fl(sqrt(x))` and `remainder = x - root*root` such that `sqrt(x-remainder) = root`.
"""
@inline function sqrt_rem(x::T) where {T}
   root = sqrt(x)
   remainder = fma(-root, root, x)
   return root, remainder
end

#=
   Those remainders are exact.
   For more informations, applications and additional references, see the introduction of
   Sylvie Boldo and Jean-Michel Muller
   Some Functions Computable with a Fused-mac
=#


text:
module Lynx

import Gtk
import Gtk.GtkBox
import Gtk.GtkButton
import Gtk.GtkCanvas
import Gtk.GtkCheckButton
import Gtk.GtkColorButton
import Gtk.GtkComboBox
import Gtk.GtkComboBoxText
import Gtk.GtkEntry
import Gtk.GtkFileChooserDialog
import Gtk.GtkFrame
import Gtk.GtkGLArea
import Gtk.GtkImage
import Gtk.GtkLabel
import Gtk.GtkLinkButton
import Gtk.GtkMessageDialog
import Gtk.GtkNotebook
import Gtk.GtkPaned
import Gtk.GtkProgressBar
import Gtk.GtkScale
import Gtk.GtkScrolledWindow
import Gtk.GtkSpinButton
import Gtk.GtkSpinner
import Gtk.GtkSwitch
import Gtk.GtkToggleButton
import Gtk.GtkTreeView
import Gtk.GtkVolumeButton
import Gtk.GtkWidget
import Gtk.GtkGrid
import Gtk.GtkWindow
import Gtk.GtkAdjustment
import Gtk.GtkExpander

import Gtk.GAccessor
import Gtk.GdkRGBA
import Gtk.GEnum

import Gtk.signal_connect
import Gtk.set_gtk_property!
import Gtk.get_gtk_property
import Gtk.@guarded
import Gtk.draw
import Gtk.width
import Gtk.height
import Gtk.GObject
import Gtk.destroy
import Gtk.signal_handler_disconnect
import Gtk.showall
import Gtk.Cairo

using Gtk: libgtk
using Gtk: libgdk
using Gtk: GConstants

using Colors
using Luxor
using MacroTools
using Observables

import Observables.ObserverFunction

const G = GAccessor

include("utils.jl")
include("extended.jl")
include("math.jl")
include("widgets.jl")
include("grid.jl")
include("props.jl")
include("canvas.jl")
include("events.jl")
include("layout.jl")
include("app.jl")
include("image.jl")

include("dialog.jl")
using .Dialog

# Submodules
export Dialog
# Constants
export Baseline
export Center
export End
export Fill
export Start
# Types
export AbstractLayout
export CanvasEvents
export CanvasOnly
export LynxApp
export MouseEvent
export ScrollEvent
export SideBar
export SymString
export Image
# Widgets
export Activable
export Box
export Button
export Canvas
export CheckBox
export ColorButton
export Container
export Dropdown
export Expander
export Frame
export GLArea
export Grid
export ImageView
export Input
export Label
export Notebook
export NullContainer
export Paned
export ProgressBar
export RangeWidget
export Scrolled
export Slider
export SpinButton
export Spinner
export Switch
export TextField
export ToggleButton
export TreeView
export Widget
export Window
export Picture
# Functions
export events
export framerate!
export getapp
export getprop
export init
export layout!
export loop!
export mapr
export middle
export onevent
export pulse!
export run!
export setapp
export setprop!
export start!
export stop!
export use!
export value
export value!
export gkey
export gwidget
export keyname
export onkeypress
export waitfor
export ondraw
export onupdate
export onmousedown
export onmouseup
export onmousemove
export onmousedrag
export onscroll
export span
export cspan
export rspan
export GridElement
export slice
export drawimage
export getpixels
export rotr
export rotl
export rotr!
export rotl!
# Macros
export @app
export @canvas
export @framerate
export @height
export @map
export @new
export @on
export @size
export @unpack
export @use
export @waitfor
export @width
export @window
export @key_str
export @showall
export @widget

# Trick to get intellisense on vs-code
if (false) include("../test/runtests.jl") end

end # module Lynx


text:
module Sherlock

# Usings and imports
using Blink
using GraphRecipes
using Interact
using InteractiveUtils
using LightGraphs
using Plots

# Exports: types.jl
export abstracttypes
export Detective
export ENTITY_TYPES
export functions
export inquire
export ismacro
export safeisabstract
export safeisfield
export safeisnotabstract
export types
export undefined

# Exports: utilities.jl
export functiontype_edges
export greet
export nothing_is_val
export trim_type
export typetype_edges
#pkgchildren, loaded_packages #irreprable error with these...

# Exports: visualizations.jl
export magnify
export MIMO
export sherlockplot
export sherlock_UI
export typetree

# Includes
include("types.jl")
include("utilities.jl")
include("visualizations.jl")

end # module


text:
import Dolo
import PyPlot
path = Dolo.pkg_path

filename = joinpath(path,"examples","models","rbc_dtcc_iid_ar1.yaml")
model = Dolo.yaml_import(filename)
@time dr = Dolo.time_iteration(model, verbose=true, maxit=10000)

n_steps=100
s0 = model.calibration[:states]
m0 = model.calibration[:exogenous]
index = findfirst(isequal(:z), model.symbols[:states])
bounds = [dr.grid_endo.min[index], dr.grid_endo.max[index]]

df = Dolo.tabulate(model, dr, :z, bounds, s0, m0)
df = Dolo.tabulate(model, dr, :z, s0, m0)

cn =[:i, :n]
cn[1]
Svalues = linspace(bounds[1], bounds[2], n_steps)
PyPlot.plot(Svalues, df[cn[2]])

for j in cn
  fig = PyPlot.figure(j,figsize=(3,3))
  fig
  PyPlot.plot(Svalues, df[j], label=j)
  PyPlot.legend()
  # PyPlot.xlabel('state = {} | mstate = {}'.format(state, i0))
end

plot_dr = Dolo.tabulate(model, dr, :k, bounds, s0, m0, cn)


text:
@test [combinations([])...] == []
@test [combinations(['a', 'b', 'c'])...] == [['a'],['b'],['c'],['a','b'],['a','c'],['b','c'],['a','b','c']]

@test [combinations("abc",3)...] == [['a','b','c']]
@test [combinations("abc",2)...] == [['a','b'],['a','c'],['b','c']]
@test [combinations("abc",1)...] == [['a'],['b'],['c']]
@test [combinations("abc",0)...] == [[]]
@test [combinations("abc",-1)...] == []

@test filter(x->iseven(x[1]),[combinations([1,2,3],2)...]) == Any[[2,3]]

# multiset_combinations
@test [multiset_combinations("aabc", 5)...] == Any[]
@test [multiset_combinations("aabc", 2)...] == Any[['a','a'],['a','b'],['a','c'],['b','c']]
@test [multiset_combinations("aabc", 1)...] == Any[['a'],['b'],['c']]
@test [multiset_combinations("aabc", 0)...] == Any[Char[]]
@test [multiset_combinations("aabc", -1)...] == Any[]
@test [multiset_combinations("", 1)...] == Any[]
@test [multiset_combinations("", 0)...] == Any[Char[]]
@test [multiset_combinations("", -1)...] == Any[]

# with_replacement_combinations
@test [with_replacement_combinations("abc", 2)...] == Any[['a','a'],['a','b'],['a','c'],
                                                              ['b','b'],['b','c'],['c','c']]
@test [with_replacement_combinations("abc", 1)...] == Any[['a'],['b'],['c']]
@test [with_replacement_combinations("abc", 0)...] == Any[Char[]]
@test [with_replacement_combinations("abc", -1)...] == Any[]
@test [with_replacement_combinations("", 1)...] == Any[]
@test [with_replacement_combinations("", 0)...] == Any[Char[]]
@test [with_replacement_combinations("", -1)...] == Any[]

#cool-lex iterator
@test_throws DomainError [CoolLexCombinations(-1, 1)...]
@test_throws DomainError [CoolLexCombinations(5, 0)...]
@test [CoolLexCombinations(4,2)...] == Vector[[1,2], [2,3], [1,3], [2,4], [3,4], [1,4]]
@test isa(iterate(CoolLexCombinations(1000, 20))[2], Combinatorics.CoolLexIterState{BigInt})

# multichoose
@test multichoose(0, 1) == 0
@test multichoose(0, 2) == 0
@test multichoose(0, 3) == 0
@test multichoose(1, 1) == 1
@test multichoose(1, 2) == 1
@test multichoose(1, 3) == 1
@test multichoose(2, 1) == 2
@test multichoose(2, 2) == 3
@test multichoose(2, 3) == 4
@test multichoose(3, 0) == 1
@test multichoose(3, 1) == 3
@test multichoose(3, 2) == 6
@test multichoose(3, 3) == 10
@test multichoose(10, 8) == 24310

# multicombinations
@test [multicombinations("abc", 0)...] == [Char[]]
@test [multicombinations("abc", 1)...] == [['a'], ['b'], ['c']]
@test [multicombinations("abc", 2)...] == [['a','a'], ['a','b'], ['a','c'], ['b','b'], ['b','c'], ['c','c']]
@test [multicombinations("abc", 3)...] == [['a','a','a'], ['a','a','b'], ['a','a','c'], ['a','b','b'],
    ['a','b','c'], ['a','c','c'], ['b','b','b'], ['b','b','c'], ['b','c','c'], ['c','c','c']]

# Power set
@test collect(powerset([])) == Any[[]]
@test collect(powerset(['a', 'b', 'c'])) == Any[[],['a'],['b'],['c'],['a','b'],['a','c'],['b','c'],['a','b','c']]
@test collect(powerset(['a', 'b', 'c'], 1)) == Any[['a'],['b'],['c'],['a','b'],['a','c'],['b','c'],['a','b','c']]
@test collect(powerset(['a', 'b', 'c'], 1, 2)) == Any[['a'],['b'],['c'],['a','b'],['a','c'],['b','c']]


text:
###############################################################################
### Regularizer - l1
###############################################################################
abstract type regularizer_l1 <: AbstractRegularizer end

g(pb::AbstractProblem{regularizer_l1}, x) = pb.λ * norm(x, 1)

softthresh(x, gamma) = sign(x) * max(0, abs(x)-gamma)
prox_αg(pb::AbstractProblem{regularizer_l1}, x, α) = softthresh.(x, α * pb.λ)

## regularization term for JuMP solve
function build_reg!(model, regtype::Type{regularizer_l1}, x, n)
    p = @variable model p[1:n]
    @constraint(model, [i = 1:n], 0 <= x[i]+p[i])
    @constraint(model, [i = 1:n], 0 <= -x[i]+p[i])
    return p
end

@inline function get_support(regtype::Type{regularizer_l1}, y; tol=1e-13)
    ind_null = SortedSet{Int}()
    for i in 1:length(y)
        abs(y[i])<tol && push!(ind_null, i)
    end
    return ind_null
end

###############################################################################
### Regularizer - TV
###############################################################################
abstract type regularizer_TV <: AbstractRegularizer end

g(pb::AbstractProblem{regularizer_TV}, x) = pb.λ* norm(x[1:end-1]-x[2:end], 1)

function prox_αg(pb::AbstractProblem{regularizer_TV}, x, α)
    res = zeros(length(x))
    TVdenoise!(res, x, α * pb.λ)
    return res
end

function build_reg!(model, regtype::Type{regularizer_TV}, x, n)
    p = @variable model p[1:n-1]
    @constraint(model, [i = 2:n], 0 <=  x[i]-x[i-1]+p[i-1])
    @constraint(model, [i = 2:n], 0 <= -x[i]+x[i-1]+p[i-1])
    return p
end

@inline function get_support(regtype::Type{regularizer_TV}, y; tol=1e-13)
    ind_null = SortedSet{Int}()
    for i in 1:length(y)-1
        abs(y[i] - y[i+1])<tol && push!(ind_null, i)
    end
    return ind_null
end

###############################################################################
### Regularizer - l1-2
###############################################################################
abstract type regularizer_l12{T} <: AbstractRegularizer end

function g(pb::AbstractProblem{regularizer_l12{T}}, x) where T
    res = 0.0

    n = size(x, 1)
    ngroups = Int(ceil(n/T))

    for i in 1:ngroups-1
        res += norm(x[T*i-3:T*i], 2)
    end
    res += norm(x[T*ngroups-3:end], 2)

    return pb.λ* res
end

function prox_αg(pb::AbstractProblem{regularizer_l12{T}}, x, α) where T
    res = copy(x)

    n = size(x, 1)
    ngroups = Int(ceil(n/T))

    for i in 1:ngroups-1
        x_gnorm = norm(x[T*i-3:T*i], 2)
        if x_gnorm >= α*pb.λ
            res[T*i-3:T*i] *= (x_gnorm - α*pb.λ) / x_gnorm
        else
            res[T*i-3:T*i] .= 0
        end
    end

    x_gnorm = norm(x[T*ngroups-3:end], 2)
    if x_gnorm >= α*pb.λ
        res[T*ngroups-3:end] *= (x_gnorm - α*pb.λ) / x_gnorm
    else
        res[T*ngroups-3:end] .= 0
    end

    return res
end

# function build_reg!(model, regtype::Type{regularizer_l12}, x, n)
#     p = @variable model p[1:n-1]
#     @constraint(model, [i = 2:n], 0 <=  x[i]-x[i-1]+p[i-1])
#     @constraint(model, [i = 2:n], 0 <= -x[i]+x[i-1]+p[i-1])
#     return p
# end

@inline function get_support(regtype::Type{regularizer_l12{T}}, y; tol=1e-13) where T
    ind_null = SortedSet{Int}()

    n = size(y, 1)
    ngroups = Int(ceil(n/T))

    for i in 1:ngroups-1
        norm(y[T*i-3:T*i], 2)<tol && push!(ind_null, i)
    end
    norm(y[T*ngroups-3:end], 2)<tol && push!(ind_null, ngroups)

    return ind_null
end

###############################################################################
### Regularizer - l_nuclear
###############################################################################
abstract type regularizer_lnuclear <: AbstractRegularizer end

function g(pb::AbstractProblem{regularizer_lnuclear}, x)
    n = Int(sqrt(size(x, 1)))
    x_mat = reshape(x, n, n)

    return pb.λ* norm(svdvals(x_mat), 1)
end

function prox_αg(pb::AbstractProblem{regularizer_lnuclear}, x, α)
    n = Int(sqrt(size(x, 1)))
    x_mat = reshape(x, n, n)

    F = svd(x_mat)
    res = vec(F.U * Diagonal(softthresh.(F.S, α*pb.λ)) * F.Vt)

    return res
end

# function build_reg!(model, regtype::Type{regularizer_lnuclear}, x, n)
#     p = @variable model p[1:n-1]
#     @constraint(model, [i = 2:n], 0 <=  x[i]-x[i-1]+p[i-1])
#     @constraint(model, [i = 2:n], 0 <= -x[i]+x[i-1]+p[i-1])
#     return p
# end

@inline function get_support(regtype::Type{regularizer_lnuclear}, y; tol=1e-13)
    n = Int(sqrt(size(y, 1)))
    y_mat = reshape(y, n, n)

    return SortedSet{Int}(1:n - rank(y_mat))
end

###############################################################################
### Regularizer - dist(⋅, B(0,1)) = max(||⋅||-1, 0)
###############################################################################
abstract type regularizer_distball{p} <: AbstractRegularizer end

g(pb::AbstractProblem{regularizer_distball{p}}, x) where p = pb.λ* max(norm(x, p)-1, 0)

function prox_αg(pb::AbstractProblem{regularizer_distball{p}}, x, α) where p
    xnorm = norm(x, p)
    if xnorm <= 1
        res = copy(x)
    elseif xnorm > 1+α*pb.λ
        res = x .* (1 - α*pb.λ / xnorm)
    else
        res = x ./ xnorm
    end
    return res
end

function build_reg!(model, regtype::Type{regularizer_distball{p}}, x, n) where p
    regval = @variable(model, regval)
    @constraint(model, 0 <= regval)
    @constraint(model, 0 <= sum(x[i]^2 for i in 1:n) - regval^2)
    return regval
end

@inline function get_support(regtype::Type{regularizer_distball{p}}, y; tol=1e-11) where p
    ind_null = SortedSet{Int}()
    abs(norm(y, p)-1) <= tol && push!(ind_null, 1)
    return ind_null
end


text:
# CUDAnative run-time library
#
# This module defines method instances that will be compiled into a device-specific image
# and will be available to the CUDAnative compiler to call after Julia has generated code.

module Runtime

using ..CUDAnative
using LLVM

## representation of a runtime method instance

struct RuntimeMethodInstance
    def::Function

    return_type::Type
    types::Tuple
    name::Symbol

    # LLVM types cannot be cached, so we can't put them in the runtime method instance.
    # the actual types are constructed upon accessing them, based on a sentinel value:
    #  - nothing: construct the LLVM type based on its Julia counterparts
    #  - function: call this generator to get the type (when more control is needed)
    llvm_return_type::Union{Nothing, Function}
    llvm_types::Union{Nothing, Function}
    llvm_name::String
end

function Base.getproperty(rt::RuntimeMethodInstance, field::Symbol)
    value = getfield(rt, field)
    if field == :llvm_types
        if value == nothing
            LLVMType[convert.(LLVMType, typ) for typ in rt.types]
        else
            value()
        end
    elseif field == :llvm_return_type
        if value == nothing
            convert.(LLVMType, rt.return_type)
        else
            value()
        end
    else
        return value
    end
end

const methods = Dict{Symbol,RuntimeMethodInstance}()
get(name::Symbol) = methods[name]

# Register a Julia function `def` as a runtime library function identified by `name`. The
# function will be compiled upon first use for argument types `types` and should return
# `return_type`. Use `Runtime.get(name)` to get a reference to this method instance.
#
# The corresponding LLVM types `llvm_types` and `llvm_return_type` will be deduced from
# their Julia counterparts. To influence that conversion, pass a callable object instead;
# this object will be evaluated at run-time and the returned value will be used instead.
#
# When generating multiple runtime functions from a single definition, make sure to specify
# different values for `name`. The LLVM function name will be deduced from that name, but
# you can always specify `llvm_name` to influence that. Never use an LLVM name that starts
# with `julia_` or the function might clash with other compiled functions.
function compile(def, return_type, types, llvm_return_type=nothing, llvm_types=nothing;
                 name=typeof(def).name.mt.name, llvm_name="ptx_$name")
    meth = RuntimeMethodInstance(def,
                                 return_type, types, name,
                                 llvm_return_type, llvm_types, llvm_name)
    if haskey(methods, name)
        error("Runtime function $name has already been registered!")
    end
    methods[name] = meth
    meth
end

## exception handling

function report_exception(ex)
    @cuprintf("""
        ERROR: a %s exception occurred during kernel execution.
               Run Julia on debug level 2 for device stack traces.
        """, ex)
    return
end

compile(report_exception, Nothing, (Ptr{Cchar},))

function report_exception_name(ex)
    @cuprintf("""
        ERROR: a %s exception occurred during kernel execution.
        Stacktrace:
        """, ex)
    return
end

function report_exception_frame(idx, func, file, line)
    @cuprintf(" [%i] %s at %s:%i\n", idx, func, file, line)
    return
end

compile(report_exception_frame, Nothing, (Cint, Ptr{Cchar}, Ptr{Cchar}, Cint))
compile(report_exception_name, Nothing, (Ptr{Cchar},))

## GC

@enum AddressSpace begin
    Generic         = 1
    Tracked         = 10
    Derived         = 11
    CalleeRooted    = 12
    Loaded          = 13
end

# LLVM type of a tracked pointer
function T_prjlvalue()
    T_pjlvalue = convert(LLVMType, Any, true)
    LLVM.PointerType(eltype(T_pjlvalue), Tracked)
end

function gc_pool_alloc(sz::Csize_t)
    ptr = malloc(sz)
    return unsafe_pointer_to_objref(ptr)
end

compile(gc_pool_alloc, Any, (Csize_t,), T_prjlvalue)

end

max_stars_repo_name: Cyvadra/CUDAnative

max_issues_repo_name: Cyvadra/CUDAnative

max_forks_repo_name: Cyvadra/CUDAnative


text:
# Everything related to transformation on sparse matrices

function scale_row!(A::SMat{T}, i::Int, c::T) where T
  for j in 1:length(A.rows[i].values)
    A.rows[i].values[j] *= c
  end
end

doc"""
  swap_rows!{T}(A::SMat{T}, i::Int, j::Int)

  swaps, inplace, the i-th row and the j-th
"""
function swap_rows!(A::SMat{T}, i::Int, j::Int) where T
  A[i], A[j] = A[j], A[i]
end

doc"""
    invert_rows!{T}(A::SMat{T})

> Inplace, inverts the rows, ie. swaps the last and the 1st, the 2nd last and the
> 2nd, ...
"""
function invert_rows!(A::SMat{T}) where T
  for i=1:div(A.r, 2)
    A[i], A[A.r+1-i] = A[A.r+1-i], A[i]
  end
end

doc"""
    swap_cols!{T}(A::SMat{T}, i::Int, j::Int)

> Swap the i-th and j-th column inplace.
"""
function swap_cols!(A::SMat{T}, i::Int, j::Int) where T
  @assert 1 <= i <= cols(A) && 1 <= j <= cols(A)

  if i == j
    return nothing
  end

  for r in A.rows
    if i in r.pos
      i_i = findfirst(r.pos, i)
      val_i = r.values[i_i]
      if j in r.pos
        i_j = findfirst(r.pos, j)
        val_j = r.values[i_j]

        r.values[i_i], r.values[i_j] = r.values[i_j], r.values[i_i]
      else
        t = r.values[i_i]
        deleteat!(r.pos, i_i)
        deleteat!(r.values, i_i)
        k = searchsortedfirst(r.pos, j)
        insert!(r.pos, k, j)
        insert!(r.values, k, t)
      end
    else
      if j in r.pos
        i_j = findfirst(r.pos, j)
        val_j = r.values[i_j]

        t = r.values[i_j]
        deleteat!(r.pos, i_j)
        deleteat!(r.values, i_j)
        k = searchsortedfirst(r.pos, i)
        insert!(r.pos, k, i)
        insert!(r.values, k, t)
      end
    end
  end
  return nothing
end

# rows j -> row i*c + row j
doc"""
  add_scaled_row!{T}(A::SMat{T}, i::Int, j::Int, c::T)

  adds, inplace, the c*i-th row to the j-th
"""
function add_scaled_row!(A::SMat{T}, i::Int, j::Int, c::T) where T
  A.nnz = A.nnz - length(A[j])
  A.rows[j] = add_scaled_row(A[i], A[j], c)
  A.nnz = A.nnz + length(A[j])
end

function add_scaled_row(Ai::SRow{T}, Aj::SRow{T}, c::T) where T
  sr = SRow{T}()
  pi = 1
  pj = 1
  @assert c != 0
  while pi <= length(Ai.pos) && pj <= length(Aj.pos)
    if Ai.pos[pi] < Aj.pos[pj]
      push!(sr.pos, Ai.pos[pi])
      push!(sr.values, c*Ai.values[pi])
      pi += 1
    elseif Ai.pos[pi] > Aj.pos[pj]
      push!(sr.pos, Aj.pos[pj])
      push!(sr.values, Aj.values[pj])
      pj += 1
    else
      n = c*Ai.values[pi] + Aj.values[pj]
      if n != 0
        push!(sr.pos, Ai.pos[pi])
        push!(sr.values, n)
      end
      pi += 1
      pj += 1
    end
  end
  while pi <= length(Ai.pos)
    push!(sr.pos, Ai.pos[pi])
    push!(sr.values, c*Ai.values[pi])
    pi += 1
  end
  while pj <= length(Aj.pos)
    push!(sr.pos, Aj.pos[pj])
    push!(sr.values, Aj.values[pj])
    pj += 1
  end
  return sr
end

# col j -> col i*c + col j
doc"""
    add_scaled_col!{T}(A::SMat{T}, i::Int, j::Int, c::T)

> Adds, inplace, the c*i-th column to the j-th column.
"""
function add_scaled_col!(A::SMat{T}, i::Int, j::Int, c::T) where T
  @assert c != 0

  @assert 1 <= i <= cols(A) && 1 <= j <= cols(A)

  for r in A.rows
    if i in r.pos
      i_i = findfirst(r.pos, i)
      val_i = r.values[i_i]
      if j in r.pos
        i_j = findfirst(r.pos, j)
        val_j = r.values[i_j]

        r.values[i_j] += c*r.values[i_i]
      else
        k = searchsortedfirst(r.pos, j)
        insert!(r.pos, k, j)
        insert!(r.values, k, c*r.values[i_i])
      end
    end
  end
  return nothing
end

# row i -> a*row i + b * row j
# row j -> c*row i + d * row j
doc"""
  transform_row!{T}(A::SMat{T}, i::Int, j::Int, a::T, b::T, c::T, d::T)

  Inplace, replaces the i-th row and the j-th row by
  [a,b; c,d] * [i-th-row ; j-th row]
"""
function transform_row!(A::SMat{T}, i::Int, j::Int, a::T, b::T, c::T, d::T) where T
  A.nnz = A.nnz - length(A[i]) - length(A[j])
  A.rows[i], A.rows[j] = transform_row(A[i], A[j], a, b, c, d)
  A.nnz = A.nnz + length(A[i]) + length(A[j])
end
function transform_row(Ai::SRow{T}, Aj::SRow{T}, a::T, b::T, c::T, d::T) where T
  sr = SRow{T}()
  tr = SRow{T}()
  pi = 1
  pj = 1
  while pi <= length(Ai) && pj <= length(Aj)
    if Ai.pos[pi] < Aj.pos[pj]
      if a != 0
        push!(sr.pos, Ai.pos[pi])
        push!(sr.values, a*Ai.values[pi])
      end
      if c != 0
        push!(tr.pos, Ai.pos[pi])
        push!(tr.values, c*Ai.values[pi])
      end
      pi += 1
    elseif Ai.pos[pi] > Aj.pos[pj]
      if b != 0
        push!(sr.pos, Aj.pos[pj])
        push!(sr.values, b*Aj.values[pj])
      end
      if d != 0
        push!(tr.pos, Aj.pos[pj])
        push!(tr.values, d*Aj.values[pj])
      end
      pj += 1
    else
      m = a*Ai.values[pi] + b*Aj.values[pj]
      n = c*Ai.values[pi] + d*Aj.values[pj]
      if m != 0
        push!(sr.pos, Ai.pos[pi])
        push!(sr.values, m)
      end
      if n != 0
        push!(tr.pos, Ai.pos[pi])
        push!(tr.values, n)
      end
      pi += 1
      pj += 1
    end
  end
  while pi <= length(Ai.pos)
    if a != 0
      push!(sr.pos, Ai.pos[pi])
      push!(sr.values, a*Ai.values[pi])
    end
    if c != 0
      push!(tr.pos, Ai.pos[pi])
      push!(tr.values, c*Ai.values[pi])
    end
    pi += 1
  end
  while pj <= length(Aj.pos)
    if b != 0
      push!(sr.pos, Aj.pos[pj])
      push!(sr.values, b*Aj.values[pj])
    end
    if d != 0
      push!(tr.pos, Aj.pos[pj])
      push!(tr.values, d*Aj.values[pj])
    end
    pj += 1
  end

  return sr, tr
end

################################################################################
#
#  Inverses of elementary matrices
#
################################################################################

function inv(t::TrafoSwap)
  return T
end

function inv(t::TrafoAddScaled)
  return TrafoAddScaled(t.i, t.j, -t.s)
end

function inv(t::TrafoPartialDense)
  return TrafoPartialDense(t.i, t.rows, t.cols, inv(t.U))
end

# TrafoParaAddScaled is missing

################################################################################
#
#  Application of a transformation on the left side of a sparse matrix
#
################################################################################

# The following function do not update the number of nonzero entries
# properly

function apply_left!(A::SMat{T}, t::TrafoScale{T}) where T
  scale_row!(A, t.i, t.c)
  return nothing
end

function apply_left!(A::SMat{T}, t::TrafoSwap{T}) where T
  swap_rows!(A, t.i, t.j)
  return nothing
end

function apply_left!(A::SMat{T}, t::TrafoAddScaled{T}) where T
  add_scaled_row!(A, t.i, t.j, t.s)
  return nothing
end

function apply_left!(A::SMat{T}, t::TrafoParaAddScaled{T}) where T
  transform_row!(A, t.i, t.j, t.a, t.b, t.c, t.d)
  return nothing
end

function apply_left!(A::SMat{T}, t::TrafoDeleteZero{T}) where T
  deleteat!(A.rows, t.i)
  A.r -= 1
end

function apply_left!(A::SMat{T}, t::TrafoPartialDense{S}) where {T, S}
  R = parent(A.rows[1].values[1])
  i = t.i
  h = sub(A, t.rows, t.cols)
  deleteat!(A.rows, t.rows)
  A.r -= length(t.rows)
  @assert length(A.rows) == A.r
  # make dense matrix
  hdense = typeof(t.U)(h)

  #println(t.U)

  hdense = t.U * hdense

  h = _SMat(hdense, R = R, zerorows = true)

  for k in 1:length(h.rows)
    j = h.rows[k]
    rw = SRow{T}()
    for e in 1:length(j.pos)
      push!(rw.pos, j.pos[e] + i-1)
      push!(rw.values, j.values[e])
    end
    insert!(A.rows, i + k - 1, rw)
    A.r += 1
  end
  return nothing
end

################################################################################
#
#  Application of a transformation on the right side of a sparse matrix
#
################################################################################

# The following function do not update the number of nonzero entries
# properly
function apply_right!(A::SMat{T}, t::TrafoSwap{T}) where T
  swap_cols!(A, t.i, t.j)
  return nothing
end

function apply_right!(A::SMat{T}, t::TrafoAddScaled{T}) where T
  add_scaled_col!(A, t.j, t.i, t.s)
  return nothing
end

function apply_right!(A::SMat{T}, t::TrafoPartialDense{S}) where {T, S}
  # this works only if there are zeros left of the block to which we apply t
  i = t.i
  h = sub(A, t.rows, t.cols)
  deleteat!(A.rows, t.rows)
  A.r -= length(t.rows)
  @assert length(A.rows) == A.r
  # make dense matrix
  hdense = typeof(t.U)(h)

  #println(t.U)

  hdense = hdense * t.U

  h = _SMat(hdense, R = parent(A.rows[1].values[1]), zerorows = true)

  for k in 1:length(h.rows)
    j = h.rows[k]
    rw = SRow{T}()
    for e in 1:length(j.pos)
      push!(rw.pos, j.pos[e] + i-1)
      push!(rw.values, j.values[e])
    end
    insert!(A.rows, i + k - 1, rw)
    A.r += 1
  end
  return nothing
end

################################################################################
#
#  Application of a transformation on the right hand side of a row
#
################################################################################

function apply_right!(x::Array{T, 1}, t::TrafoAddScaled{T}) where T
  x[t.i] = x[t.i] + x[t.j]*t.s
  return nothing
end

function apply_right!(x::Array{T, 1}, t::TrafoScale{T}) where T
  x[t.i] = x[t.i] * t.c
end

function apply_right!(x::Array{T, 1}, t::TrafoSwap{T}) where T
  r = x[t.i]
  x[t.i] = x[t.j]
  x[t.j] = r
  return nothing
end

function apply_right!(x::Array{T, 1}, t::TrafoParaAddScaled{T}) where T
  r = t.a * x[t.i] + t.c * x[t.j]
  s = t.b * x[t.i] + t.d * x[t.j]
  x[t.i] = r
  x[t.j] = s
  return nothing
end

function apply_right!(x::Array{T, 1}, t::TrafoDeleteZero{T}) where T
  # move ith position to the back
  for j in length(x):-1:t.i+1
    r = x[j]
    x[j] = x[j - 1]
    x[j - 1] = r
  end
end

function apply_right!(x::Array{T, 1}, t::TrafoPartialDense) where T
  s = matrix(parent(x[1]), 1, rows(t.U), x[t.rows])
  #println("s :$s")
  s = s*t.U
  for (i, j) in zip(t.rows,1:rows(t.U))
    x[i] = s[1, j]
  end
  return nothing
end

################################################################################
#
#  Application to the left of an Array of ideals
#
################################################################################

function apply_left!(x::Vector{NfOrdFracIdl}, y::TrafoSwap)
  x[y.i], x[y.j] = x[y.j], x[y.i]
end

function apply_left!(x::Vector{NfOrdFracIdl}, y::TrafoAddScaled)
  x[y.j] = x[y.j] * x[y.i]^Int(y.s)
end

function apply_left!(x::Vector{NfOrdFracIdl}, y::TrafoPartialDense)
  z = view(deepcopy(x), y.cols)
  xx = view(x, y.cols)
  for i in 1:rows(y.U)  ## use power product instead
    xx[i] = z[1]^Int(y.U[i, 1])
    for j in 2:cols(y.U)
      xx[i] *= z[j]^Int(y.U[i, j])
    end
  end
end


text:
module Assignment09

export question1,
       question2,
       question3,
       question4,
       question5

include("assignment.jl")

end


text:
@doc """
Calculate a causal graph that is consistent with the data:

* `G` : graph of potential connections
* `df` : data frame of observed values for nodes in `G`
""" ->
function pc_algo(G::GenericGraph, df::DataFrame)
end


text:
using LowRankModels
import StatsBase: sample, Weights

# tests OvALoss

srand(1);
m,n,k = 200,50,2;
kfit = k+1
K = 4; # number of categories
d = n*K;
# matrix to encode
X_real, Y_real = randn(m,k), randn(k,d);
XY = X_real*Y_real;
# subtract the mean so we can compare the truth with the fit;
# the loss function is invariant under shifts
losses = fill(OvALoss(K),n)
yidxs = get_yidxs(losses)
for i=1:m
	for j=1:n
		mef = mean(XY[i,yidxs[j]])
		XY[i,yidxs[j]] = XY[i,yidxs[j]] - mef
	end
end

A = zeros(Int, (m, n))
for i=1:m
	for j=1:n
		wv = Weights(Float64[exp(-XY[i, K*(j-1) + l]) for l in 1:K])
		l = sample(wv)
		A[i,j] = l
	end
end

# and the model
losses = fill(OvALoss(K),n)
rx, ry = QuadReg(), QuadReg();
glrm = GLRM(A,losses,rx,ry,kfit, scale=false, offset=false, X=randn(kfit,m), Y=randn(kfit,d));

# fit w/o initialization
@time X,Y,ch = fit!(glrm);
XYh = X'*Y;
@show ch.objective
println("After fitting, parameters differ from true parameters by $(vecnorm(XY - XYh)/sqrt(prod(size(XY)))) in RMSE")
A_imputed = impute(glrm)
println("After fitting, $(sum(A_imputed .!= A) / prod(size(A))*100)\% of imputed entries are wrong")
println("(Picking randomly, $((K-1)/K*100)\% of entries would be wrong.)\n")

# initialize
init_svd!(glrm)
XYh = glrm.X' * glrm.Y
println("After initialization with the svd, parameters differ from true parameters by $(vecnorm(XY - XYh)/sqrt(prod(size(XY)))) in RMSE")
A_imputed = impute(glrm)
println("After initialization with the svd, $(sum(A_imputed .!= A) / prod(size(A))*100)\% of imputed entries are wrong")
println("(Picking randomly, $((K-1)/K*100)\% of entries would be wrong.)\n")

# fit w/ initialization
@time X,Y,ch = fit!(glrm);
XYh = X'*Y;
@show ch.objective
println("After fitting, parameters differ from true parameters by $(vecnorm(XY - XYh)/sqrt(prod(size(XY)))) in RMSE")
A_imputed = impute(glrm)
println("After fitting, $(sum(A_imputed .!= A) / prod(size(A))*100)\% of imputed entries are wrong")
println("(Picking randomly, $((K-1)/K*100)\% of entries would be wrong.)\n")


text:
using Documenter, CoordinatedSupplyChains

makedocs(sitename="CoordinatedSupplyChains.jl Documentation",
	pages = [
        	"Home" => "index.md",
        	"Tutorials" => "tutorial.md"
        	]
)

deploydocs(
    repo = "github.com/Tominapa/CoordinatedSupplyChains.jl.git",
)


text:
# Autogenerated wrapper script for Antlr4CppRuntime_jll for aarch64-apple-darwin
export libantlr4_cpp_runtime

JLLWrappers.@generate_wrapper_header("Antlr4CppRuntime")
JLLWrappers.@declare_library_product(libantlr4_cpp_runtime, "@rpath/libantlr4-runtime.4.9.3.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libantlr4_cpp_runtime,
        "lib/libantlr4-runtime.4.9.3.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
struct ContourLike end

contourlevels(spacing=1.7, n=12) = (spacing^i for i=0:(n-1))

# 1D plot
@recipe function f(A::NMRData{T,1}; normalize=true) where T
    Afwd = DimensionalData.forwardorder(A) # make sure data axes are in forwards order
    dim = dims(Afwd, 1)

    # force 1D to be a line plot
    seriestype := :path
    markershape --> :none

    # set default title
    title --> DimensionalData.refdims_title(Afwd)
    legend --> false

    xguide --> "$(Afwd[dim,:label]) chemical shift / ppm"
    xflip --> true
    xgrid --> false
    xtick_direction --> :out

    yguide --> ""
    yshowaxis --> false
    yticks --> nothing

    delete!(plotattributes, :normalize)
    val(dim), parent(Afwd)./(normalize ? scale(A) : 1)
end

# multiple 1D plots
@recipe function f(v::Vector{<:NMRData{T,1} where T}; normalize=true, vstack=false)
    # force 1D to be a line plot
    seriestype := :path
    markershape --> :none

    # get the first entry to determine axis label
    Afwd = DimensionalData.forwardorder(v[1])
    dim = dims(Afwd, 1)

    xguide --> "$(Afwd[dim,:label]) chemical shift / ppm"
    xflip --> true
    xgrid --> false
    xtick_direction --> :out

    yguide --> ""
    yshowaxis --> false
    yticks --> nothing

    delete!(plotattributes, :vstack)
    delete!(plotattributes, :normalize)

    voffset = 0
    vdelta = maximum([
            maximum(abs.(A)) / (normalize ? A[:ns]*A[:rg] : 1)
            for A in v]) / length(v)

    # TODO add guide lines
    # if vstack
    #     yticks --> voffset .+ (0:length(v)-1)*vdelta
    # else
    #     yticks --> [0,]
    # end

    for A in v
        @series begin
            seriestype := :path
            Afwd = DimensionalData.forwardorder(A) # make sure data axes are in forwards order
            dim = dims(Afwd, 1)
            label --> label(A)
            val(dim), parent(Afwd)./(normalize ? scale(A) : 1) .+ voffset
        end
        vstack && (voffset += vdelta)
    end
end

# 2D plot
@recipe f(d::D) where {D<:NMRData{T,2} where T} = SimpleTraits.trait(HasPseudoDimension{D}), d

@recipe function f(::Type{Not{HasPseudoDimension{D}}}, d::D) where {D<:NMRData{T, 2} where T}
    dfwd = DimensionalData.forwardorder(d) # make sure data axes are in forwards order
    dfwd = DimensionalData.maybe_permute(dfwd, (YDim, XDim))
    y, x = dims(dfwd)

    # set default title
    title --> DimensionalData.refdims_title(dfwd)
    legend --> false
    framestyle --> :box

    xguide --> "$(dfwd[x,:label]) chemical shift / ppm"
    xflip --> true
    xgrid --> false
    xtick_direction --> :out

    yguide --> "$(dfwd[y,:label]) chemical shift / ppm"
    yflip --> true
    ygrid --> false
    ytick_direction --> :out

    # generate light and dark colours for plot contours, based on supplied colour
    # - create a 5-tone palette with the same hue as the passed colour, and select the
    # fourth and second entries to provide dark and light shades
    basecolor = get(plotattributes, :linecolor, :blue)
    colors = sequential_palette(hue(convert(HSV,parse(Colorant, basecolor))),5)[[4,2]]

    #delete!(plotattributes, :normalize)
    #scale = normalize ? A[:ns]*A[:rg] : 1
    #val(dim), parent(A) ./ scale
    @series begin
        levels --> 5*dfwd[:noise].*contourlevels()
        linecolor := colors[1]
        primary := true
        val(x), val(y), data(dfwd)
    end
    @series begin
        levels --> -5*dfwd[:noise].*contourlevels()
        linecolor := colors[2]
        primary := false
        val(x), val(y), data(dfwd)
    end
end

@recipe function f(::Type{HasPseudoDimension{D}}, d::D) where {D<:NMRData{T, 2} where T}
    # TODO define recipe for pseudo2D data
    dfwd = DimensionalData.forwardorder(d) # make sure data axes are in forwards order
    dfwd = DimensionalData.maybe_permute(dfwd, (YDim, XDim))
    y, x = dims(dfwd)
    @info "plot recipe for pseudo2D data not yet defined"
    val(x), val(y), data(dfwd)
end

# multiple 2D plots
@recipe f(v::Vector{D}; normalize=true) where {D<:NMRData{T,2}} where {T} = SimpleTraits.trait(HasPseudoDimension{D}), v

@recipe function f(::Type{Not{HasPseudoDimension{D}}}, v::Vector{D}; normalize=true) where {D<:NMRData{T,2}} where {T}
    @info "plotting vector of 2D NMR data"
    n = length(v)
    hues = map(h->HSV(h,0.5,0.5), (0:n-1) .* (360/n))

    # force 1D to be a line plot
    seriestype := :contour

    # get the first entry to determine axis label
    dfwd = DimensionalData.forwardorder(v[1]) # make sure data axes are in forwards order
    dfwd = DimensionalData.maybe_permute(dfwd, (YDim, XDim))
    y, x = dims(dfwd)

    # set default title
    title --> DimensionalData.refdims_title(dfwd)
    legend --> false
    framestyle --> :box

    xguide --> "$(dfwd[x,:label]) chemical shift / ppm"
    xflip --> true
    xgrid --> false
    xtick_direction --> :out

    yguide --> "$(dfwd[y,:label]) chemical shift / ppm"
    yflip --> true
    ygrid --> false
    ytick_direction --> :out

    delete!(plotattributes, :normalize)

    h = 0.0
    for d in v
        dfwd = DimensionalData.forwardorder(d) # make sure data axes are in forwards order
        dfwd = DimensionalData.maybe_permute(dfwd, (YDim, XDim))
        y, x = dims(dfwd)
        label --> label(d)
        # generate light and dark colours for plot contours, based on supplied colour
        # - create a 5-tone palette with the same hue as the passed colour, and select the
        # fourth and second entries to provide dark and light shades
        colors = sequential_palette(h, 5)[[4,2]]

        @series begin
            levels --> 5*dfwd[:noise].*contourlevels()
            linecolor := colors[1]
            primary := true
            val(x), val(y), data(dfwd)
        end
        @series begin
            levels --> -5*dfwd[:noise].*contourlevels()
            linecolor := colors[2]
            primary := false
            val(x), val(y), data(dfwd)
        end
        h += 360.0/n
    end
end

@recipe function f(::Type{HasPseudoDimension{D}}, v::Vector{D}; normalize=true) where {D<:NMRData{T,2}} where {T}
    @info "plot recipe for vector of pseudo2D NMR data not yet implemented"
    delete!(plotattributes, :normalize)
    # TODO just make repeat calls to single plot recipe
    for d in v
        @series begin
            HasPseudoDimension{D}, d
        end
    end
end


text:
#enum cusparseStatus_t
#error messages from CUSPARSE

"""
Status messages from CUSPARSE's C API.
"""
const cusparseStatus_t = UInt32
const CUSPARSE_STATUS_SUCCESS                   = 0
const CUSPARSE_STATUS_NOT_INITIALIZED           = 1
const CUSPARSE_STATUS_ALLOC_FAILED              = 2
const CUSPARSE_STATUS_INVALID_VALUE             = 3
const CUSPARSE_STATUS_ARCH_MISMATCH             = 4
const CUSPARSE_STATUS_MAPPING_ERROR             = 5
const CUSPARSE_STATUS_EXECUTION_FAILED          = 6
const CUSPARSE_STATUS_INTERNAL_ERROR            = 7
const CUSPARSE_STATUS_MATRIX_TYPE_NOT_SUPPORTED = 8

#enum cusparseAction_t
"""
Perform operation on indices only (`CUSPARSE_ACTION_SYMBOLIC`) or
on both data and indices (`CUSPARSE_ACTION_NUMERIC`). Used in
conversion routines.
"""
const cusparseAction_t = UInt32
const CUSPARSE_ACTION_SYMBOLIC = 0
const CUSPARSE_ACTION_NUMERIC  = 1

#enum cusparseDirection_t
"""
Parse dense matrix by rows (`CUSPARSE_DIRECTION_ROW`) or columns
(`CUSPARSE_DIRECTION_COL`) to compute its number of non-zeros.
"""
const cusparseDirection_t = UInt32
const CUSPARSE_DIRECTION_ROW = 0
const CUSPARSE_DIRECTION_COL = 1

#enum cusparseHybPartition_t
"""
How to partition the HYB matrix in a [`CudaSparseMatrixHYB`](@ref).
There are three choices:
* `CUSPARSE_HYB_PARTITION_AUTO` - let CUSPARSE decide internally for best performance.
* `CUSPARSE_HYB_PARTITION_USER` - set the partition manually in the conversion function.
* `CUSPARSE_HYB_PARTITION_MAX` - use the maximum partition, putting the matrix in ELL format.
"""
const cusparseHybPartition_t = UInt32
const CUSPARSE_HYB_PARTITION_AUTO = 0
const CUSPARSE_HYB_PARTITION_USER = 1
const CUSPARSE_HYB_PARTITION_MAX  = 2

#enum cusparseFillMode_t
"""
Determines if a symmetric/Hermitian/triangular matrix has its upper
(`CUSPARSE_FILL_MODE_UPPER`) or lower (`CUSPARSE_FILL_MODE_LOWER`)
triangle filled.
"""
const cusparseFillMode_t = UInt32
const CUSPARSE_FILL_MODE_LOWER = 0
const CUSPARSE_FILL_MODE_UPPER = 1

#enum cusparseDiagType_t
"""
Determines if the diagonal of a matrix is all ones (`CUSPARSE_DIAG_TYPE_UNIT`)
or not all ones (`CUSPARSE_DIAG_TYPE_NON_UNIT`).
"""
const cusparseDiagType_t = UInt32
const CUSPARSE_DIAG_TYPE_NON_UNIT = 0
const CUSPARSE_DIAG_TYPE_UNIT     = 1

#enum cusparsePointerMode_t
"""
Determines if scalar arguments to a function are present on the host CPU
(`CUSPARSE_POINTER_MODE_HOST`) or on the GPU (`CUSPARSE_POINTER_MODE_DEVICE`).
"""
const cusparsePointerMode_t = UInt32
const CUSPARSE_POINTER_MODE_HOST   = 0
const CUSPARSE_POINTER_MODE_DEVICE = 1

#enum cusparseOperation_t
"""
Determines whether to perform an operation, such as a matrix multiplication
or solve, on the matrix as-is (`CUSPARSE_OPERATION_NON_TRANSPOSE`), on the
matrix's transpose (`CUSPARSE_OPERATION_TRANSPOSE`), or on its conjugate
transpose (`CUSPARSE_OPERATION_CONJUGATE_TRANSPOSE`).
"""
const cusparseOperation_t = UInt32
const CUSPARSE_OPERATION_NON_TRANSPOSE       = 0
const CUSPARSE_OPERATION_TRANSPOSE           = 1
const CUSPARSE_OPERATION_CONJUGATE_TRANSPOSE = 2

#enum cusparseMatrixType_t
"""
Indicates whether a matrix is a general matrix (`CUSPARSE_MATRIX_TYPE_GENERAL`),
symmetric (`CUSPARSE_MATRIX_TYPE_SYMMETRIC`), Hermitian
(`CUSPARSE_MATRIX_TYPE_HERMITIAN`), or triangular
(`CUSPARSE_MATRIX_TYPE_TRIANGULAR`). Note that for some matrix types
(those in [`CompressedSparse`](@ref)), this can be inferred for some function
calls.
"""
const cusparseMatrixType_t = UInt32
const CUSPARSE_MATRIX_TYPE_GENERAL    = 0
const CUSPARSE_MATRIX_TYPE_SYMMETRIC  = 1
const CUSPARSE_MATRIX_TYPE_HERMITIAN  = 2
const CUSPARSE_MATRIX_TYPE_TRIANGULAR = 3

#enum cusparseSolvePolicy_t
"""
Indicates whether to keep level info in solvers (`CUSPARSE_SOLVE_POLICY_USE_LEVEL`)
or whether to not use it (`CUSPARSE_SOLVE_POLICY_NO_LEVEL`).
"""
const cusparseSolvePolicy_t = UInt32
const CUSPARSE_SOLVE_POLICY_NO_LEVEL  = 0
const CUSPARSE_SOLVE_POLICY_USE_LEVEL = 1

#enum cusparseIndexBase_t
"""
Indicates whether a sparse object is zero-indexed (`CUSPARSE_INDEX_BASE_ZERO`)
or one-indexed (`CUSPARSE_INDEX_BASE_ONE`). CUSPARSE.jl supports both. Julia
sparse matrices are one-indexed, but you may wish to pass matrices from other
libraries which use zero-indexing (e.g. C language ODE solvers).
"""
const cusparseIndexBase_t = UInt32
const CUSPARSE_INDEX_BASE_ZERO = 0
const CUSPARSE_INDEX_BASE_ONE  = 1

#struct cusparseMatDescr_t
"""
Describes shape and properties of a CUSPARSE matrix. A convenience wrapper.

Contains:
* `MatrixType` - a [`cusparseMatrixType_t`](@ref)
* `FillMode` - a [`cusparseFillMode_t`](@ref)
* `DiagType` - a [`cusparseDiagType_t`](@ref)
* `IndexBase` - a [`cusparseIndexBase_t`](@ref)
"""
struct cusparseMatDescr_t
    MatrixType::cusparseMatrixType_t
    FillMode::cusparseFillMode_t
    DiagType::cusparseDiagType_t
    IndexBase::cusparseIndexBase_t
    function cusparseMatDescr_t(MatrixType,FillMode,DiagType,IndexBase)
        new(MatrixType,FillMode,DiagType,IndexBase)
    end
end

"""
An opaque struct containing information about the solution approach
CUSPARSE will take. Generated by [`sv_analysis`](@ref) or
[`sm_analysis`](@ref) and passed to [`sv_solve!`](@ref), [`sm_solve`](@ref),
[`ic0!`](@ref), or [`ilu0!`](@ref).
"""
const cusparseSolveAnalysisInfo_t = Ptr{Cvoid}
const bsrsm2Info_t = Ptr{Cvoid}
const bsrsv2Info_t = Ptr{Cvoid}
const csrsv2Info_t = Ptr{Cvoid}
const csric02Info_t = Ptr{Cvoid}
const csrilu02Info_t = Ptr{Cvoid}
const bsric02Info_t = Ptr{Cvoid}
const bsrilu02Info_t = Ptr{Cvoid}

const cusparseContext = Cvoid
const cusparseHandle_t = Ptr{cusparseContext}

#complex numbers

const cuComplex = Complex{Float32}
const cuDoubleComplex = Complex{Float64}

const CusparseFloat = Union{Float64,Float32,ComplexF64,ComplexF32}
const CusparseReal = Union{Float64,Float32}
const CusparseComplex = Union{ComplexF64,ComplexF32}

const cusparseHybMat_t = Ptr{Cvoid}


text:
struct DiffEqNLSolveTag end

struct DiffCache{T<:AbstractArray, S<:AbstractArray}
    du::T
    dual_du::S
end

Base.@pure function DiffCache(T, size, ::Type{Val{chunk_size}}) where chunk_size
    DiffCache(fill(zero(T), size...), fill(zero(Dual{typeof(ForwardDiff.Tag(DiffEqNLSolveTag(),T)),T,chunk_size}), size...))
end

Base.@pure DiffCache(u::AbstractArray) = DiffCache(eltype(u),size(u),Val{ForwardDiff.pickchunksize(length(u))})
Base.@pure DiffCache(u::AbstractArray,nlsolve) = DiffCache(eltype(u),size(u),Val{get_chunksize(nlsolve)})
Base.@pure DiffCache(u::AbstractArray,T::Type{Val{CS}}) where {CS} = DiffCache(eltype(u),size(u),T)

get_du(dc::DiffCache, ::Type{T}) where {T<:Dual} = dc.dual_du
get_du(dc::DiffCache, T) = dc.du

# Default nlsolve behavior, should be removed

Base.@pure determine_chunksize(u,alg::DEAlgorithm) = determine_chunksize(u,get_chunksize(alg))
Base.@pure function determine_chunksize(u,CS)
  if CS != 0
    return CS
  else
    return ForwardDiff.pickchunksize(length(u))
  end
end

struct NLSOLVEJL_SETUP{CS,AD} end
Base.@pure NLSOLVEJL_SETUP(;chunk_size=0,autodiff=true) = NLSOLVEJL_SETUP{chunk_size,autodiff}()
(::NLSOLVEJL_SETUP)(f,u0; kwargs...) = (res=NLsolve.nlsolve(f,u0; kwargs...); res.zero)
function (p::NLSOLVEJL_SETUP{CS,AD})(::Type{Val{:init}},f,u0_prototype) where {CS,AD}
  AD ? autodiff = :forward : autodiff = :central
  OnceDifferentiable(f, u0_prototype, u0_prototype, autodiff,
                     ForwardDiff.Chunk(determine_chunksize(u0_prototype,CS)))
end

get_chunksize(x) = 0
get_chunksize(x::NLSOLVEJL_SETUP{CS,AD}) where {CS,AD} = CS

macro cache(expr)
  name = expr.args[2].args[1].args[1]
  fields = expr.args[3].args[2:2:end]
  cache_vars = Expr[]
  rand_vars = Expr[]
  jac_vars = Pair{Symbol,Expr}[]
  ratenoise_vars = Expr[]
  for x in fields
    if x.args[2] == :uType || x.args[2] == :rateType ||
       x.args[2] == :kType || x.args[2] == :uNoUnitsType #|| x.args[2] == :possibleRateType
      push!(cache_vars,:(c.$(x.args[1])))
    elseif x.args[2] == :JCType
      push!(cache_vars,:(c.$(x.args[1]).duals...))
    elseif x.args[2] == :GCType
      push!(cache_vars,:(c.$(x.args[1]).duals))
    elseif x.args[2] == :DiffCacheType
      push!(cache_vars,:(c.$(x.args[1]).du))
      push!(cache_vars,:(c.$(x.args[1]).dual_du))
    elseif x.args[2] == :JType || x.args[2] == :WType
      push!(jac_vars,x.args[1] => :(c.$(x.args[1])))
    elseif x.args[2] == :randType
      push!(rand_vars,:(c.$(x.args[1])))
    elseif x.args[2] == :rateNoiseType || x.args[2] == :rateNoiseCollectionType
      # Should be a pair for handling non-diagonal
      push!(ratenoise_vars,:(c.$(x.args[1])))
    end
  end
  quote
    $expr
    $(esc(:full_cache))(c::$name) = tuple($(cache_vars...))
    $(esc(:jac_iter))($(esc(:c))::$name) = tuple($(jac_vars...))
    $(esc(:rand_cache))($(esc(:c))::$name) = tuple($(rand_vars...))
    $(esc(:ratenoise_cache))($(esc(:c))::$name) = tuple($(ratenoise_vars...))
  end
end


text:
abstract type MarkovProcess end

# This type sould define generator(), which returns a transition matrix 𝕋 such that
# 𝕋f = lim_{t→0} E[f(x_t)|x_0=x]/t

"""
 computes the stationary distribution corresponding to the MarkovProcess X
"""
function stationary_distribution(X::MarkovProcess; δ = 0.0, ψ = Zeros(length(X.x)))
    δ >= 0 ||  throw(ArgumentError("δ needs to be positive"))
    if δ > 0
        g = abs.((δ * I - generator(X)') \ (δ * ψ))
    else
        η, g = principal_eigenvalue(generator(X)')
        abs(η) <= 1e-5 || @warn "Principal Eigenvalue does not seem to be zero"
    end
    g ./ sum(g)
end

"""
    Returns the Diffusion Process `x_t` with SDE

        dx_t = μ(x_t) dt + σ(x_t) dZ_t

"""
mutable struct DiffusionProcess <: MarkovProcess
    x::AbstractVector{<:Real}
    μx::AbstractVector{<:Real}
    σx::AbstractVector{<:Real}
    function DiffusionProcess(x::AbstractVector{<:Real}, μx::AbstractVector{<:Real}, σx::AbstractVector{<:Real})
        length(x) == length(μx) == length(σx) || throw(ArgumentError("Vector for grid, drift, and volatility should have the same size"))
        new(x, μx, σx)
    end
end

state_space(X::DiffusionProcess) = X.x

"""
    Returns the discretized version of the infinitesimal generator of the Diffusion Process

        𝕋: f ⭌ lim 1/t * E[f(x_t)|x_0=x]
                 = μx * ∂f + 0.5 * σx^2 * ∂^2f

    defined on the set of functions f such that

        ∂f(x) = 0

    at the border of the state space

    The transpose of this operator corresponds to

        𝕋': g ⭌ v * g - ∂(μx * g) + 0.5 * ∂^2(σx^2 * g)

    defined on the set of functions g such that

        -μx * g(x) + 0.5 * ∂(σx^2 * g) = 0

    at the border of state space
"""
generator(X::DiffusionProcess) = generator(X.x, X.μx, X.σx)

function generator(x::AbstractVector, μx::AbstractVector, σx::AbstractVector)
    n = length(x)
    𝕋 = Tridiagonal(zeros(n-1), zeros(n), zeros(n-1))
    @inbounds for i in 1:n
        Δxp = x[min(i, n-1)+1] - x[min(i, n-1)]
        Δxm = x[max(i-1, 1) + 1] - x[max(i-1, 1)]
        Δx = (Δxm + Δxp) / 2
        # upwinding to ensure off diagonals are posititive
        if μx[i] >= 0
            𝕋[i, min(i + 1, n)] += μx[i] / Δxp
            𝕋[i, i] -= μx[i] / Δxp
        else
            𝕋[i, i] += μx[i] / Δxm
            𝕋[i, max(i - 1, 1)] -= μx[i] / Δxm
        end
        𝕋[i, max(i - 1, 1)] += 0.5 * σx[i]^2 / (Δxm * Δx)
        𝕋[i, i] -= 0.5 * σx[i]^2 * 2 / (Δxm * Δxp)
        𝕋[i, min(i + 1, n)] += 0.5 * σx[i]^2 / (Δxp * Δx)
    end
    # ensure rows sum to zero with machine precision
    c = sum(𝕋, dims = 2)
    for i in 1:n
        𝕋[i, i] -= c[i]
    end
    return 𝕋
end

"""
    Returns the discretized version of the operator ∂

        δ: f ⭌ ∂f

"""
function ∂(X::DiffusionProcess)
    Diagonal(X.μx) \ generator(X.x, X.μx, Zeros(length(X.x)))
end

"""
    Returns the Ornstein Uhlenbeck process defined by the SDE

        dx_t = -κ * (x_t - xbar) * dt + σ * dZ_t

"""
function OrnsteinUhlenbeck(; xbar = 0.0, κ = 0.1, σ = 1.0, p = 1e-10, length = 100,
    xmin = quantile(Normal(xbar, σ / sqrt(2 * κ)), p), xmax = quantile(Normal(xbar, σ / sqrt(2 * κ)), 1 - p))
    # it's important to take low p to have the right tail index of Additive functional
    if xmin > 0
        x = range(xmin^(1/pow), stop = xmax^(1/pow), length = length).^pow
    else
        x = range(xmin, stop = xmax, length = length)
    end
    DiffusionProcess(x, κ .* (xbar .- x), σ * Ones(Base.length(x)))
end

"""
    Returns the Cox Ingersoll Ross process defined by the SDE

        dx_t = -κ * (x - xbar) * dt + σ * sqrt(x) * dZ_t

"""
function CoxIngersollRoss(; xbar = 0.1, κ = 0.1, σ = 1.0, p = 1e-10, length = 100, α = 2 * κ * xbar / σ^2, β = σ^2 / (2 * κ), xmin = quantile(Gamma(α, β), p), xmax = quantile(Gamma(α, β), 1 - p), pow = 2)
    # check 0 is not attainable
    @assert (2 * κ * xbar) / σ^2 > 1
    x = range(xmin^(1/pow), stop = xmax^(1/pow), length = length).^pow
    DiffusionProcess(x, κ .* (xbar .- x), σ .* sqrt.(x))
end


text:
const EPS = 1e-1

function agc_step(input, inputindexoffset, output, outputindexoffset,
        state, stateindexoffset, epsilon, target, n)

    StateLimit = 1 - EPS
    OneMinusEpsOverThree = (1.0 - epsilon) / 3.0
    EpsOverTarget = epsilon / target

    prevState = state[1 + stateindexoffset]

    f = 0.0

    for i in 1:n-1
        oi = outputindexoffset + i
        ii = inputindexoffset + i
        si = stateindexoffset + i

        output[oi] = abs(input[ii] * (1.0 - state[si]))
        f = output[oi] * EpsOverTarget + OneMinusEpsOverThree *
            (prevState + state[si] + state[si + 1])

        f = min(f, StateLimit)
        prevState = state[si]
        state[si] = f
    end

    i = n
    oi = outputindexoffset + i
    ii = inputindexoffset + i
    si = stateindexoffset + i

    output[oi] = abs(input[ii] * (1.0 - state[si]))
    f = output[oi] * EpsOverTarget + OneMinusEpsOverThree *
        (prevState + state[si] + state[si])
    f = min(f, StateLimit)
    state[si] = f

    nothing
end

function agc(input, nChannels, nSamples, nStages, agcParams, state, output)

    for j in 0:nSamples-1
        agc_step(input, j * nChannels, output, j * nChannels, state, 0, agcParams[2], agcParams[1], nChannels)

        for i in 1:nStages-1
            agc_step(output, j * nChannels, output, j * nChannels, state, i * nChannels, agcParams[2 + 2i], agcParams[1 + 2i], nChannels)
        end

    end

    nothing

end

function agc(input::AbstractMatrix, agcParams, state = nothing)

    nChannels, nSamples = size(input)
    nStages = size(agcParams)[2]

    if isnothing(state)
        state = zeros(nChannels, nStages)
    end

    output = zeros(nChannels, nSamples)

    agc(input, nChannels, nSamples, nStages, agcParams, state, output)

    output, state

end


text:
using StatsBase

mutable struct LanguageFind
    lang::AbstractString
    function LanguageFind(text::AbstractString, ngram::Integer = 0)
        length(text) < 1 && throw(DomainError("Text is empty"))
        (ngram < 0 || ngram > 4) && throw(DomainError("Ngram is not valid"))
        if(ngram == 0)
            lang = sort(collect(countmap(collect([i for i in [check_language(text, 2),check_language(text, 3),check_language(text, 4)] if i != ""]))), by=x->x[2], rev=true)[1][1]
        else
            lang = check_language(text, ngram) end
        new(lang)
    end
end

function read_corpus(PATH::String)
    s = ""
    open(PATH, "r") do f
        s=read(f, String) end
    return s
end

function process_text(TEXT::String)
    temp_text = replace(TEXT, r"[\n\r\t]" => " ")
    temp_text = replace(temp_text, r"[^\p{L}]" => " ")
    temp_text = replace(temp_text, r"\s\s+" => " ")
    return strip(lowercase(temp_text))
end

function count_letter_map(TEXT::SubString)
    counts = Dict{Any,Int64}()
    counts = countmap(collect(TEXT))
    return counts
end

function count_n_grams(TEXT::SubString, n::Integer)
    arr = []
    collection = collect(TEXT)
    for i in 1:length(collection)-n+1
        temp_gram = ""
        for j in 1:n
            temp_gram *= collection[i+j-1] end
            push!(arr, temp_gram) end
    counts = countmap(arr)
    return counts
end

function add_count_maps(DICT1::Dict, DICT2::Dict)
    return merge!(+, DICT1, DICT2)
end

function density_calculate(ARR::Array, CORPUS::Array)
    arr_test = []
    arr_corpus = []
    for i in CORPUS
        pos = findfirst(x->x[1]==i[1], ARR)
        if(pos != nothing)
            push!(arr_corpus, i[2])
            push!(arr_test, ARR[pos][2])
        else
            push!(arr_corpus, i[2])
            push!(arr_test, 0) end end
    return arr_test, arr_corpus
end

function density_normalize(ARR::Array)
    arr = []
    temp_sum = 0
    for i in 1:length(ARR)
        temp_sum += ARR[i] end
    for i in 1:length(ARR)
        push!(arr, ARR[i]/temp_sum) end
    return arr
end

function density_calculate_sum(ARR::Array, CORPUS::Array)
    arr_corpus = []
    for i in CORPUS push!(arr_corpus, i[2]) end
    normal_corpus = density_normalize(arr_corpus)
    sum_weight = 0
    for i in 1:length(CORPUS)
        pos = findfirst(x->x[1]==CORPUS[i][1], ARR)
        if(pos != nothing) sum_weight += (ARR[pos][2]*normal_corpus[i]) end end
    return sum_weight
end

function distance_delta(TEXT::String, CORPUS::Array, n::Integer)
    temp_sort = sort(collect(count_n_grams(process_text(TEXT), n)), by=x->x[2], rev=true)
    return density_calculate_sum(temp_sort, CORPUS)
end

function check_language(text::String, ngram::Integer = 3)
    best_score = 0
    lang_code = ""
    ngram_files = readdir(joinpath(dirname(pathof(LanguageFinder)), "Wikipedia", "ngrams", string(ngram)))
    for i in ngram_files
        f = open(joinpath(dirname(pathof(LanguageFinder)), "Wikipedia", "ngrams", string(ngram), string(i)), "r")
        raw_text = readlines(f)
        close(f)
        dictionary = Dict{Any,Int64}()
        for j in raw_text
            push!(dictionary,split(j, ",")[1] => parse(Int, split(j, ",")[2])) end
        sorted_array = sort(collect(dictionary), by=x->x[2], rev=true)
        temp_distance = distance_delta(text, sorted_array, ngram)
        if(temp_distance > best_score)
            lang_code = i[1:2]
            best_score = distance_delta(text, sorted_array, ngram) end end
    return lang_code
end

max_stars_repo_name: nusretipek/LanguageFinder

max_issues_repo_name: nusretipek/LanguageFinder

max_forks_repo_name: nusretipek/LanguageFinder


text:
struct SimulationExpection <: Exception
    message::String
end

get_organism_ids(snapshot) = [get_id(snapshot, o) for o in get_organisms(snapshot)]
get_ids(snapshot, organisms::Vector) = [get_id(snapshot, o) for o in organisms]

function get_fitness(snapshot, abstracted_organism_to_replace::Organism, new_genotype)
    get_fitness(snapshot, abstracted_organism_to_replace.id, new_genotype)
end
function get_fitness(snapshot, key_to_replace::UInt64, new_genotype)
    new_organism = simulate_organism!(snapshot, key_to_replace, new_genotype)
    fitness = get_fitness(snapshot, new_organism)

    return fitness
end
function get_fitness(snapshot, organism)
    num_reproducing_daughters = length([
        d for d in get_daughters(snapshot, organism)
        if get_genotype(snapshot, organism) == get_genotype(snapshot, d) || !isempty(get_daughters(snapshot, d))
    ])
    age = get_age(snapshot, organism)
    fitness = num_reproducing_daughters / age

    return age == 0 ? 0 : fitness
end

function simulate_snapshot!(termination_predicate, snapshot, logger = DoNothingLogger())
    snapshot = deepcopy(snapshot)
    set_logger!(snapshot, logger)
    run_until!(termination_predicate, snapshot)
    return snapshot
end

function simulate_organism!(snapshot, abstracted_organism_to_replace::Organism, new_genotype)
    simulate_organism!(snapshot, abstracted_organism_to_replace.id, new_genotype)
end
function simulate_organism!(snapshot, key_to_replace::UInt64, new_genotype)
    snapshot = deepcopy(snapshot)

    new_organism = replace_organism!(snapshot, key_to_replace, new_genotype)

    # run simulation until death of the new organism and its (non-identical) children
    # (we assume that identical children are able to reproduce themselves)

    run_until!(snapshot) do snapshot
        living_organisms = get_organisms(snapshot)

        return !(new_organism in living_organisms) &&
                all([
                    get_genotype(snapshot, new_organism) == get_genotype(snapshot, d) ||
                    !(d in get_organisms(snapshot))

                    for d in get_daughters(snapshot, new_organism)
                ])
    end

    return new_organism
end

function run_n_timesteps!(snapshot, n_timesteps)
    executed_timesteps = 0

    run_until!(snapshot) do _
        executed_timesteps += 1
        return executed_timesteps == n_timesteps
    end
end

"""
To implement by each model:

get_time(snapshot)
get_organisms(snapshot)

get_logger(snapshot)
set_logger!(snapshot)

get_id(snapshot, organism)

get_genotype(snapshot, organism)
get_parent_genotype(snapshot, organism)

get_position(snapshot, organism)
get_environment(snapshot, organism)

get_time_birth(snapshot, organism)
get_time_death(snapshot, organism)
get_age(snapshot, organism)

get_daughters(snapshot, organism)

get_abstracted_organism(snapshot, organism, parent_id)

replace_organism!(snapshot, key_to_replace, new_genotype)
run_until!(predicate, snapshot)
"""

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
using IterTools
d = readdlm("input.txt", Int)
println(sum(maximum(d, 2) - minimum(d, 2)))
println(sum(mapslices(r -> [Int(b/a) for (a, b) in IterTools.product(r, r) if a != b && isinteger(b/a)], d, 2)))


text:
using MPSKit,TensorKit,Plots

let
    #defining the hamiltonian
    th = nonsym_ising_ham(lambda = 0.3)
    (sxt,syt,szt) = nonsym_spintensors(1//2);

    #initilizing a random mps
    ts = InfiniteMPS([ℂ^2],[ℂ^12]);

    #Finding the groundstate
    (ts,pars,_)=find_groundstate(ts,th,Vumps(maxiter=400));

    len=20;deltat=0.05;totaltime=3.0

    #apply a single spinflip at the middle site
    mpco = MPSComoving(ts,copy.([ts.AC[1];ts.AR[2:len]]),ts)
    @tensor mpco.AC[Int(round(len/2))][-1 -2;-3]:=mpco.AC[Int(round(len/2))][-1,1,-3]*sxt[-2,1]
    normalize!(mpco);

    pars = params(mpco,th)

    szdat = [expectation_value(mpco,szt)]

    for i in 1:(totaltime/deltat)
        (mpco,pars) = changebonds(mpco,th,OptimalExpand()&SvdCut(trschemes=[truncdim(20)])) # grow the bond dimension by 1, and truncate at bond dimension 20
        (mpco,pars) = timestep(mpco,th,deltat,Tdvp(),pars)
        push!(szdat,expectation_value(mpco,szt))
    end

    display(heatmap(real.(reduce((a,b)->[a b],szdat))))

    println("Enter to continue ...")
    readline();
end