text:
mutable struct MapGenerateState{T,U}
    score::Float64
    noise::Float64
    weight::Float64
    subtraces::Vector{U}
    retval::Vector{T}
    num_nonempty::Int
end

function process!(gen_fn::Map{T,U}, args::Tuple, choices::ChoiceMap,
                  key::Int, state::MapGenerateState{T,U}) where {T,U}
    local subtrace::U
    local retval::T
    kernel_args = get_args_for_key(args, key)
    submap = get_submap(choices, key)
    (subtrace, weight) = generate(gen_fn.kernel, kernel_args, submap)
    state.weight += weight
    state.noise += project(subtrace, EmptySelection())
    state.num_nonempty += (isempty(get_choices(subtrace)) ? 0 : 1)
    state.score += get_score(subtrace)
    state.subtraces[key] = subtrace
    retval = get_retval(subtrace)
    state.retval[key] = retval
end

function generate(gen_fn::Map{T,U}, args::Tuple, choices::ChoiceMap) where {T,U}
    len = length(args[1])
    state = MapGenerateState{T,U}(0., 0., 0., Vector{U}(undef,len), Vector{T}(undef,len), 0)
    # TODO check for keys that aren't valid constraints
    for key=1:len
        process!(gen_fn, args, choices, key, state)
    end
    trace = VectorTrace{MapType,T,U}(gen_fn,
        PersistentVector{U}(state.subtraces), PersistentVector{T}(state.retval),
        args, state.score, state.noise, len, state.num_nonempty)
    (trace, state.weight)
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct JobDef_TasksEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function JobDef_TasksEntry(; kwargs...)
        obj = new(meta(JobDef_TasksEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct JobDef_TasksEntry (mapentry)
const __meta_JobDef_TasksEntry = Ref{ProtoMeta}()
function meta(::Type{JobDef_TasksEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_JobDef_TasksEntry)
            __meta_JobDef_TasksEntry[] = target = ProtoMeta(JobDef_TasksEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => Int32, :value => AbstractString]
            meta(target, JobDef_TasksEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_JobDef_TasksEntry[]
    end
end
function Base.getproperty(obj::JobDef_TasksEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct JobDef <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function JobDef(; kwargs...)
        obj = new(meta(JobDef), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct JobDef
const __meta_JobDef = Ref{ProtoMeta}()
function meta(::Type{JobDef})
    ProtoBuf.metalock() do
        if !isassigned(__meta_JobDef)
            __meta_JobDef[] = target = ProtoMeta(JobDef)
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :tasks => Base.Dict{Int32,AbstractString}]
            meta(target, JobDef, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_JobDef[]
    end
end
function Base.getproperty(obj::JobDef, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :tasks
        return (obj.__protobuf_jl_internal_values[name])::Base.Dict{Int32,AbstractString}
    else
        getfield(obj, name)
    end
end

mutable struct ClusterDef <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ClusterDef(; kwargs...)
        obj = new(meta(ClusterDef), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct ClusterDef
const __meta_ClusterDef = Ref{ProtoMeta}()
function meta(::Type{ClusterDef})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ClusterDef)
            __meta_ClusterDef[] = target = ProtoMeta(ClusterDef)
            allflds = Pair{Symbol,Union{Type,String}}[:job => Base.Vector{JobDef}]
            meta(target, ClusterDef, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_ClusterDef[]
    end
end
function Base.getproperty(obj::ClusterDef, name::Symbol)
    if name === :job
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{JobDef}
    else
        getfield(obj, name)
    end
end

export JobDef_TasksEntry, JobDef, ClusterDef
# mapentries: "JobDef_TasksEntry" => ("Int32", "AbstractString")

ext: jl


text:
using NCDatasets, ColorSchemes, RPRMakie
using ImageShow, FileIO

# Taken from https://lazarusa.github.io/BeautifulMakie/GeoPlots/topography/
cmap = dataset = Dataset(joinpath(@__DIR__, "ETOPO1_halfdegree.nc"))
lon = dataset["lon"][:]
lat = dataset["lat"][:]
data = Float32.(dataset["ETOPO1avg"][:, :])

function glow_material(data_normed)
    emission_weight = map(data_normed) do i
        return Float32(i < 0.7 ? 0.0 : i)
    end
    emission_color = map(data_normed) do i
        em = i * 2
        return RGBf(em * 2.0, em * 0.4, em * 0.3)
    end

    return (
        reflection_weight = 1,
        reflection_color = RGBf(0.5, 0.5, 1.0),
        reflection_metalness = 0,
        reflection_ior = 1.4,
        diffuse_weight = 1,
        emission_weight = emission_weight',
        emission_color = emission_color',
    )
end

RPRMakie.activate!(iterations=32, plugin=RPR.Northstar)
fig = Figure(; resolution=(2000, 800))
radiance = 30000
lights = [EnvironmentLight(1.0, load(RPR.assetpath("studio026.exr"))),
            PointLight(Vec3f(0, 100, 100), RGBf(radiance, radiance, radiance))]

ax = LScene(fig[1, 1]; show_axis=false, scenekw=(lights=lights,))

mini, maxi = extrema(data)
data_normed = ((data .- mini) ./ (maxi - mini))

material = glow_material(data_normed)

pltobj = surface!(ax, lon, lat, data_normed .* 20;
                    material=material, colormap=[:black, :white, :brown],
                    colorrange=(0.2, 0.8) .* 20)
# Set the camera to a nice angle
cam = cameracontrols(ax.scene)
cam.eyeposition[] = Vec3f(3, -300, 300)
cam.lookat[] = Vec3f(0)
cam.upvector[] = Vec3f(0, 0, 1)
cam.fov[] = 23

save("earth.png", ax.scene)


text:
# --- Binary managment

# ----------------------------------------------------
# --- Modules
# ----------------------------------------------------
using Random

# ----------------------------------------------------
# --- Source code
# ----------------------------------------------------

"""
---
Create a binary sequence and populate input buffer with bits
The array is of type UInt8 with x00 or x01)
If stated, randSeed controls the seed of the random generator
# --- Syntax
	  genBitsequence!(buffer,nbBits,randSeed=-1);
# --- Input parameters
- buffer	: Buffer to populate [Array{UInt8,nbBits}]
- nbBits	: Number of bits to generate [Int]
- randSeed	: Seed of random process (default -> -1, no seed is used)
# --- Output parameters
- buffer	: Populated buffer [Array{UInt8}]
"""
function genBitSequence!(buffer::Array{UInt8},randSeed=-1)
	# --- Setting the random seed if given
	if randSeed != -1
		# --- Seed as the second parameter
		Random.seed!(randSeed);
	end
	# --- Generating binary sequence
	#buffer .= Int8.(rand([0,1],Int(nbBits)));
	Random.rand!(buffer,[0,1]);
	# --- Switch to "pure" random system
	if randSeed != -1
		RandomDevice();
	end
end

"""
---
Create a binary sequence and return a buffer with nbBits bits
The array is of type UInt8 with x00 or x01)
If stated, randSeed controls the seed of the random generator
# --- Syntax
	  genBitsequence(nbBits,randSeed=-1);
# --- Input parameters
- nbBits	: Number of bits to generate [Int]
- randSeed	: Seed of random process (default -> -1) [Int]
# --- Output parameters
- buffer	: Populated buffer [Array{UInt8}]
"""
function genBitSequence(nbBits,randSeed=-1)
	# --- Setting the random seed if given
	if randSeed != -1
		# --- Seed as the second parameter
		Random.seed!(randSeed);
	end
	# --- Create the input buffer
	buffer	= zeros(UInt8,nbBits);
	# --- Call the bang method
	genBitSequence!(buffer,randSeed);
	# --- Return the buffer
	return buffer;
end

"""
---
Create a byte sequence and populate input buffer with nbytes bytes
The array is of type UInt8 with x00 to xff
If stated, randSeed controls the seed of the random generator
# --- Syntax
	  genByteSequence!(buffer,nbBytes,randSeed=-1);
# --- Input parameters
- buffer	: Buffer to populate [Array{UInt8,nbByte}]
- randSeed	: Seed of random process (default -> -1)
# --- Output parameters
- buffer	: Populated buffer [Array{UInt8}]
# ---
# v 1.0 - Robin Gerzaguet.
"""
function genByteSequence!(buffer::Array{UInt8},randSeed=-1)
	# --- Setting the random seed if given
	if randSeed != -1
		# --- Seed as the second parameter
		Random.seed!(randSeed);
	end
	# --- Generating byte sequence
	Random.rand!(buffer, 0x00:0x01:0xff);
	# --- Switch to "pure" random system
	if randSeed != -1
		RandomDevice();
	end
	# --- Export buffer
	return buffer;
end

"""
---
Create a byte sequence and return a populated  buffer with nbytes bytes
The array
# --- Syntax
	  genByteSequence(nbBytes,randSeed=-1);
# --- Input parameters
- nbBytes	: Number of byte to generate [Int]
- randSeed	: Seed of random process (default -> -1)
# --- Output parameters
- buffer	: Populated buffer [Array{UInt8}]
# ---
# v 1.0 - Robin Gerzaguet.
"""
function genByteSequence(nbBytes,randSeed=-1)
	# --- Create buffer
	buffer = zeros(UInt8,nbBytes);
	# --- Call !
	genByteSequence!(buffer,randSeed);
end


text:
using Interact
using Plots
include("analysis_functions.jl")
include("plot_fcns.jl")

dosfiles = string.("DOSes/", readdir("./DOSes/"))
electrode_surfaces = [dosfile[7:end-7] for dosfile in dosfiles]

@manipulate for Ef_eq=-2.0:0.05:1.0, λ=0.05:0.01:0.5, T=100:20:400, electrode_surface=electrode_surfaces, show_electrolyte_DOS=true
    kT = 8.617e-5*T
    dosfile = string("DOSes/", electrode_surface, "dos.txt")
    dos_f, average_dos, min_E, max_E = get_dos(dosfile; Ef=Ef_eq)
    if occursin("Li1",electrode_surface)
        E_range=1.5
    else
        E_range=4
    end
    E_vals = range(-E_range, E_range, length=500)
    dos_vals=dos_f.(E_vals)
    fd_vals_e = fermi_dirac.(E_vals; kT=kT)
    fd_vals_h = 1 .- fd_vals_e
    dos_vals_e = fd_vals_e .* dos_vals
    dos_vals_h = fd_vals_h .* dos_vals
    if show_electrolyte_DOS
        ox_dos_vals = 0.5.*maximum(dos_vals).*marcus_integrand.(E_vals, λ, 0, true; kT=kT)
        red_dos_vals = 0.5.*maximum(dos_vals).*marcus_integrand.(E_vals, λ, 0, false; kT=kT)
        p1 = plot(E_vals, hcat(dos_vals_e, dos_vals_h, ox_dos_vals, red_dos_vals), label=["electrons" "holes" "reduction" "oxidation"], title="Electrode state occupations", xlabel="E-E_f", ylabel="# (arb.)", color=[:blue :orange :blue :orange], linestyle=[:solid :solid :dash :dash])
    else
        p1 = plot(E_vals, hcat(dos_vals_e, dos_vals_h), label=["electrons" "holes"], title="Electrode state occupations", xlabel="E-E_f", ylabel="# (arb.)")
    end
    p3 = plot_comparison(dos_f, E_range, min_E, max_E, average_dos; kT=kT, λ=λ, plot_title="Rate constants")
    plot(p1, p3, layout=(2,1), size=(500,500))
end

#w = Window()
#body!(w, ui);


text:
"""
Estimate Stochastic Differential Equation

This module uses @debug to report estimation error messages.
"""
module EstimateSDE

export isdiagonalizable,TimeSeriesModel,VectorSde,
        statedim,measdim,logpdf,testprocess,
        autocovz,autocov,statcovz,statcov,randtimes, convert
export Sde
export checkstable, isstable, isvalid
export posterior_samples
export logpdf_reg_refprior
export logpdf2pdf
export describe_model
export powerspectrum
export invscode
export unique_index
export modelvariation
export MlEstimator
export Problem,ProblemProb,ProblemLike, transformation
export PmeanEstimator, State, State, statstate
export SdeProcess, testprocess, SdeProcess, order, kldiv, WhiteNoise
export detstate, lowerorder, powerspectrum_ar
export ArModelStruct, predpar, interpolate,initial_estimate
export SdeRoots,diag_vectorsde
export SdeEngine,SdeRootsEngine,ElimQ,MatFrac,VectorSdeEngine
export VectorSdeProcess,logpdf_ml_σϵ
export logpdf_state, paramstuple
export transformation, transform_ss, sdetype
export AddRoot, InitialModel, ResampleAr
export sortroots, logpdf_ml_white_noise
export Obs, semiposdef_mvrand, combineroots, splitroots
export sde,ml_estimate_init, engine
export parint, AddComplexRoot,RemoveRoots,default_init_tors
export ScalarSde
export estimate,mean,rand
export SquaredExp,defaultengine, datacov, ReferencePrior
export std, EstimateInfo, EstimationStage, label
export systemstatus, tracetable, reporttrace
export issorted_roots, issorted_rootstuple
export in
export safe_kldiv
export maxlike_est, bringin
export BatchDataSettings,batch_data
export BatchMlEstimation,batch_mle_trace_single,batch_mle_trace
export default_random_init_tors
export batchsettings
export batch_posterior_samples
export trace_estimates
export settings_resultsdir
export batch_model_evaluation
export derived_estimates
export copy_icmla_figures

using Distributions
using Optim
using Polynomials
using Zygote
using Parameters
using TransformVariables
using LogDensityProblems
using LogDensityProblems: logdensity_and_gradient
using DynamicHMC
using DocStringExtensions: TYPEDEF, FIELDS, SIGNATURES
using ToeplitzMatrices
using LinearAlgebra
using Random
using Printf
using DataStructures
using Requires
using .Threads
using DSP
using QuickTypes
using Statistics
using StaticArrays
using Intervals
using MCMCDiagnostics
using Plots
using Dates
using InteractiveUtils
using Serialization
using DataFrames
using Pipe
using ProgressMeter
using Pkg.TOML
using REPL.TerminalMenus

basedir = abspath(joinpath(@__DIR__,".."))

"Results directory"
function resultsdir(parts...)
  d = joinpath(homedir(),"results","EstimateSDE",parts...)
  if !ispath(d) mkpath(d) end
  return d
end

function __init__()
  @require MD5="6ac74813-4b46-53a4-afec-0b5dc9d7885c" begin
      @require PyCall="438e738f-606a-5dbb-bf0a-cddfbfd45ab0" begin
        include("stan.jl")
      end
  end
end

include("armodels.jl")
using .ArModels
include("utilities.jl")
include("state.jl")
include("vector-sde.jl")
include("scalar-sde.jl")
include("white-noise.jl")
include("squared-exp.jl")
include("roots.jl")
include("likelihood.jl")
include("problem.jl")
include("random.jl")
include("transformations.jl")
include("estimators-init.jl")
include("estimators.jl")
include("kullback-leibler.jl")
include("report.jl")
include("batch.jl")

end # module


text:
function matrix_construct(N::Int64, q::Complex{Float64}, alphas::Array{Float64, 1}; flag::String="False", bc::String = "Neumann")
    l = ones(N);
    if bc == "Neumann"
        d = [4r^2 for r in 0:N]; # N+1 array
        if flag == "True"  # This is Mathieu's canonical form. Or if first element is non-zero
            l[1] = sqrt(2) * l[1];
        end
    elseif bc == "Dirichlet"
        d = [4r^2 for r in 1:N+1]; # N+1 array
    else
        error("BC not recognized. Try `Neumann` or `Dirichlet`")
    end
    A = spdiagm(0 => d);
    for i in 0:length(alphas[2:end])
        if alphas[i+1] != 0
            as = alphas[i+1]*l[1:N-i];
            A = A + spdiagm((i+1) => q.*as, -(i+1)=> q.*as);
        end
    end
    return A
end

function eigval_iterator(N::Int64, q::Array{Complex{Float64}, 1}, alphas::Array{Float64, 1}; flag::String="False", bc::String="Neumann")
    if q[1] == 0
        q[1] = 0.5*(q[1]+q[2]);
    end
    A = matrix_construct(N, q[1], alphas; flag, bc);
    l, = eigs(A, nev= N-1, which=:SM);
    l = round.(l, digits=2);
    I = sortperm(l, lt = (x,y) -> real(x)==real(y) ? imag(x)<imag(y) : real(x)<real(y));
    a = l[I];
    for qf in q[2:end]
        A = matrix_construct(N, qf, alphas; flag, bc);
        l, = eigs(A, nev=N-1, which=:SM);
        l = round.(l, digits=2);
        I = sortperm(l, lt = (x,y) -> real(x)==real(y) ? imag(x)<imag(y) : real(x)<real(y));
        a = hcat(a, l[I]);
    end
    if bc == "Neumann"
        var = "a"  # eigenvalue a_2n
        m = 0  # sets lowest index in a
    elseif bc == "Dirichlet"
        var = "b"  # eigenvalue b_2n+2
        m = 2  # sets lowest index in b
    end
    vals = Dict(var*string(m) => a[1, :], "q"=> q);  # initialized Dict
    for n in 1:length(l)-1
        vals[var*string(2n+m)] = a[n+1, :];
    end
    return vals
end

max_stars_repo_name: Mikejmnez/HillsEigenvalues

max_issues_repo_name: Mikejmnez/HillsEigenvalues

max_forks_repo_name: Mikejmnez/HillsEigenvalues


text:
using UnPack: @unpack
using MuladdMacro: @muladd
using DiffEqBase: DiffEqBase
using FastBroadcast: @..
using StochasticDiffEq: StochasticDiffEq
using LinearAlgebra: LAPACK, diagm, diag, mul!, diagind
using RecursiveArrayTools: ArrayPartition

using NQCDynamics: get_temperature

export MDEF_BAOAB

struct MDEF_BAOAB <: StochasticDiffEq.StochasticDiffEqAlgorithm end

StochasticDiffEq.alg_compatible(::DiffEqBase.AbstractSDEProblem,::MDEF_BAOAB) = true

struct MDEF_BAOABConstantCache{uType,uEltypeNoUnits} <: StochasticDiffEq.StochasticDiffEqConstantCache
    k::uType
    half::uEltypeNoUnits
end
struct MDEF_BAOABCache{uType,rType,vType,uTypeFlat,uEltypeNoUnits,rateNoiseType,compoundType} <: StochasticDiffEq.StochasticDiffEqMutableCache
    tmp::uType
    utmp::rType
    dutmp::vType
    k::vType
    flatdutmp::uTypeFlat
    tmp1::uTypeFlat
    tmp2::uTypeFlat
    gtmp::compoundType
    noise::rateNoiseType
    half::uEltypeNoUnits
    c1::Matrix{uEltypeNoUnits}
    c2::Matrix{uEltypeNoUnits}
end

function StochasticDiffEq.alg_cache(::MDEF_BAOAB,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
    k = zero(rate_prototype.x[1])
    MDEF_BAOABConstantCache(k, uEltypeNoUnits(1//2))
end

function StochasticDiffEq.alg_cache(::MDEF_BAOAB,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
    tmp = zero(u)
    dutmp = zero(u.x[1])
    utmp = zero(u.x[2])
    k = zero(rate_prototype.x[1])
    flatdutmp = zero(dutmp[:])
    tmp1 = zero(flatdutmp)
    tmp2 = zero(flatdutmp)

    gtmp = zeros(length(utmp), length(utmp))
    noise = zero(rate_prototype.x[1][:])

    half = uEltypeNoUnits(1//2)
    c1 = zeros(length(utmp), length(utmp))
    c2 = zero(c1)

    MDEF_BAOABCache(tmp, utmp, dutmp, k, flatdutmp, tmp1, tmp2, gtmp, noise, half, c1, c2)
end

function StochasticDiffEq.initialize!(integrator, cache::MDEF_BAOABConstantCache)
    @unpack t,uprev,p = integrator
    du1 = uprev.x[1]
    u1 = uprev.x[2]

    cache.k .= integrator.f.f1(du1,u1,p,t)
end

function StochasticDiffEq.initialize!(integrator, cache::MDEF_BAOABCache)
    @unpack t,uprev,p = integrator
    du1 = integrator.uprev.x[1]
    u1 = integrator.uprev.x[2]

    integrator.f.f1(cache.k,du1,u1,p,t)
end

@muladd function StochasticDiffEq.perform_step!(integrator,cache::MDEF_BAOABConstantCache,f=integrator.f)
    @unpack t,dt,sqdt,uprev,p,W = integrator
    @unpack k, half = cache
    du1 = uprev.x[1]
    u1 = uprev.x[2]

    # B
    du2 = du1 + half*dt*k

    # A
    u2 = u1 + half*dt*du2

    # O
    Λ = integrator.g(u2,p,t+dt*half)
    σ = sqrt(get_temperature(p, t+dt*half)) ./ sqrt.(repeat(p.atoms.masses; inner=ndofs(p)))
    γ, c = LAPACK.syev!('V', 'U', Λ) # symmetric eigen
    clamp!(γ, 0, Inf)
    c1 = diagm(exp.(-γ.*dt))
    c2 = diagm(sqrt.(1 .- diag(c1).^2))
    noise = σ.*W.dW[:] / sqdt
    du3 = c*c1*c'*du2[:] + c*c2*c'*noise
    du3 = reshape(du3, size(du2))

    # A
    u = u2 + half*dt*du3

    # B
    k .= f.f1(du3,u,p,t+dt)
    du = du3 + half*dt*k

    integrator.u = ArrayPartition((du, u))
end

@muladd function StochasticDiffEq.perform_step!(integrator,cache::MDEF_BAOABCache,f=integrator.f)
    @unpack t,dt,sqdt,uprev,u,p,W = integrator
    @unpack utmp, dutmp, k, half, gtmp = cache
    du1 = uprev.x[1]
    u1 = uprev.x[2]

    step_B!(dutmp, du1, half*dt, k)

    step_A!(utmp, u1, half*dt, dutmp)

    integrator.g(gtmp,utmp,p,t+dt*half)
    step_O!(cache, integrator)

    step_A!(u.x[2], utmp, half*dt, dutmp)

    integrator.f.f1(k,dutmp,u.x[2],p,t+dt)
    step_B!(u.x[1], dutmp, half*dt, k)
end

step_B!(v2, v1, dt, k) = @.. v2 = v1 + dt*k
step_A!(r2, r1, dt, v)  = @.. r2 = r1 + dt*v

function step_O!(cache, integrator)
    @unpack t, dt, W, p, sqdt = integrator
    @unpack dutmp, flatdutmp, tmp1, tmp2, gtmp, noise, half, c1, c2 = cache

    Λ = gtmp
    σ = sqrt(get_temperature(p, t+dt*half)) ./ sqrt.(repeat(p.atoms.masses; inner=ndofs(p)))

    @.. noise = σ*W.dW[:] / sqdt

    γ, c = LAPACK.syev!('V', 'U', Λ) # symmetric eigen
    clamp!(γ, 0, Inf)
    for (j,i) in enumerate(diagind(c1))
        c1[i] = exp(-γ[j]*dt)
        c2[i] = sqrt(1 - c1[i]^2)
    end

    # equivalent to: flatdutmp .= c*c1*c'dutmp[:] + c*c2*c'noise
    copyto!(flatdutmp, dutmp)
    mul!(tmp1, transpose(c), flatdutmp)
    mul!(tmp2, c1, tmp1)
    mul!(flatdutmp, c, tmp2)

    mul!(tmp1, transpose(c), noise)
    mul!(tmp2, c2, tmp1)
    mul!(tmp1, c, tmp2)

    @.. flatdutmp += tmp1
    copyto!(dutmp, flatdutmp)
end

DynamicsMethods.select_algorithm(::Simulation{<:DynamicsMethods.ClassicalMethods.AbstractMDEF}) = MDEF_BAOAB()


text:
struct Travelers end

# number of agents in Travelers is 2
n_agents(simpleGame::Travelers) = 2

# set of actions of specific agent in Travelers
ordered_actions(simpleGame::Travelers, i::Int) = 2:100
# set of actions of all agents in Travelers: [2:100] x [2:100]
ordered_joint_actions(simpleGame::Travelers) = vec(collect(Iterators.product([ordered_actions(simpleGame, i) for i in 1:n_agents(simpleGame)]...)))

# number of element of set of actions of all agents
n_joint_actions(simpleGame::Travelers) = length(ordered_joint_actions(simpleGame))
# number of element of set of actions of specific agent
n_actions(simpleGame::Travelers, i::Int) = length(ordered_actions(simpleGame, i))

# function to calculate how much money Agent "i" is rewarded
function reward(simpleGame::Travelers, i::Int, a)
    if i == 1
        noti = 2
    else
        noti = 1
    end
    if a[i] == a[noti]
        r = a[i]
    elseif a[i] < a[noti]
        r = a[i] + 2
    else
        r = a[noti] - 1
    end
    return r
end

# function to return how much money all agents are rewarded
# [a + 2, a - 2]
function joint_reward(simpleGame::Travelers, a)
    return [reward(simpleGame, i, a) for i in 1:n_agents(simpleGame)]
end

# recursive function
function SimpleGame(simpleGame::Travelers)
    return SimpleGame(
        0.9, # discount factor
        vec(collect(1:n_agents(simpleGame))), # vector of all agents [1, 2]
        [ordered_actions(simpleGame, i) for i in 1:n_agents(simpleGame)], # [2:100 2:100] # joint action space
        (a) -> joint_reward(simpleGame, a) # an arrow function
    )
end


text:
using JSON

# --- Utility ---

flatten(x::Array{<:Array, 1}) = Iterators.flatten(x)|> collect|> flatten
flatten(x::Array{<:Number, 1}) = x

shape(x::Array{<:Array, 1}) = vcat(shape(first(x)), [length(x)])
shape(x::Array{<:Number, 1}) = [length(x)]

mdim(x::Array{<:Array, 1}) = reshape(flatten(x), shape(x)...)
mdim(x::Array{<:Number, 1}) = x

transform(x::Array, ::Type{Array{T, N}}) where T <: Number where N = mdim(x)
transform(x::Array, ::Type{Array{T, 1}}) where T <: Array = x
transform(x::Number, ::Type{T}) where T <: Number = x

function convert_type(::Type{Array{T, N}}) where T <: Number where N
    reduce((t,_) -> Array{t, 1}, 1:N; init=T)
end
convert_type(t::Type{Array{T, 1}}) where T <: Array = t
convert_type(t::Type{T}) where T <: Number = t

# -- IO ---

"""Saves any JSON serializable `object` into JSON file in `filepath`."""
function save_json(object::Any, filepath::AbstractString)
    open(filepath, "w") do io
        JSON.print(io, object)
    end
end

"""Loads values from JSON file in `filepath` to DataType supplied by `dtype`."""
function load_json(dtype::DataType, filepath::AbstractString)
    objects = JSON.parsefile(filepath)
    fields = (
        objects[string(s)] |>
        v -> convert(convert_type(t), v) |>
        v -> transform(v, t)
        for (s, t) in zip(fieldnames(dtype), fieldtypes(dtype)))
    dtype(fields...)
end

max_stars_repo_name: nnhjy/JuMP-project-template

max_issues_repo_name: nnhjy/JuMP-project-template

max_forks_repo_name: nnhjy/JuMP-project-template


text:
include("Graph.jl")
include("SimMetric.jl")

function topK(G, S)
    dl = Dict{Int32, Int32}()
    foreach(x -> dl[x] = haskey(dl, x) ? dl[x]+1 : 1, G.L)
    maxk = maximum(values(dl))-1
    tot = zeros(maxk)
    totn = zeros(maxk)
    for i = 1 : G.n
        tp = Vector{Tuple{Int32, Float64}}()
        foreach(j -> (i != j) ? push!(tp, (j, S[i, j])) : nothing, 1 : G.n)
        sort!(tp, by=x->-x[2])
        s = 0.0
        for j = 1 : maxk
            s += (G.L[i] == G.L[tp[j][1]]) ? 1 : 0
            totn[j] += 1
            tot[j] += (s / j)
        end
    end
    ret = zeros(maxk)
    foreach(i -> ret[i] = tot[i] / totn[i], 1 : maxk)
    return ret
end

function topKofForestSimApprox(G)
    dl = Dict{Int32, Int32}()
    foreach(x -> dl[x] = haskey(dl, x) ? dl[x]+1 : 1, G.L)
    maxk = maximum(values(dl))-1
    tot = zeros(maxk)
    totn = zeros(maxk)
    S = approxForest(G)
    A = Vector{Tuple{Int32, Float64}}()
    foreach(i -> push!(A, (i, S[i])), 1 : G.n)
    sort!(A, by=x->x[2])
    pos = zeros(Int32, G.n)
    foreach(i -> pos[A[i][1]] = i, 1 : G.n)

    sim(x, y) = min(S[x], S[y]) / max(S[x], S[y])

    predict(x, k) = begin
        l = pos[x]-1
        r = pos[x]+1
        ret = Vector{Int32}()
        for rep = 1 : k
            if l == 0
                push!(ret, A[r][1])
                r += 1
            elseif r > G.n
                push!(ret, A[l][1])
                l -= 1
            elseif sim(x, A[l][1]) > sim(x, A[r][1])
                push!(ret, A[l][1])
                l -= 1
            else
                push!(ret, A[r][1])
                r += 1
            end
        end
        return ret
    end

    for i = 1 : G.n
        cntk = maxk
        tp = predict(i, cntk)
        s = 0.0
        for j = 1 : cntk
            s += (G.L[i] == G.L[tp[j]]) ? 1 : 0
            totn[j] += 1
            tot[j] += (s / j)
        end
    end

    ret = zeros(maxk)
    foreach(i -> ret[i] = tot[i] / totn[i], 1 : maxk)
    return ret
end


text:
module LVServer

using ZMQ, JSON3, ImageCore, Colors, PkgVersion

# version of the communication protocoll: First, and, presumably, the last as well
const PROTOC_V = 0x01

# update these constants on update of any of the packages!
const LabVIEW0_LOCAL_CURRENT = v"0.1.1"
const LabVIEW0_LOCAL_MIN = v"0.1.0"

include("./Julia/ZMQ-server.jl")

# see setglobals() and get_script_path() using these globals
scriptexists = false
scriptOK = false
scriptexcep = nothing
# see setreorderbytes()
reorderbytes = :undefined

export server_0mq4lv, get_script_path, setglobals, setreorderbytes, get_LVlib_path # functions
export scriptexists, scriptOK, scriptexcep # global variables

end


text:
#= DSP is removed from test dependencies in this version.
import MechanicalSketch
import MechanicalSketch: empty_figure, PALETTE, O, HE, WI, EM, finish, ∙, Point
import MechanicalSketch: @import_expand, set_scale_sketch, scale_pt_to_unit, settext
import MechanicalSketch: noise, normalize_datarange, place_image, @layer
import MechanicalSketch: poly, dimension_aligned, sethue, arrow, circle, prettypoly
import MechanicalSketch: Length, noise_between_wavelengths, Greys_9, BinLegend
import DSP.Periodograms: spectrogram, Spectrogram
import DSP.Util:         nextfastfft
import DSP.Windows:      tukey
import MechanicalSketch: MechanicalUnits.dimension, NoDims, upreferred, ComplexQuantity, Quantity
import Statistics:       mean

let
if !@isdefined m²
    @import_expand(~m, s)
end
empty_figure(filename = joinpath(@__DIR__, "test_30.png"));

global curpoint = O + (-WI / 2 + EM, -HE / 2 + 3.5EM)
v_b = (0.05m/s, 0.5m/s)
settext("
From test 21 to 27, streamline velocity bounds are
    vb = $(collect(v_b))
", curpoint, markup=true)

curpoint += (0, 4EM)
Δt = 2.0s
wave_per_streamline = 2
λ_b = Δt ∙ v_b  / wave_per_streamline
λ_calibration = collect(range(λ_b[1], λ_b[2], length = 4))

settext("
Tracing a particle over Δt = $Δt forward ($(Δt / 2)) and back ($(Δt / 2)) we want to sample $wave_per_streamline harmonic noise waves with length λ.
    λb = Δt ∙ vb / $wave_per_streamline = $λ_b
Between bounds, we want a linearly increasing noise function amplitude.
", curpoint, markup = true)

curpoint += (0, 3.5EM)
n_pixels = round(Int, 0.9 * WI)
physwidth = 3.0m
set_scale_sketch(physwidth, n_pixels)
length_one_pixel = physwidth / n_pixels
include("test_functions_29.jl")
include("test_functions_30.jl")
no = normalize_datarange([noise_between_wavelengths(λ_b..., x) for x in (1:n_pixels) * length_one_pixel])
greylegend = BinLegend(;maxlegend = 1.0, noofbins = 128, colorscheme = reverse(Greys_9))
@layer begin
    place_image(curpoint, greylegend.(no); centered = false)
    place_image(curpoint + (0,1), greylegend.(no); centered = false)
    place_image(curpoint + (0,2), greylegend.(no); centered = false)
    sethue(PALETTE[3])
    dimension_aligned(curpoint, curpoint + (λ_calibration[1], 0.0m), offset = -EM,
        fromextension = (0, EM), toextension = (0, 0.5EM))
    dimension_aligned(curpoint, curpoint + (λ_calibration[4], 0.0m), offset = -3EM,
        fromextension = (EM, EM), toextension = (EM, 5EM))
end

curpoint += (0, 3EM)
settext("
The shape looks clearer shown as a curve:", curpoint)

curpoint += (0, EM)
points = map(1:n_pixels, no) do xul, yul
    curpoint + (xul, -yul * EM)
end
poly(points, :stroke)

curpoint += (0, 3EM)
settext("
Spectrum from sampling the curve over $(100 * maximum(λ_b)):
", curpoint)

curpoint += (0, 2EM)
longnoise = normalize_datarange([noise_between_wavelengths(λ_b..., x) for x in (1:100n_pixels)*length_one_pixel])
spectrum = sampledspectrum_29(longnoise, length_one_pixel, maximumwavelength = maximum(λ_b))
draw_spectrum_29(curpoint, spectrum)

curpoint += (0, 5EM)
(λ_1, a_1), (λ_2, a_2) = two_largest_maxima_amplitude_wavelength(spectrum)
a_calibration = map(λ-> amplitude_interpolated(spectrum, λ), λ_calibration)
a_calibration_normalized = round.(a_calibration / maximum(a_calibration), digits = 2)
settext("
    Local amplitude maxima are (λ, amplitude) = [$((λ_1, a_1)), $((λ_2, a_2))]
    Normalized amplitudes at λcalibration = $λ_calibration :
                                                                      $a_calibration_normalized, which is acceptably linear.
    ",
    curpoint, markup=true)

curpoint += (0, 2EM)
binwidth = 1 / (n_pixels - 1)
binno(x) = 1 + round(Int, x / (binwidth ), RoundDown)
binends = range(binwidth, 1.0, length = n_pixels)
bincounts = [0 for bin in binends]
for x in longnoise
    bincounts[binno(x)] += 1
end
relative_frequency = normalize_datarange(bincounts / sum(bincounts))
histogrampoints = Point.(collect(0.9 * WI * binends), -EM * relative_frequency )
poly(curpoint .+ histogrampoints, :stroke)

@layer begin
    sethue(PALETTE[3])
    arrow(curpoint, curpoint + (0.0,  -EM))
    settext("Relative frequency", curpoint + (0.0,  -EM ), markup=true)
    arrow(curpoint, curpoint + (0.95 * WI, 0.0))
    settext("Noise values, max 1.0", curpoint + ( 0.8 * scale_pt_to_unit(m) * WI, 0.0m), markup=true)
end

finish()
end
=#


text:
@testset "fftshift alternatives" begin

    @testset "Test fftshift_view and ifftshift_view" begin
        Random.seed!(42)
        x = randn((2,1,4,1,6,7,4,7))
        dims = (4,6,7)
        @test fftshift(x,dims) == FourierTools.fftshift_view(x, dims)
        @test ifftshift(x,dims) == FourierTools.ifftshift_view(x, dims)
        @test x === FourierTools.optional_collect(ifftshift_view(fftshift_view(x)))
        @test x === FourierTools.optional_collect(fftshift_view(ifftshift_view(x)))
        @test x === FourierTools.optional_collect(ifftshift_view(fftshift_view(x, dims), dims))
        @test x === FourierTools.optional_collect(fftshift_view(ifftshift_view(x, dims), dims))

        x = randn((13, 13, 14))
        @test fftshift(x) == FourierTools.fftshift_view(x)
        @test ifftshift(x) == FourierTools.ifftshift_view(x)
        @test fftshift(x, (2,3)) == FourierTools.fftshift_view(x, (2,3))
        @test ifftshift(x, (2,3) ) == FourierTools.ifftshift_view(x, (2,3))

    end
end


text:
const USE_GPU = haskey(ENV, "USE_GPU") ? parse(Bool, ENV["USE_GPU"]) : false
const do_viz  = haskey(ENV, "DO_VIZ")  ? parse(Bool, ENV["DO_VIZ"])  : true
const do_save = haskey(ENV, "DO_SAVE") ? parse(Bool, ENV["DO_SAVE"]) : false
const nx = haskey(ENV, "NX") ? parse(Int, ENV["NX"]) : 512 - 1
const ny = haskey(ENV, "NY") ? parse(Int, ENV["NY"]) : 512 - 1
###
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 2)
else
    @init_parallel_stencil(Threads, Float64, 2)
end
using Plots, Printf, Statistics, LinearAlgebra

@parallel function smooth!(A2::Data.Array, A::Data.Array, fact::Data.Number)
    @inn(A2) = @inn(A) + 1.0/4.1/fact*(@d2_xi(A) + @d2_yi(A))
    return
end

@parallel function compute_maxloc!(Musτ2::Data.Array, Musτ::Data.Array)
    @inn(Musτ2) = @maxloc(Musτ)
    return
end

@parallel function compute_timesteps!(dτVx::Data.Array, dτVy::Data.Array, dτPt::Data.Array, Musτ::Data.Array, Vsc::Data.Number, Ptsc::Data.Number, min_dxy2::Data.Number, β_n::Data.Number, max_nxy::Int)
    @all(dτVx) = 1.0/4.1/Vsc*min_dxy2/@av_xi(Musτ)/(1.0+β_n)
    @all(dτVy) = 1.0/4.1/Vsc*min_dxy2/@av_yi(Musτ)/(1.0+β_n)
    @all(dτPt) = 4.1*(1.0+β_n)/Ptsc/max_nxy*@all(Musτ)
    return
end

@parallel function compute_P!(∇V::Data.Array, Pt::Data.Array, Vx::Data.Array, Vy::Data.Array, dτPt::Data.Array, dx::Data.Number, dy::Data.Number)
    @all(∇V)  = @d_xa(Vx)/dx + @d_ya(Vy)/dy
    @all(Pt)  = @all(Pt) - @all(dτPt)*@all(∇V)
    return
end

@parallel function compute_τ!(∇V::Data.Array, τxx::Data.Array, τyy::Data.Array, τxy::Data.Array, Vx::Data.Array, Vy::Data.Array, Mus::Data.Array, β_n::Data.Number, dx::Data.Number, dy::Data.Number)
    @all(τxx) = 2.0*@all(Mus)*(@d_xa(Vx)/dx - 1.0/3.0*@all(∇V) + β_n*@all(∇V))
    @all(τyy) = 2.0*@all(Mus)*(@d_ya(Vy)/dy - 1.0/3.0*@all(∇V) + β_n*@all(∇V))
    @all(τxy) = 2.0*@av(Mus)*(0.5*(@d_yi(Vx)/dy + @d_xi(Vy)/dx))
    return
end

@parallel function compute_dV!(Rx::Data.Array, Ry::Data.Array, dVxdτ::Data.Array, dVydτ::Data.Array, Pt::Data.Array, τxx::Data.Array, τyy::Data.Array, τxy::Data.Array, dampX::Data.Number, dampY::Data.Number, dx::Data.Number, dy::Data.Number)
    @all(Rx)    = @d_xi(τxx)/dx + @d_ya(τxy)/dy - @d_xi(Pt)/dx
    @all(Ry)    = @d_yi(τyy)/dy + @d_xa(τxy)/dx - @d_yi(Pt)/dy
    @all(dVxdτ) = dampX*@all(dVxdτ) + @all(Rx)
    @all(dVydτ) = dampY*@all(dVydτ) + @all(Ry)
    return
end

@parallel function compute_V!(Vx::Data.Array, Vy::Data.Array, dVxdτ::Data.Array, dVydτ::Data.Array, dτVx::Data.Array, dτVy::Data.Array)
    @inn(Vx) = @inn(Vx) + @all(dτVx)*@all(dVxdτ)
    @inn(Vy) = @inn(Vy) + @all(dτVy)*@all(dVydτ)
    return
end

@parallel_indices (iy) function bc_x!(A::Data.Array)
    A[1  , iy] = A[2    , iy]
    A[end, iy] = A[end-1, iy]
    return
end

@parallel_indices (ix) function bc_y!(A::Data.Array)
    A[ix, 1  ] = A[ix, 2    ]
    A[ix, end] = A[ix, end-1]
    return
end

@views function Stokes2D()
    # Physics
    lx, ly    = 10.0, 10.0  # domain extends
    μs0       = 1.0         # matrix viscosity
    μsi       = 1e-3        # inclusion viscosity
    εbg       = 1.0         # background strain-rate
    # Numerics
    iterMax   = 1e5         # maximum number of pseudo-transient iterations
    nout      = 500         # error checking frequency
    Vdmp      = 4.0         # damping paramter for the momentum equations
    β_n       = 2.0         # numerical compressibility
    Vsc       = 1.0         # relaxation paramter for the momentum equations pseudo-timesteps limiters
    Ptsc      = 3.0         # relaxation paramter for the pressure equation pseudo-timestep limiter
    ε         = 1e-8        # nonlinear absolute tolerence
    # nx, ny    = 1*128-1, 1*128-1    # numerical grid resolution; should be a mulitple of 32-1 for optimal GPU perf
    # Derived numerics
    dx, dy    = lx/nx, ly/ny # cell sizes
    min_dxy2  = min(dx,dy)^2
    max_nxy   = max(nx,ny)
    dampX     = 1.0-Vdmp/nx # damping term for the x-momentum equation
    dampY     = 1.0-Vdmp/ny # damping term for the y-momentum equation
    # Array allocations
    Pt        = @zeros(nx  ,ny  )
    dτPt      = @zeros(nx  ,ny  )
    ∇V        = @zeros(nx  ,ny  )
    τxx       = @zeros(nx  ,ny  )
    τyy       = @zeros(nx  ,ny  )
    τxy       = @zeros(nx-1,ny-1)
    Rx        = @zeros(nx-1,ny-2)
    Ry        = @zeros(nx-2,ny-1)
    dVxdτ     = @zeros(nx-1,ny-2)
    dVydτ     = @zeros(nx-2,ny-1)
    dτVx      = @zeros(nx-1,ny-2)
    dτVy      = @zeros(nx-2,ny-1)
    Mus2      = @zeros(nx  ,ny  )
    Musτ      = @zeros(nx  ,ny  )
    # Initial conditions
    Rad2      =  zeros(nx  ,ny  )
    Vx        =  zeros(nx+1,ny  )
    Vy        =  zeros(nx  ,ny+1)
    Rad2     .= [((ix-1)*dx +0.5*dx -0.5*lx)^2 + ((iy-1)*dy +0.5*dy -0.5*ly)^2 for ix=1:size(Rad2,1), iy=1:size(Rad2,2)]
    Vx        = Data.Array( -εbg.*[((ix-1)*dx -0.5*lx) for ix=1:size(Vx,1), iy=1:size(Vx,2)] )
    Vy        = Data.Array(  εbg.*[((iy-1)*dy -0.5*ly) for ix=1:size(Vy,1), iy=1:size(Vy,2)] )
    Mus       = μs0*ones(nx,ny)
    Mus[Rad2.<1.0] .= μsi
    Mus       = Data.Array( Mus )
    Mus2     .= Mus
    for ism=1:10
        @parallel smooth!(Mus2, Mus, 1.0)
        Mus, Mus2 = Mus2, Mus
    end
    Musτ     .= Mus
    @parallel compute_maxloc!(Musτ, Mus)
    @parallel (1:size(Musτ,2)) bc_x!(Musτ)
    @parallel (1:size(Musτ,1)) bc_y!(Musτ)
    # Time loop
    @parallel compute_timesteps!(dτVx, dτVy, dτPt, Musτ, Vsc, Ptsc, min_dxy2, β_n, max_nxy)
    err=2*ε; iter=0; err_evo1=[]; err_evo2=[]
    while err > ε && iter <= iterMax
        if (iter==11)  global wtime0 = Base.time()  end
        @parallel compute_P!(∇V, Pt, Vx, Vy, dτPt, dx, dy)
        @parallel compute_τ!(∇V, τxx, τyy, τxy, Vx, Vy, Mus, β_n, dx, dy)
        @parallel compute_dV!(Rx, Ry, dVxdτ, dVydτ, Pt, τxx, τyy, τxy, dampX, dampY, dx, dy)
        @parallel compute_V!(Vx, Vy, dVxdτ, dVydτ, dτVx, dτVy)
        @parallel (1:size(Vx,1)) bc_y!(Vx)
        @parallel (1:size(Vy,2)) bc_x!(Vy)
        iter += 1
        if iter % nout == 0
            Vmin, Vmax = minimum(Vx), maximum(Vx)
            Pmin, Pmax = minimum(Pt), maximum(Pt)
            norm_Rx    = norm(Rx)/(Pmax-Pmin)*lx/sqrt(length(Rx))
            norm_Ry    = norm(Ry)/(Pmax-Pmin)*lx/sqrt(length(Ry))
            norm_∇V    = norm(∇V)/(Vmax-Vmin)*lx/sqrt(length(∇V))
            # norm_Rx = norm(Rx)/length(Rx); norm_Ry = norm(Ry)/length(Ry); norm_∇V = norm(∇V)/length(∇V)
            err = maximum([norm_Rx, norm_Ry, norm_∇V])
            push!(err_evo1, maximum([norm_Rx, norm_Ry, norm_∇V])); push!(err_evo2,iter)
            @printf("Total steps = %d, err = %1.3e [norm_Rx=%1.3e, norm_Ry=%1.3e, norm_∇V=%1.3e] \n", iter, err, norm_Rx, norm_Ry, norm_∇V)
        end
    end
    # Performance
    wtime    = Base.time() - wtime0
    A_eff    = (3*2)/1e9*nx*ny*sizeof(Data.Number)  # Effective main memory access per iteration [GB] (Lower bound of required memory access: Te has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
    wtime_it = wtime/(iter-10)                      # Execution time per iteration [s]
    T_eff    = A_eff/wtime_it                       # Effective memory throughput [GB/s]
    @printf("Total steps = %d, err = %1.3e, time = %1.3e sec (@ T_eff = %1.2f GB/s) \n", iter, err, wtime, round(T_eff, sigdigits=2))
    # Visualisation
    if do_viz
        X, Y, Yv  = dx/2:dx:lx-dx/2, dy/2:dy:ly-dy/2, 0:dy:ly
        p1 = heatmap(X,  Y, Array(Musτ)', aspect_ratio=1, xlims=(X[1],X[end]), ylims=(Y[1],Y[end]), c=:viridis, title="Pressure")
        p2 = heatmap(X, Yv, Array(Vy)', aspect_ratio=1, xlims=(X[1],X[end]), ylims=(Yv[1],Yv[end]), c=:viridis, title="Vy")
        p4 = heatmap(X[2:end-1], Yv[2:end-1], log10.(abs.(Array(Ry)')), aspect_ratio=1, xlims=(X[2],X[end-1]), ylims=(Yv[2],Yv[end-1]), c=:viridis, title="log10(Ry)")
        p5 = plot(err_evo2,err_evo1, legend=false, xlabel="# iterations", ylabel="log10(error)", linewidth=2, markershape=:circle, markersize=3, labels="max(error)", yaxis=:log10)
        display(plot(p1, p2, p4, p5))
    end
    if do_save
        !ispath("../output") && mkdir("../output")
        open("../output/out_Stokes2D0.txt","a") do io
            println(io, "$(nx) $(ny) $(iter)")
        end
    end
    return
end

Stokes2D()


text:
# Autogenerated wrapper script for libsvm_jll for aarch64-linux-gnu
export libsvm, svm_predict, svm_scale, svm_train

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("libsvm")
JLLWrappers.@declare_library_product(libsvm, "libsvm.so")
JLLWrappers.@declare_executable_product(svm_predict)
JLLWrappers.@declare_executable_product(svm_scale)
JLLWrappers.@declare_executable_product(svm_train)
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libsvm,
        "lib/libsvm.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        svm_predict,
        "bin/svm-predict",
    )

    JLLWrappers.@init_executable_product(
        svm_scale,
        "bin/svm-scale",
    )

    JLLWrappers.@init_executable_product(
        svm_train,
        "bin/svm-train",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "712.minimum-ascii-delete-sum-for-two-strings.jl" begin
    @test minimum_delete_sum("sea", "eat") == 231
    @test minimum_delete_sum("delete", "leet") == 403
end


text:
using Test

include("atbash-cipher.jl")

@testset "encoding from English to atbash" begin
    @test encode("yes") == "bvh"
    @test encode("no") == "ml"
    @test encode("OMG") == "lnt"
    @test encode("O M G") == "lnt"
    @test encode("mindblowingly") == "nrmwy oldrm tob"
    @test encode("Testing,1 2 3, testing.") == "gvhgr mt123 gvhgr mt"
    @test encode("Truth is fiction.") == "gifgs rhurx grlm"
    @test encode("The quick brown fox jumps over the lazy dog.") == "gsvjf rxpyi ldmul cqfnk hlevi gsvoz abwlt"
end

@testset "decoding from atbash to English" begin
    @test decode("vcvix rhn") == "exercism"
    @test decode("zmlyh gzxov rhlug vmzhg vkkrm thglm v") == "anobstacleisoftenasteppingstone"
    @test decode("gvhgr mt123 gvhgr mt") == "testing123testing"
    @test decode("gsvjf rxpyi ldmul cqfnk hlevi gsvoz abwlt") == "thequickbrownfoxjumpsoverthelazydog"
    @test decode("vc vix    r hn") == "exercism"
    @test decode("zmlyhgzxovrhlugvmzhgvkkrmthglmv") == "anobstacleisoftenasteppingstone"
end


text:
"""
Subtype of `LDScoreRegression`.
"""

mutable struct Hsq <: LDScoreRegression
    y; x; w; N; M;
    n_blocks; intercept; slow; step1_ii; old_weights;
    __null_intercept__

    constrain_intercept; n_annot; intercept_se; twostep_filtered;

    function Hsq(
        y, x, w, N, M; n_blocks = 200, slow = false, old_weights = false,
        intercept = nothing, step1_ii = nothing,
    )
        hsq = new(
            y, x, w, N, M, n_blocks, intercept, slow, step1_ii, old_weights, 1,
        )
        return ld_score_regression(
            hsq,
            y, x, w, N, M, n_blocks, intercept, slow, step1_ii, old_weights,
        )
    end
end

function hsq_weights(ld, w_ld, N, M, hsq; intercept=nothing)
    if intercept == nothing intercept = 1.0 end
    M = M[1][1]
    hsq = max(maximum(hsq), 0.0)
    hsq = min(hsq, 1.0)
    ld = fmax(ld, 1.0)
    w_ld = fmax(w_ld, 1.0)
    c = hsq .* N ./ M

    het_w = 1.0 ./ (2 .* (c .* ld .+ intercept).^2)
    oc_w = 1.0 ./ w_ld
    w = het_w .* oc_w
    return w
end

function update_weights(reg::Hsq, ld, w_ld, N, M, hsq, intercept)
    if intercept == nothing intercept = reg.__null_intercept__ end

    return hsq_weights(ld, w_ld, N, M, hsq; intercept=intercept)
end

function overlap_output(
    hsq::Hsq,
    category_names, overlap_matrix, M_annot, M_tot, print_coefficients,
)
    overlap_matrix_prop = zeros([hsq.n_annot, hsq.n_annot])
    for i in range(hsq.n_annot)
        overlap_matrix_prop[i, :] = overlap_matrix[i, :] / M_annot
    end

    # TODO missing code here

    #=
    df = DataFrames.DataFrame({
        "Category": category_names,
        "Prop._SNPs": one_d_convert(prop_M_overlap),
        "Prop._h2": one_d_convert(prop_hsq_overlap),
        "Prop._h2_std_error": one_d_convert(prop_hsq_overlap_se),
        "Enrichment": one_d_convert(enrichment),
        "Enrichment_std_error": one_d_convert(enrichment_se),
        "Enrichment_p":diff_p,
        "Coefficient": one_d_convert(hsq.coef),
        "Coefficient_std_error": hsq.coef_se,
        "Coefficient_z-score": one_d_convert(hsq.coef) / one_d_convert(hsq.coef_se)
    })
    if print_coefficients
        return df[
            !,
            [
                "Category", "Prop._SNPs", "Prop._h2", "Prop._h2_std_error",
                "Enrichment","Enrichment_std_error", "Enrichment_p",
                "Coefficient", "Coefficient_std_error","Coefficient_z-score",
            ],
        ]
    else
        return df[
            !,
            [
                "Category", "Prop._SNPs", "Prop._h2", "Prop._h2_std_error",
                "Enrichment","Enrichment_std_error", "Enrichment_p",
            ],
        ]
    end
    =#
end


text:
include("OpCodes.jl")
using .OpCodes
using OffsetArrays

ic = IntComputer(readline("data/Dec17_data.txt"))
asc = getoutput(compute!(ic))
field = string(Char.(asc)...)

spfield = split(field, '\n')[1:end-2]
bf = OffsetArray([spfield[i][j] != '.' for i in 1:length(spfield),
    j in 1:length(spfield[1])], (-1, -1))

isintersection(bf, i, j) = bf[i,j] && bf[i-1,j] && bf[i+1,j] && bf[i, j-1] && bf[i, j+1]

ints = [x for x in Iterators.product(1:45, 1:47) if isintersection(bf, x...)]
sum(prod,ints) # answer to problem 1

# problem 2
corner = Tuple(findfirst([spfield[i][j] == '^' for i in 1:length(spfield), j in 1:length(spfield[1])])) .- (1,1)

function hasscaffold(bf, (i, j))
    (i < 0 || i > size(bf, 1)-1) && return false
    (j < 0 || j > size(bf, 2)-1) && return false
    bf[i,j]
end

function followpath(bf, corner)
    pos = corner
    right =  [0  1;
             -1  0]
    left = -right
    dir = [-1, 0]
    path = []
    stepcounter = 0
    while true
        @show pos
        if hasscaffold(bf, pos .+ dir)
            pos = pos .+ dir
            stepcounter += 1
        else
            push!(path, stepcounter)
            stepcounter = 0
            if hasscaffold(bf, pos .+ left * dir)
                push!(path, 'L')
                dir = left * dir
            elseif hasscaffold(bf, pos .+ right * dir)
                push!(path, 'R')
                dir = right * dir
            else
                return path
            end
        end
    end
end

findrepeats(v) =

p = followpath(bf, corner)[2:end]

str = "A,A,B,C,B,C,B,C\n"
[Int(Char(x[1])) for x in str[1:2:end]]


text:
##MDP Types
#Convenience abstract types to remove dependence on external packages for hte purposes of testing
function test_behavior_model_creation()
	println("\t\tTesting IDMMOBILBehavior creation")
	ps = IDMMOBILBehavior[IDMMOBILBehavior(x[1],x[2],x[3],idx) for (idx,x) in enumerate(Iterators.product(["cautious","normal","aggressive"],[27.,31.,35.],[4.]))][:]
end

function test_behavior_model_equality()
	ps = IDMMOBILBehavior[IDMMOBILBehavior(x[1],x[2],x[3],idx) for (idx,x) in enumerate(Iterators.product(["cautious","normal","aggressive"],[27.,31.,35.],[4.]))][:]
	test_equality("IDMMOBILBehavior",ps)
end

function test_behavior_model_hashing()
	ps = IDMMOBILBehavior[IDMMOBILBehavior(x[1],x[2],x[3],idx) for (idx,x) in enumerate(Iterators.product(["cautious","normal","aggressive"],[27.,31.,35.],[4.]))][:]
	test_hashing("IDMMOBILBehavior",ps)
end

function test_carstate_creation()
	println("\t\tTesting CarState creation")
	#nothing to do here--no constructors to test
end

function test_carstate_equality()
	#println("\t\tTesting CarState equality")
	bs = IDMMOBILBehavior[IDMMOBILBehavior(x[1],x[2],x[3],idx) for (idx,x) in enumerate(Iterators.product(["cautious","normal","aggressive"],[27.,31.,35.],[4.]))][:]
	ps = CarState[]
	push!(ps,CarState(1.,1,3.,0,bs[1],0))
	push!(ps,CarState(1.,1,3.,0,bs[2],0))
	push!(ps,CarState(1.,1,5.,0,bs[1],0))
	push!(ps,CarState(1.,1,3.,1,bs[1],0))
	push!(ps,CarState(1.,2,3.,0,bs[1],0))
	push!(ps,CarState(2.,1,3.,0,bs[1],0))

	test_equality("CarState",ps)

end

function test_carstate_hashing()
	#println("\t\tTesting CarState hashing")
	bs = IDMMOBILBehavior[IDMMOBILBehavior(x[1],x[2],x[3],idx) for (idx,x) in enumerate(Iterators.product(["cautious","normal","aggressive"],[27.,31.,35.],[4.]))][:]
	ps = CarState[]
	push!(ps,CarState(1.,1,3.,0,bs[1],0))
	push!(ps,CarState(1.,1,3.,0,bs[2],0))
	push!(ps,CarState(1.,1,5.,0,bs[1],0))
	push!(ps,CarState(1.,1,3.,1,bs[1],0))
	push!(ps,CarState(1.,2,3.,0,bs[1],0))
	push!(ps,CarState(2.,1,3.,0,bs[1],0))

	test_hashing("CarState",ps)
end

function test_MLState_creation()
	println("\t\tTesting MLState creation")
	#nothing to do here--no constructors to test
end

function test_MLState_equality()
	#println("\t\tTesting MLState equality")
	bs = IDMMOBILBehavior[IDMMOBILBehavior(x[1],x[2],x[3],idx) for (idx,x) in enumerate(Iterators.product(["cautious","normal","aggressive"],[27.,31.,35.],[4.]))][:]
	cs = CarState[]
	push!(cs,CarState(1.,1,3.,0,bs[1],0))
	push!(cs,CarState(1.,1,3.,0,bs[2],0))
	push!(cs,CarState(1.,1,5.,0,bs[1],0))
	push!(cs,CarState(1.,1,3.,1,bs[1],0))
	push!(cs,CarState(1.,2,3.,0,bs[1],0))
	push!(cs,CarState(2.,1,3.,0,bs[1],0))

	ps = MLState[]
	push!(ps,MLState(false,0.0,0.0,2,5.,cs[1:1]))
	push!(ps,MLState(false,0.0,0.0,3,5.,cs[1:1]))
	push!(ps,MLState(false,0.0,0.0,2,4.,cs[1:1]))
	push!(ps,MLState(false,0.0,0.0,2,5.,cs[1:2]))
	push!(ps,MLState(true,0.0,0.0,2,5.,cs[1:2]))

	test_equality("MLState",ps)
end

function test_MLState_hashing()
	#println("\t\tTesting MLState hashing")
	bs = IDMMOBILBehavior[IDMMOBILBehavior(x[1],x[2],x[3],idx) for (idx,x) in enumerate(Iterators.product(["cautious","normal","aggressive"],[27.,31.,35.],[4.]))][:]
	cs = CarState[]
	push!(cs,CarState(1.,1,3.,0,bs[1],0))
	push!(cs,CarState(1.,1,3.,0,bs[2],0))
	push!(cs,CarState(1.,1,5.,0,bs[1],0))
	push!(cs,CarState(1.,1,3.,1,bs[1],0))
	push!(cs,CarState(1.,2,3.,0,bs[1],0))
	push!(cs,CarState(2.,1,3.,0,bs[1],0))

	ps = MLState[]
	push!(ps,MLState(false,0.0,0.0,2,5.,cs[1:1]))
	push!(ps,MLState(false,0.0,0.0,3,5.,cs[1:1]))
	push!(ps,MLState(false,0.0,0.0,2,4.,cs[1:1]))
	push!(ps,MLState(false,0.0,0.0,2,5.,cs[1:2]))
	push!(ps,MLState(true,0.0,0.0,2,5.,cs[1:2]))

	test_hashing("MLState",ps)

    # all crashed states should be the same
    @assert(MLState(true,0.0,0.0,3,5.,cs[1:2]) == MLState(true,0.0,0.0,3,4.,cs[1:2]))
end

function test_MLAction_creation()
	println("\t\tTesting MLAction creation")
	#nothing to do here--no constructors to test
	#as = [MLAction(x[1],x[2]) for x in Iterators.product([-1;0;1],[-1;0;1])][:]
end

function test_MLAction_equality()
	as = [MLAction(x[1],x[2]) for x in Iterators.product([-1;0;1],[-1;0;1])][:]
	test_equality("MLAction",as)
end

function test_MLAction_hashing()
	as = [MLAction(x[1],x[2]) for x in Iterators.product([-1;0;1],[-1;0;1])][:]
	test_hashing("MLAction",as)
end

function test_MLObs_creation()
	println("\t\tTesting MLObs creation")
	#nothing to do here--no constructors to test
end

function test_MLMDP_creation()
	println("\t\tTesting OriginalMDP type creation")
	p = OriginalMDP()
end

function test_n_state()
	println("\t\tTesting n_states")
	pp = PhysicalParam(2,nb_vel_bins=8,lane_length=2.5) #2.5=>10
	p = OriginalMDP(nb_cars=1,phys_param=pp)
	p.nb_col = 4
	@assert(n_states(p) == (4*8)*(4*10*8*3*9+1))
	#idk if its worht it to do more
end

function test_repr()
    println("\t\tTesting the repr() function")
    bs = IDMMOBILBehavior[IDMMOBILBehavior(x[1],x[2],x[3],idx) for (idx,x) in enumerate(Iterators.product(["cautious","normal","aggressive"],[27.,31.,35.],[4.]))][:]
	c = CarState(1.,1,3.,0,bs[1],0)
    c2 = eval(Meta.parse(repr(c)))
    @test isa(c2, CarState)
end

#don't think you need tests for n_actions, the *Space types, and the associated domain() and length() functions...

function test_mdp_types()
	println("\tTesting MDP Type Units...")
	test_behavior_model_creation()
	test_behavior_model_equality()
	test_behavior_model_hashing()
	test_carstate_creation()
	test_carstate_equality()
	test_carstate_hashing()
	test_MLState_creation()
	test_MLState_equality()
	test_MLState_hashing()
	test_MLAction_creation()
	test_MLAction_equality()
	test_MLAction_hashing()
    test_repr()
end


text:
""" Functions for the types of random variate draws often used in individual-based epidemiological models.
"""
module Sampling

export bernoulli_sample, choose, delay_geom_sample

using Distributions: Exponential, Geometric, cdf
using Random: randsubseq
using StatsBase: sample

function rate_to_prob(x::AbstractFloat)
    0.0 <= x || throw(ArgumentError("rate $x not in [0,Inf)"))
    cdf(Exponential(), x)
end

"""
    delay_geom_sample(n::Integer, rate::AbstractFloat, dt::AbstractFloat)

Sample time steps until an event fires given a `rate` and `dt` where the waiting time
follows a Geometric distribution. The minimum delay is `1`, otherwise there could be
"instantaneous" events, which are better (more efficiently) simulated as _processes_.
"""
function delay_geom_sample(n::Integer, rate::AbstractFloat, dt::AbstractFloat)
    prob = rate_to_prob(rate * dt)
    rand(Geometric(prob), n) .+ 1
end

function delay_geom_sample(n::Integer, rate::AbstractVector{T}, dt::AbstractFloat) where {T <: AbstractFloat}
    n == length(rate) || throw(ArgumentError("number of draws must be equal to the length of 'rate'"))
    out = Vector{Int64}(undef, n)
    for i = 1:n
        prob = rate_to_prob(rate[i] * dt)
        out[i] = rand(Geometric(prob)) + 1
    end
    return out
end

"""
    bernoulli_sample(target::AbstractVector, prob::AbstractFloat)

Sample without replacement from `target` with success probability `prob`.
"""
function bernoulli_sample(target::T, prob::AbstractFloat) where {T <: Integer}
    0.0 <= prob <= 1.0 || throw(ArgumentError("probability $prob not in [0,1]"))
    if rand() < prob
        return [target]
    else
        return T[]
    end
end

function bernoulli_sample(target::AbstractVector, prob::AbstractFloat)
    randsubseq(target, prob)
end

"""
    bernoulli_sample(target::AbstractVector, rate::AbstractFloat, dt::AbstractFloat)

Sample without replacement from `target` with success probability calculated
from ``1 - e^{-\\mathrm{rate} * \\mathrm{dt}}``
"""
function bernoulli_sample(target::T, rate::AbstractFloat, dt::AbstractFloat) where {T <: Integer}
    prob = rate_to_prob(rate * dt)
    if rand() < prob
        return [target]
    else
        return T[]
    end
end

function bernoulli_sample(target::AbstractVector, rate::AbstractFloat, dt::AbstractFloat)
    randsubseq(target, rate_to_prob(rate * dt))
end

"""
    bernoulli_sample(target::AbstractVector{T}, prob::Vector{R})

Sample without replacement from `target` where each element has a unique success probability
given in the vector `prob`.
"""
function bernoulli_sample(target::AbstractVector{T}, prob::Vector{R}) where {T <: Integer, R <: AbstractFloat}
    length(target) == length(prob) || throw(ArgumentError("target and prob not of equal length"))
    samp = Vector{T}(undef, length(target))
    runif = rand(length(target))
    j = 0
    for i = 1:length(target)
        0.0 <= prob[i] <= 1.0 || throw(ErrorException("probability $(prob[i]) not in [0,1]"))
        if runif[i] < prob[i]
            j += 1
            samp[j] = target[i]
        end
    end

    # return sampled set
    if j > 0
        return samp[1:j]
    else
        return T[]
    end
end

"""
    bernoulli_sample(target::AbstractVector{T}, rate::Vector{R}, dt::AbstractFloat)

Sample without replacement from `target` where each element's success probability is calculated
from ``1 - e^{-\\mathrm{rate} * \\mathrm{dt}}``.
"""
function bernoulli_sample(target::AbstractVector{T}, rate::Vector{R}, dt::AbstractFloat) where {T <: Integer, R <: AbstractFloat}
    length(target) == length(rate) || throw(ArgumentError("target and rate not of equal length"))
    prob = map((x) -> rate_to_prob(x), rate * dt)
    samp = Vector{T}(undef, length(target))
    runif = rand(length(target))
    j = 0
    for i = 1:length(target)
        if runif[i] < prob[i]
            j += 1
            samp[j] = target[i]
        end
    end

    # return sampled set
    if j > 0
        return samp[1:j]
    else
        return T[]
    end
end

"""
    choose(target::T, K::Integer) where {T <: Integer}
    choose(target::AbstractVector{T}, K::Integer) where {T <: Integer}

Return a vector of size `K` with that number of random elements selected
without replacement from `target`.
"""
function choose(target::T, K::Integer) where {T <: Integer}
    if K > 0 # should K > 1 be an error?
        return target
    else
        return T[]
    end
end

function choose(target::AbstractVector{T}, K::Integer) where {T <: Integer}
    if K > 0
        target[sample(1:length(target), K, replace = false)]
    else
        return T[]
    end
end

end # end module


text:
struct Branch
    id::String
    #src::Int # bus id from which emits the branch
    #dst::Int # bus it to which goes the branch

    # Metier
    limit::Float64
end

function get_limit(branch::Branch)
    return branch.limit
end

################
## INFO / LOG ##
################

function get_info(branch::Branch)::String
    info::String =
        string(branch.id)
        # * ":" *
        # string(branch.src) *
        # "->" *
        # string(branch.dst)
    return info
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
#= -------------------------------------------------------------------------
# @file image_processing.jl
#
# @date 07/01/15 17:51:05
# @author Martin Noblia
# @email martin.noblia@openmailbox.org
#
# @brief
# Images processing functions in Julia language
# @detail
#
  Licence:
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License

---------------------------------------------------------------------------=#
module ImageProcessing
#-------------------------------------------------------------------------
# imports
#-------------------------------------------------------------------------
using Images, Color, Docile

#-------------------------------------------------------------------------
# functions
#-------------------------------------------------------------------------
"""
Conpute the threshold of a RGB Image

Input:
-----

img: Image{RGB} (from Image.jl)

c: RGB center value(like a center of mass)

s: Array with the absolute distances which form a desire cube

example:
-------
```julia
using Images, ImageView, Color
c = RGB(117, 66, 38) /255 # normalized
s = [113, 80, 77] / 255
bin = RGB_threshold(img, c, s)

```
"""
function RGB_threshold(img::Image, c::RGB, s::Array)
    x, y = size(img)
    r = red(img)
    g = green(img)
    b = blue(img)
    r_thresh = (c.r - s[1]/2) .< r .< (c.r + s[1]/2)
    g_thresh = (c.g - s[2]/2) .< g .< (c.g + s[2]/2)
    b_thresh = (c.b - s[3]/2) .< b .< (c.b + s[3]/2)
    return r_thresh & g_thresh & b_thresh
end

"""
Colored a label components Array

Input:
-----
label: Array{Int64, 2}

Output:
------

label_color: Array{RGB{Float64}}
"""
function colored_labels(label::Array{Int, 2})
    x, y = size(label)
    nums_labels = maximum(label)
    r = [RGB(a, 0, 0) for a in rand(nums_labels)]
    g = [RGB(0, b, 0) for b in rand(nums_labels)]
    b = [RGB(0, 0, c) for c in rand(nums_labels)]
    colors = r + g + b
    label_color = zeros(RGB{Float64}, x, y)
    for j in 1:y
        for i in 1:x
            for k in 1:nums_labels
                if(label[i,j] == k)
                    label_color[i,j] = colors[k]
                end
            end
        end
    end
    return label_color

end

"""
Transform a image in chromatics coordinates

Input:
-----

img: Image{RGB} (from Images.jl)

Output:
------

out: Image{RGB}
"""
function chromatics_coord(img::Image)

    x, y = size(img)
    a = zeros(img)
    for i in 1:x
        for j in 1:y
            r = img[i, j].r
            g = img[i, j].g
            b = img[i, j].b
            s = r + g + b
            if s == 0
                continue
            end

            a[i, j] = RGB(r / s, g / s, b / s)
        end
    end
    return a
end

"""
Conpute the white patch algorithm to a RGB image

Input:
-----

img: Image{RGB} (from Images.jl)

Output:
------

out: Image{RGB}
"""
function white_patch(img::AbstractArray)
    x, y = size(img)
    out = zeros(img)
    r_m = maximum(red(img))
    g_m = maximum(green(img))
    b_m = maximum(blue(img))

    for i in 1:x
        for j in 1:y
            r = img[i, j].r
            g = img[i, j].g
            b = img[i, j].b
            out[i, j] = RGB(r / r_m, g / g_m, b / b_m)
        end
    end
    return out
end

"""
Compute the histogram of a grayscale image or Vector

Input:
-----

img: Image{Gray} - Array{Int, 2} - Array{Float64, 2}

Output:
------

counts: Vector{Float64,1}
"""
function histogram(img)
    _, counts = hist(img[:], -1/256:1/256:1)

    s₁,s₂ = size(img_gray)
    return counts / (s₁ * s₂)
end

"""
Compute and applied the histogram equalization of a image

Input:
-----

img_gray: Image{Gray}, Array{Float64, 2}

Output:
------

img_eq: Image{Gray}
"""
function eq_hist(img_gray)
    h = histogram(img_gray) # calculate the histogram
    h_eq = cumsum(h) # calculate the cumulative sum
    gr = [Gray(i) for i in h_eq] # generated the equalized grayscale map
    A = float(data(img_gray)) # convert to float
    dataint = iround(Uint8, 254*A + 1 ) # convert 1-254
    img_eq = ImageCmap(dataint, gr)

    return img_eq
end

"""
Compute the false color map of a Image using a look-up table

Input:
-----

img: Image{Gray}

Output:
------

img_fc: Image{RGB}

"""
function fc(img)
    # LUT colors
    red = zeros(256)
    green = zeros(256)
    blue = zeros(256)

    # red LUT
    red[1:43] = 1.0
    red[43:85] = linspace(1, 0, 43)
    red[172:214] = linspace(0, 1, 43)
    red[214:end] = 1.0
    red = red * RGB(1, 0, 0)
    # green LUT
    green[1:43] = linspace(0, 1, 43)
    green[43:129] = 1.0
    green[129:171] = linspace(1, 0, 43)
    green = green * RGB(0, 1, 0)
    # blue LUT
    blue[86:128] = linspace(0, 1, 43)
    blue[129:214] = 1.0
    blue[214:end] = linspace(1, 0, 43)
    blue = blue * RGB(0, 0, 1)

    A = float(data(img)) # convert to float
    dataint = iround(Uint8, 254*A + 1 ) # convert
    false_color = red + green + blue # blend the colors
    img_fc = ImageCmap(dataint, false_color)

    return img_fc
end

"""
Compute the optimal segmentation value of a image(Gonzalez-Woods)

Input:
-----

img_gray: Image{Gray}

Output:
------

T: Optimal segmentation value: Float64
"""
function segment_Gonz(img_gray)
    T = 0.5 * (minimum(img_gray) + maximum(img_gray))
    flag = false
    while ~flag
        g = img_gray .< T
        T_next = 0.5 * (mean(img_gray[g]) + mean(img_gray[~g]) )
        flag = abs(T - T_next) < 0.5
        T = T_next
    end
    return T
end

end


text:
@testset "Dynamics constraints" begin
model = Cartpole()
n,m = size(model)
N = 11
dt = 0.1
Z = Traj([rand(n) for k = 1:N], [rand(m) for k = 1:N-1], fill(dt,N))

dyn = TO.DynamicsConstraint(model, N)
@test TO.evaluate(dyn, Z[1], Z[2]) ≈ discrete_dynamics(RK3, model, Z[1]) - state(Z[2])
F = TO.gen_jacobian(dyn)
F0 = zero(F)
@test TO.jacobian!(F, dyn, Z[1], Z[2]) == false
discrete_jacobian!(RK3, F0, model, Z[1])
@test F0 ≈ F
F .= 0
@test TO.jacobian!(F, dyn, Z[1], Z[2], 2) == true
@test F ≈ -Matrix(I,n,n+m)

dyn = TO.DynamicsConstraint{RK4}(model, N)
@test TO.evaluate(dyn, Z[1], Z[2]) ≈ discrete_dynamics(RK4, model, Z[1]) - state(Z[2])
F = TO.gen_jacobian(dyn)
F0 = zero(F)
@test TO.jacobian!(F, dyn, Z[1], Z[2]) == false
discrete_jacobian!(RK4, F0, model, Z[1])
@test F0 ≈ F
F .= 0
@test TO.jacobian!(F, dyn, Z[1], Z[2], 2) == true
@test F ≈ -Matrix(I,n,n+m)

@test length(dyn) == n
@test TO.widths(dyn) == (n+m,n)

dyn = TO.DynamicsConstraint{HermiteSimpson}(model, N)
conval = TO.ConVal(n, m, dyn, 1:N-1)
TO.evaluate!(conval, Z)
X = states(Z)
U = controls(Z)
Xm = dyn.xMid
fVal = dyn.fVal
fValm = dynamics(model, Xm[1], 0.5*(U[1] + U[2]))
@test conval.vals[1] ≈ X[1] - X[2] + dt * (fVal[1] + 4*fValm + fVal[2])/6

@test length(dyn) == n
@test TO.widths(dyn) == (n+m,n+m)
@test size(conval.jac[2]) == (n,n+m)
@test TO.upper_bound(dyn) == zeros(n)
@test TO.lower_bound(dyn) == zeros(n)
end


text:
using Revise, Pkg

#TODO: Check daily data, debug whatever is wrong w/ the table
#if this is hard maybe do the summary table and walk from there
#other idea is to try valuewieghting

if  pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

Pkg.activate(pwd())
using LeverageChanges

#LeverageChanges.recompressfiles()
#LeverageChanges.recompressfile("CRSP-D")
#LeverageChanges.recompressfile("CRSP-M")

#refreshcomp (PreprocessCOMP): refreshes the compustat data from the csv.lz4 file
#refreshcrsp (PreprocessCRSP): refreshes the crsp data from the csv.lz4 file
#refreshmerge[1] (MergeCRSPCOMP): refreshes the compustat-crsp merge.
#refreshfactors (FormFactors): refreshes the compustat-crsp merge.
#[1]True if refreshcomp, refreshcrsp, or validatemerged are true
#[2]True if refreshcomp, refreshcrsp, refreshmerge, or validatemerged are true
#hello
#usequarterlydates

function leveragechangesscript()
  #LeverageChanges.recompressfile("CCM-A")
  #LeverageChanges.recompressfile("CRSP-D-I")
  #LeverageChanges.recompressfile("CRSP-D-II")
  #LeverageChanges.recompressfile("CRSP-M")
  #LeverageChanges.recompressfile("CCM")
  #LeverageChanges.recompressfile("COMP-A")

  LeverageChanges.replicate(
    refreshcomp = false,
    refreshtvolmcapcrsp = false,
    refreshtretcrsp = false,
    refreshdataseries=false,
    refreshportfolios=true,
    parallel = true,
    trialrunonly = true,
    table1 = true)
end

#LeverageChanges.archivefile("lev-daily.jls.lz4")
#LeverageChanges.unarchivefile("lev-daily.jls.lz4")
@time leveragechangesscript()

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using OrdinaryDiffEq, Zygote, DiffEqSensitivity, Test

p_model = [1f0]

u0 = Float32.([0.0])

function dudt(du, u, p, t)
    du[1] = p[1]
end

prob = ODEProblem(dudt,u0,(0f0,99.9f0))

function predict_neuralode(p)
    _prob = remake(prob,p=p)
    Array(solve(_prob,Tsit5(), saveat=0.1))
end

loss(p) = sum(abs2,predict_neuralode(p))/length(p)

p_model_ini = copy(p_model)

@test !iszero(Zygote.gradient(loss,p_model_ini)[1])

## https://github.com/SciML/DiffEqSensitivity.jl/issues/675

u0 = Float32[2.0; 0.0] # Initial condition
p = [-0.1 2.0; -2.0 -0.1]

datasize = 30 # Number of data points
tspan = (0.0f0, 1.5f0) # Time range
# tsteps = range(tspan[1], tspan[2], length = datasize) # Split time range into equal steps for each data point
tsteps = (rand(datasize) .* (tspan[2] - tspan[1])  .+ tspan[1]) |> sort

function f(du,u,p,t)
    du .= p*u
end

function loss(p)
    prob = ODEProblem(f,u0,tspan,p)
    sol = solve(prob,Tsit5(),saveat=tsteps,sensealg=InterpolatingAdjoint())
    sum(sol)
end

Zygote.gradient(loss, p)[1]

@test !(Zygote.gradient(loss, p)[1] .|> iszero |> all)


text:
include(find_in_path("PatternDispatch.jl"))

module TestAmbWarning
using PatternDispatch

@pattern f(x::Int, y) = 1
@pattern f(x, y::Int) = 2

end


text:
const CSlice = Union{Colon,Slice}
const IterIndex = Union{Int,UnitRange{Int},CSlice}

# isindex == true  => want the indexes (keys) of the array
# isindex == false => want the values of the array
# isstored == true  => visit only stored entries
# isstored == false => visit all indexes
struct ArrayIndexingWrapper{A, I<:Tuple{Vararg{IterIndex}}, isindex, isstored}
    data::A
    indexes::I
    function ArrayIndexingWrapper{A,I,IND,STO}(data::A, indexes::I) where {A,I,IND,STO}
        ind = to_slices(indexes, data)
        new{A,typeof(ind),IND,STO}(data, ind)
    end
end

to_slice(::Colon, data, d) = Slice(OneTo(size(data, d)))
to_slice(a, data, d) = a
to_slices(ind::Tuple, data) = tuple(to_slice.(ind, Ref(data), 1:length(ind))...)

# Base.keys(W::ArrayIndexingWrapper) = CartesianIndices(W.indexes)

# Internal type for storing instantiated index iterators but returning
# array values
struct ValueIterator{A<:AbstractArray,I}
    data::A
    iter::I
end

struct SyncedIterator{I,O<:Tuple,F<:Tuple{Vararg{Function}}}
    iter::I
    items::O
    itemfuns::F
end

# declare that an array/iterhint should not control which index
# positions are visited, but only follow the lead of other objects
struct Follower{T}
    value::T
end

const ArrayOrWrapper = Union{AbstractArray,ArrayIndexingWrapper}
const AllElements{A,I,isindex} = Union{AbstractArray,ArrayIndexingWrapper{A,I,isindex,false}}
const StoredElements{A,I,isindex} = ArrayIndexingWrapper{A,I,isindex,true}

# storageorder has to be type-stable because it controls the output of
# sync, which is used in iteration
abstract type StorageOrder end
struct FirstToLast <: StorageOrder end
struct OtherOrder{p} <: StorageOrder end
struct NoOrder <: StorageOrder end  # combination of reshape+permutedims=>undefined

# For iterating over the *values* of an array in column-major order
struct FirstToLastIterator{N,AA}
    parent::AA
    itr::CartesianIndices{N}
end

# Contiguous ranges
abstract type Contiguity end
struct Contiguous <: Contiguity end
struct NonContiguous <: Contiguity end
struct MaybeContiguous <: Contiguity end  # intermediate type used in assessing contiguity

# Contiguous cartesian ranges. Sometimes needed for IndexCartesian arrays.
struct ContigCartIterator{N}
    arrayrange::CartesianIndices{N}
    columnrange::CartesianIndices{N}
end

# Note: sparse types are in sparse.jl


text:
function sample!(fin_el::Tetrahedron, s::Matrix{Float64} , wt::Vector{Float64})
  #
  # This subroutine returns the local coordinates and weighting coefficients
  # of the integrating points.
  #

  root3 = 1.0 / sqrt(3.0)
  r15 = 0.2 * sqrt(15.0)
  nip = size(s,1)

  w = [5.0/9.0, 8.0/9.0, 5.0/9.0]
  v = [5.0/9.0*w; 8.0/9.0*w; 5.0/9.0*w]

  if nip == 1
    s[1,1] = 0.25
    s[1,2] = 0.25
    s[1,3] = 0.25
    wt[1] = 1//6
  elseif nip == 4
    s[1,1] = 0.58541020
    s[1,2] = 0.13819660
    s[1,3] = s[1,2]
    s[2,2] = s[1,1]
    s[2,3] = s[1,2]
    s[2,1] = s[1,2]
    s[3,3] = s[1,1]
    s[3,1] = s[1,2]
    s[3,2] = s[1,2]
    s[4,1] = s[1,2]
    s[4,2] = s[1,2]
    s[4,3] = s[1,2]
    wt[1:4] = 0.25/6.0
  elseif nip == 5
    s[1,1:3] = 0.25
    s[2,1] = 0.5
    s[2,2] = 1 // 6
    s[2,3] = s[2,2]
    s[3,2] = 0.5
    s[3,3] = 1 // 6
    s[3,1] = s[3,3]
    s[4,3] = 0.5
    s[4,1] = 1 // 6
    s[4,2] = s[4,1]
    s[5,1] = 1 // 6
    s[5,2] = s[5,1]
    s[5,3] = s[5,1]
    wt[1] = -0.8/6.0
    wt[2] = 9 // 120
    wt[3:5] = wt[2]
    wt[:] = wt[:]/6.0
  else
    println("Wrong number of integrating points for a tetrahedron.")
  end
end


text:
## ----------------------------------------------------------------------------
function plot_to_img(p::AbstractPlot; ext = "png", clear = true)
    fn = string(tempname(), ".", ext)
    try
        Plots.savefig(p, fn)
        return FileIO.load(fn)
    finally
        clear && rm(fn; force = true)
    end
end

## ----------------------------------------------------------------------------
function _make_mat(imgs)
    isempty(imgs) && error("imgs is empty")
    length(unique!(size.(imgs))) != 1 && error("All imgs must have the same size")

    w, h = size(first(imgs))
    d = length(imgs)
    img_mat = Array{RGB{N0f8}}(undef, w, h, d)
    @views for (i, img) in enumerate(imgs)
        img_mat[:, :, i] = img[:, :]
    end
    img_mat
end

## ----------------------------------------------------------------------------
function __sgif(imgs::Vector, arg, args...;
        fps = 10.0
    )
    file = dfname(arg, args...)
    !endswith(file, ".gif") && error("filename must end with .gif")
    mat = _make_mat(imgs)
    imgs = nothing # save memory ?
    FileIO.save(file, mat; fps)
    return file
end

function _sgif(ps::Vector{T}, arg, args...;
        kwargs...
    ) where {T<:AbstractPlot}
    imgs = plot_to_img.(ps)
    __sgif(imgs, arg, args...; kwargs...)
end

function _sgif(sourcepaths::Vector{String}, arg, args...;
        kwargs...
    )
    imgs = FileIO.load.(sourcepaths)
    __sgif(imgs, arg, args...; kwargs...)
end

sgif(ps, arg, args...; kwargs...) = _sgif(ps, arg, args...; kwargs...)
sgif(ps, fn::String = string(tempname(), ".gif"); kwargs...) = _sgif(ps, fn; kwargs...)

## ----------------------------------------------------------------------------
function make_group_gif(freedim, sourcedir::String;
        filter = (filename) -> true,
        sortby = identity,
        destdir = sourcedir,
        verbose = true,
        fps = 3.0
    )

    dfg = group_files(freedim, sourcedir; filter)
    gifs = []
    for ((k, head, params, ext), elms) in dfg
        gifname = dfname(head..., params, "gif")
        giffile = joinpath(destdir, gifname)

        selms = sort!(collect(elms); by = (p) -> sortby(first(p)))
        paths = [file for (kv, file) in selms]

        sgif(paths, giffile; fps)
        verbose && @info("Gif produced", gifname)
        push!(gifs, giffile)
    end
    gifs
end


text:
#RSDeltaSigmaPort: Implement missing core functionality using ControlSystems.jl
#-------------------------------------------------------------------------------

#=TODO (Goal)
Want to not depend on ControlSystems.jl (slows down compile time significantly).

UPDATE: Much better in Julia 1.6 - especially using Revise.jl.
Maybe ok to depend on ControlSystems.jl.
=#

_c2d(sys, Ts, method=:zoh) = ControlSystems.c2d(sys, Ts, method)

_ss(A,B,C,D) = ControlSystems.ss(A,B,C,D) #Ts=0.0 default: continuous-time system
_ss(A,B,C,D,Ts) = ControlSystems.ss(A,B,C,D,Ts)

function _zpk(sys::ControlSystems.StateSpace)
	H=ControlSystems.zpk(sys)
	z, p, k = ControlSystems.zpkdata(H)

	#Cannot access sys.Ts for continuous time:
	#(use sys because Ts does not get ported over in zpk() for some reason)
	Ts = isdiscrete(sys) ? sys.Ts : 0.0
	return _zpk(z, p, k, Ts)
end

function _zp2ss(z, p, k)
	H = ControlSystems.zpk(z, p, k)
	return ControlSystems.ss(H)
end

function _zp2tf(z, p, k)
	H = ControlSystems.zpk(z, p, k)
	_num = ControlSystems.num(H)
	_den = ControlSystems.den(H)
	return (_num, _den)
end

function _freqz(num, den, w)
	H = ControlSystems.tf(num, den, 1)
	return ControlSystems.freqresp(H, w)
end

_tf(num, den, ts) = ControlSystems.tf(num, den, ts)

_minreal(sys, tol) = ControlSystems.minreal(sys, tol)

function _impulse(sys, t)
	y, t, x = ControlSystems.impulse(sys, t)
	return (y, t, x)
end

function _lsim(sys, u)
	if !isdiscrete(sys)
		throw("_lsim: must specify t-vector")
	end
	Ts = sys.Ts
	t = range(0, step=Ts, length=length(u))
	y, t, x = ControlSystems.lsim(sys, u, t)
	return (y, t, x)
end

#Last line


text:
###################################################
## network/stats/statistics.jl
## Statistic computations on results (errors...)
###################################################

"""
    `tripsLogError`: compute log timings error
    - `timings`: timing reference
    - `trips` : list of trip data
"""
function tripsLogError(timings::NetworkTimings, trips::Vector{NetworkTrip})
    tt = getPathTimes(timings)
    times = timings.times
    error = 0.
    weight = 0.
    for t in trips
        if t.roadProj
            error += t.weight * log(t.time / (tt[t.orig[2], t.dest[1]] + t.orig[3] * times[t.orig[1], t.orig[2]] + t.dest[3] * times[t.dest[1], t.dest[2]]))^2
        else
            err = t.weight * log(t.time / tt[t.orig[2],t.dest[1]])^2
            if err < Inf
                error += err
                weight += t.weight
            end
        end
    end
    return error/weight
end

"""
    `tripsLogBias`: compute log bias between trimings and trips
    - `timings`: timing reference
    - `trips` : list of trip data
"""
function tripsLogBias(timings::NetworkTimings, trips::Vector{NetworkTrip})
    tt = getPathTimes(timings)
    times = timings.times
    bias = 0.
    weight = 0.
    for t in trips
        if t.roadProj
            bias += t.weight * (log(tt[t.orig[2], t.dest[1]] + t.orig[3] * times[t.orig[1], t.orig[2]] + t.dest[3] * times[t.dest[1], t.dest[2]]) - log(t.time))
        else
            err = t.weight * log(tt[t.orig[2], t.dest[1]]/t.time)
            if abs(err) < Inf
                bias += err
                weight += t.weight
            end
        end
    end
    return bias/weight
end

"""
    `allPathsLogError`: compute log timings error for all paths
    - `timingsRef`: reference times
    - `timingsNew`: times to compare with
"""
function allPathsLogError(timingsRef::NetworkTimings, timingsNew::NetworkTimings)
    tt1 = getPathTimes(timingsRef)
    tt2 = getPathTimes(timingsNew)

    res = log.(tt1./tt2).^2
    #remove NAN where o==d
    for i in 1:size(tt1)[1]
        res[i,i] = 0.
    end
    return sum(res)/(length(tt1)-size(tt1)[1])
end

"""
    `allPathsLogBias`: compute log bias of timings for all paths
    - `timingsRef`: reference times
    - `timingsNew`: times to compare with
"""
function allPathsLogBias(timingsRef::NetworkTimings, timingsNew::NetworkTimings)
    tt1 = getPathTimes(timingsRef)
    tt2 = getPathTimes(timingsNew)

    res = log.(tt2./tt1)
    #remove NAN where o==d
    for i in 1:size(tt1)[1]
        res[i,i] = 0.
    end
    return sum(res)/(length(tt1)-size(tt1)[1])
end

"""
    `roadTimeLogError`: compute log road time error
    - `timingsRef`: reference times
    - `timingsNew`: times to compare with
"""
function roadTimeLogError(timingsRef::NetworkTimings, timingsNew::NetworkTimings)
    g = timingsRef.network.graph
    t1 = timingsRef.times
    t2 = timingsNew.times
    error = 0.
    for o in vertices(g), d in outneighbors(g,o)
        error += log(t1[o,d]/t2[o,d])^2
    end
    return error/ne(g)
end

"""
    `roadTimeLogBias`: compute log bias of road times
    - `timingsRef`: reference times
    - `timingsNew`: times to compare with
"""
function roadTimeLogBias(timingsRef::NetworkTimings, timingsNew::NetworkTimings)
    g = timingsRef.network.graph
    t1 = timingsRef.times
    t2 = timingsNew.times
    error = 0.
    for o in vertices(g), d in outneighbors(g,o)
        error += log(t2[o,d]/t1[o,d])
    end
    return error/ne(g)
end

"""
    `nwTripsLogError`: compute log error on NetworkTrips inside NetworkData
"""
nwTripsLogError(timings::NetworkTimings, nd::NetworkData) = tripsLogError(timings, nd.trips)

"""
    `nwTripsLogBias`: compute log bias on NetworkTrips inside NetworkData
"""
nwTripsLogBias(timings::NetworkTimings, nd::NetworkData) = tripsLogBias(timings, nd.trips)

"""
    `tripsRealLogError` : computes log error on trips using their real times as computed from timingsRef
    Used to see if method denoises input data.
"""
function tripsRealLogError(timingsRef::NetworkTimings, timingsNew::NetworkTimings, trips::Vector{NetworkTrip})
    tt1 = getPathTimes(timingsRef)
    times1 = timingsRef.times
    tt2 = getPathTimes(timingsNew)
    times2 = timingsNew.times
    error = 0.
    count = 0
    for t in trips
        if t.roadProj
            error += log((tt1[t.orig[2], t.dest[1]] + t.orig[3] * times1[t.orig[1], t.orig[2]] + t.dest[3] * times1[t.dest[1], t.dest[2]])/(tt2[t.orig[2], t.dest[1]] + t.orig[3] * times2[t.orig[1], t.orig[2]] + t.dest[3] * times2[t.dest[1], t.dest[2]]))^2
        else
            err = log(tt1[t.orig[2],t.dest[1]] / tt2[t.orig[2],t.dest[1]])^2
            if err < Inf
                error += err
            else
                count += 1
            end
        end
    end
    return error/(length(trips)-count)
end

"""
    `tripsRealLogBias` : computes log bias on trips using their real times as computed from timingsRef
    Used to see if method denoises input data. Bias is (predicted - real).
"""
function tripsRealLogBias(timingsRef::NetworkTimings, timingsNew::NetworkTimings, trips::Vector{NetworkTrip})
    tt1 = getPathTimes(timingsRef)
    times1 = timingsRef.times
    tt2 = getPathTimes(timingsNew)
    times2 = timingsNew.times
    error = 0.
    count = 0
    for t in trips
        if t.roadProj
            error += log(tt2[t.orig[2], t.dest[1]] + t.orig[3] * times2[t.orig[1], t.orig[2]] + t.dest[3] * times2[t.dest[1], t.dest[2]]) - log(tt1[t.orig[2], t.dest[1]] + t.orig[3] * times1[t.orig[1], t.orig[2]] + t.dest[3] * times1[t.dest[1], t.dest[2]])
        else
            err = log(tt1[t.orig[2],t.dest[1]] / tt2[t.orig[2],t.dest[1]])
            if abs(err) < Inf
                error += err
            else
                count += 1
            end
        end
    end
    return error/(length(trips)-count)
end

"""
    `networkTripsRealLogError`: compute log error on NetworkTrips inside NetworkData
"""
nwTripsRealLogError(timingsRef::NetworkTimings, timingsNew::NetworkTimings, nd::NetworkData) = tripsRealLogError(timingsRef, timingsNew, nd.trips)

"""
    `networkTripsRealLogBias`: compute log bias on NetworkTrips inside NetworkData
"""
nwTripsRealLogBias(timingsRef::NetworkTimings, timingsNew::NetworkTimings, nd::NetworkData) = tripsRealLogBias(timingsRef, timingsNew, nd.trips)


text:
using KitBase

cd(@__DIR__)
ks, ctr, a1face, a2face, simTime = KitBase.initialize("cavity.txt")
simTime = KitBase.solve!(ks, ctr, a1face, a2face, simTime)

# equivalent low-level procedures
using KitBase.ProgressMeter: @showprogress
simTime = 0.0
res = zeros(4)
dt = KitBase.timestep(ks, ctr, simTime)
nt = floor(ks.set.maxTime / dt) |> Int
@showprogress for iter = 1:nt
    KitBase.reconstruct!(ks, ctr)
    KitBase.evolve!(ks, ctr, a1face, a2face, dt; mode = :kfvs, bc = :maxwell)
    KitBase.update!(ks, ctr, a1face, a2face, dt, res; coll = :bgk, bc = :maxwell)
end

# lower-level backend
@showprogress for iter = 1:nt
    @inbounds Threads.@threads for j = 1:ks.pSpace.ny
        for i = 2:ks.pSpace.nx
            KitBase.flux_kfvs!(
                a1face[i, j].fw,
                a1face[i, j].fh,
                a1face[i, j].fb,
                ctr[i-1, j].h,
                ctr[i-1, j].b,
                ctr[i, j].h,
                ctr[i, j].b,
                ks.vSpace.u,
                ks.vSpace.v,
                ks.vSpace.weights,
                dt,
                ks.ps.dy[i, j],
            )
        end
    end

    # vertical flux
    vn = ks.vSpace.v
    vt = -ks.vSpace.u
    @inbounds Threads.@threads for j = 2:ks.pSpace.ny
        for i = 1:ks.pSpace.nx
            KitBase.flux_kfvs!(
                a2face[i, j].fw,
                a2face[i, j].fh,
                a2face[i, j].fb,
                ctr[i, j-1].h,
                ctr[i, j-1].b,
                ctr[i, j].h,
                ctr[i, j].b,
                vn,
                vt,
                ks.vSpace.weights,
                dt,
                ks.ps.dy[i, j],
            )
            a2face[i, j].fw .= KitBase.global_frame(a2face[i, j].fw, 0.0, 1.0)
        end
    end

    # boundary flux
    @inbounds Threads.@threads for j = 1:ks.pSpace.ny
        KitBase.flux_boundary_maxwell!(
            a1face[1, j].fw,
            a1face[1, j].fh,
            a1face[1, j].fb,
            [1.0, 0.0, 0.0, 1.0],
            ctr[1, j].h,
            ctr[1, j].b,
            ks.vSpace.u,
            ks.vSpace.v,
            ks.vSpace.weights,
            ks.gas.K,
            dt,
            ks.ps.dy[1, j],
            1.0,
        )

        KitBase.flux_boundary_maxwell!(
            a1face[ks.pSpace.nx+1, j].fw,
            a1face[ks.pSpace.nx+1, j].fh,
            a1face[ks.pSpace.nx+1, j].fb,
            [1.0, 0.0, 0.0, 1.0],
            ctr[ks.pSpace.nx, j].h,
            ctr[ks.pSpace.nx, j].b,
            ks.vSpace.u,
            ks.vSpace.v,
            ks.vSpace.weights,
            ks.gas.K,
            dt,
            ks.ps.dy[ks.pSpace.nx, j],
            -1.0,
        )
    end

    @inbounds Threads.@threads for i = 1:ks.pSpace.nx
        KitBase.flux_boundary_maxwell!(
            a2face[i, 1].fw,
            a2face[i, 1].fh,
            a2face[i, 1].fb,
            [1.0, 0.0, 0.0, 1.0],
            ctr[i, 1].h,
            ctr[i, 1].b,
            vn,
            vt,
            ks.vSpace.weights,
            ks.gas.K,
            dt,
            ks.ps.dx[i, 1],
            1,
        )
        a2face[i, 1].fw .= KitBase.global_frame(a2face[i, 1].fw, 0.0, 1.0)

        KitBase.flux_boundary_maxwell!(
            a2face[i, ks.pSpace.ny+1].fw,
            a2face[i, ks.pSpace.ny+1].fh,
            a2face[i, ks.pSpace.ny+1].fb,
            [1.0, 0.0, -0.15, 1.0],
            ctr[i, ks.pSpace.ny].h,
            ctr[i, ks.pSpace.ny].b,
            vn,
            vt,
            ks.vSpace.weights,
            ks.gas.K,
            dt,
            ks.ps.dy[i, ks.pSpace.ny],
            -1,
        )
        a2face[i, ks.pSpace.ny+1].fw .=
            KitBase.global_frame(a2face[i, ks.pSpace.ny+1].fw, 0.0, 1.0)
    end

    # update
    @inbounds Threads.@threads for j = 1:ks.pSpace.ny
        for i = 1:ks.pSpace.nx
            KitBase.step!(
                ctr[i, j].w,
                ctr[i, j].prim,
                ctr[i, j].h,
                ctr[i, j].b,
                a1face[i, j].fw,
                a1face[i, j].fh,
                a1face[i, j].fb,
                a1face[i+1, j].fw,
                a1face[i+1, j].fh,
                a1face[i+1, j].fb,
                a2face[i, j].fw,
                a2face[i, j].fh,
                a2face[i, j].fb,
                a2face[i, j+1].fw,
                a2face[i, j+1].fh,
                a2face[i, j+1].fb,
                ks.vSpace.u,
                ks.vSpace.v,
                ks.vSpace.weights,
                ks.gas.K,
                ks.gas.γ,
                ks.gas.μᵣ,
                ks.gas.ω,
                ks.gas.Pr,
                ks.ps.dx[i, j] * ks.ps.dy[i, j],
                dt,
                zeros(4),
                zeros(4),
                :bgk,
            )
        end
    end
end

# visulization
using Plots
plot(ks, ctr)
#savefig("cavity.png")

# low-level backend
begin
    using Plots
    sol = zeros(4, ks.pSpace.nx, ks.pSpace.ny)
    for i in axes(sol, 2)
        for j in axes(sol, 3)
            sol[1:3, i, j] .= ctr[i, j].prim[1:3]
            sol[4, i, j] = 1.0 / ctr[i, j].prim[4]
        end
    end
    contourf(ks.pSpace.x[1:ks.pSpace.nx, 1], ks.pSpace.y[1, 1:ks.pSpace.ny], sol[4, :, :]')
end


text:
mutable struct CustomVariance
    meanratio::Float64
    sn::Float64 # sum nominator
    sd::Float64 # sum denominator
    n::Int # count
    m2::Float64 # var*(n-1)
    var::Float64 # weigthedVariance
    function CustomVariance(num, denom)
        n = length(num)
        @assert n == length(denom) # length(w)
        sn = sum(num)
        sd = sum(denom)
        meanratio = sn / sd
        m2 = zero(eltype(num))
        for i = 1:n
            @inbounds m2 += denom[i] * (num[i] / denom[i] - meanratio)^2
        end
        if n > 1
            var = m2 / (n - 1)
        else
            var = zero(eltype(num))
        end
        return new(meanratio, sn, sd, n, m2, var)
    end
    function CustomVariance()
        return new(0.0, 0.0, 0.0, 0, 0.0, 0.0)
    end
end

function removeSinglePoint!(a::CustomVariance, num, denom)
    sn = a.sn - num
    sd = a.sd - denom
    meanratio = sn / sd
    # merge: m2 = a.m2 + b.m2 + 2*(a.sn*(a.meanratio-meanratio)+§*(num/denom-meanratio)) + a.sd*(meanratio^2-a.meanratio^2) + denom*(meanratio^2-num/denom^2)
    # unmerge:
    a.m2 = a.m2 + 2 * (num * (a.meanratio - num / denom) + sn * (a.meanratio - meanratio)) + denom * ((num / denom)^2 - a.meanratio^2) + sd * (meanratio^2 - a.meanratio^2)
    a.sn = sn
    a.sd = sd
    a.n -= 1
    a.meanratio = meanratio
    a.var = a.m2 / (a.n - 1)
    return nothing
end

function addSinglePoint!(a::CustomVariance, num, denom)
    sn = a.sn + num
    sd = a.sd + denom
    a.n = a.n + 1
    meanratio = sn / sd
    a.m2 = a.m2 + 2 * (a.sn * (a.meanratio - meanratio) + num * (num / denom - meanratio)) + a.sd * (meanratio^2 - a.meanratio^2) + denom * (meanratio^2 - (num / denom)^2)
    a.sn = sn
    a.sd = sd
    a.meanratio = meanratio
    a.var = a.m2 / (a.n - 1)
    return nothing
end

function Base.merge!(a::CustomVariance, b::CustomVariance)
    # merges b into a
    # i.e. the points form b are 'added' to a
    sn = a.sn + b.sn
    sd = a.sd + b.sd
    n = a.n + b.n
    meanratio = sn / sd
    m2 = a.m2 + b.m2 + 2 * (a.sn * (a.meanratio - meanratio) + b.sn * (b.meanratio - meanratio)) + a.sd * (meanratio^2 - a.meanratio^2) + b.sd * (meanratio^2 - b.meanratio^2)
    var = m2 / (n - 1)
    a.sn = sn
    a.sd = sd
    a.n = n
    a.meanratio = meanratio
    a.m2 = m2
    a.var = var
    return nothing
end

function Base.merge(a::CustomVariance, b::CustomVariance)
    c = deepcopy(a)
    merge!(c, b)
    return c
end

"""
function unmerge!(a::CustomVariance,b::CustomVariance)
this function updates the mean and variance of a
The assumption is that the underlying data of b is a subset of a.
If this is not the case, the result of this function is most likely meaningless
That is to say "subtract" the series b from the series a and update the mean and variance of a accordingly
"""
function unmerge!(a::CustomVariance, b::CustomVariance)
    sn = a.sn - b.sn
    sd = a.sd - b.sd
    n = a.n - b.n
    meanratio = sn / sd
    # merge: m2 = a.m2 + b.m2 + 2*(a.sn*(a.meanratio-meanratio)+b.sn*(b.meanratio-meanratio)) + a.sd*(meanratio^2-a.meanratio^2) + b.sd*(meanratio^2-b.meanratio^2)
    # unmerge:
    m2 = a.m2 - b.m2 + 2 * (b.sn * (a.meanratio - b.meanratio) + sn * (a.meanratio - meanratio)) + b.sd * (b.meanratio^2 - a.meanratio^2) + sd * (meanratio^2 - a.meanratio^2)
    var = m2 / (n - 1)
    a.sn = sn
    a.sd = sd
    a.n = n
    a.meanratio = meanratio
    a.m2 = m2
    a.var = var
    return nothing
end


text:
module MeshIO

using GeometryBasics
using ColorTypes
using Printf

using GeometryBasics: raw, value, decompose_normals, convert_simplex
using FileIO: FileIO, @format_str, Stream, File, stream, skipmagic

import Base.show

include("io/off.jl")
include("io/ply.jl")
include("io/stl.jl")
include("io/obj.jl")
include("io/2dm.jl")
include("io/msh.jl")
include("io/gts.jl")
include("io/ifs.jl")

"""
    load(fn::File{MeshFormat}; pointtype=Point3f, uvtype=Vec2f,
         facetype=GLTriangleFace, normaltype=Vec3f)

"""
function load(fn::File{format}; element_types...) where {format}
    open(fn) do s
        skipmagic(s)
        load(s; element_types...)
    end
end

function save(fn::File{format}, msh::AbstractMesh) where {format}
    open(fn, "w") do s
        save(s, msh)
    end
end

if Base.VERSION >= v"1.4.2"
    include("precompile.jl")
    _precompile_()
end

# `filter(f, ::Tuple)` is not available on Julia 1.3
# https://github.com/JuliaLang/julia/pull/29259
function filtertuple(f, xs::Tuple)
    return @static if VERSION < v"1.4.0-DEV.551"
        Base.afoldl((ys, x) -> f(x) ? (ys..., x) : ys, (), xs...)
    else
        filter(f, xs)
    end
end

end # module

max_stars_repo_name: JuliaIO/MeshIO

max_issues_repo_name: JuliaIO/MeshIO

max_forks_repo_name: JuliaIO/MeshIO