text:
abstract type RestartStrategy end

# Tabu restart
mutable struct TabuRestart <: RestartStrategy
    index::Int
    tenure::Int
    limit::Int
    reset_percentage::Float64
end

function restart(tabu_strat,::Val{:tabu}; rp = 1.0, index = 1)
    limit = tenure(tabu_strat, :tabu) - tenure(tabu_strat, :pick)
    return TabuRestart(index, tenure(tabu_strat, :tabu), limit, rp)
end

function check_restart!(rs::TabuRestart; tabu_length)
    a = rs.index * (tabu_length + rs.limit - rs.tenure)
    b = (rs.index + 1) * rs.limit
    # a = tabu_length + rs.limit - rs.tenure
    # b = rs.limit
    if rand() ≤ a / b
        rs.index += 1
        return true
    end
    return false
end

# Restart sequences
mutable struct RestartSequence{F <: Function} <: RestartStrategy
    index::Int
    current::Int
    last_restart::Int
    next::F

    RestartSequence(seq) = new{typeof(seq)}(1, seq(1), 1, seq)
end

current(r) = r.current

function next!(r)
    r.index += 1
    r.current = r.next(r.index)
    r.last_restart = 1
    return r.current
end

inc_last!(rs) = rs.last_restart += 1

function check_restart!(rs::RestartSequence; tabu_length = nothing)
    proceed = rs.current > rs.last_restart
    proceed ? inc_last!(rs) : next!(rs)
    return !proceed
end

## Universal restart sequence

function oeis(n, b, ::Val{:A082850})
    m = log(b,n+1)
    return isinteger(m) ? Int(m) : oeis(n - (b^floor(m) - 1), :A082850)
end
oeis(n, b, ::Val{:A182105}) = b^(oeis(n, :A082850)-1)
oeis(n, ref::Symbol, b = 2) = oeis(n, b, Val(ref))
restart(::Any, ::Val{:universal}) = RestartSequence(n -> oeis(n, :A182105))

# Generic restart constructor
restart(tabu, strategy::Symbol) = restart(tabu, Val(strategy))


text:
using MatrixCompletion.MathLib
import LinearAlgebra

function create_symmetric_matrix(n)
    a = rand(n,n)*5
    return a+a'
end

function correct_output_sparseeigen(input,k)
    eigen_dcp = LinearAlgebra.eigen(input);
    eigen_val = eigen_dcp.values;
    eigen_vec = eigen_dcp.vectors;
    first_k_idx = Base.sortperm(eigen_val,rev=true)[1:k];
    return eigen_val[first_k_idx],eigen_vec[:,first_k_idx];
end

@testset "$(format("Math Library: Projection[SemidefiniteCone]"))" begin
  let
    for i in 1:50
      data = create_symmetric_matrix(100)
      rk = rand(2:30)
      tc = project(SemidefiniteCone(rank = rk), data)
      λ0, X0 = correct_output_sparseeigen(data, rk)
      tc_comp = X0 * LinearAlgebra.Diagonal(λ0) * X0'
      @test LinearAlgebra.norm(tc_comp - tc)^2 / LinearAlgebra.norm(tc_comp) ^ 2 < 0.02
    end
  end
end

ext: jl


text:
module Continuation
include("definitions.jl")
include("she.jl") # swift-hohenberg
include("she_parse.jl")  # parse swift-hohenberg data
include("spiral.jl") # Spiral chimeras
include("spiral_parse.jl")
include("newton.jl")
end


text:
const BigFloats = Union{BigFloat,Complex{BigFloat}}

function fft(x::AbstractVector{F}) where F<:Fun
    n,T = length(x),mapreduce(eltype,promote_type,x)
    if ispow2(n) return fft_pow2(x) end
    ks = range(zero(real(T)),stop=n-one(real(T)),length=n)
    Wks = exp(-im*convert(T,π)*ks.^2/n)
    xq,wq = x.*Wks,conj([exp(-im*convert(T,π)*n);reverse(Wks);Wks[2:end]])
    return Wks.*conv(xq,wq)[n+1:2n]
end

ifft(x::AbstractVector{F}) where {F<:Fun} = conj(fft(conj(x)))/length(x)
function ifft!(x::AbstractVector{F}) where F<:Fun
    y = conj(fft(conj(x)))/length(x)
    x[:] = y
    return x
end

function conv(u::StridedVector{F}, v::StridedVector) where F<:Fun
    nu,nv = length(u),length(v)
    n = nu + nv - 1
    np2 = nextpow2(n)
    pad!(u,np2),pad!(v,np2)
    y = ifft_pow2(fft_pow2(u).*fft_pow2(v))
    #TODO This would not handle Dual/ComplexDual numbers correctly
    T = promote_type(mapreduce(eltype,promote_type,u),mapreduce(eltype,promote_type,v))
    y = T<:Real ? real(y[1:n]) : y[1:n]
end

######################################################################
# TO BE DEPRECATED FOR V0.4 UPGRADE
######################################################################

# plan_fft for BigFloats (covers Laurent svfft)

plan_fft(x::Vector{F}) where {F<:Fun} = fft
plan_ifft(x::Vector{F}) where {F<:Fun} = ifft
plan_ifft!(x::Vector{F}) where {F<:Fun} = ifft

# Chebyshev transforms and plans for BigFloats
# no plan exists at the moment, so we make a dummy plan
plan_chebyshevtransform!(x::Vector{T}, kind) where {T<:BigFloats} =
    error("In-place variant not implemented for BigFloat")

plan_ichebyshevtransform!(x::Vector{T}, kind) where {T<:BigFloats} =
    error("In-place variant not implemented for BigFloat")

plan_chebyshevtransform(x::AbstractVector{T}, ::Val{kind}) where {T<:BigFloats,kind} =
    ChebyshevTransformPlan{T,kind,false,Nothing}()
plan_ichebyshevtransform(x::AbstractVector{T}, ::Val{kind}) where {T<:BigFloats,kind} =
    IChebyshevTransformPlan{T,kind,false,Nothing}()

#following Chebfun's @Chebtech1/vals2coeffs.m and @Chebtech2/vals2coeffs.m
function *(P::ChebyshevTransformPlan{T,1,false,Nothing}, x::AbstractVector{T}) where T<:BigFloats
    n = length(x)
    if n == 1
        x
    else
        w = [2exp(im*convert(T,π)*k/2n) for k=0:n-1]
        ret = w.*ifft([x;reverse(x)])[1:n]
        ret = T<:Real ? real(ret) : ret
        ret[1] /= 2
        ret
    end
end

function *(P::ChebyshevTransformPlan{T,2,false}, x::AbstractVector{T}) where T<:BigFloats
    n = length(x)
    if n == 1
        x
    else
        ret = ifft([x;x[end:-1:2]])[1:n]
        ret = T<:Real ? real(ret) : ret
        ret[2:n-1] *= 2
        ret
    end
end

#following Chebfun's @Chebtech1/vals2coeffs.m and @Chebtech2/vals2coeffs.m
function *(P::IChebyshevTransformPlan{T,1,false,Nothing}, x::AbstractVector{T}) where T<:BigFloats
    n = length(x)
    if n == 1
        x
    else
        w = [exp(-im*convert(T,π)*k/2n)/2 for k=0:2n-1]
        w[1] *= 2;w[n+1] *= 0;w[n+2:end] *= -1
        ret = fft(w.*[x;one(T);x[end:-1:2]])
        ret = T<:Real ? real(ret) : ret
    end
end
function *(P::IChebyshevTransformPlan{T,2,false}, x::AbstractVector{T}) where T<:BigFloats
    n = length(x)
    if n == 1
        x
    else
        ##TODO: make thread safe
        x[1] *= 2;x[end] *= 2
        ret = chebyshevtransform(x, Val(2))
        x[1] /=2;x[end] /=2
        ret[1] *= 2;ret[end] *= 2
        ret *= .5*(n-1)
        ret
    end
end

# Fourier space plans for BigFloat

plan_transform(sp::Fourier{D,R},x::AbstractVector{T}) where {T<:BigFloat,D,R} =
    TransformPlan{T,typeof(sp),false,Nothing}(sp,nothing)
plan_itransform(sp::Fourier{D,R},x::AbstractVector{T}) where {T<:BigFloat,D,R} =
    ITransformPlan{T,typeof(sp),false,Nothing}(sp,nothing)

function *(::TransformPlan{T,Fourier{D,R},false},x::AbstractVector{T}) where {T<:BigFloat,D,R}
    l = length(x); n = div(l+1,2)
    v = fft(x)
    rmul!(v,convert(T,2)/l)
    v[1] /= 2
    mod(l,2) == 1 ? ApproxFunBase.interlace(real(v[1:n]),-imag(v[2:n])) :
      [ApproxFunBase.interlace(real(v[1:n]),-imag(v[2:n]));-real(v[n+1])/2]
end

function *(::ITransformPlan{T,Fourier{D,R},false},x::AbstractVector{T}) where {T<:BigFloat,D,R}
    l = length(x); n = div(l+1,2)
    # v = complex([x[1:n];x[n-1:-1:2]],[0;-x[2n-2:-1:n+1];0;x[n+1:2n-2]])
    v = mod(l,2) == 1 ?
        complex([2x[1];x[3:2:end];x[end:-2:3]],[0;-x[2:2:end];x[end-1:-2:2]]) :
        complex([2x[1];x[3:2:end-1];-2x[end];x[end-1:-2:3]],[0;-x[2:2:end];x[end-2:-2:2]])
    rmul!(v,convert(T,l)/2)
     real(ifft(v))
end

# SinSpace plans for BigFloat

plan_transform(sp::SinSpace{D,R},x::AbstractVector{T}) where {T<:BigFloat,D,R} =
    TransformPlan{T,typeof(sp),false,Nothing}(sp,nothing)
plan_itransform(sp::SinSpace{D,R},x::AbstractVector{T}) where {T<:BigFloat,D,R} =
    ITransformPlan{T,typeof(sp),false,Nothing}(sp,nothing)

function *(::TransformPlan{T,SinSpace{D,R},false},x::AbstractVector{T}) where {T<:BigFloat,D,R}
    v=imag(fft([0;-x;0;reverse(x)]))[2:length(x)+1]
    rmul!(v,convert(T,1)/(length(x)+1))
    v
end

*(::ITransformPlan{T,SinSpace{D,R},false},x::AbstractVector{T}) where {T<:BigFloat,D,R} =
    imag(fft([0;-x;0;reverse(x)]))[2:length(x)+1]/2

# Fourier space & SinSpace plans for Complex{BigFloat}

for SP in (:Fourier,:SinSpace), (pl,TransPlan) in ((:plan_transform,:TransformPlan),
                                                    (:plan_itransform,:ITransformPlan))
    @eval begin
        $pl(sp::$SP{D,R},x::AbstractVector{T}) where {T<:Complex{BigFloat},D,R} =
                $TransPlan(sp,$pl(sp,Array{T}(undef, length(x))),Val{false})
        *(P::$TransPlan{T,$SP{D,R},false},x::Vector{T}) where {T<:Complex{BigFloat},D,R} =
            complex(P.plan*real(x),P.plan*imag(x))
    end
end


text:
# Covariance Matrix Adaptation Evolution Strategy
# ==================
#
# Implementation: (μ/μ_I,λ)-CMA-ES
#
# μ is the number of parents
# λ is the number of offspring.
#
@with_kw struct CMAES <: Optimizer
    N::Int
    initPopulation::Individual = ones(N)
    initStrategy::Strategy = strategy(τ = sqrt(N), τ_c = N^2, τ_σ = sqrt(N))
    μ::Integer = 1
    λ::Integer = 1
end

function optimize(objfun::Function, opt::CMAES;
                    iterations::Integer = 1_000,
                    tol::Float64 = 1e-10,
                    verbose = false)
    @unpack N,initPopulation,initStrategy,μ,λ = opt
    @assert μ < λ "Offspring population must be larger then parent population"

    # Initialize parent population
    individual = getIndividual(initPopulation, N)
    population = fill(individual, μ)
    offspring = Array{typeof(individual)}(undef, λ)
    fitpop = fill(Inf, μ)
    fitoff = fill(Inf, λ)

    parent = copy(individual)
    C = Diagonal{Float64}(I, N)
    s = zeros(N)
    s_σ = zeros(N)
    σ = 1.0
    E = zeros(N, λ)
    W = zeros(N, λ)

    # Generation cycle
    count = 0
    while true
        SqrtC = (C + C') / 2.0
        try
            SqrtC = cholesky(SqrtC).U
        catch
            println("Break on Cholesky")
            break
        end

        for i in 1:λ
            # offspring are generated by transforming standard normally distributed random vectors using a transformation matrix
            E[:,i] = randn(N)
            W[:,i] = σ * (SqrtC * E[:,i])
            offspring[i] = parent + W[:,i]   # (L1)
            fitoff[i] = objfun(offspring[i]) # Evaluate fitness
        end

        # Select new parent population
        idx = sortperm(fitoff)[1:μ]
        for i in 1:μ
            population[i] = offspring[idx[i]]
            fitpop[i] = fitoff[idx[i]]
        end

        w = vec(mean(W[:,idx], dims=2))
        ɛ = vec(mean(E[:,idx], dims=2))
        parent += w            #  forming recombinant perent for next generation (L2)
        s = (1.0 - 1.0/initStrategy[:τ])*s +
            (sqrt(μ/initStrategy[:τ] * (2.0 - 1.0/initStrategy[:τ]))/σ)*w      # (L3)
        C = (1.0 - 1.0/initStrategy[:τ_c]).*C + (s./initStrategy[:τ_c])*s'     # (L4)
        s_σ = (1.0 - 1.0/initStrategy[:τ_σ])*s_σ +
            sqrt(μ/initStrategy[:τ_σ]*(2.0 - 1.0/initStrategy[:τ_σ]))*ɛ        # (L5)
        σ = σ*exp(((s_σ'*s_σ)[1] - N)/(2*N*sqrt(N)))                           # (L6)

        # termination condition
        count += 1
        if count == iterations || σ <= tol
            break
        end
        verbose && println("BEST: $(fitpop[1]): $(σ)")
    end

    return population[1], fitpop[1], count
end


text:
using FWI
using jInv.Mesh;
using jInv.LinearSolvers
using jInv.InverseSolve
using jInv.Utils
using EikonalInv
using ForwardHelmholtz
using Multigrid

include("getAnalyticalMediums.jl");

function checkDerivativeFWI(f::Function,x0,param,Iact;out::Bool=true,tol::Number=1.9,nSuccess::Int=3)
	if out
		println(@sprintf("%9s\t%9s\t%9s\t%9s\t%9s\t%5s","h","E0","E1","O1","O2","OK?"))
	end
	v = 0.1*randn(size(x0))
	v = Iact*Iact'*v;
	# if eltype(x0)<:Complex
		# v += 1im*randn(size(v))
	# end
	f0,dvf = f(x0,v,param)

	# dvf   = real(dvf)
	Error = zeros(10,2)
	Order = zeros(10,2)
	Success = zeros(10)
	for j=1:5
		ft = f(x0+10.0^(-j)*v,zeros(size(x0)),param)[1]            # function value
		Error[j,1] = norm(f0-ft)/norm(f0)           # Error TaylorPoly 0
		Error[j,2] = norm(f0 .+10.0^(-j)*dvf .- ft)/norm(f0) # Error TaylorPoly 1
		if j>1
			Order[j,:] = log10.(Error[j-1,:]./Error[j,:]);
		end
		if (Order[j,2]>tol) || (Error[j,1]/Error[j,2] > 100); Success[j]=1; end
		if out
			println(@sprintf("%1.3e\t%1.3e\t%1.3e\t%1.3e\t%1.3e\t%5d",
							10.0^(-j), Error[j,1],Error[j,2], Order[j,1],Order[j,2],Success[j]))
		end
	end
	pass = sum(Success) > nSuccess
	return  pass,Error,Order
end

println("########### 2D TEST SETTING ##############################");
n = [65,129]
Omega = [0.0,4.0,0.0,8.0];
Minv = getRegularMesh(Omega,n-1);
println(Minv.n)
m = getAnalyticalConstGradInv2D(Minv.n+1,Minv.h)[1];
n = [65,129]
Omega = [0.0,8.0,0.0,4.0];
Minv = getRegularMesh(Omega,n-1);
m = m';

pad = 16;
jump = 15;
offset = 100;

# omega = [0.5./maximum(Minv.h), 0.2./maximum(Minv.h)];
omega = [0.5./maximum(Minv.h)];

WaveletCoef = rand(Complex128,size(omega));
println(omega)

RCVfile = string("rcvMap.dat");
SRCfile = string("srcMap.dat");
writeSrcRcvLocFile(SRCfile,Minv,pad,jump);
writeSrcRcvLocFile(RCVfile,Minv,pad,2);
srcNodeMap = readSrcRcvLocationFile(SRCfile,Minv);
rcvNodeMap = readSrcRcvLocationFile(RCVfile,Minv);
Q = generateSrcRcvProjOperators(Minv.n+1,srcNodeMap);
P = generateSrcRcvProjOperators(Minv.n+1,rcvNodeMap);

ABLPad = pad+5;
gamma = getABL(Minv,true,ones(Int64,Minv.dim)*ABLPad,omega[1]);
attenuation = 0.01*omega[1];
gamma += attenuation; # adding Attenuation.

# Ainv = getMUMPSsolver([],0,0,0);
# Ainv = getJuliaSolver();

###################################################################################################
levels      = 2;
numCores 	= 4;
BLAS.set_num_threads(numCores);
maxIter     = 50;
relativeTol = 1e-4;
relaxType   = "SPAI";
relaxParam  = 1.0;
relaxPre 	= 2;
relaxPost   = 2;
cycleType   ='W';
coarseSolveType = "NoMUMPS";
MG = getMGparam(levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType,0.0,0.0);
shift = 0.15;
Hparam = HelmholtzParam(Minv,zeros(0),zeros(0),0.0,true,true);
Ainv = getShiftedLaplacianMultigridSolver(Hparam, MG,shift);

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

println("Number of Sources: ",size(Q,2));

(pFor,contDiv,SourcesSubInd) = getFWIparam(omega,WaveletCoef,vec(gamma),Q,P,Minv,Ainv,workers(),16);

# setup active cells - this is important so that Sommerfeld BC will not interfere with the sensitivities.

N = prod(Minv.n+1);
Iact = speye(N);
mask = zeros(N);
mask[gamma[:] .>= 0.2*maximum(gamma)] = 1;
Iact = Iact[:,mask .== 0.0];

function f(m,v=zeros(size(m)),pFor = [])
	(D,pFor) = getData(m[:],pFor);
	D = fetch(D[1]);
	pf = fetch(pFor[1]);
	Jv = 0;
	if vecnorm(v)>0.0
		Jv = getSensMatVec(v[:],m[:],pf);
	end
	return D[:],Jv;
end
checkDerivativeFWI(f,m[:],pFor,Iact);

pf = fetch(pFor[1]);
v = rand(N);
Jv = getSensMatVec(v,m[:],pf);
Z = rand(Complex128,size(Jv));
JtZ = getSensTMatVec(Z,m[:],pf);

println(dot(Jv[:],Z[:]))
println(dot(JtZ[:],v))

if real(dot(JtZ[:],v) - dot(Jv[:],Z[:])) > 1e-5
	warn("Checking that (J^T)^T = J for FWI falied: ",real(dot(JtZ[:],v) - dot(Jv[:],Z[:])) );
end

############################## STOCHASTIC #####################################################
# println("####################### STOCHASTIC SELECTION ################################")
# sourcesSelection = randperm(size(Q,2))[1:div(size(Q,2),2)];
# setSourceSelection(pFor[1], sourcesSelection)

# function f(m,v=zeros(size(m)),pFor = [])
	# (D,pFor) = getData(m[:],pFor,Mesh2Mesh);
	# D = fetch(D[1]);
	# pf = fetch(pFor[1]);
	# Jv = 0;
	# if vecnorm(v)>0.0
		# Jv = getSensMatVec(v[:],m[:],pf);
	# end
	# return D[:],Jv;
# end
# checkDerivativeFWI(f,m[:],pFor,Iact);

# pf = fetch(pFor[1]);
# v = rand(N);
# Jv = getSensMatVec(v,m[:],pf);
# Z = rand(Complex128,size(Jv));
# JtZ = getSensTMatVec(Z,m[:],pf);

# println(dot(Jv[:],Z[:]))
# println(dot(JtZ[:],v))

# if real(dot(JtZ[:],v) - dot(Jv[:],Z[:])) > 1e-8
	# warn("Checking that (J^T)^T = J for FWI falied: ",real(dot(JtZ[:],v) - dot(Jv[:],Z[:])) );
# end

rm("srcMap.dat");
rm("rcvMap.dat");


text:
using SparseArrays

abstract type constraint end

# Constraint Types and settings

include("affineConstraints.jl")
include("socpConstraints.jl")
include("socpNewConstraints.jl")
include("socpAngleConstraints.jl")

#= Augmented Lagrangian Constraints
These are presented in the order:
- AffineEquality
- AffineInequality
- (Old) SOCP Constraints
- (Old) Functional SOCP Constraints with Slack Variables
- Simple SOCP Constraints
- Simple SOCP Constraints for N-constraints
- SOCP Constraints for Angles
- SOCP Constraints for Angles for N-constraints
=#

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
module Core
  include("token.jl")
  include("pos_tag.jl")
  include("sentence.jl")
  include("sequence.jl")
  include("cache.jl")
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
"""
    smc!(model::SMCModel, smcio::SMCIO)
Run the SMC algorithm for the given model and input/output arguments.
"""
function smc!(model::SMCModel, smcio::SMCIO{Particle}) where Particle
  zetas = smcio.zetas
  zetaAncs = smcio.internal.zetaAncs
  lws = smcio.internal.lws
  ws = smcio.ws
  as = smcio.internal.as
  engine = getRNG()
  pScratch = smcio.internal.particleScratch

  logZhats = smcio.logZhats
  lZ = 0.0

  for p = 1:smcio.n
    p > 1 && _copyParticles!(zetaAncs, zetas, as)

    _mutateParticles!(zetas, engine, p, model.M!, zetaAncs, pScratch)
    _logWeightParticles!(lws, p, model.lG, zetas, pScratch)

    maxlw::Float64 = maximum(lws)
    ws .= exp.(lws .- maxlw)
    mws::Float64 = mean(ws)
    lZ += maxlw + log(mws)
    ws ./= mws
    @inbounds logZhats[p] = lZ

    p < smcio.n && _resample!(smcio)

    _intermediateOutput!(smcio, p)
  end
end

function _csmc!(model::SMCModel, smcio::SMCIO{Particle},
  ref::Vector{Particle}, refout::Vector{Particle}) where Particle
  zetas = smcio.zetas
  zetaAncs = smcio.internal.zetaAncs
  lws = smcio.internal.lws
  ws = smcio.ws
  as = smcio.internal.as
  engine = getRNG()
  pScratch = smcio.internal.particleScratch

  logZhats = smcio.logZhats
  lZ = 0.0

  for p = 1:smcio.n
    p > 1 && _copyParticles!(zetaAncs, zetas, as)

    _mutateParticles!(zetas, engine, p, model.M!, zetaAncs, pScratch, ref[p])
    _logWeightParticles!(lws, p, model.lG, zetas, pScratch)

    maxlw::Float64 = maximum(lws)
    ws .= exp.(lws .- maxlw)
    mws::Float64 = mean(ws)
    lZ += maxlw + log(mws)
    ws ./= mws
    @inbounds logZhats[p] = lZ

    p < smcio.n && _cresample!(smcio)

    _intermediateOutput!(smcio, p)
  end
end

"""
    csmc!(model::SMCModel, smcio::SMCIO, ref::Vector{Particle}, refout::Vector{Particle})

Run the conditional SMC algorithm for the given model, input/output arguments,
reference path and output path.

It is permitted for ref and refout to be the same.
"""
function csmc!(model::SMCModel, smcio::SMCIO{Particle},
  ref::Vector{Particle}, refout::Vector{Particle}) where Particle

  _csmc!(model, smcio, ref, refout)

  pickParticle!(refout, smcio)
end

"""
    csmc!(model::SMCModel, lM::F, smcio::SMCIO, ref::Vector{Particle}, refout::Vector{Particle})

Run the conditional SMC algorithm for the given model, input/output arguments,
reference path and output path. Backward samping is used to choose the new
reference path, which requires the provision of the lM function.

It is permitted for ref and refout to be the same.
"""
function csmc!(model::SMCModel, lM::F, smcio::SMCIO{Particle},
  ref::Vector{Particle}, refout::Vector{Particle}) where {F<:Function, Particle}

  _csmc!(model, smcio, ref, refout)

  pickParticleBS!(refout, smcio, lM)
end


text:
# Gaver-Stehfest method

# 1.
# Stehfest, Harald.
# Algorithm 368: Numerical inversion of Laplace transforms [D5].
# Communications of the ACM 13.1 (1970): 47-49.

# 2.
# Kuhlman, Kristopher L.
# Review of inverse Laplace transform algorithms for Laplace-space numerical approaches.
# Numerical Algorithms 63.2 (2013): 339-355.

"""
    stehfest_coeffs(N::Integer)

Computes the Stehfest coefficients (summation weights) Vk for N terms. Coefficients are computed using BigFloats.
Only depends on the approximation order (N) and the precision. N must be even but can be precomputed and used only once.
Coefficients get very large in magnitude and oscillate rapidly. Must be careful to avoid catastrophic cancellation for lower precision.
"""
function stehfest_coeffs(N::Integer)
    if isodd(N)
        N += 1
        @warn "N must be even... increment N += 1"
    end
    v = zeros(BigFloat, N)
    aux = zero(eltype(v))
    for k in 1:N
        for j in fld(k + 1, 2):minimum(Int, [k, div(N, 2)])
            aux = big(j)^(div(N, 2)) * factorial(big(2 * j))
            aux /= factorial(big(div(N, 2) - j)) * factorial(big(j)) * factorial(big(j - 1))
            aux /= factorial(big(k - j)) * factorial(big(2 * j - k))
            v[k] += aux
        end
        v[k] *= (-1)^(k + div(N, 2))
    end
    return v
end

"""
    gaverstehfest(func::Function, t::AbstractFloat, v = stehfest_coeffs(N))

Evaluate the inverse Laplace transform of `func` at the point `t` using the Gaver-Stehfest algorithm.
`v` are the stehfest coefficients computed with `stehfest_coeffs` that only depend on number of terms.
N (which must be even) defaults to 36 should depend on the precision used and desired accuracy.
The precision of the stehfest coefficients is used in the computation.

In double precision (Float64), N should be <= 18 to provide best accuracy.
Usually only moderate accuracy can be achieved in double precision ~rtol ≈ 1e-4.
For higher accuracy, BigFloats should be used with N = 36.

Increasing precision should be accompanied by an increase in the number of coefficients used.
Increasing precision without increasing number of coefficients will not yield better accuracy. The inverse is generally true as well.

This method is not robust to oscillating F(t) and must be smooth.

# Examples
```
julia> F(s) = 1 / (s + 1) # where analytical inverse is f(t) = exp(-t)

julia> InverseLaplace.gaverstehfest(F, 2.0) # computes with default 36 coefficients
0.1353352832366128315426471959891170863784520272858265151734040491181311503059561

julia> InverseLaplace.gaverstehfest(F, 2.0, v = stehfest_coeffs(20)) # computes with custom number
0.1353353114073885136885007645878189977740624364169818512599342310325432753817199

# to calculate in double precision convert stehfest coefficients to Float64
julia> InverseLaplace.gaverstehfest(F, 2.0, v = convert(Vector{Float64}, stehfest_coeffs(18)))
0.13533650985980258
```
"""
function gaverstehfest(func::Function, t::AbstractFloat; v = stehfest_coeffs(36))
    N = length(v)
    a = zero(eltype(v))
    for k in 1:N
        bk = convert(eltype(v), k)
        a += v[k] * func(bk * log(convert(eltype(a), 2)) / t)
    end
    return a * log(convert(eltype(a), 2)) / t
end

"""
    gaverstehfest(func::Function, t::AbstractArray, v = stehfest_coeffs(N))

Evaluate the inverse Laplace transform of `func` over an array of `t` using the Gaver-Stehfest algorithm.
Computes coefficients once and calculates f(t) across available threads.

# Example
```
julia> gaverstehfest(s -> 1/(s + 1), 2.0:4.0)
0.1353355048178631463198819259249043857320738467297190227379428405365703748720082
0.04978728177016550841951683309410878070827215175940039571536203126773133604403101
0.01831383641619355549133471411401755204406266755309342902744499924574072011058623
```
"""
function gaverstehfest(f::Function, t::AbstractArray; v = stehfest_coeffs(18))
    N = length(v)
    a = zeros(eltype(v), length(t))
    Threads.@threads for ind in eachindex(t)
        a[ind] = gaverstehfest(f, t[ind], v = v)
    end
    return a
end


text:
using PyCall
export donthcomp

"""
This function was adapted by ADT from the subroutine donthcomp in
donthcomp.f, distributed with XSpec.
Nthcomp documentation:
https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XSmodelNthcomp.html
Refs:
Zdziarski, Johnson & Magdziarz 1996, MNRAS, 283, 193,
as extended by Zycki, Done & Smith 1999, MNRAS 309, 561
Note that the subroutine has been modified so that parameter 4
is ignored, and the seed spectrum is always a blackbody.
ear: Energy vector, listing "Energy At Right" of bins (keV)
param: list of parameters; see the 5 parameters listed below.
The original fortran documentation for this subroutine is included below:
Driver for the Comptonization code solving Kompaneets equation
seed photons  -  (disk) blackbody
reflection + Fe line with smearing

Model parameters:
1: photon spectral index
2: plasma temperature in keV
3: (disk)blackbody temperature in keV
4: type of seed spectrum (0 - blackbody, 1 - diskbb)
5: redshift
"""
function donthcomp(ear, param)
    ne = length(ear)

    zfactor = 1 + param[5]

    xth, nth, spt = _thcompton(param[3] / 511, param[2] / 511, param[1])

    xninv = 511 / zfactor
    ih = 2
    xx = 1 / xninv
    while (ih < nth) && (xx > xth[ih])
        ih = ih + 1
    end
    il = ih - 1
    spp = spt[il] + (spt[ih] - spt[il]) * (xx - xth[il]) / (xth[ih] - xth[il])
    normfac = 1 / spp

    photar = zeros(ne)
    prim = zeros(ne)
    j = 1
    for i = 1:ne
        while (j <= nth) && (511.0 * xth[j] < ear[i] * zfactor)
            j = j + 1
        end
        if (j <= nth)
            if (j > 1)
                jl = j - 1
                prim[i] =
                    spt[jl] + (
                        (ear[i] / 511.0 * zfactor - xth[jl]) * (spt[jl + 1] - spt[jl]) /
                        (xth[jl + 1] - xth[jl])
                    )
            else
                prim[i] = spt[1]
            end
        end
    end
    for i = 2:ne
        photar[i] = (
            0.5 *
            (prim[i] / ear[i]^2 + prim[i - 1] / ear[i - 1]^2) *
            (ear[i] - ear[i - 1]) *
            normfac
        )
    end
    return photar
end

function _thcompton(tempbb, theta, gamma)
    """
    This function was adapted by ADT from the subroutine thcompton in
    donthcomp.f, distributed with XSpec.
    Nthcomp documentation:
    https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XSmodelNthcomp.html
    Refs:
    Zdziarski, Johnson & Magdziarz 1996, MNRAS, 283, 193,
    as extended by Zycki, Done & Smith 1999, MNRAS 309, 561
    The original fortran documentation for this subroutine is included below:
    Thermal Comptonization; solves Kompaneets eq. with some
    relativistic corrections. See Lightman Zdziarski (1987), ApJ
    The seed spectrum is a blackbody.
    version: January 96
    #c  input parameters:
    #real * 8 tempbb,theta,gamma
    """
    #c use internally Thomson optical depth
    tautom = sqrt(2.250 + 3.0 / (theta * ((gamma + 0.50)^2 - 2.250))) - 1.50

    # Initialise arrays
    dphdot = zeros(900)
    rel = zeros(900)
    c2 = zeros(900)
    sptot = zeros(900)
    bet = zeros(900)
    x = zeros(900)

    #c JMAX  -  # OF PHOTON ENERGIES
    #c delta is the 10 - log interval of the photon array.
    delta = 0.02
    deltal = delta * log(10.0)
    xmin = 1e-4 * tempbb
    xmax = 40.0 * theta
    jmax = min(899, Int(round(log10(xmax / xmin) / delta)) + 1)

    #c X  -  ARRAY FOR PHOTON ENERGIES
    # Energy array is normalized by 511 keV, the rest energy of an electron
    for j=1:jmax+1
        x[j] = xmin * 10 ^ ((j-1) * delta)
    end

    #c compute c2(x), and rel(x) arrays
    #c c2(x) is the relativistic correction to Kompaneets equation
    #c rel(x) is the Klein - Nishina cross section divided by the
    #c Thomson crossection
    for j = 1:jmax
        w = x[j]
        #c c2 is the Cooper's coefficient calculated at w1
        #c w1 is x(j + 1 / 2) (x(i) defined up to jmax + 1)
        w1 = sqrt(x[j] * x[j + 1])
        c2[j] = (w1^4 / (1.0 + 4.60 * w1 + 1.1 * w1 * w1))
        if (w <= 0.05)
            #c use asymptotic limit for rel(x) for x less than 0.05
            rel[j] = (1.0 - 2.0 * w + 26.0 * w * w * 0.2)
        else
            z1 = (1.0 + w) / w^3
            z2 = 1.0 + 2.0 * w
            z3 = log(z2)
            z4 = 2.0 * w * (1.0 + w) / z2
            z5 = z3 / 2.0 / w
            z6 = (1.0 + 3.0 * w) / z2 / z2
            rel[j] = (0.75 * (z1 * (z4 - z3) + z5 - z6))
        end
    end

    #c the thermal emission spectrum
    jmaxth = min(900, Int(round(log10(50 * tempbb / xmin) / delta)))
    if jmaxth > jmax
        jmaxth = jmax
    end
    planck = 15.0 / (π * tempbb)^4
    for j=1:jmaxth
        dphdot[j] = planck * x[j]^2 / (exp(x[j] / tempbb)-1)
    end

    #c compute beta array, the probability of escape per Thomson time.
    #c bet evaluated for spherical geometry and nearly uniform sources.
    #c Between x = 0.1 and 1.0, a function flz modifies beta to allow
    #c the increasingly large energy change per scattering to gradually
    #c eliminate spatial diffusion
    jnr = Int(round(log10(0.10 / xmin) / delta + 1))
    jnr = min(jnr, jmax - 1)
    jrel = Int(round(log10(1 / xmin) / delta + 1))
    jrel = min(jrel, jmax)
    xnr = x[jnr]
    xr = x[jrel]
    for j = 1:(jnr - 1)
        taukn = tautom * rel[j]
        bet[j] = 1.0 / tautom / (1.0 + taukn / 3.0)
    end
    for j = jnr:jrel
        taukn = tautom * rel[j]
        arg = (x[j] - xnr) / (xr - xnr)
        flz = 1 - arg
        bet[j] = 1.0 / tautom / (1.0 + taukn / 3.0 * flz)
    end
    for j = jrel+1:jmax
        bet[j] = 1.0 / tautom
    end

    dphesc = _thermlc(tautom, theta, deltal, x, jmax, dphdot, bet, c2)

    #c     the spectrum in E F_E
    for j = 1:(jmax - 1)
        sptot[j] = dphesc[j] * x[j]^2
    end

    return x, jmax, sptot
end

function _thermlc(tautom, theta, deltal, x, jmax, dphdot, bet, c2)
    """
    This function was adapted by ADT from the subroutine thermlc in
    donthcomp.f, distributed with XSpec.
    Nthcomp documentation:
    https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XSmodelNthcomp.html
    Refs:
    Zdziarski, Johnson & Magdziarz 1996, MNRAS, 283, 193,
    as extended by Zycki, Done & Smith 1999, MNRAS 309, 561
    The original fortran documentation for this subroutine is included below:
    This program computes the effects of Comptonization by
    nonrelativistic thermal electrons in a sphere including escape, and
    relativistic corrections up to photon energies of 1 MeV.
    the dimensionless photon energy is x = hv / (m * c * c)
    The input parameters and functions are:
    dphdot(x), the photon production rate
    tautom, the Thomson scattering depth
    theta, the temperature in units of m*c*c
    c2(x), and bet(x), the coefficients in the K - equation and the
      probability of photon escape per Thomson time, respectively,
      including Klein - Nishina corrections
    The output parameters and functions are:
    dphesc(x), the escaping photon density
    """
    dphesc = zeros(900)  # Initialise the output
    a = zeros(900)
    b = zeros(900)
    c = zeros(900)
    d = zeros(900)
    alp = zeros(900)
    u = zeros(900)
    g = zeros(900)
    gam = zeros(900)

    #c u(x) is the dimensionless photon occupation number
    c20 = tautom / deltal

    #c determine u
    #c define coefficients going into equation
    #c a(j) * u(j + 1) + b(j) * u(j) + c(j) * u(j - 1) = d(j)
    for j = 2:(jmax - 1)
        w1 = sqrt(x[j] * x[j + 1])
        w2 = sqrt(x[j - 1] * x[j])
        #c  w1 is x(j + 1 / 2)
        #c  w2 is x(j - 1 / 2)
        a[j] = -c20 * c2[j] * (theta / deltal / w1 + 0.5)
        t1 = -c20 * c2[j] * (0.5 - theta / deltal / w1)
        t2 = c20 * c2[j - 1] * (theta / deltal / w2 + 0.5)
        t3 = x[j]^3 * (tautom * bet[j])
        b[j] = t1 + t2 + t3
        c[j] = c20 * c2[j - 1] * (0.5 - theta / deltal / w2)
        d[j] = x[j] * dphdot[j]
    end

    #c define constants going into boundary terms
    #c u(1) = aa * u(2) (zero flux at lowest energy)
    #c u(jx2) given from region 2 above
    x32 = sqrt(x[1] * x[2])
    aa = (theta / deltal / x32 + 0.5) / (theta / deltal / x32 - 0.5)

    #c zero flux at the highest energy
    u[jmax] = 0.0

    #c invert tridiagonal matrix
    alp[2] = b[2] + c[2] * aa
    gam[2] = a[2] / alp[2]
    for j = 3:(jmax - 1)
        alp[j] = b[j] - c[j] * gam[j - 1]
        gam[j] = a[j] / alp[j]
    end
    g[2] = d[2] / alp[2]
    for j = 3:(jmax - 2)
        g[j] = (d[j] - c[j] * g[j - 1]) / alp[j]
    end
    g[jmax - 1] =
        (d[jmax - 1] - a[jmax - 1] * u[jmax] - c[jmax-1] * g[jmax - 2]) /
        alp[jmax - 1]
    u[jmax - 1] = g[jmax - 1]
    for j = 3:jmax-1
        jj = jmax + 1 - j
        u[jj] = g[jj] - gam[jj] * u[jj + 1]
    end
    u[1] = aa * u[2]
    #c compute new value of dph(x) and new value of dphesc(x)
    dphesc[1:jmax] = @. x[1:jmax] * x[1:jmax] * u[1:jmax] * bet[1:jmax] * tautom

    return dphesc
end


text:
module DIO
using MLStyle
import Libdl
import Parameters
const DEBUG = false
export @exportapi, @apisetup, PyOType

include("juliainfo.jl")
include("python.jl")
include("utils.jl")
include("static.jl")

const Setup_AutoBound_Symbols = Symbol[]
const PyAPI_FieldSymbols = Symbol[]
const PyAPI_FieldTypes = Symbol[]
const PyAPI_Fields = Any[]
const PyAPI_Construct = Expr[]

"""
Automatically bind the loaded Python APIs(`PyAPI_Struct`) to the first argument,
and enter a function `MyFunc(...) = <callee module>.MyFunc(PyAPI_Struct, ...)`
to the caller module:
    ```
    @apisetup begin
        PyObject_Call = PySym(:PyObject_Call)
    end
    @exportapi MyFunc
    function MyFunc(apis, ...)
        ccall(apis.PyObject_Call, ...)
    end
    ```
"""
macro exportapi(sym::Symbol)
    _exportapi(sym)
end
function _exportapi(sym::Symbol)
    push!(Setup_AutoBound_Symbols, sym)
end
macro apisetup(define::Expr)
    _apisetup(define)
end

function _apisetup(define::Expr)
    Meta.isexpr(define, :block) ||
        error("malformed use of @apisetup require a block of assignments but got $(define).")

    for each in define.args
        @switch each begin
        @case ::LineNumberNode
            nothing
        @case :($(a :: Symbol) = $expr)
            if a in PyAPI_FieldSymbols
                error("duplicate symbol reference from Python DLL: $a")
            end
            t = Symbol("_Type_Of_$a")
            push!(PyAPI_FieldTypes, t)
            push!(PyAPI_FieldSymbols, a)
            push!(PyAPI_Fields, :($a :: $t))
            push!(PyAPI_Construct, :($a = $expr))
            nothing
        @case _
            # TODO
            error("malformed use of @apisetup: expect assignments but got: $each")
        end
    end
end

Parameters.@with_kw struct PyOType
    PY_VERSION :: Tuple{Int, Int, Int, String, Int}
    builtins :: PyPtr
    print :: PyPtr
    bool :: PyPtr
    int :: PyPtr
    float :: PyPtr
    str :: PyPtr
    type :: PyPtr
    None :: PyPtr
    True :: PyPtr
    False :: PyPtr
    complex :: PyPtr
    tuple :: PyPtr
    list :: PyPtr
    set :: PyPtr
    dict :: PyPtr
    import_module :: PyPtr
end

const IsInitialized = Ref(false)
const LibdlFlags = Libdl.RTLD_GLOBAL
_initialize() = IsInitialized[] = true
macro setup(libpy_addr :: Addr)
    if IsInitialized[]
        Base.@warn "Re-setup is not allowed. Restart your runtime and import this package again."
    else
        py_dll_struct_cons_expr =
            @q let $(PyAPI_Construct...)
                ($(PyAPI_FieldSymbols...), )
            end

        unhygienic_part1 = @q begin
            $__source__
            const PyDLL = reinterpret(Ptr{Nothing}, $libpy_addr)
            PySym(sym::Symbol) = $Libdl.dlsym(PyDLL, sym)
            function PySymMaybe(sym::Symbol)
                p = $Libdl.dlsym_e(PyDLL, sym)
                p == C_NULL ? nothing : p
            end
        end

        # compute all required symbols from libpython.so/dll, and other related values
        hygienic_part = :($(esc(:PyAPI_Args)) = $py_dll_struct_cons_expr)

        unhygienic_part2 = @q begin
            # compute types of all required symbols and related values
            $(Expr(:tuple, PyAPI_FieldTypes...)) =  map(typeof, PyAPI_Args)

            # declare the type of 'PyAPI_Struct'
            struct PyAPI_Type
                PyO :: PyOType
                $(PyAPI_Fields...)
            end

            # construct 'PyAPI_Struct'
            const PyAPI_Struct = PyAPI_Type(PyO, PyAPI_Args...)
            PyAPI_Args = nothing

            $( [:($each(args...) = $(@__MODULE__).$each(PyAPI_Struct, args...))
                for each in Setup_AutoBound_Symbols]...)
            # set up exception code for API functions
            $( [:($(@__MODULE__).DIO_ExceptCode(::typeof($each)) = DIO_ExceptCode($(@__MODULE__).$each))
                for each in Setup_AutoBound_Symbols]...)
            $_initialize()
        end

        r = @q begin
            # 'PyO' is computed when starting Julia from Python
            PyO = $(esc(:PyO))
            $(esc(unhygienic_part1))
            PySym = $(esc(:PySym))
            PySymMaybe = $(esc(:PySymMaybe))
            $hygienic_part
            $(esc(unhygienic_part2))
        end
        DEBUG && @info r
        return r
    end
end
include("macros.jl")
include("symbols.jl")
include("support.jl")
include("dynamic.jl")

end


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

using Markdown
using InteractiveUtils

# ╔═╡ bc0083ae-2c24-11ec-27d1-d3151362feba
using VegaLite, DataFrames, CSV, RDatasets

# ╔═╡ af8af33a-f75c-4ac1-9c3e-7afa932fd4c1
md"## Plotting with VegaLite
Vegalite is a modern grammar-of-graphics programming language. See information and gallery here https://vega.github.io/vega-lite/
"

# ╔═╡ 5b856367-3bf0-4f0c-b82f-6ae3de31e3b5
md"Load the relevant libraries, and get the iris dataset"

# ╔═╡ 9742e502-f9df-4576-8783-a3b81c40a685

# ╔═╡ 2c955350-cd80-4429-8215-a13ada4b9efd
iris = dataset("datasets", "iris")

# ╔═╡ 08a306fe-684b-42ae-8586-a0e61fc9d889
md"The most basic plots involve a DataFrame that you `pipe` (using `|>`) to the `@vlplot` macro. You choose a `mark` (or `geom`) and define `x` and `y`."

# ╔═╡ e563a075-5fd9-4541-99ee-07a84080287b
iris |> @vlplot(:point, x = :SepalLength, y = :SepalWidth)

# ╔═╡ 16b538b2-1d30-4b91-988f-59448783545e
md"By default, scales start at 0. To control the scales for x and y, pass a specification to x and y"

# ╔═╡ 4a59b6e3-1619-485e-a862-a8ab2c863e5a
iris |> @vlplot(:point,
    x = {:SepalLength, scale={zero=false}},
    y = {:SepalWidth, scale={zero=false}}
)

# ╔═╡ 0a84be01-7907-4482-9f35-f92e81c7ef5d
md"Everything inside the `@vlplot`  first bracket counts as a `mapping` or `aesthetic`. Specify `color` to have different colors for different species"

# ╔═╡ 506b7d99-cce0-4f0d-afc4-4e151038e0d8
iris |> @vlplot(:point,
    x = {:SepalLength, scale={zero=false}},
    y = {:SepalWidth, scale={zero=false}},
    color = :Species
)

# ╔═╡ f97a6cdd-a871-42f3-9ba6-ac8be800fbb2
md"You can stack `@vlplot`s on top of each other with`+`. Unspecified arguments will get passed over from the previous. Here to add both a line and points"

# ╔═╡ 6d9354ee-ed0a-49af-9c40-835c1f20d35e
iris |> @vlplot(
    x = {:SepalLength, scale={zero=false}},
    y = {:SepalWidth, scale={zero=false}},
    color = :Species) +
    @vlplot(:point) + @vlplot(:line)

# ╔═╡ a5bb11a0-b215-487e-ae9d-b72ef4cd6465
md"In practice what you probably want is a trendline. To add something to the plot that is a derived function of the data, use a `transform`. Note that `groupby` in the regression transform requires a vector"

# ╔═╡ 2a2bed9a-7f97-4975-aea8-62e17bf37923
iris |> @vlplot(
    x = {:SepalLength, scale={zero=false}},
    y = {:SepalWidth, scale={zero=false}},
    color = :Species) +
    @vlplot(:point) +
	@vlplot(:line,
		transform = [{
			regression=:SepalWidth,
			on = :SepalLength,
			groupby = [:Species]
		}])

# ╔═╡ 66073b7c-e972-451a-ab77-63dcb3dbe7f1
md"Add a `count()` transformation function to use a bar plot to plot  a histogram"

# ╔═╡ c08554d6-c6d9-408e-bc96-1fa9eaedb4d3
iris |> @vlplot(:bar,
    x={
        :SepalLength,
        bin={maxbins=10},
    },
    y="count()")

# ╔═╡ 983c4ab3-969a-4421-b9ce-6b6cec0e8438
md"You can easily create faceted plots by mapping a variable to the `row` or `col` argument"

# ╔═╡ 7e81c820-c212-49a1-a21f-46161cb4716c
iris |> @vlplot(:bar,
    x={
        :SepalLength,
        bin={maxbins=10},
    },
    y="count()",
    row = :Species
)

# ╔═╡ c06bcfcd-34f9-4254-93f1-206271982843
md"Let's load a different dataset - a version of the sleep dataset"

# ╔═╡ d508a53d-015c-483c-aae6-2e65c858d18e
allison = CSV.read("data/allison.csv", DataFrame, missingstring = "NA")

# ╔═╡ b1fa7bf2-3a96-4f44-8ca1-a07094467fdf
md"A basic plot reveals a very skewed distribution of x values"

# ╔═╡ effc71fb-261a-406f-b305-604bf8a2207a
allison |> @vlplot(:point, x = :BodyWt, y = :Gestation)

# ╔═╡ 5e76f741-56a8-4b0b-8fb8-a7a0aaf6763e
md"We can address that by applying a different `scale` to `x`"

# ╔═╡ 8cd88777-a081-40ec-83cc-63451043c71a
allison |> @vlplot(:point, x = {:BodyWt, scale = {type = :log}},
        y = :Gestation, color = :Predation)

# ╔═╡ 9731b8e1-a464-4280-82eb-ff4e3831227b
md"Apply and `:o` to the `Predation` variable to signal that we are dealing with an ordinal variable`"

# ╔═╡ b4e8e702-bd06-4d64-b05f-c2a901c453f6
allison |> @vlplot(:point, x = {:BodyWt, scale = {type = :log}},
        y = :Gestation, color = "Predation:o")

# ╔═╡ dc120eee-0308-425b-a22b-0725b0a18f53
md"""
### Exercise 1
Create different histograms for total sleep for different predation categories
"""

# ╔═╡ ad5bcb43-7806-40a4-8210-601ee67f8bb2

# ╔═╡ f73627ad-a8fb-40c2-810c-6052e1a4216a
md"""
### Exercise 2
Create one or more diagrams investigating the relationship between predation risk and gestation time
"""

# ╔═╡ bbadec0d-6242-481a-bc9c-f7ef132dbea5

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
RDatasets = "ce6b1742-4840-55fa-b093-852dadbb1d8b"
VegaLite = "112f6efa-9a02-5b7d-90c0-432ed331239a"

[compat]
CSV = "~0.8.5"
DataFrames = "~1.2.2"
RDatasets = "~0.7.5"
VegaLite = "~2.6.0"
"""

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

julia_version = "1.7.0-rc1"
manifest_format = "2.0"

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

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

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

[[deps.CSV]]
deps = ["Dates", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode"]
git-tree-sha1 = "b83aa3f513be680454437a0eee21001607e5d983"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.8.5"

[[deps.CategoricalArrays]]
deps = ["DataAPI", "Future", "Missings", "Printf", "Requires", "Statistics", "Unicode"]
git-tree-sha1 = "fbc5c413a005abdeeb50ad0e54d85d000a1ca667"
uuid = "324d7699-5711-5eae-9e2f-1d82baa6b597"
version = "0.10.1"

[[deps.CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "31d0151f5716b655421d9d75b7fa74cc4e744df2"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.39.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.ConstructionBase]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f74e9d5388b8620b4cee35d4c5a618dd4dc547f4"
uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
version = "1.3.0"

[[deps.Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "d785f42445b63fc86caa08bb9a9351008be9b765"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.2.2"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[deps.DataValues]]
deps = ["DataValueInterfaces", "Dates"]
git-tree-sha1 = "d88a19299eba280a6d062e135a43f00323ae70bf"
uuid = "e7dc6d0d-1eca-5fa6-8ad6-5aecde8b7ea5"
version = "0.4.13"

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

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

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

[[deps.ExprTools]]
git-tree-sha1 = "b7e3d17636b348f005f11040025ae8c6f645fe92"
uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
version = "0.1.6"

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

[[deps.FilePaths]]
deps = ["FilePathsBase", "MacroTools", "Reexport", "Requires"]
git-tree-sha1 = "919d9412dbf53a2e6fe74af62a73ceed0bce0629"
uuid = "8fc22ac5-c921-52a6-82fd-178b2807b824"
version = "0.8.3"

[[deps.FilePathsBase]]
deps = ["Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "7fb0eaac190a7a68a56d2407a6beff1142daf844"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.12"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[deps.Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[deps.HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "14eece7a3308b4d8be910e265c724a6ba51a9798"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.16"

[[deps.IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

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

[[deps.InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

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

[[deps.JSONSchema]]
deps = ["HTTP", "JSON", "URIs"]
git-tree-sha1 = "2f49f7f86762a0fbbeef84912265a1ae61c4ef80"
uuid = "7d188eb4-7ad8-530c-ae41-71a32a6d4692"
version = "0.3.4"

[[deps.LazyArtifacts]]
deps = ["Artifacts", "Pkg"]
uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"

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

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

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

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

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

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

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

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

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

[[deps.MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

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

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

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

[[deps.Mocking]]
deps = ["Compat", "ExprTools"]
git-tree-sha1 = "29714d0a7a8083bba8427a4fbfb00a540c681ce7"
uuid = "78c3b35d-d492-501b-9361-3d52fe80e533"
version = "0.7.3"

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

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

[[deps.NodeJS]]
deps = ["Pkg"]
git-tree-sha1 = "905224bbdd4b555c69bb964514cfa387616f0d3a"
uuid = "2bd173c7-0d6d-553b-b6af-13a54713934c"
version = "1.3.0"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "bfd7d8c7fd87f04543810d9cbd3995972236ba1b"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "1.1.2"

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

[[deps.PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "a193d6ad9c45ada72c14b731a318bedd3c2f00cf"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.3.0"

[[deps.PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "69fd065725ee69950f3f58eceb6d144ce32d627d"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.2.2"

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

[[deps.RData]]
deps = ["CategoricalArrays", "CodecZlib", "DataFrames", "Dates", "FileIO", "Requires", "TimeZones", "Unicode"]
git-tree-sha1 = "19e47a495dfb7240eb44dc6971d660f7e4244a72"
uuid = "df47a6cb-8c03-5eed-afd8-b6050d6c41da"
version = "0.8.3"

[[deps.RDatasets]]
deps = ["CSV", "CodecZlib", "DataFrames", "FileIO", "Printf", "RData", "Reexport"]
git-tree-sha1 = "06d4da8e540edb0314e88235b2e8f0429404fdb7"
uuid = "ce6b1742-4840-55fa-b093-852dadbb1d8b"
version = "0.7.5"

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

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

[[deps.RecipesBase]]
git-tree-sha1 = "44a75aa7a527910ee3d1751d1f0e4148698add9e"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.1.2"

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

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

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

[[deps.SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "54f37736d8934a12a200edea2f9206b03bdf3159"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.7"

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

[[deps.Setfield]]
deps = ["ConstructionBase", "Future", "MacroTools", "Requires"]
git-tree-sha1 = "fca29e68c5062722b5b4435594c3d1ba557072a3"
uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46"
version = "0.7.1"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

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

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

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

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.TableTraitsUtils]]
deps = ["DataValues", "IteratorInterfaceExtensions", "Missings", "TableTraits"]
git-tree-sha1 = "78fecfe140d7abb480b53a44f3f85b6aa373c293"
uuid = "382cd787-c1b6-5bf2-a167-d5b971a19bda"
version = "1.0.2"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "fed34d0e71b91734bf0a7e10eb1bb05296ddbcd0"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.0"

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

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.TimeZones]]
deps = ["Dates", "Future", "LazyArtifacts", "Mocking", "Pkg", "Printf", "RecipesBase", "Serialization", "Unicode"]
git-tree-sha1 = "a5688ffdbd849a98503c6650effe79fe89a41252"
uuid = "f269a46b-ccf7-5d73-abea-4c690281aa53"
version = "1.5.9"

[[deps.TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[deps.URIParser]]
deps = ["Unicode"]
git-tree-sha1 = "53a9f49546b8d2dd2e688d216421d050c9a31d0d"
uuid = "30578b45-9adc-5946-b283-645ec420af67"
version = "0.4.1"

[[deps.URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

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

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

[[deps.Vega]]
deps = ["DataStructures", "DataValues", "Dates", "FileIO", "FilePaths", "IteratorInterfaceExtensions", "JSON", "JSONSchema", "MacroTools", "NodeJS", "Pkg", "REPL", "Random", "Setfield", "TableTraits", "TableTraitsUtils", "URIParser"]
git-tree-sha1 = "43f83d3119a868874d18da6bca0f4b5b6aae53f7"
uuid = "239c3e63-733f-47ad-beb7-a12fde22c578"
version = "2.3.0"

[[deps.VegaLite]]
deps = ["Base64", "DataStructures", "DataValues", "Dates", "FileIO", "FilePaths", "IteratorInterfaceExtensions", "JSON", "MacroTools", "NodeJS", "Pkg", "REPL", "Random", "TableTraits", "TableTraitsUtils", "URIParser", "Vega"]
git-tree-sha1 = "3e23f28af36da21bfb4acef08b144f92ad205660"
uuid = "112f6efa-9a02-5b7d-90c0-432ed331239a"
version = "2.6.0"

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

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

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

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

# ╔═╡ Cell order:
# ╟─af8af33a-f75c-4ac1-9c3e-7afa932fd4c1
# ╟─5b856367-3bf0-4f0c-b82f-6ae3de31e3b5
# ╠═9742e502-f9df-4576-8783-a3b81c40a685
# ╠═bc0083ae-2c24-11ec-27d1-d3151362feba
# ╠═2c955350-cd80-4429-8215-a13ada4b9efd
# ╟─08a306fe-684b-42ae-8586-a0e61fc9d889
# ╠═e563a075-5fd9-4541-99ee-07a84080287b
# ╟─16b538b2-1d30-4b91-988f-59448783545e
# ╠═4a59b6e3-1619-485e-a862-a8ab2c863e5a
# ╟─0a84be01-7907-4482-9f35-f92e81c7ef5d
# ╠═506b7d99-cce0-4f0d-afc4-4e151038e0d8
# ╟─f97a6cdd-a871-42f3-9ba6-ac8be800fbb2
# ╠═6d9354ee-ed0a-49af-9c40-835c1f20d35e
# ╟─a5bb11a0-b215-487e-ae9d-b72ef4cd6465
# ╠═2a2bed9a-7f97-4975-aea8-62e17bf37923
# ╟─66073b7c-e972-451a-ab77-63dcb3dbe7f1
# ╠═c08554d6-c6d9-408e-bc96-1fa9eaedb4d3
# ╟─983c4ab3-969a-4421-b9ce-6b6cec0e8438
# ╠═7e81c820-c212-49a1-a21f-46161cb4716c
# ╟─c06bcfcd-34f9-4254-93f1-206271982843
# ╠═d508a53d-015c-483c-aae6-2e65c858d18e
# ╟─b1fa7bf2-3a96-4f44-8ca1-a07094467fdf
# ╠═effc71fb-261a-406f-b305-604bf8a2207a
# ╟─5e76f741-56a8-4b0b-8fb8-a7a0aaf6763e
# ╠═8cd88777-a081-40ec-83cc-63451043c71a
# ╟─9731b8e1-a464-4280-82eb-ff4e3831227b
# ╠═b4e8e702-bd06-4d64-b05f-c2a901c453f6
# ╟─dc120eee-0308-425b-a22b-0725b0a18f53
# ╠═ad5bcb43-7806-40a4-8210-601ee67f8bb2
# ╟─f73627ad-a8fb-40c2-810c-6052e1a4216a
# ╠═bbadec0d-6242-481a-bc9c-f7ef132dbea5
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_path: Materials Tuesday morning/Vegalite.jl

max_stars_repo_name: EcoJulia/Julia-PhD-course-Copenhagen

max_issues_repo_path: Materials Tuesday morning/Vegalite.jl

max_issues_repo_name: EcoJulia/Julia-PhD-course-Copenhagen

max_forks_repo_path: Materials Tuesday morning/Vegalite.jl

max_forks_repo_name: EcoJulia/Julia-PhD-course-Copenhagen


text:
##################
# Sine series
##################

"Sine series on the interval `[0,1]`."
struct SineSeries{T} <: Dictionary{T,T}
    n   ::  Int
end

show(io::IO, b::SineSeries{Float64}) = print(io, "SineSeries($(length(b)))")
show(io::IO, b::SineSeries{T}) where T = print(io, "SineSeries{$(T)}($(length(b)))")

SineSeries(n::Int) = SineSeries{Float64}(n)

similar(b::SineSeries, ::Type{T}, n::Int) where {T} = SineSeries{T}(n)

isbasis(b::SineSeries) = true
isorthogonal(b::SineSeries, ::FourierWeight) = true

hasinterpolationgrid(b::SineSeries) = false
hasderivative(b::SineSeries) = false #for now
hasantiderivative(b::SineSeries) = false #for now
hastransform(b::SineSeries, d::GridBasis{T,G}) where {T,G <: PeriodicEquispacedGrid} = false #for now
hasextension(b::SineSeries) = true

extension(::Type{T}, src::SineSeries, dest::SineSeries; options...) where {T} = IndexExtension{T}(src, dest, 1:length(src))
restriction(::Type{T}, src::SineSeries, dest::SineSeries; options...) where {T} = IndexRestriction{T}(src, dest, 1:length(dest))

size(b::SineSeries) = (b.n,)

period(b::SineSeries{T}, idx) where {T} = T(2)

interpolation_grid(b::SineSeries{T}) where {T} = EquispacedGrid(b.n, T(0), T(1))

##################
# Native indices
##################

const SineFrequency = NativeIndex{:sine}

frequency(idxn::SineFrequency) = value(idxn)

Base.show(io::IO, idx::SineFrequency) =
	print(io, "Sine frequency $(value(idx))")

ordering(b::SineSeries) = NativeIndexList{:sine}(length(b))

##################
# Evaluation
##################

support(b::SineSeries{T}) where {T} = UnitInterval{T}()

unsafe_eval_element(b::SineSeries{T}, idx::SineFrequency, x) where {T} =
    sinpi(T(x) * frequency(idx))

function unsafe_eval_element_derivative(b::SineSeries{T}, idx::SineFrequency, x, order) where {T}
	@assert order == 1
    arg = T(pi) * frequency(idx)
    arg * cos(arg * x)
end

##################
# Differentiation
##################

derivative_dict(Φ::CosineSeries{T}, order::Int) where {T} =
	iseven(order) ? Φ : SineSeries{T}(length(Φ)-1)

diff_scaling_function(Φ::CosineSeries{T}, idx::Int, symbol) where {T} = symbol(T(π)*idx)

function differentiation(::Type{T}, src::CosineSeries, dest::CosineSeries, order::Int; options...) where {T}
	if orderiszero(order)
		@assert src==dest
		IdentityOperator{T}(src)
	else
		@assert iseven(order)
		sign = (-1)^(order>>1)
		_pseudodifferential_operator(T, src, dest, x->sign*x^order; options...)
	end
end

function differentiation(::Type{T}, src::CosineSeries, dest::SineSeries, order::Int; options...) where {T}
	@assert isodd(order)
	sign = (-1)^((order-1)>>1)
	_pseudodifferential_operator(T, src, dest, x->sign*x^order; options...)
end

## Inner products

hasmeasure(dict::SineSeries) = true

measure(dict::SineSeries{T}) where {T} = FourierWeight{T}()

gram(::Type{T}, dict::SineSeries, ::FourierWeight; options...) where {T} = ScalingOperator(one(T)/2, dict)

function innerproduct_native(b1::SineSeries, i::SineFrequency, b2::SineSeries, j::SineFrequency, m::FourierWeight;
			T = coefficienttype(b1), quad = :analytic, options...)
	if quad == :analytic
		innerproduct_sine_full(i, j, T)
	else
		innerproduct1(b1, i, b2, j, m)
	end
end

function innerproduct_sine_full(i, j, T)
	if i == j
		one(T)/2
	else
		zero(T)
	end
end


text:
module MatrixEquations
# Release V1.0

const BlasFloat = Union{Float64,Float32,ComplexF64,ComplexF32}
const BlasReal = Union{Float64,Float32}
const BlasComplex = Union{ComplexF64,ComplexF32}

using LinearAlgebra
using LinearOperators

include("lapackutil.jl")
using .LapackUtil: tgsyl!, lanv2, ladiv, lag2, lacn2!

export utqu, utqu!, qrupdate!, rqupdate!, isschur, triu2vec, vec2triu
export lanv2, ladiv, lag2, lacn2!
export lyapc, lyapd, lyapcs!, lyapds!
export plyapc, plyaps, plyapcs!, plyapd, plyapds!, plyap2, pglyap2
export arec, ared, garec, gared
export sylvc, sylvd, gsylv, sylvcs!, sylvds!, gsylvs!
export sylvsys, dsylvsys, sylvsyss!, dsylvsyss!, tgsyl!
export sylvckr, sylvdkr, gsylvkr, sylvsyskr, dsylvsyskr
export opnorm1, opnorm1est, oprcondest, opsepest, trmatop
export lyapop, invlyapop, invlyapsop
export sylvop, invsylvop, invsylvsop
export sylvsysop, invsylvsysop, invsylvsyssop

include("meutil.jl")
include("sylvester.jl")
include("lyapunov.jl")
include("riccati.jl")
include("sylvkr.jl")
include("plyapunov.jl")
include("meoperators.jl")
include("condest.jl")
end


text:
using Flux
using WindMixing
using JLD2
using FileIO
using OceanParameterizations
using OrdinaryDiffEq
using Random
using GalacticOptim
using LinearAlgebra

BLAS.set_num_threads(1)

T_fraction_str = "0.8"
T_fraction = parse(Float32, T_fraction_str)
NN_type = ARGS[1]
# N_sims = parse(Int, ARGS[2])
N_sims = 18
rate_str = ARGS[2]
rate = parse(Float64, rate_str)
params_type = ARGS[3]

# T_fraction = parse(Float32, "0.8")
# NN_type = "relu"

# train_files_all = [
#   "wind_-5e-4_diurnal_5e-8",
#   "wind_-3.5e-4_diurnal_5e-8",
#   "wind_-2e-4_diurnal_5e-8",

#   "wind_-5e-4_diurnal_1e-8",
#   "wind_-3.5e-4_diurnal_1e-8",
#   "wind_-2e-4_diurnal_1e-8",

#   "wind_-5e-4_diurnal_3e-8",
#   "wind_-3.5e-4_diurnal_3e-8",
#   "wind_-2e-4_diurnal_3e-8",
# ]

train_files_all = [
  "wind_-5e-4_cooling_3e-8_new",
  "wind_-5e-4_cooling_1e-8_new",
  "wind_-2e-4_cooling_3e-8_new",
  "wind_-2e-4_cooling_1e-8_new",
  "wind_-5e-4_heating_-3e-8_new",
  "wind_-2e-4_heating_-1e-8_new",
  "wind_-2e-4_heating_-3e-8_new",
  "wind_-5e-4_heating_-1e-8_new",

  "wind_-3.5e-4_cooling_2e-8_new",
  "wind_-3.5e-4_heating_-2e-8_new",

  "wind_-5e-4_cooling_2e-8_new",
  "wind_-3.5e-4_cooling_3e-8_new",
  "wind_-3.5e-4_cooling_1e-8_new",
  "wind_-2e-4_cooling_2e-8_new",
  "wind_-3.5e-4_heating_-3e-8_new",
  "wind_-3.5e-4_heating_-1e-8_new",
  "wind_-2e-4_heating_-2e-8_new",
  "wind_-5e-4_heating_-2e-8_new",
]

train_files = train_files_all[1:N_sims]

#
PATH = pwd()
OUTPUT_PATH = joinpath(PATH, "training_output")
# OUTPUT_PATH = "D:\\University Matters\\MIT\\CLiMA Project\\OceanParameterizations.jl\\training_output"

VIDEO_PATH = joinpath(PATH, "Output")

EXTRACTED_OUTPUT_PATH = joinpath(PATH, "extracted_training_output")

FILE_NAME = "NDE_$(N_sims)sim_windcooling_windheating_$(params_type)_divide1f5_gradient_smallNN_$(NN_type)_rate_2e-4_T$(T_fraction_str)_$(rate_str)"
# FILE_NAME = "test_$(params_type)_$(T_fraction)_$(NN_type)"
FILE_PATH = joinpath(OUTPUT_PATH, "$(FILE_NAME).jld2")
@assert !isfile(FILE_PATH)

EXTRACTED_FILE_PATH = joinpath(EXTRACTED_OUTPUT_PATH, "$(FILE_NAME)_extracted.jld2")

# if params_type == "old"
#   ν₀ = 1f-4
#   ν₋ = 1f-1
#   ΔRi = 1f-1
#   Riᶜ = 0.25f0
#   Pr = 1f0
# elseif params_type == "18simBFGST0.8grad"
#   PARAMETERS_PATH = joinpath(EXTRACTED_OUTPUT_PATH, "parameter_optimisation_18sim_windcooling_windheating_5params_BFGS_T0.8_grad_extracted.jld2")

#   parameters_file = jldopen(PARAMETERS_PATH)
#   mpp_parameters = parameters_file["parameters"]
#   close(parameters_file)

#   ν₀, ν₋, ΔRi, Riᶜ, Pr = mpp_parameters
# elseif params_type == "18simBFGST0.8nograd"
#   PARAMETERS_PATH = joinpath(EXTRACTED_OUTPUT_PATH, "parameter_optimisation_18sim_windcooling_windheating_5params_BFGS_T0.8_nograd_extracted.jld2")

#   parameters_file = jldopen(PARAMETERS_PATH)
#   mpp_parameters = parameters_file["parameters"]
#   close(parameters_file)

#   ν₀, ν₋, ΔRi, Riᶜ, Pr = mpp_parameters
# end

# PARAMETERS_PATH = joinpath(EXTRACTED_OUTPUT_PATH, "parameter_optimisation_$(N_sims)sim_windcooling_windheating_5params_ADAM1e-2_T$(T_fraction_str)_nograd_extracted.jld2")

# parameters_file = jldopen(PARAMETERS_PATH)
# mpp_parameters = parameters_file["parameters"]
# close(parameters_file)

# ν₀, ν₋, ΔRi, Riᶜ, Pr = mpp_parameters

# FILE_PATH_uw = joinpath(PATH, "extracted_training_output", "uw_NN_training_1sim_-1e-3_extracted.jld2")
# FILE_PATH_vw = joinpath(PATH, "extracted_training_output", "vw_NN_training_1sim_-1e-3_extracted.jld2")
# FILE_PATH_wT = joinpath(PATH, "extracted_training_output", "wT_NN_training_1sim_-1e-3_extracted.jld2")

# uw_file = jldopen(FILE_PATH_uw, "r")
# vw_file = jldopen(FILE_PATH_vw, "r")
# wT_file = jldopen(FILE_PATH_wT, "r")

# uw_NN = uw_file["neural_network"]
# vw_NN = vw_file["neural_network"]
# wT_NN = wT_file["neural_network"]

FILE_PATH_NN = joinpath(PATH, "extracted_training_output",
                        "NDE_$(N_sims)sim_windcooling_windheating_$(params_type)_divide1f5_gradient_smallNN_$(NN_type)_rate_2e-4_T$(T_fraction_str)_extracted.jld2")

@assert isfile(FILE_PATH_NN)

file = jldopen(FILE_PATH_NN, "r")

uw_NN = file["neural_network/uw"]
vw_NN = file["neural_network/vw"]
wT_NN = file["neural_network/wT"]

train_parameters = file["training_info/parameters"]

ν₀ = train_parameters["ν₀"]
ν₋ = train_parameters["ν₋"]
ΔRi = train_parameters["ΔRi"]
Riᶜ = train_parameters["Riᶜ"]
Pr = train_parameters["Pr"]

train_optimizers = [[ADAM(rate)]]
train_optimizers[1][1].beta = file["optimizer/β"]
train_optimizers[1][1].state = file["optimizer/state"]

close(file)

# N_inputs = 96
# hidden_units = 400
# N_outputs = 31

# if NN_type == "mish"
#   activation = mish
# elseif NN_type == "swish"
#   activation = swish
# elseif NN_type == "leakyrelu"
#   activation = leakyrelu
# elseif NN_type == "relu"
#   activation = relu
# else
#   activation = tanh
# end

# weights, re = Flux.destructure(Chain(Dense(N_inputs, hidden_units, activation), Dense(hidden_units, N_outputs)))

# uw_NN = re(weights ./ 1f5)
# vw_NN = re(weights ./ 1f5)
# wT_NN = re(weights ./ 1f5)

# gradient_scaling = 5f-3
diurnal = occursin("diurnal", train_files[1])
training_fractions = (T=T_fraction, ∂T∂z=T_fraction, profile=0.5f0)
train_parameters = Dict(
                               "ν₀" => ν₀,
                               "ν₋" => ν₋,
                              "ΔRi" => ΔRi,
                              "Riᶜ" => Riᶜ,
                               "Pr" => Pr,
                                "κ" => 1f0,
    "modified_pacanowski_philander" => true,
            "convective_adjustment" => false,
                   "smooth_profile" => false,
                        "smooth_NN" => false,
                        "smooth_Ri" => false,
                   "train_gradient" => true,
                     "zero_weights" => true,
                #  "gradient_scaling" => gradient_scaling,
               "training_fractions" => training_fractions,
                          "diurnal" => diurnal,
    )

train_epochs = [1]
train_tranges = [1:9:1153]
train_iterations = [130]

# train_epochs = [1]
# train_tranges = [1:20:200]
# train_iterations = [5]
# train_optimizers = [[ADAM(2e-4)]]

timestepper = ROCK4()

function train(FILE_PATH, train_files, train_epochs, train_tranges, train_parameters, train_optimizers, train_iterations, uw_NN, vw_NN, wT_NN, timestepper)
    @info "Writing metadata"
    write_metadata_NDE_training(FILE_PATH, train_files, train_epochs, train_tranges, train_parameters, train_optimizers, uw_NN, vw_NN, wT_NN)
    for i in 1:length(train_epochs)
        @info "iteration $i/$(length(train_epochs)), time range $(train_tranges[i])"
        if train_parameters["modified_pacanowski_philander"]
            uw_NN, vw_NN, wT_NN = train_NDE(uw_NN, vw_NN, wT_NN, train_files, train_tranges[i], timestepper, train_optimizers[i], train_epochs[i], FILE_PATH, i,
                                                         maxiters = train_iterations[i],
                                    modified_pacanowski_philander = train_parameters["modified_pacanowski_philander"],
                                            convective_adjustment = train_parameters["convective_adjustment"],
                                                               ν₀ = train_parameters["ν₀"],
                                                               ν₋ = train_parameters["ν₋"],
                                                              ΔRi = train_parameters["ΔRi"],
                                                              Riᶜ = train_parameters["Riᶜ"],
                                                               Pr = train_parameters["Pr"],
                                                                κ = train_parameters["κ"],
                                                   smooth_profile = train_parameters["smooth_profile"],
                                                        smooth_NN = train_parameters["smooth_NN"],
                                                        smooth_Ri = train_parameters["smooth_Ri"],
                                                   train_gradient = train_parameters["train_gradient"],
                                                     zero_weights = train_parameters["zero_weights"],
                                                #  gradient_scaling = train_parameters["gradient_scaling"],
                                               training_fractions = train_parameters["training_fractions"],
                                                          diurnal = train_parameters["diurnal"],
                                    )
        else
            uw_NN, vw_NN, wT_NN = train_NDE(uw_NN, vw_NN, wT_NN, train_files, train_tranges[i], timestepper, train_optimizers[i], train_epochs[i], FILE_PATH, i,
                                                         maxiters = train_iterations[i],
                                    modified_pacanowski_philander = train_parameters["modified_pacanowski_philander"],
                                            convective_adjustment = train_parameters["convective_adjustment"],
                                                                κ = train_parameters["κ"],
                                                   smooth_profile = train_parameters["smooth_profile"],
                                                        smooth_NN = train_parameters["smooth_NN"],
                                                        smooth_Ri = train_parameters["smooth_Ri"],
                                                   train_gradient = train_parameters["train_gradient"],
                                                     zero_weights = train_parameters["zero_weights"],
                                                #  gradient_scaling = train_parameters["gradient_scaling"],
                                               training_fractions = train_parameters["training_fractions"],
                                                          diurnal = train_parameters["diurnal"],
                                                )
        end
    end
    return uw_NN, vw_NN, wT_NN
end

uw_NN_res, vw_NN_res, wT_NN_res = train(FILE_PATH, train_files, train_epochs, train_tranges, train_parameters, train_optimizers, train_iterations,
                                    uw_NN, vw_NN, wT_NN, timestepper)

extract_NN(FILE_PATH, EXTRACTED_FILE_PATH, "NDE")

# test_files = [
#   "wind_-5e-4_diurnal_5e-8",
#   "wind_-5e-4_diurnal_3e-8",
#   "wind_-5e-4_diurnal_1e-8",

#   "wind_-3.5e-4_diurnal_5e-8",
#   "wind_-3.5e-4_diurnal_3e-8",
#   "wind_-3.5e-4_diurnal_1e-8",

#   "wind_-2e-4_diurnal_5e-8",
#   "wind_-2e-4_diurnal_3e-8",

#   "wind_-2e-4_diurnal_1e-8",

#   "wind_-2e-4_diurnal_2e-8",
#   "wind_-2e-4_diurnal_3.5e-8",
#   "wind_-3.5e-4_diurnal_2e-8",
#   "wind_-3.5e-4_diurnal_3.5e-8",
#   "wind_-5e-4_diurnal_2e-8",
#   "wind_-5e-4_diurnal_3.5e-8",
# ]

test_files = [
  "wind_-5e-4_cooling_3e-8_new",
  "wind_-5e-4_cooling_1e-8_new",
  "wind_-2e-4_cooling_3e-8_new",
  "wind_-2e-4_cooling_1e-8_new",
  "wind_-5e-4_heating_-3e-8_new",
  "wind_-2e-4_heating_-1e-8_new",
  "wind_-2e-4_heating_-3e-8_new",
  "wind_-5e-4_heating_-1e-8_new",

  "wind_-3.5e-4_cooling_2e-8_new",
  "wind_-3.5e-4_heating_-2e-8_new",

  "wind_-5e-4_cooling_2e-8_new",
  "wind_-3.5e-4_cooling_3e-8_new",
  "wind_-3.5e-4_cooling_1e-8_new",
  "wind_-2e-4_cooling_2e-8_new",
  "wind_-3.5e-4_heating_-3e-8_new",
  "wind_-3.5e-4_heating_-1e-8_new",
  "wind_-2e-4_heating_-2e-8_new",
  "wind_-5e-4_heating_-2e-8_new",
]

animate_training_results(test_files, FILE_NAME,
                         EXTRACTED_DATA_DIR=EXTRACTED_OUTPUT_PATH, OUTPUT_DIR=VIDEO_PATH)


text:
using StaticArrays

#Quadrature rulesfrom Shunn&Ham Paper

const triangle1a = [
    SVector{3}([0.333333333333333 0.333333333333333 0.333333333333333])
]

const triangle1w = [
  1.0,
]

const triangle3a = [
    SVector{3}([0.666666666666667 0.166666666666667 0.166666666666667]),
    SVector{3}([0.166666666666667 0.666666666666667 0.166666666666667]),
    SVector{3}([0.166666666666667 0.166666666666667 0.666666666666667])
]

const triangle3w = [
    0.333333333333333,
    0.333333333333333,
    0.333333333333333
]

const triangle6a = [
    SVector{3}([0.816847572980440 0.091576213509780 0.091576213509780]),
    SVector{3}([0.091576213509780 0.816847572980440 0.091576213509780]),
    SVector{3}([0.091576213509780 0.091576213509780 0.816847572980440]),
    SVector{3}([0.445948490915964 0.445948490915964 0.108103018168071]),
    SVector{3}([0.445948490915964 0.108103018168071 0.445948490915964]),
    SVector{3}([0.108103018168071 0.445948490915964 0.445948490915964])
]

const triangle6w = [
    0.109951743655333,
    0.109951743655333,
    0.109951743655333,
    0.223381589678000,
    0.223381589678000,
    0.223381589678000
]

const triangle10a = [
    SVector{3}([0.888871894660413 0.055564052669793 0.055564052669793]),
    SVector{3}([0.055564052669793 0.888871894660413 0.055564052669793]),
    SVector{3}([0.055564052669793 0.055564052669793 0.888871894660413]),
    SVector{3}([0.295533711735893 0.634210747745723 0.070255540518384]),
    SVector{3}([0.295533711735893 0.070255540518384 0.634210747745723]),
    SVector{3}([0.070255540518384 0.295533711735893 0.634210747745723]),
    SVector{3}([0.634210747745723 0.295533711735893 0.070255540518384]),
    SVector{3}([0.634210747745723 0.070255540518384 0.295533711735893]),
    SVector{3}([0.070255540518384 0.634210747745723 0.295533711735893]),
    SVector{3}([0.333333333333333 0.333333333333333 0.333333333333333])
]

const triangle10w = [
    0.041955512996649,
    0.041955512996649,
    0.041955512996649,
    0.112098412070887,
    0.112098412070887,
    0.112098412070887,
    0.112098412070887,
    0.112098412070887,
    0.112098412070887,
    0.201542988584730
]

const triangle15a = [
    SVector{3}([0.928258244608533 0.035870877695734 0.035870877695734]),
    SVector{3}([0.035870877695734 0.928258244608533 0.035870877695734]),
    SVector{3}([0.035870877695734 0.035870877695734 0.928258244608533]),
    SVector{3}([0.516541208464066 0.241729395767967 0.241729395767967]),
    SVector{3}([0.241729395767967 0.516541208464066 0.241729395767967]),
    SVector{3}([0.241729395767967 0.241729395767967 0.516541208464066]),
    SVector{3}([0.474308787777079 0.474308787777079 0.051382424445843]),
    SVector{3}([0.474308787777079 0.051382424445843 0.474308787777079]),
    SVector{3}([0.051382424445843 0.474308787777079 0.474308787777079]),
    SVector{3}([0.201503881881800 0.751183631106484 0.047312487011716]),
    SVector{3}([0.201503881881800 0.047312487011716 0.751183631106484]),
    SVector{3}([0.047312487011716 0.201503881881800 0.751183631106484]),
    SVector{3}([0.751183631106484 0.201503881881800 0.047312487011716]),
    SVector{3}([0.751183631106484 0.047312487011716 0.201503881881800]),
    SVector{3}([0.047312487011716 0.751183631106484 0.201503881881800])
]

const triangle15w = [
    0.017915455012303,
    0.017915455012303,
    0.017915455012303,
    0.127712195881265,
    0.127712195881265,
    0.127712195881265,
    0.076206062385535,
    0.076206062385535,
    0.076206062385535,
    0.055749810027115,
    0.055749810027115,
    0.055749810027115,
    0.055749810027115,
    0.055749810027115,
    0.055749810027115
]

const triangle21a = [
    SVector{3}([0.943774095634672 0.028112952182664 0.028112952182664]),
    SVector{3}([0.028112952182664 0.943774095634672 0.028112952182664]),
    SVector{3}([0.028112952182664 0.028112952182664 0.943774095634672]),
    SVector{3}([0.645721803061365 0.177139098469317 0.177139098469317]),
    SVector{3}([0.177139098469317 0.645721803061365 0.177139098469317]),
    SVector{3}([0.177139098469317 0.177139098469317 0.645721803061365]),
    SVector{3}([0.405508595867433 0.405508595867433 0.188982808265134]),
    SVector{3}([0.405508595867433 0.188982808265134 0.405508595867433]),
    SVector{3}([0.188982808265134 0.405508595867433 0.405508595867433]),
    SVector{3}([0.148565812270887 0.817900980028499 0.033533207700614]),
    SVector{3}([0.148565812270887 0.033533207700614 0.817900980028499]),
    SVector{3}([0.033533207700614 0.148565812270887 0.817900980028499]),
    SVector{3}([0.817900980028499 0.148565812270887 0.033533207700614]),
    SVector{3}([0.817900980028499 0.033533207700614 0.148565812270887]),
    SVector{3}([0.033533207700614 0.817900980028499 0.148565812270887]),
    SVector{3}([0.357196298615681 0.604978911775132 0.037824789609186]),
    SVector{3}([0.357196298615681 0.037824789609186 0.604978911775132]),
    SVector{3}([0.037824789609186 0.357196298615681 0.604978911775132]),
    SVector{3}([0.604978911775132 0.357196298615681 0.037824789609186]),
    SVector{3}([0.604978911775132 0.037824789609186 0.357196298615681]),
    SVector{3}([0.037824789609186 0.604978911775132 0.357196298615681])
]

const triangle21w = [
    0.010359374696538,
    0.010359374696538,
    0.010359374696538,
    0.075394884326738,
    0.075394884326738,
    0.075394884326738,
    0.097547802373242,
    0.097547802373242,
    0.097547802373242,
    0.028969269372473,
    0.028969269372473,
    0.028969269372473,
    0.028969269372473,
    0.028969269372473,
    0.028969269372473,
    0.046046366595935,
    0.046046366595935,
    0.046046366595935,
    0.046046366595935,
    0.046046366595935,
    0.046046366595935
]

const triangle28a = [
    SVector{3}([0.960045625755613 0.019977187122193 0.019977187122193]),
    SVector{3}([0.019977187122193 0.960045625755613 0.019977187122193]),
    SVector{3}([0.019977187122193 0.019977187122193 0.960045625755613]),
    SVector{3}([0.736556464940005 0.131721767529998 0.131721767529998]),
    SVector{3}([0.131721767529998 0.736556464940005 0.131721767529998]),
    SVector{3}([0.131721767529998 0.131721767529998 0.736556464940005]),
    SVector{3}([0.333333333333333 0.333333333333333 0.333333333333333]),
    SVector{3}([0.485135346793461 0.485135346793461 0.029729306413079]),
    SVector{3}([0.485135346793461 0.029729306413079 0.485135346793461]),
    SVector{3}([0.029729306413079 0.485135346793461 0.485135346793461]),
    SVector{3}([0.107951981846011 0.867911210117951 0.024136808036039]),
    SVector{3}([0.107951981846011 0.024136808036039 0.867911210117951]),
    SVector{3}([0.024136808036039 0.107951981846011 0.867911210117951]),
    SVector{3}([0.867911210117951 0.107951981846011 0.024136808036039]),
    SVector{3}([0.867911210117951 0.024136808036039 0.107951981846011]),
    SVector{3}([0.024136808036039 0.867911210117951 0.107951981846011]),
    SVector{3}([0.270840772921567 0.700872570380723 0.028286656697710]),
    SVector{3}([0.270840772921567 0.028286656697710 0.700872570380723]),
    SVector{3}([0.028286656697710 0.270840772921567 0.700872570380723]),
    SVector{3}([0.700872570380723 0.270840772921567 0.028286656697710]),
    SVector{3}([0.700872570380723 0.028286656697710 0.270840772921567]),
    SVector{3}([0.028286656697710 0.700872570380723 0.270840772921567]),
    SVector{3}([0.316549598844617 0.536654684206138 0.146795716949245]),
    SVector{3}([0.316549598844617 0.146795716949245 0.536654684206138]),
    SVector{3}([0.146795716949245 0.316549598844617 0.536654684206138]),
    SVector{3}([0.536654684206138 0.316549598844617 0.146795716949245]),
    SVector{3}([0.536654684206138 0.146795716949245 0.316549598844617]),
    SVector{3}([0.146795716949245 0.536654684206138 0.316549598844617])
]

const triangle28w = [
    0.005272170280495,
    0.005272170280495,
    0.005272170280495,
    0.044552936679504,
    0.044552936679504,
    0.044552936679504,
    0.083608212215637,
    0.033815712804198,
    0.033815712804198,
    0.033815712804198,
    0.015710461340183,
    0.015710461340183,
    0.015710461340183,
    0.015710461340183,
    0.015710461340183,
    0.015710461340183,
    0.028205136280616,
    0.028205136280616,
    0.028205136280616,
    0.028205136280616,
    0.028205136280616,
    0.028205136280616,
    0.066995957127830,
    0.066995957127830,
    0.066995957127830,
    0.066995957127830,
    0.066995957127830,
    0.066995957127830
]

const triangle36a = [
    SVector{3}([0.957657154441070 0.021171422779465 0.021171422779465]),
    SVector{3}([0.021171422779465 0.957657154441070 0.021171422779465]),
    SVector{3}([0.021171422779465 0.021171422779465 0.957657154441070]),
    SVector{3}([0.798831205208225 0.100584397395888 0.100584397395888]),
    SVector{3}([0.100584397395888 0.798831205208225 0.100584397395888]),
    SVector{3}([0.100584397395888 0.100584397395888 0.798831205208225]),
    SVector{3}([0.457923384576135 0.271038307711932 0.271038307711932]),
    SVector{3}([0.271038307711932 0.457923384576135 0.271038307711932]),
    SVector{3}([0.271038307711932 0.271038307711932 0.457923384576135]),
    SVector{3}([0.440191258403832 0.440191258403832 0.119617483192335]),
    SVector{3}([0.440191258403832 0.119617483192335 0.440191258403832]),
    SVector{3}([0.119617483192335 0.440191258403832 0.440191258403832]),
    SVector{3}([0.101763679498021 0.879979641427232 0.018256679074748]),
    SVector{3}([0.101763679498021 0.018256679074748 0.879979641427232]),
    SVector{3}([0.018256679074748 0.101763679498021 0.879979641427232]),
    SVector{3}([0.879979641427232 0.101763679498021 0.018256679074748]),
    SVector{3}([0.879979641427232 0.018256679074748 0.101763679498021]),
    SVector{3}([0.018256679074748 0.879979641427232 0.101763679498021]),
    SVector{3}([0.394033271669987 0.582562022863673 0.023404705466341]),
    SVector{3}([0.394033271669987 0.023404705466341 0.582562022863673]),
    SVector{3}([0.023404705466341 0.394033271669987 0.582562022863673]),
    SVector{3}([0.582562022863673 0.394033271669987 0.023404705466341]),
    SVector{3}([0.582562022863673 0.023404705466341 0.394033271669987]),
    SVector{3}([0.023404705466341 0.582562022863673 0.394033271669987]),
    SVector{3}([0.226245530909229 0.751530614542782 0.022223854547989]),
    SVector{3}([0.226245530909229 0.022223854547989 0.751530614542782]),
    SVector{3}([0.022223854547989 0.226245530909229 0.751530614542782]),
    SVector{3}([0.751530614542782 0.226245530909229 0.022223854547989]),
    SVector{3}([0.751530614542782 0.022223854547989 0.226245530909229]),
    SVector{3}([0.022223854547989 0.751530614542782 0.226245530909229]),
    SVector{3}([0.635737183263105 0.249079227621332 0.115183589115563]),
    SVector{3}([0.635737183263105 0.115183589115563 0.249079227621332]),
    SVector{3}([0.115183589115563 0.635737183263105 0.249079227621332]),
    SVector{3}([0.249079227621332 0.635737183263105 0.115183589115563]),
    SVector{3}([0.249079227621332 0.115183589115563 0.635737183263105]),
    SVector{3}([0.115183589115563 0.249079227621332 0.635737183263105])

]

const triangle36w = [
    0.005639123786910,
    0.005639123786910,
    0.005639123786910,
    0.027148968192278,
    0.027148968192278,
    0.027148968192278,
    0.063100912533359,
    0.063100912533359,
    0.063100912533359,
    0.051752795679899,
    0.051752795679899,
    0.051752795679899,
    0.009866753574646,
    0.009866753574646,
    0.009866753574646,
    0.009866753574646,
    0.009866753574646,
    0.009866753574646,
    0.022008204800147,
    0.022008204800147,
    0.022008204800147,
    0.022008204800147,
    0.022008204800147,
    0.022008204800147,
    0.016644570076736,
    0.016644570076736,
    0.016644570076736,
    0.016644570076736,
    0.016644570076736,
    0.016644570076736,
    0.044326238118914,
    0.044326238118914,
    0.044326238118914,
    0.044326238118914,
    0.044326238118914,
    0.044326238118914
]

const triangleA = [
    triangle1a,
    triangle3a,
    triangle6a,
    triangle10a,
    triangle15a,
    triangle21a,
    triangle28a,
    triangle36a,
]

const triangleW = [
    triangle1w,
    triangle3w,
    triangle6w,
    triangle10w,
    triangle15w,
    triangle21w,
    triangle28w,
    triangle36w
]


text:
using DynamicalSystemsBase, SimpleDiffEq, DiffEqBase
using Statistics, Test

ds = Systems.lorenz()
ALG = SimpleATsit5()

# %%

i1 = tangent_integrator(ds, 3; alg = ALG)
i2 = tangent_integrator(ds, 3; alg = ALG,
internalnorm = DiffEqBase.ODE_DEFAULT_NORM)

step!(i1)
step!(i2)
dts = zeros(2, 500)
for i in 1:500
    step!(i1); step!(i2)
    dts[:, i] .= (i1.dt, i2.dt)
end

dt1 = mean(dts[1, :])
dt2 = mean(dts[2, :])

@test dt1 > dt2

# %%
s = [get_state(ds), get_state(ds) .+ rand(3)]
i1 = parallel_integrator(ds, deepcopy(s);
alg = ALG, internalnorm = DiffEqBase.ODE_DEFAULT_NORM)
i2 = parallel_integrator(ds, deepcopy(s);
alg = ALG)

step!(i1)
step!(i2)
dts = zeros(2, 1000)
for i in 1:1000
    step!(i1); step!(i2)
    dts[:, i] .= (i1.dt, i2.dt)
end

dt1 = mean(dts[1, :])
dt2 = mean(dts[2, :])

@test dt2 > dt1


text:
using Decimals
using Base.Test

# Rounding
@test round(decimal(7.123456), 0) == decimal(7)
@test round(decimal(7.123456), 2) == decimal(7.12)
@test round(decimal(7.123456), 3) == decimal(7.123)
@test round(decimal(7.123456), 5) == decimal(7.12346)
@test round(decimal(7.123456), 6) == decimal(7.123456)

@test round(0.123, 1) == 0.1
@test round([0.1111 0.2222 0.8888], 2) == [0.11 0.22 0.89]

function tet()
    a = decimal("1.0000001")
    for i = 1:27
        a *= a
    end
    return a
end

# set DIGITS = 20 (aaljuffali's example)
@test tet() == Decimal(0, 67453047074102193157641340, -20)


text:
# Update packages
using Pkg
Pkg.update()
Pkg.build()


text:
module Unmarked

using StatsBase, StatsModels, StatsFuns, Distributions, Random, Printf
using LinearAlgebra, DataFrames, Optim, NLSolversBase, ForwardDiff
using CategoricalArrays
using Combinatorics: combinations
using PrettyTables: pretty_table, ft_printf
using ProgressMeter: Progress, next!
using Gadfly: plot, style, Guide, Geom, pt, vstack, hstack

import Base: show, getindex, length
import StatsBase: aic, aicc, bic, coef, coefnames, coeftable, deviance, dof,
                  loglikelihood, modelmatrix, predict, nobs, stderror, vcov

export UmData, Occu, Nmix, RN
export @formula
export aic, aicc, bic, coef, coefnames, coeftable, deviance, dof, gof,
       loglikelihood, modelmatrix, nobs, predict, simulate, stderror, vcov

#Fitting functions
export occu, nmix, rn

#Submodel extractors
export detection, occupancy, abundance

#Utility functions
export allsub, whiskerplot, effectsplot

include("utils.jl")
include("links.jl")
include("formula.jl")
include("data.jl")
include("design.jl")
include("fit.jl")
include("gof.jl")
include("predict.jl")
include("occupancy.jl")
include("nmix.jl")
include("roylenichols.jl")
include("plots.jl")

end

max_stars_repo_name: kenkellner/julia-hme

max_issues_repo_name: kenkellner/julia-hme

max_forks_repo_name: kenkellner/julia-hme


text:
function get_project_deps(api::GitHub.GitHubAPI,
                          clone_hostname::HostnameForClones,
                          repo::GitHub.Repo;
                          auth::GitHub.Authorization,
                          master_branch::Union{DefaultBranch, AbstractString},
                          subdir::AbstractString,
                          include_jll::Bool)
    original_directory = pwd()
    tmp_dir = mktempdir(; cleanup = true)
    url_with_auth = "https://x-access-token:$(auth.token)@$(clone_hostname.hostname)/$(repo.full_name).git"
    cd(tmp_dir)
    my_retry(() -> run(`git clone $(url_with_auth) REPO`))
    cd(joinpath(tmp_dir, "REPO"))
    default_branch = git_get_current_branch()
    master_branch_name = git_decide_master_branch(master_branch, default_branch)
    run(`git checkout $(master_branch_name)`)
    project_file = joinpath(tmp_dir, "REPO", subdir, "Project.toml")
    dep_to_current_compat_entry, dep_to_current_compat_entry_verbatim,
                                 dep_to_latest_version,
                                 deps_with_missing_compat_entry = get_project_deps(project_file; include_jll)
    cd(original_directory)
    rm(tmp_dir; force = true, recursive = true)
    result = dep_to_current_compat_entry, dep_to_current_compat_entry_verbatim,
                                          dep_to_latest_version,
                                          deps_with_missing_compat_entry
    return result
end

function get_project_deps(project_file::String; include_jll::Bool)
    dep_to_current_compat_entry = Dict{Package, Union{Pkg.Types.VersionSpec, Nothing}}()
    dep_to_current_compat_entry_verbatim = Dict{Package, Union{String, Nothing}}()
    dep_to_latest_version = Dict{Package, Union{VersionNumber, Nothing}}()
    deps_with_missing_compat_entry = Set{Package}()
    get_project_deps!(dep_to_current_compat_entry,
                      dep_to_current_compat_entry_verbatim,
                      dep_to_latest_version,
                      deps_with_missing_compat_entry,
                      project_file;
                      include_jll)
    result = dep_to_current_compat_entry, dep_to_current_compat_entry_verbatim,
                                          dep_to_latest_version,
                                          deps_with_missing_compat_entry
    return result
end

function get_project_deps!(dep_to_current_compat_entry::Dict{Package, Union{Pkg.Types.VersionSpec, Nothing}},
                           dep_to_current_compat_entry_verbatim::Dict{Package, Union{String, Nothing}},
                           dep_to_latest_version::Dict{Package, Union{VersionNumber, Nothing}},
                           deps_with_missing_compat_entry::Set{Package},
                           project_file::String;
                           include_jll::Bool)
    project = TOML.parsefile(project_file)
    if haskey(project, "deps")
        deps = project["deps"]
        add_compat_section!(project)
        compat = project["compat"]
        stdlib_uuids = gather_stdlib_uuids()
        for d in deps
            name = d[1]
            uuid = UUIDs.UUID(d[2])
            if uuid in stdlib_uuids
                @debug("Skipping stdlib: $(uuid)")
            elseif (endswith(lowercase(strip(name)), lowercase(strip("_jll")))) && (!include_jll)
                @debug("Skipping JLL package: $(name)")
            else
                package = Package(name, uuid)
                dep_to_latest_version[package] = nothing
                compat_entry = convert(String, strip(get(compat, name, "")))::String
                if length(compat_entry) > 0
                    dep_to_current_compat_entry_verbatim[package] = compat_entry
                    dep_to_current_compat_entry[package] = try
                        Pkg.Types.semver_spec(compat_entry)
                    catch
                        nothing
                    end
                else
                    push!(deps_with_missing_compat_entry, package)
                    dep_to_current_compat_entry[package] = nothing
                    dep_to_current_compat_entry_verbatim[package] = nothing
                end
            end
        end
    else
        @info("This project has no dependencies.")
    end
    result = dep_to_current_compat_entry, dep_to_current_compat_entry_verbatim,
                                          dep_to_latest_version,
                                          deps_with_missing_compat_entry
    return result
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Definition of the types and structures.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export Design_Variable_MGEO_Canonical
export Design_Variable_MGEO_Var
export Design_Variable_MGEO_Real
export MGEO_Structure, Pareto_Point
export MGEO_Canonical, MGEO_Var, MGEO_Real

MGEO_Canonical() = Val{:MGEO_Canonical}
MGEO_Var()       = Val{:MGEO_Var}
MGEO_Real()      = Val{:MGEO_Real}

################################################################################
#                                  Structures
################################################################################

#                              Public Structures
# ==============================================================================

"""
    struct Design_Variable_MGEO_Canonical

Structure that defines the design variables for MGEO Canonical.

# Fields

* `bits`: Number of bits.
* `min`: Minimum allowed value for the design variable.
* `max`: Maximum allowed value for the design variable.
* `full_scale`: Full scale of the design variable.
* `index`: Index in the string.
* `name`: Name of the variable.

"""
struct Design_Variable_MGEO_Canonical
	bits::Int64
	min::Float64
	max::Float64
	full_scale::Int64
	index::Int64
	name::AbstractString
end

"""
    struct Design_Variable_MGEO_Var

Structure that defines the design variables for MGEO Var.

# Fields

* `bits`: Number of bits.
* `min`: Minimum allowed value for the design variable.
* `max`: Maximum allowed value for the design variable.
* `full_scale`: Full scale of the design variable.
* `index`: Index in the string.
* `name`: Name of the variable.

"""
struct Design_Variable_MGEO_Var
	bits::Int64
	min::Float64
	max::Float64
	full_scale::Int64
	index::Int64
	name::AbstractString
end

"""
    struct Design_Variable_MGEO_Real

Structure that defines the design variables for MGEO Real.

# Fields

* `min`: Minimum allowed value for the design variable.
* `max`: Maximum allowed value for the design variable.
* `σ`: Standard deviation used when modifying the variable.
* `name`: Name of the variable.

"""
struct Design_Variable_MGEO_Real
	min::Float64
	max::Float64
    σ::Float64
	name::AbstractString
end

"""
    struct MGEO_Structure{Nv, Nf, T}

Structure to store the configuration of MGEO.

`Nv` is the number of design variables, `Nf` is the number of objective
functions, and `T` is the MGEO algorithm that will be used.

# Fields

* `tau`: Parameter to set the search determinism.
* `ngen_max`: Maximum number of generations.
* `run_max`: Maximum number of independent runs (reinitializations).
* `design_vars`: Structure that store the configuration for each design variable.
* `mgeo_eps`: Epsilon to compare objective functions.

"""
struct MGEO_Structure{Nv, Nf, T}
    τ::Float64
    ngen_max::Int64
    run_max::Int64
    design_vars::SVector{Nv, T}
    mgeo_eps::Float64
end

"""
    struct Pareto_Point{Nv, Nf}

Structure that defines a point in the Pareto frontier.

`Nv` is the number of design variables and `Nf` is the number of objective
functions.

# Fields

* `vars`: Design variables.
* `f`: Objective functions.

"""
struct Pareto_Point{Nv, Nf}
    vars::SVector{Nv, Float64}
    f::SVector{Nf, Float64}
end

#                              Private Structures
# ==============================================================================

"""
    struct sRank

Structure used to sort the candidate points in MGEO Canonical and MGEO Var.

# Fields

* `valid`: Is the value valid?
* `index`: Index of the flipped bit.
* `f`: Objective function value after flipping the bit.

"""
struct sRank
    valid::Bool
    index::Int64
    f::Float64
end

"""
    struct sRank_Real

Structure used to sort the candidate points in MGEO Real.

# Fields

* `valid`: Is the value valid?
* `index`: Index of the modified variable.
* `var_value`: Value of the modified variable.
* `f`: Objective function value after flipping the bit.

"""
struct sRank_Real
    valid::Bool
    index::Int64
    var_value::Float64
    f::Float64
end