text:
# Dispersion models for NLO-relevant materials adapted from
# https://github.com/doddgray/optics_modeling/blob/master/nlo/NLO_tools.py

"""
################################################################################
#																			   #
#							Physical Material Models						   #
#																			   #
################################################################################
"""

# ng(fₙ::Function) = λ -> ( (n,n_pb) = Zygote.pullback(fₙ,λ); ( n - λ * n_pb(1)[1] ) )
export λ, Material, materials, n, ng, gvd, ε, ε⁻¹, unique_axes, plot_data, nn̂g, nĝvd, rotate #, fn, fng, fgvd, fε, fε⁻¹,

c = Unitful.c0      # Unitful.jl speed of light
@parameters λ, T
Dλ = Differential(λ)
DT = Differential(T)
ng(n_sym::Num) = n_sym - λ * expand_derivatives(Dλ(n_sym))
gvd(n_sym::Num) = λ^3 * expand_derivatives(Dλ(Dλ(n_sym))) # gvd = uconvert( ( 1 / ( 2π * c^2) ) * _gvd(lm_um,T_C)u"μm", u"fs^2 / mm" )

struct Material{T}
	ε::T
	fε::Function
	fng::Function
	fgvd::Function
	# fχ⁽²⁾::Function
end

# Material(x::AbstractArray) = Material(ε_tensor(x))
# Material(x::Number) = Material(ε_tensor(x))
# Material(x::Tuple) = Material(ε_tensor(x))

# Material(x,fe,fng) = Material(ε_tensor(x),fe,fng)
# Material(x) = Material(ε_tensor(x),lm->ε_tensor(x),lm->sqrt.(ε_tensor(x)))
# Material(mat::Material) = mat	# maybe

Material(x) = Material(lm->ε_tensor(x),lm->sqrt.(ε_tensor(x)))
Material(mat::Material) = mat	# maybe

# import Base.convert
# convert(::Type{Material}, x) = Material(x)

materials(shapes::Vector{<:GeometryPrimitives.Shape}) = Zygote.@ignore(unique(Material.(getfield.(shapes,:data)))) # unique!(getfield.(shapes,:data))

#### Fixed index Material constructors and methods ####

# Material(n::T where T<:Real) = Material{T}(ε_tensor(n))							# isotropic fixed index Material
# Material(ns::Tuple{T,T,T} where T<:Real) = Material{T}(ε_tensor(ns))			# anisotropic fixed index Material
# Material(ε::AbstractMatrix{T} where T<:Real) = Material{T}(ε_tensor(ε))			# fixed dielectric tensor material

# ε(mat::Material) = mat.ε
# ε⁻¹(mat::Material) = inv(mat.ε)
# n(mat::Material) = sqrt.(diag(mat.ε))
# n(mat::Material,axind::Int) = sqrt(mat.ε[axind,axind])
nn̂g(mat::Material,lm::Real) = getfield(mat,:fng)(lm) * sqrt.(getfield(mat,:fε)(lm))
nĝvd(mat::Material,lm::Real) = getfield(mat,:fgvd)(lm) * sqrt.(getfield(mat,:fε)(lm))

#### Symbolic Material model methods ####
# ng(mat::Material) = ng.(n(mat))
# ng(mat::Material,axind::Int) = ng(n(mat,axind))
# gvd(mat::Material) = gvd.(n(mat))
# gvd(mat::Material,axind::Int) = gvd(n(mat,axind))

# ### methods to generate numerical functions for symbolic material dispersions ###
# # non-mutating numerical function generators #
# fε(mat::Material{Num}) = eval(build_function(ε(mat),λ)[1])
# fε⁻¹(mat::Material{Num}) = eval(build_function(ε⁻¹(mat),λ)[1])
# fn(mat::Material{Num}) = eval(build_function(n(mat),λ)[1])
# fn(mat::Material{Num},axind::Int) = eval(build_function(n(mat,axind),λ)[1])
# fng(mat::Material{Num}) = eval(build_function(ng(mat),λ)[1])
# fng(mat::Material{Num},axind::Int) = eval(build_function(ng(mat,axind),λ)[1])
# fgvd(mat::Material{Num}) = eval(build_function(gvd(mat),λ)[1])
# fgvd(mat::Material{Num},axind::Int) = eval(build_function(gvd(mat,axind),λ)[1])
# # mutating numerical function generators #
# fε!(mat::Material{Num}) = eval(build_function(ε(mat),λ)[2])
# fε⁻¹!(mat::Material{Num}) = eval(build_function(ε⁻¹(mat),λ)[2])
# fn!(mat::Material{Num}) = eval(build_function(n(mat),λ)[2])
# fn!(mat::Material{Num},axind::Int) = eval(build_function(n(mat,axind),λ)[2])
# fng!(mat::Material{Num}) = eval(build_function(ng(mat),λ)[2])
# fng!(mat::Material{Num},axind::Int) = eval(build_function(ng(mat,axind),λ)[2])
# fgvd!(mat::Material{Num}) = eval(build_function(gvd(mat),λ)[2])
# fgvd!(mat::Material{Num},axind::Int) = eval(build_function(gvd(mat,axind),λ)[2])

# # Fallback methods for fixed-index Materials
# ng(mat::Material{<:Real}) = n(mat)
# ng(mat::Material{<:Real},axind::Int) = n(mat,axind)
# function gvd(mat::Material{T})  where T<:Real
# 	zeros(T,3)
# end
# function gvd(mat::Material{T},axind::Int) where T<:Real
# 	zero(T)
# end

# fε(mat::Material) = x->ε(mat)
# fε⁻¹(mat::Material) = x->ε⁻¹(mat)
# fn(mat::Material) = x->n(mat)
# fn(mat::Material,axind::Int) = x->n(mat,axind)
# fng(mat::Material) = x->ng(mat)
# fng(mat::Material,axind::Int) = x->ng(mat,axind)
# fgvd(mat::Material) = x->gvd(mat)
# fgvd(mat::Material,axind::Int) = x->gvd(mat,axind)

"""
################################################################################
#																			   #
#							   	  Rotations							   		   #
#																			   #
################################################################################
"""

function rotate(χ::AbstractArray{T,2},𝓡::TR) where {T<:Real, TR<:StaticMatrix{3,3}}
	@tullio χᵣ[i,j] := 𝓡[a,i] * 𝓡[b,j] * χ[a,b]  fastmath=true
end

function rotate(χ::AbstractArray{T,3},𝓡::TR) where {T<:Real, TR<:StaticMatrix{3,3}}
	@tullio χᵣ[i,j,k] := 𝓡[a,i] * 𝓡[b,j] * 𝓡[c,k] * χ[a,b,c]  fastmath=true
end

function rotate(χ::AbstractArray{T,4},𝓡::TR) where {T<:Real, TR<:StaticMatrix{3,3}}
	@tullio χᵣ[i,j,k,l] := 𝓡[a,i] * 𝓡[b,j] * 𝓡[c,k] * 𝓡[d,l] * χ[a,b,c,d]  fastmath=true
end

# struct RotatedMaterial{T,TR}
# 	parent::Material{T}
# 	fε::Function
# 	fng::Function
# 	fgvd::Function
# end

"""
################################################################################
#																			   #
#							    Utility methods					   			   #
#																			   #
################################################################################
"""

function unique_axes(mat::Material)
	e11,e22,e33 = diag(getfield(mat,:ε))
	if isequal(e11,e22)
		isequal(e11,e33) ? (return ( [1,], [""] )) : (return ( [1,3], ["₁,₂","₃"] )) # 1 == 2 == 3 (isotropic) : 1 == 2 != 3 (uniaxial)
	elseif isequal(e22,e33)
		return ( [1,2], ["₁","₂,₃"] )	# 1 != 2 == 3 (uniaxial)
	else
		isequal(e11,e33) ? (return ( [1,2], ["₁,₃","₂"] )) : (return ( [1,2,3], ["₁","₂","₃"] )) # 1 == 3 != 2 (uniaxial) : 1 != 2 != 3 (biaxial)
	end
end

"""
################################################################################
#																			   #
#							   Plotting methods					   			   #
#																			   #
################################################################################
"""

function plot_data(mats::AbstractVector{<:Material})
	fes = getfield.(mats,:fε)
	axind_axstr_unq = unique_axes.(mats)
	axind_unq = getindex.(axind_axstr_unq,1)
	axstr_unq = getindex.(axind_axstr_unq,2)
	fns = vcat(map((ff,as)->[(x->sqrt(ff(x)[a,a])) for a in as ], fes, axind_unq)...)
	mat_names = chop.(String.(nameof.(fes)),head=2,tail=0)	# remove "ε_" from function names
	names = "n" .* vcat([.*(axstr_unq[i], " (", mat_names[i],")") for i=1:length(mats)]...) # "n, n_i or n_i,j (Material)" for all unique axes and materials
	return fns, names
end
plot_data(mat::Material) = plot_data([mat,])
plot_data(mats::NTuple{N,<:Material} where N) = plot_data([mats...])

function uplot(x::Union{Material, AbstractVector{<:Material}, NTuple{N,<:Material} };
		xlim=[0.5,1.8], xlabel="λ [μm]", ylabel="n", kwargs...)  where N
	fns, name = plot_data(x)
	UnicodePlots.lineplot(fns, xlim[1], xlim[2];
	 	xlim,
		ylim=map((a,b)->a(b,digits=1),(floor,ceil),ylims(fns;xlims=xlim)),
		name,
		xlabel,
		ylabel,
		kwargs...
		)
end

function uplot!(plt::UnicodePlots.Plot,x::Union{Material, AbstractVector{<:Material}, NTuple{N,<:Material} };
		xlim=[0.5,1.8], xlabel="λ [μm]", ylabel="n")  where N
	fns, name = plot_data(x)
	UnicodePlots.lineplot!(plt, fns; name ) #, xlim[1], xlim[2];
	 	# xlim,
		# ylim=round.( ylims(plt,ylims(fns;xlims=xlim)) ,digits=1),
		# name,
		# xlabel,
		# ylabel,
		# )
end
#
import Base: show
Base.show(io::IO, ::MIME"text/plain", mat::Material) = uplot(mat;title="MIME version") #print(io, "Examplary instance of Material\n", m.x, " ± ", m.y)
Base.show(io::IO, mat::Material) = uplot(mat) #print(io, m.x, '(', m.y, ')')

################################################################################
#                                Load Materials                                #
################################################################################
include("material_lib/MgO_LiNbO3.jl")
include("material_lib/SiO2.jl")
include("material_lib/Si3N4.jl")
include("material_lib/αAl2O3.jl")
# include("material_lib/LiB3O5.jl")
# include("material_lib/silicon.jl")
# include("material_lib/GaAs.jl")
# include("material_lib/MgF2.jl")
# include("material_lib/HfO2.jl")


text:
function [J]=Jac(U,CG,Param)
nz=Param.Grid.nz;
dz=Param.Grid.dz;
nCol=size(U,1);
nJ=nCol*nz;
nz=Param.Grid.nz;
RhoPos=Param.RhoPos;
uPos=Param.uPos;
vPos=Param.vPos;
wPos=Param.wPos;
ThPos=Param.ThPos;
% W=[zeros(size(U,1),1)...
%    ,U(:,:,wPos)];

RhoCol=permute(U(:,:,RhoPos),[2 1]);
ThCol=permute(U(:,:,ThPos),[2 1]);
if strcmp(Param.Thermo,'Energy')
  W=[zeros(size(U,1),1)...
    ,U(:,:,wPos)];
  WC=0.5*(W(:,1:nz)+W(:,2:nz+1));
  KE=0.5*(U(:,:,Param.uPos).*U(:,:,Param.uPos)+U(:,:,Param.vPos).*U(:,:,Param.vPos)+WC.*WC);
  Pres=Pressure(U(:,:,ThPos),U(:,:,RhoPos),KE,Param);
  PresCol=permute(Pres,[2 1]);
  PhiCol=Param.Grav*repmat(Param.Grid.zP,1,size(U,1));
else
  Pres=Pressure(U(:,:,ThPos),U(:,:,ThPos),U(:,:,ThPos),Param);
  PresCol=permute(Pres,[2 1]);
end

if strcmp(Param.ModelType,'Curl')
  D=zeros(nz,nCol);
  D(1:nz-1,:)=-0.5*(RhoCol(1:nz-1,:)+RhoCol(2:nz,:))/dz;
  JRhoW=spdiags([reshape(-D,nJ,1) reshape(D,nJ,1)],[-1 0]...
    ,nJ,nJ);
  if strcmp(Param.Thermo,'Energy')
    dPdTh=(Param.Rd/Param.Cvd)./RhoCol;
  else
    dPdTh=dPresdTh(ThCol,Param);
  end
  Dp=zeros(nz,nCol);
  Dm=zeros(nz,nCol);
  Dp(1:nz-1,:)=dPdTh(2:nz,:)...
    ./(0.5*(RhoCol(1:nz-1,:)+RhoCol(2:nz,:)))...
    /dz;
  Dm(1:nz-1,:)=dPdTh(1:nz-1,:)...
    ./(0.5*(RhoCol(1:nz-1,:)+RhoCol(2:nz,:)))...
    /dz;
  JWTh=spdiags([[reshape(Dm,nJ,1);0] [0;reshape(-Dp,nJ,1)]],[0 1]...
    ,nJ,nJ);

  if strcmp(Param.Thermo,'Energy')
    Dp=zeros(nz,nCol);
    Dm=zeros(nz,nCol);
    Dp(1:nz-1,:)=(Param.Rd/Param.Cvd)...
      ./(0.5*(RhoCol(1:nz-1,:)+RhoCol(2:nz,:)))...
      .*(PhiCol(2:nz,:))/dz;
    Dm(1:nz-1,:)=(Param.Rd/Param.Cvd)...
      ./(0.5*(RhoCol(1:nz-1,:)+RhoCol(2:nz,:)))...
      .*(PhiCol(1:nz-1,:))/dz;
    D(1:nz-1,:)=-0.5*Param.Grav...
      ./(0.5*(RhoCol(1:nz-1,:)+RhoCol(2:nz,:)));
    JWRho=spdiags([[reshape(-Dm+D,nJ,1);0] [0;reshape(Dp+D,nJ,1)]],[0 1]...
      ,nJ,nJ);
  else
    D(1:nz-1,:)=0.5*(PresCol(2:nz,:)-PresCol(1:nz-1,:))/Param.Grid.dz...
      ./(0.5*(RhoCol(1:nz-1,:)+RhoCol(2:nz,:))).^2;
    JWRho=spdiags([[reshape(D,nJ,1);0] [0;reshape(D,nJ,1)]],[0 1]...
      ,nJ,nJ);
  end

  if strcmp(Param.Thermo,'Energy')
    D(1:nz-1,:)=-0.5*(ThCol(1:nz-1,:)+ThCol(2:nz,:)...
      +PresCol(1:nz-1,:)+PresCol(2:nz,:))/dz;
  else
    D(1:nz-1,:)=-0.5*(ThCol(1:nz-1,:)+ThCol(2:nz,:))/dz;
  end
  JThW=spdiags([reshape(-D,nJ,1) reshape(D,nJ,1)],[-1 0]...
    ,nJ,nJ);
else
  % Conservative
  D=zeros(nz,nCol);
  D(1:nz-1,:)=-1/dz;
  JRhoW=spdiags([reshape(-D,nJ,1) reshape(D,nJ,1)],[-1 0]...
    ,nJ,nJ);

  dPdTh=dPresdTh(ThCol,Param);
  Dp=zeros(nz,nCol);
  Dm=zeros(nz,nCol);
  Dp(1:nz-1,:)=dPdTh(2:nz,:)/dz;
  Dm(1:nz-1,:)=dPdTh(1:nz-1,:)/dz;
  JWTh=spdiags([[reshape(Dm,nJ,1);0] [0;reshape(-Dp,nJ,1)]],[0 1]...
    ,nJ,nJ);
  D(1:nz-1,:)=-0.5*Param.Grav;
  JWRho=spdiags([[reshape(D,nJ,1);0] [0;reshape(D,nJ,1)]],[0 1]...
    ,nJ,nJ);
  Temp=ThCol./RhoCol;
  D(1:nz-1,:)=-0.5*(Temp(1:nz-1,:)+Temp(2:nz,:))/dz;
  JThW=spdiags([reshape(-D,nJ,1) reshape(D,nJ,1)],[-1 0]...
    ,nJ,nJ);
end
J=[sparse(nJ,nJ) sparse(nJ,nJ) sparse(nJ,nJ) JRhoW sparse(nJ,nJ)
  sparse(nJ,nJ) sparse(nJ,nJ) sparse(nJ,nJ) sparse(nJ,nJ) sparse(nJ,nJ)
  sparse(nJ,nJ) sparse(nJ,nJ) sparse(nJ,nJ) sparse(nJ,nJ) sparse(nJ,nJ)
  JWRho sparse(nJ,nJ)  sparse(nJ,nJ) sparse(nJ,nJ) JWTh
  sparse(nJ,nJ) sparse(nJ,nJ) sparse(nJ,nJ) JThW sparse(nJ,nJ)];
end

ext: jl


text:
"""
    brute_force(fla::Formula)::SatSolverResult

Tries all possible variable assignments and stops at the first satisfying assignment.
"""
function brute_force(fla::Formula)::SatSolverResult
  pa = PartialAssignment(fill(-1, fla.variable_count))

  while true
    if all((cl) -> is_clause_satisfied(cl, pa), fla)
      return SatisfiableFormula(pa)
    end

    last_negative = findlast(pa.assignment) do x
      -1 == x
    end

    if nothing == last_negative
      break
    end

    pa[last_negative] = 1
    for i in (last_negative + 1):fla.variable_count
      pa[i] = -1
    end
  end

  UnsatisfiableFormula()
end


text:
module MaxHelpingHandOneWayCameraTrigger

using ..Ahorn, Maple

@mapdef Trigger "MaxHelpingHand/OneWayCameraTrigger" OneWayCameraTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight,
    left::Bool=true, right::Bool=true, up::Bool=true, down::Bool=true, flag::String="", blockPlayer::Bool=false)

const placements = Ahorn.PlacementDict(
    "One-Way Camera Trigger (max480's Helping Hand)" => Ahorn.EntityPlacement(
        OneWayCameraTrigger,
        "rectangle"
    )
)

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text: println("This is a stub file that does nothing and is probably obsolete..? But still recommended by the current julia docs.")

max_stars_repo_name: Skenvy/Collatz

max_issues_repo_name: Skenvy/Collatz

max_forks_repo_name: Skenvy/Collatz


text:
# Monte Carlo solution with primitive thermodynamic estimators.

"""
    get_sample_pt(sys::System{S,M}, pseudosp::SamplingParameters{S,M,P}, sp::SamplingParameters{S_,M,P})

Compute a sample for `sys` using `pseudosp` and `sp`.
"""
function get_sample_pt(sys::System{S,M}, pseudosp::SamplingParameters{S,M,P}, sp::SamplingParameters{S_,M,P}) where {S,S_,M,P}
    # Choose a surface.
    s_ = sample(1:S_, sp.weights)

    # Sample coordinates.
    qs = zeros(P, M)
    for m in 1:M
        qs[:, m] .= rand(sp.mvns[m, s_])
    end

    # Calculate the numerator, energy, heat capacity, and denominator matrices.
    num = Matrix{Float64}(I, S, S)
    num_Es = [Matrix{Float64}(I, S, S) for _ in 1:P]
    num_Cvs = [Matrix{Float64}(I, S, S) for _ in 1:P]
    denom = Matrix{Float64}(I, S_, S_)

    for i1 in 1:P
        # Next index in cyclic path.
        i2 = i1%P+1

        # Numerator.
        FM_num, scaling = sampling_matrix_free_particle(pseudosp, qs[i1, :], qs[i2, :])
        EM1, IM = sampling_matrix_interaction(sys, pseudosp, qs[i1, :])
        EM2 = sampling_matrix_energy(pseudosp, qs[i1, :], qs[i2, :])
        EM = EM1 - EM2
        if i1 == 1
            CM1 = sampling_matrix_heat_capacity(pseudosp, qs[i1, :], qs[i2, :])
            CM = EM1 * EM + CM1 - EM * EM2
        end
        IFM = IM * FM_num
        IEFM = IM * EM * FM_num
        num *= IM * FM_num
        for i in 1:P
            if i == i1
                num_Es[i] *= IEFM
                if i1 == 1
                    num_Cvs[i] *= IM * CM * FM_num
                else
                    num_Cvs[i] *= IEFM
                end
            else
                num_Es[i] *= IFM
                if i1 == 1
                    num_Cvs[i] *= IEFM
                else
                    num_Cvs[i] *= IFM
                end
            end
        end

        # Denominator.
        FM_denom, _ = sampling_matrix_free_particle(sp, qs[i1, :], qs[i2, :]; scaling=scaling)
        denom *= FM_denom
    end

    trace_num = tr(num)
    trace_num_E = mean(tr(num_E) for num_E in num_Es)
    trace_num_Cv = mean(tr(num_Cv) for num_Cv in num_Cvs)
    trace_denom = tr(denom)

    [trace_num/trace_denom,
     trace_num_E/trace_denom,
     trace_num_Cv/trace_denom]
end

"""
Monte Carlo solution for an arbitrary system using primitive thermodynamic
estimators.
"""
struct SamplingPrimitiveThermodynamic <: Sampling
    "Partition function."
    Z::Float64
    "Partition function standard error."
    Z_err::Float64
    "Energy."
    E::Float64
    "Energy standard error."
    E_err::Float64
    "Heat capacity."
    Cv::Float64
    "Heat capacity standard error."
    Cv_err::Float64
end

"""
    SamplingPrimitiveThermodynamic(sys::System, beta::Float64, P::Int, num_samples::Int; sampling_sys::Maybe{System}=nothing, progress_output::IO=stderr)

Calculate the solution for `sys` at `beta` with `P` links and `num_samples`
random samples.

If `sampling_sys` is provided, it is used for sampling. Otherwise, sampling
defaults to the simplified diagonal subsystem of `sys`.

The progress meter is written to `progress_output`.
"""
function SamplingPrimitiveThermodynamic(sys::System, beta::Float64, P::Int, num_samples::Int; sampling_sys::Maybe{System}=nothing, progress_output::IO=stderr)
    sys_diag = diag(sys)
    sys_diag_simple = simplify(sys_diag)
    pseudosp = SamplingParameters(sys_diag_simple, beta, P)

    sampling_sys === nothing && (sampling_sys = sys_diag_simple)
    sp = SamplingParameters(sampling_sys, beta, P)

    samples = zeros(Float64, 3, num_samples)
    meter = Progress(num_samples, output=progress_output)
    for n in ProgressWrapper(1:num_samples, meter)
        samples[:, n] .= get_sample_pt(sys, pseudosp, sp)
    end

    simple = Analytical(sampling_sys, beta)
    normalization = simple.Z

    f_E(samples, samples_E, samples_Cv) =
        samples_E ./ samples
    f_Cv(samples, samples_E, samples_Cv) =
        (samples_Cv ./ samples - (samples_E ./ samples).^2) * beta^2

    Z = mean(samples[1, :]) * normalization
    Z_err = std(samples[1, :]) / sqrt(num_samples) * normalization
    E, E_err = jackknife(f_E, [samples[n, :] for n in 1:size(samples, 1)]...)
    Cv, Cv_err = jackknife(f_Cv, [samples[n, :] for n in 1:size(samples, 1)]...)

    SamplingPrimitiveThermodynamic(Z, Z_err, E, E_err, Cv, Cv_err)
end


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module SeqTests

using Test, OEISUtils, SeqUtils, Triangles
export SeqTest, GenerateAllTest, data_installed

const ShowTest = true

"""

Alias for is_oeis_installed.
"""
data_installed() = is_oeis_installed()

"""

Perform tests for an array of sequences of given type assuming the given offset.
"""
function SeqTest(seqarray, kind, offset = 0)
    if kind == 'V'
        return SeqVTest(seqarray, offset)
    end
    if kind == 'B'
        return SeqBTest(seqarray)
    end
    if kind == 'L'
        return SeqLTest(seqarray)
    end
    if kind == 'T'
        return SeqTTest(seqarray)
    end
    if kind == 'P'
        return SeqPTest(seqarray)
    end

    error("Test function not found!")
end

function SeqVTest(seq, offset = 0)
    name = SeqName(seq)
    O = oeis_local(name, 10)
    S = ZArray(10, seq, offset)
    if ShowTest
        println("V --> ", name)
        println(O)
        println(S)
    end
    @test all(S[1:10] .== O[1:10])
end

function SeqBTest(seqarray)
    for seq ∈ seqarray
        name = SeqName(seq)
        # the parameter is not 'length' but 'search bound'.
        O = oeis_local(name, 12)
        S = seq(300)
        if ShowTest
            println("B --> ", name)
            println(O)
            println(S)
        end
        @test all(S[0:11] .== O[0:11])
    end
end

function SeqLTest(seqarray)
    for seq ∈ seqarray
        name = SeqName(seq)
        O = oeis_local(name, 12)
        S = seq(12)
        if ShowTest
            println("L --> ", name)
            println(O)
            println(S)
        end
        @test all(S .== O)
    end
end

function SeqTTest(seqarray)
    for seq ∈ seqarray
        name = SeqName(seq)
        O = oeis_local(name, 21)
        S = seq(6)
        if ShowTest
            println("T --> ", name)
            println(O)
            println(S)
            # ShowAsΔ(O); ShowAsΔ(S)
        end
        @test all(S .== O)
    end
end

function SeqPTest(seqarray)
    for seq ∈ seqarray
        name = SeqName(seq)
        O = oeis_local(name, 28)
        S = seq(7)
        if ShowTest
            println("P --> ", name)
            #ShowAsΔ(O); ShowAsΔ(S)
        end
    end
end

"""

Test to generate sequences from a given list of sequences.
"""
function GenerateAllTest(A::Array{Function,1})
    err = String[]
    for a ∈ A
        stra = string(a)
        stra * " " |> print
        try
            if startswith(stra, 'L')
                a(10) |> println
            elseif startswith(stra, "TL")
                if applicable(a, 1)
                    for k ∈ 0:4 a(k) |> println end
                else
                    for k ∈ 0:4 a(0, k) |> println end
                    for k ∈ 0:4 a(2, k) |> println end
                end
            elseif startswith(stra, 'T')
                reduce(vcat, [row for row in a(5)]) |> println
            elseif startswith(stra, "is")
                [a(k) for k ∈ 0:12] |> println
                [k for k ∈ 0:12 if a(k)] |> println
            elseif startswith(stra, 'I')
                collect(a(6)) |> println
            elseif startswith(stra, 'F')
                collect(IterTools.takewhile(k -> k <= 20, a(21))) |> println
            elseif startswith(stra, 'P')
                for k ∈ 0:4 a(k) |> println end
            elseif startswith(stra, 'R')
                [a(Float64(k)) for k ∈ 0:4] |> println
            elseif startswith(stra, 'G')
                Coefficients(a, 6) |> println
            elseif startswith(stra, 'M')
                a(6) |> println
            elseif startswith(stra, 'V')
                if applicable(a, 1)
                    [a(k) for k ∈ 0:6] |> println
                else
                    [a(0, k) for k ∈ 0:6] |> println
                    [a(2, k) for k ∈ 0:6] |> println
                end
            else
                if applicable(a, 1)
                    [a(k) for k ∈ 0:6] |> println
                else
                    [a(0, k) for k ∈ 0:6] |> println
                    [a(2, k) for k ∈ 0:6] |> println
                end
            end
        catch
            "ERROR!" |> println
            push!(err, stra)
        end
    end
    numerr = length(err)
    if numerr > 0
        @warn(numerr, " errors found:")
        for e ∈ err println(e) end
    end
    return numerr
end

#START-TEST-########################################################

using Test

function test()
    @testset "GenerateAll" begin
        #global SEQUENCES
        #@test GenerateAllTest(SEQUENCES) == 0
        #GenerateAllTest(SEQUENCES)
    end
end

function demo()
end

function perf()
end

function main()
    #test()
    demo()
    perf()
end

main()

end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
"""
    struct ProfilePoint
        value::Float64               # x value of profile point
        loss::Float64                # y value of profile point (loss function at value)
        params::Array{Float64, 1}    # vector of optimal values of loss_func arguments
        ret::Symbol                  # return value from NLOpt.optimize()
        counter::Union{Int, Nothing} # number of loss_func() calls to calculate the value
    end

Structure storing one point from profile function.
ret values - :FORCED_STOP, :MAXEVAL_REACHED, :FTOL_REACHED

"""
struct ProfilePoint
    value::Float64
    loss::Float64
    params::Array{Float64, 1}
    ret::Symbol
    counter::Union{Int, Nothing}
end

"""
    struct EndPoint
        value::Float64                        # value of endpoint or nothing
        profilePoints::Array{ProfilePoint, 1} # vector of profile points
        status::Symbol                        # result of analysis
        direction::Symbol                     # :right or :left
        counter::Int                          # number of loss_func() calls to calculate the endpoint
        supreme::Union{Float64, Nothing}      # maximal value inside profile interval
    end

Structure storing one endpoint for confidence interval.

status values - :BORDER_FOUND_BY_SCAN_TOL, :BORDER_FOUND_LOSS_TOL,
 :SCAN_BOUND_REACHED, :MAX_ITER_STOP, :LOSS_ERROR_STOP
"""
struct EndPoint
    value::Union{Float64, Nothing}
    profilePoints::Array{ProfilePoint, 1}
    status::Symbol
    direction::Symbol
    counter::Int
    supreme::Union{Float64, Nothing}
end


text:
# file for utility functions
function flattenArray(A::AbstractArray{AbstractArray{T}, 1}) where T
# copies array-of-array A into a flattened version B

  n = length(A)
  num_entries = 0
  for i=1:n  # count number of entries total
    num_entries += length(A[i])
  end

  B = zeros(T, num_entries)

  B_index = 1
  for i=1:n
    for j = 1:length(A[i])
      B[B_index] = A[i][j]
      B_index += 1
    end
  end

  return B

end

function getElIndex(a::Array, b)

  for i = 1:length(a)
    if a[i] == b
      return i
    end
  end

  return 0
end

function getMinandMax(arr::AbstractArray{T}) where T
# need to check this function for type stability

min_entry = typemax(T)
max_entry = typemin(T)

@inbounds for i=1:length(arr)
  entry_i = arr[i]

  if entry_i < min_entry
    min_entry = entry_i
  end

  if entry_i > max_entry
    max_entry = entry_i
  end
end

return min_entry, max_entry

end

"""
  Count the number of instances of val in arr

  Inputs:
    val: the value to find
    arr: the array to search in

  Outputs:
    cnt: the number of instances of val in arr
"""
function countin(val, arr::AbstractArray)

  cnt = 0
  for i=1:length(arr)
    if arr[i] == val
      cnt += 1
    end
  end

  return cnt
end

function common(A::AbstractArray, B::AbstractArray)
# find common element in the two sets

  # get the larger array
  if length(A) > length(B)
    arr1 = A
    arr2 = B
  else
    arr1 = B
    arr2 = A
  end

  # the output array is no larger than the smaller array
  arr_ret = zeros(eltype(arr2), size(arr2))
  num_found = 0
  for i=1:length(arr1)  # loop over larger array
    # check every entry in arr2 for the current entry in arr1
    el_i = arr1[i]
    for j=1:length(arr2)
      if el_i == arr2[j]
	num_found += 1
	arr_ret[num_found] = el_i
	break
      end
    end
  end

    return arr_ret, num_found
end

"""
  Gets the first element in A that is also in B.  If none found, returns
  first element of A
"""
function first_common(A::AbstractArray, B::AbstractArray)

  for i=1:length(A)
    A_i = A[i]
    for j=1:length(B)
      if A_i == B[j]
        return A_i
      end
    end
  end

  throw(ErrorException("No common element found"))
  return A[1]
end

function calcNewNode(i, offset_pumi, offset_orient)
# this function calculates the new node index on the entity
# using offset_pumi, the offset
# that maps the SBP node index i to the pumi node index, and offset_orient,
# which maps the pumi node index to the pumi node index that takes into account
# the orientation of the mesh entity
# the returned node
  tmp = abs(offset_pumi - i)
  tmp2 = abs(offset_orient - tmp) - 1
  return tmp2
end

"""
  Calculate the volume of the mesh as the sum of the volume of the elements

  Methods are available for 2D and 3D.
"""
function calcVolume(mesh::PumiMesh2DG)

  area = 0.0
  for i=1:mesh.numEl
    r1x = mesh.vert_coords[1, 2, i] - mesh.vert_coords[1, 1, i]
    r1y = mesh.vert_coords[2, 2, i] - mesh.vert_coords[2, 1, i]

    r2x = mesh.vert_coords[1, 3, i] - mesh.vert_coords[1, 1, i]
    r2y = mesh.vert_coords[2, 3, i] - mesh.vert_coords[2, 1, i]

    area_i = abs(0.5*(r1x*r2y - r1y*r2x))

    area +=  area_i
  end

  area = MPI.Allreduce(area, MPI.SUM, mesh.comm)

  return area
end

function calcVolume(mesh::PumiMesh3DG)

  volume = 0.0
  for i=1:mesh.numEl
    r1x = mesh.vert_coords[1, 2, i] - mesh.vert_coords[1, 1, i]
    r1y = mesh.vert_coords[2, 2, i] - mesh.vert_coords[2, 1, i]
    r1z = mesh.vert_coords[3, 2, i] - mesh.vert_coords[3, 1, i]

    r2x = mesh.vert_coords[1, 3, i] - mesh.vert_coords[1, 1, i]
    r2y = mesh.vert_coords[2, 3, i] - mesh.vert_coords[2, 1, i]
    r2z = mesh.vert_coords[3, 3, i] - mesh.vert_coords[3, 1, i]

    r3x = mesh.vert_coords[1, 4, i] - mesh.vert_coords[1, 1, i]
    r3y = mesh.vert_coords[2, 4, i] - mesh.vert_coords[2, 1, i]
    r3z = mesh.vert_coords[3, 4, i] - mesh.vert_coords[3, 1, i]

    # cross product
    c1 = r2y*r3z - r2z*r3y
    c2 = -(r2x*r3z - r2z*r3x)
    c3 = r2x*r3y - r2y*r3x

    volume_i = (r1x*c1 + r1y*c2 + r1z*c3)/6

    # dot product
    volume += volume_i
  end

  volume = MPI.Allreduce(volume, MPI.SUM, mesh.comm)
  println(mesh.f, "volume = ", volume)

  return volume
end

"""
  Calculate the volume of the mesh by integrating the 1/|dxidx|.
"""
function calcVolumeIntegral(mesh::PumiMesh, sbp)

  volume = 0.0
  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      volume += sbp.w[j]*real(1./mesh.jac[j, i])
    end
  end

  return volume
end

function accumulateAtVerts(mesh::PumiMeshDG, u_volume::Abstract3DArray, u_verts::AbstractMatrix)
# takes u_volume which is n x numVertsPerElement x numEl and does an addition
# reduction to u_verts, which is n x numVerts

  @assert size(u_volume, 1) == size(u_verts, 1)

  n = size(u_volume, 1)
  down_verts = Array{Ptr{Void}}(12)

  fill!(u_verts, 0.0)
  for i=1:mesh.numEl
    el_i = mesh.elements[i]
    nverts = apf.getDownward(mesh.m_ptr, el_i, 0, down_verts)

    for j=1:nverts
      vert_j = down_verts[j]
      vertnum_j = apf.getNumberJ(mesh.vert_Nptr, vert_j, 0, 0) + 1

      for k=1:n
        u_verts[k, vertnum_j] +=  u_volume[k, j, i]
      end
    end
  end

  return nothing
end

"""
  This function adds a new boundary condition for the specified geometric
  entity, and removes that geometric entity from any other boundary conditions.

  **Inputs**

   * new_geo: geometric entity number for the BC to be added
   * name: name of new boundary condition

  **Inputs/Outputs**

   * opts: options dictionary

   This function uses the following keys:

    * numBC
    * BC*
"""
function updateBCs(opts::Dict, new_geo::Integer, name::AbstractString)

  numBC = opts["numBC"]
  for i=1:numBC
    bndry_geo = opts["BC$i"]
    idx = findfirst(bndry_geo, new_geo)
    if idx != 0
      new_arr = zeros(eltype(bndry_geo), length(bndry_geo)-1)
      new_arr[1:(idx-1)] = bndry_geo[1:(idx-1)]
      new_arr[idx:end] = bndry_geo[(idx+1):end]
      opts["BC$i"] = new_arr
    end
  end

  numBC = numBC + 1
  opts["numBC"]
  opts["numBC"] = numBC
  opts["BC$numBC"] = [new_geo]
  opts["BC$(numBC)_name"] = name

  return nothing
end

function checkMeshPtr(mesh::PumiMesh)

  if mesh.m_ptr == C_NULL
    error("mesh has already been finalized")
  end

  return nothing
end


text:
"""
    waitmove(dev::StepperSystem)

## Description
Wait for stepper complete movement.
## Arguments
- dev: object of StepperSystem type with the serial connection (dev.con)
"""
function waitmove(dev::StepperSystem)

    status = 0
    while status == 0
        status = bytesavailable(dev.con)
    end

    sleep(0.2)
end

"""
    stepper_zero!(dev::StepperSystem)

## Description
Makes current position the origin.
## Arguments
- dev: object of StepperSystem type
"""
function stepper_zero!(dev::StepperSystem)
    dev.pos = dev.pos * 0.0
end

"""
    getpos(dev::StepperSystem)

## Description
Function to get the current position. Returns a vector of floats.
## Arguments
- dev: object of StepperSystem type
"""
getpos(dev::StepperSystem) = dev.pos[1:end]

"""
    manhattan_msg(dev::StepperSystem, id, steps)

## Description
Function to pass one instruction per line
## Arguments
- dev: object of StepperSystem type
- id: stepper motor id
- steps: number of steps
"""
function manhattan_msg(dev::StepperSystem, id, steps)

    n =length(id)

    msg = []
    for j in 1:n
        msg_temp = id[j]*";"*string(steps[j])*";"
        push!(msg, msg_temp)
        if isopen(dev.con)
            write(dev.con, msg_temp)
            waitmove(dev)
        else
            throw(AssertionError("Please open the connection!"))
        end
    end

    return msg
end

"""
    oneline_msg(dev::StepperSystem, id, steps)

## Description
Function to pass all instructions in one line
## Arguments
- dev: object of StepperSystem type
- id: stepper motor id
- steps: number of steps
"""
function oneline_msg(dev::StepperSystem, id, steps)

    msg = prod(id.*";".*string.(steps).*";")
    if isopen(dev.con)
        write(dev.con, msg)
        waitmove(dev)
    else
        throw(AssertionError("Please open the connection!"))
    end

    return msg
end

"""
    test_msg(dev::StepperSystem, id, steps)

## Description
Function test message for a stepper motor system
## Arguments
- dev: object of StepperSystem type
- id: stepper motor id
- steps: number of steps
"""
function test_msg(dev::StepperSystem, id, steps)
    msg = prod(id.*";".*string.(steps).*";")
end

"""
    stepper_move!(dev::StepperSystem, new_coords; relat=true, order=1:length(dev.id), method="manhattan")

## Description
Move system
## Arguments
- dev: object of StepperSystem type
- new_coords: array of the new coordinates
- relat: if true the relative movement is performed
- order: trigger order
- method: there are three methods available. For "manhattan_msg" steps are written one at a time. For "oneline_msg" everything is passed right through. The "test_msg" method don't need any connections (dev.con = missing).
## Examples
```jldoctest
julia> dev = stepper_open(2);

julia> stepper_config!(dev, motorID=["x","y"]);

julia> stepper_move!(dev, [10.0, 20.0])
2-element Array{Any,1}:
 "x;10;"
 "y;20;"

 julia> getpos(dev)
 2-element Array{Float64,1}:
 10.0
 20.0

 julia> stepper_move!(dev, [10.0, 20.0], order=[2,1], method=StepperControl.oneline_msg)
 "y;20;x;10;"

 julia> getpos(dev)
 2-element Array{Float64,1}:
 20.0
 40.0
```
"""
function stepper_move!(dev::StepperSystem, new_coords; relat=true, order=1:length(dev.id), method=manhattan_msg)

    n = length(new_coords)

    new_coords = float(new_coords)
    stepper_id = dev.id[order]

#    if !relat
#        new_coords = deltacoord(dev, new_coords)
#    end
    steps = [dev.coord2step[i](new_coords[dev.depend[i]]) for i in 1:n]
    coords = [dev.step2coord[k](steps[dev.depend[k]]) for k in 1:n]

    if !relat
        steps_init = [dev.coord2step[i](getpos(dev)[dev.depend[i]]) for i in 1:n]
        coords_init = [dev.step2coord[k](steps_init[dev.depend[k]]) for k in 1:n]
        steps = steps - steps_init
        coords = coords - coords_init
    end

    steps_msg = steps[order]

    msg = method(dev, stepper_id, steps_msg)

    #coords = [dev.step2coord[k](steps[dev.depend[k]]) for k in 1:n]
    dev.pos[order] += coords[order]

    return msg
end


text:
function getGaussianKernel(s::Int, width::Float64)
    σ = width / 4.3 # full width at tenth maximum = 4.3 * σ
    k = 1./(σ .* sqrt(2*π)).*exp.(-0.5.*(   (0.5:s)   .-s/2).^2/(σ.^2))
    #k = g(0.5:s)
    k = k ./ sum(k)
    return k
end

function smoothArray(a::Array{Float64,1}, kernelSize::Int, width::Float64; kernel::Array{Float64,1}=Array{Float64,1}(0))
    if length(kernel)==0
        kernel = getGaussianKernel(kernelSize, width)
    end

    halfKernelSize = Int((kernelSize-1)/2)
    return conv(a, kernel)[halfKernelSize+1:end-halfKernelSize]
end

function smoothArray(a::Array{Float64,2}, kernelSize::Int, width::Float64; kernel::Array{Float64,1}=Array{Float64,1}(0))
    if length(kernel)==0
        kernel = getGaussianKernel(kernelSize, width)
    end
    #kernel=[1/6, 2/3, 1/6]
    a_smooth = zeros(size(a))
    for i = 1:size(a,1)
        a_smooth[i,:] = smoothArray(a[i,:], kernelSize, width, kernel=kernel)
    end
    for i = 1:size(a,2)
        a_smooth[:,i] = smoothArray(a_smooth[:,i], kernelSize, width, kernel=kernel)
    end

    return a_smooth
end

function smoothArray(a::scaledArray, kernelSize::Int, width::Float64; kernel::Array{Float64,1}=Array{Float64,1}(0))
    if length(kernel)==0
        kernel = getGaussianKernel(kernelSize, width)
    else
        x = reshape(a.data[1:Int(end/2)], a.dimensions)
        xs = smoothArray(x, kernelSize, width)[:]
        y = reshape(a.data[1+Int(end/2):end], a.dimensions)
        ys = smoothArray(y, kernelSize, width)[:]
    end
    return scaledArray(vcat(xs,ys), a.dimensions, a.voxelsize, a.shift)
end


text:
# For all types that pack their type into their own structure (e.g. PPE)
const TYPEKEY = "_type"

##==============================================================================
## Variable Packing and unpacking
##==============================================================================
function packVariable(dfg::G, v::DFGVariable)::Dict{String, Any} where G <: AbstractDFG
    props = Dict{String, Any}()
    props["label"] = string(v.label)
    props["timestamp"] = string(v.timestamp)
    props["nstime"] = string(v.nstime.value)
    props["tags"] = JSON2.write(v.tags)
    props["ppeDict"] = JSON2.write(v.ppeDict)
    props["solverDataDict"] = JSON2.write(Dict(keys(v.solverDataDict) .=> map(vnd -> packVariableNodeData(dfg, vnd), values(v.solverDataDict))))
    props["smallData"] = JSON2.write(v.smallData)
    props["solvable"] = v.solvable
    props["softtype"] = string(typeof(getSofttype(v)))
    props["bigData"] = JSON2.write(Dict(keys(v.bigData) .=> map(bde -> JSON2.write(bde), values(v.bigData))))
    props["bigDataElemType"] = JSON2.write(Dict(keys(v.bigData) .=> map(bde -> typeof(bde), values(v.bigData))))
    return props
end

function unpackVariable(dfg::G, packedProps::Dict{String, Any})::DFGVariable where G <: AbstractDFG
    label = Symbol(packedProps["label"])
    timestamp = DateTime(packedProps["timestamp"])
    nstime = Nanosecond(get(packedProps, "nstime", 0))
    tags =  JSON2.read(packedProps["tags"], Vector{Symbol})
    #TODO this will work for some time, but unpacking in an <: AbstractPointParametricEst would be lekker.
    ppeDict = JSON2.read(packedProps["ppeDict"], Dict{Symbol, MeanMaxPPE})
    smallData = JSON2.read(packedProps["smallData"], Dict{String, String})
    bigDataElemTypes = JSON2.read(packedProps["bigDataElemType"], Dict{Symbol, Symbol})
    bigDataIntermed = JSON2.read(packedProps["bigData"], Dict{Symbol, String})
    softtypeString = packedProps["softtype"]
    softtype = getTypeFromSerializationModule(dfg, Symbol(softtypeString))
    softtype == nothing && error("Cannot deserialize softtype '$softtypeString' in variable '$label'")

    packed = JSON2.read(packedProps["solverDataDict"], Dict{String, PackedVariableNodeData})
    solverData = Dict(Symbol.(keys(packed)) .=> map(p -> unpackVariableNodeData(dfg, p), values(packed)))

    # Rebuild DFGVariable using the first solver softtype in solverData
    variable = DFGVariable{softtype}(Symbol(packedProps["label"]), timestamp, nstime, Set(tags), ppeDict, solverData,  smallData, Dict{Symbol,AbstractBigDataEntry}(), Ref(packedProps["solvable"]))

    # Now rehydrate complete bigData type.
    for (k,bdeInter) in bigDataIntermed
        fullVal = JSON2.read(bdeInter, getfield(DistributedFactorGraphs, bigDataElemTypes[k]))
        variable.bigData[k] = fullVal
    end

    return variable
end

function packVariableNodeData(dfg::G, d::VariableNodeData)::PackedVariableNodeData where G <: AbstractDFG
  @debug "Dispatching conversion variable -> packed variable for type $(string(d.softtype))"
  return PackedVariableNodeData(d.val[:],size(d.val,1),
                                d.bw[:], size(d.bw,1),
                                d.BayesNetOutVertIDs,
                                d.dimIDs, d.dims, d.eliminated,
                                d.BayesNetVertID, d.separator,
                                d.softtype != nothing ? string(d.softtype) : nothing,
                                d.initialized,
                                d.inferdim,
                                d.ismargin,
                                d.dontmargin,
                                d.solveInProgress,
                                d.solvedCount)
end

function unpackVariableNodeData(dfg::G, d::PackedVariableNodeData)::VariableNodeData where G <: AbstractDFG
  r3 = d.dimval
  c3 = r3 > 0 ? floor(Int,length(d.vecval)/r3) : 0
  M3 = reshape(d.vecval,r3,c3)

  r4 = d.dimbw
  c4 = r4 > 0 ? floor(Int,length(d.vecbw)/r4) : 0
  M4 = reshape(d.vecbw,r4,c4)

  # TODO -- allow out of module type allocation (future feature, not currently in use)
  @debug "Dispatching conversion packed variable -> variable for type $(string(d.softtype))"
  # Figuring out the softtype
  unpackedTypeName = split(d.softtype, "(")[1]
  unpackedTypeName = split(unpackedTypeName, '.')[end]
  @debug "DECODING Softtype = $unpackedTypeName"
  st = getTypeFromSerializationModule(dfg, Symbol(unpackedTypeName))
  st == nothing && error("The variable doesn't seem to have a softtype. It needs to set up with an InferenceVariable from IIF. This will happen if you use DFG to add serialized variables directly and try use them. Please use IncrementalInference.addVariable().")

  return VariableNodeData{st}(M3,M4, d.BayesNetOutVertIDs,
    d.dimIDs, d.dims, d.eliminated, d.BayesNetVertID, d.separator,
    st(), d.initialized, d.inferdim, d.ismargin, d.dontmargin, d.solveInProgress, d.solvedCount)
end

##==============================================================================
## Factor Packing and unpacking
##==============================================================================

function packFactor(dfg::G, f::DFGFactor)::Dict{String, Any} where G <: AbstractDFG
    # Construct the properties to save
    props = Dict{String, Any}()
    props["label"] = string(f.label)
    props["timestamp"] = string(f.timestamp)
    props["nstime"] = string(f.nstime.value)
    props["tags"] = JSON2.write(f.tags)
    # Pack the node data
    fnctype = getSolverData(f).fnc.usrfnc!
    try
        packtype = getfield(_getmodule(fnctype), Symbol("Packed$(_getname(fnctype))"))
        packed = convert(PackedFunctionNodeData{packtype}, getSolverData(f))
        props["data"] = JSON2.write(packed)
    catch ex
        io = IOBuffer()
        showerror(io, ex, catch_backtrace())
        err = String(take!(io))
        msg = "Error while packing '$(f.label)' as '$fnctype', please check the unpacking/packing converters for this factor - \r\n$err"
        error(msg)
    end
    # Include the type
    props["fnctype"] = String(_getname(fnctype))
    props["_variableOrderSymbols"] = JSON2.write(f._variableOrderSymbols)
    props["solvable"] = getSolvable(f)

    return props
end

function decodePackedType(::Type{T}, packeddata::GenericFunctionNodeData{PT}) where {T<:FactorOperationalMemory, PT}
  # usrtyp = convert(FunctorInferenceType, packeddata.fnc)
  # Also look at parentmodule
  usrtyp = getfield(PT.name.module, Symbol(string(PT.name.name)[7:end]))
  fulltype = DFG.FunctionNodeData{T{usrtyp}}
  factordata = convert(fulltype, packeddata)
  return factordata
end

function unpackFactor(dfg::G, packedProps::Dict{String, Any})::DFGFactor where G <: AbstractDFG
    label = packedProps["label"]
    timestamp = DateTime(packedProps["timestamp"])
    nstime = Nanosecond(get(packedProps, "nstime", 0))
    tags = JSON2.read(packedProps["tags"], Vector{Symbol})

    data = packedProps["data"]
    datatype = packedProps["fnctype"]
    @debug "DECODING Softtype = '$(datatype)' for factor '$label'"
    packtype = getTypeFromSerializationModule(dfg, Symbol("Packed"*datatype))

    packed = nothing
    fullFactorData = nothing
    try
        packed = JSON2.read(data, GenericFunctionNodeData{packtype})
        decodeType = getFactorOperationalMemoryType(dfg)
        fullFactorData = decodePackedType(decodeType, packed)
    catch ex
        io = IOBuffer()
        showerror(io, ex, catch_backtrace())
        err = String(take!(io))
        msg = "Error while unpacking '$label' as '$datatype', please check the unpacking/packing converters for this factor - \r\n$err"
        error(msg)
    end

    # Include the type
    _variableOrderSymbols = JSON2.read(packedProps["_variableOrderSymbols"], Vector{Symbol})
    solvable = packedProps["solvable"]

    # Rebuild DFGFactor
    #TODO use constuctor to create factor
    factor = DFGFactor(Symbol(label),
                       timestamp,
                       nstime,
                       Set(tags),
                       fullFactorData,
                       solvable,
                       Tuple(_variableOrderSymbols))

    # Note, once inserted, you still need to call rebuildFactorMetadata!
    return factor
end

##==============================================================================
## Serialization
##==============================================================================

"""
  $(SIGNATURES)
Get a type from the serialization module inside DFG.
"""
function getTypeFromSerializationModule(dfg::G, moduleType::Symbol) where G <: AbstractDFG
    st = nothing
    try
        st = getfield(Main, Symbol(moduleType))
    catch ex
        @error "Unable to deserialize soft type $(d.softtype)"
        io = IOBuffer()
        showerror(io, ex, catch_backtrace())
        err = String(take!(io))
        @error err
    end
    return st
end

"""
$(SIGNATURES)
Pack a PPE into a Dict{String, Any}.
"""
function packPPE(dfg::G, ppe::P)::Dict{String, Any} where {G <: AbstractDFG, P <: AbstractPointParametricEst}
    packedPPE = JSON.parse(JSON.json(ppe)) #TODO: Maybe better way to serialize as dictionary?
    packedPPE[TYPEKEY] = string(typeof(ppe)) # Append the type
    return packedPPE
end

"""
$(SIGNATURES)
Unpack a Dict{String, Any} into a PPE.
"""
function unpackPPE(dfg::G, packedPPE::Dict{String, Any})::AbstractPointParametricEst where G <: AbstractDFG
    !haskey(packedPPE, TYPEKEY) && error("Cannot find type key '$TYPEKEY' in packed PPE data")
    type = pop!(packedPPE, TYPEKEY)
    (type == nothing || type == "") && error("Cannot deserialize PPE, type key is empty")
    ppe = Unmarshal.unmarshal(
            DistributedFactorGraphs.getTypeFromSerializationModule(dfg, Symbol(type)),
            packedPPE)
    return ppe
end


text:
module DragWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const AbstractDragForceModel = @jimport org.orekit.forces.drag.AbstractDragForceModel
const AbstractForceModel = @jimport org.orekit.forces.AbstractForceModel
const Atmosphere = @jimport org.orekit.models.earth.atmosphere.Atmosphere
const Class = @jimport java.lang.Class
const DragForce = @jimport org.orekit.forces.drag.DragForce
const DragSensitive = @jimport org.orekit.forces.drag.DragSensitive
const Field = @jimport org.hipparchus.Field
const FieldAbsoluteDate = @jimport org.orekit.time.FieldAbsoluteDate
const FieldRotation = @jimport org.hipparchus.geometry.euclidean.threed.FieldRotation
const FieldSpacecraftState = @jimport org.orekit.propagation.FieldSpacecraftState
const FieldTimeDerivativesEquations = @jimport org.orekit.propagation.numerical.FieldTimeDerivativesEquations
const FieldVector3D = @jimport org.hipparchus.geometry.euclidean.threed.FieldVector3D
const ForceModel = @jimport org.orekit.forces.ForceModel
const Frame = @jimport org.orekit.frames.Frame
const IsotropicDrag = @jimport org.orekit.forces.drag.IsotropicDrag
const JString = @jimport java.lang.String
const NavigableSet = @jimport java.util.NavigableSet
const Object = @jimport java.lang.Object
const ParameterDriver = @jimport org.orekit.utils.ParameterDriver
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const Rotation = @jimport org.hipparchus.geometry.euclidean.threed.Rotation
const SpacecraftState = @jimport org.orekit.propagation.SpacecraftState
const Stream = @jimport java.util.stream.Stream
const TimeDerivativesEquations = @jimport org.orekit.propagation.numerical.TimeDerivativesEquations
const TimeScale = @jimport org.orekit.time.TimeScale
const TimeSpanDragForce = @jimport org.orekit.forces.drag.TimeSpanDragForce
const TimeSpanMap = @jimport org.orekit.utils.TimeSpanMap
const TimeSpanMap_Span = @jimport org.orekit.utils.TimeSpanMap$Span
const Vector3D = @jimport org.hipparchus.geometry.euclidean.threed.Vector3D

export AbstractDragForceModel
export DragForce
export DragSensitive
export IsotropicDrag
export TimeSpanDragForce
export acceleration
export add_contribution
export add_drag_sensitive_valid_after
export add_drag_sensitive_valid_before
export depends_on_position_only
export drag_acceleration
export equals
export extract_drag_sensitive_range
export extract_parameters
export get_atmosphere
export get_class
export get_drag_parameters_drivers
export get_drag_sensitive
export get_drag_sensitive_span
export get_events_detectors
export get_field_events_detectors
export get_parameter_driver
export get_parameters
export get_parameters_drivers
export get_spacecraft
export get_transitions
export hash_code
export init
export is_supported
export notify
export notify_all
export to_string
export wait

include("abstract_drag_force_model.jl")
include("drag_force.jl")
include("drag_sensitive.jl")
include("isotropic_drag.jl")
include("time_span_drag_force.jl")

end


text:
using Pkg
Pkg.update()
Pkg.add(PackageSpec(url="https://github.com/plotly/DashBase.jl.git"))
Pkg.add(PackageSpec(url="https://github.com/plotly/dash-html-components.git", rev="master"))
Pkg.add(PackageSpec(url="https://github.com/plotly/dash-core-components.git", rev="master"))
Pkg.add(PackageSpec(url="https://github.com/plotly/dash-table.git", rev="master"))
Pkg.add(PackageSpec(url="https://github.com/plotly/Dash.jl.git", rev=ENV["CIRCLE_BRANCH"]))
Pkg.add("HTTP")
Pkg.build("Dash")
Pkg.build("DashHtmlComponents")
Pkg.build("DashCoreComponents")
Pkg.build("HTTP")
Pkg.test("Dash", coverage=true)
function precompile_pkgs()
     for pkg in collect(keys(Pkg.installed()))
        if !isdefined(Symbol(pkg), :Symbol) && pkg != "Compat.jl"
            @info("Importing $(pkg)...")
            try (@eval import $(Symbol(pkg))) catch end
        end
    end
end
precompile_pkgs()


text:
export add_negate, add_round_off

function add_negate(mgr::Ptr{Manager}, f::Ptr{Node})
    result = ccall((:Cudd_addNegate, libcudd),
        Ptr{Node}, (Ptr{Manager}, Ptr{Node}), mgr, f)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end

function add_round_off(mgr::Ptr{Manager}, f::Ptr{Node}, n::Int)
    result = ccall((:Cudd_addRoundOff, libcudd),
        Ptr{Node}, (Ptr{Manager}, Ptr{Node}, Cuint), mgr, f, n)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end


text:
using MonkeyLang
using Documenter

DocMeta.setdocmeta!(MonkeyLang, :DocTestSetup, :(using MonkeyLang); recursive = true)

makedocs(;
    modules = [MonkeyLang],
    authors = "Gabriel Wu <wuzihua@pku.edu.cn> and contributors",
    repo = "https://github.com/lucifer1004/MonkeyLang.jl/blob/{commit}{path}#{line}",
    sitename = "MonkeyLang.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://lucifer1004.github.io/MonkeyLang.jl",
        assets = String[],
    ),
    pages = ["Home" => "index.md"],
)

deploydocs(; repo = "github.com/lucifer1004/MonkeyLang.jl", devbranch = "main")


text: @time @testset "Properties.tune" begin include("unit_tune.jl") end


text:
using PowerModelsWildfire
using Test

import PowerModels
import Cbc
import JuMP
import Memento
import InfrastructureModels

Memento.setlevel!(Memento.getlogger(InfrastructureModels), "error")
PowerModels.logger_config!("error")

mip_optimizer = JuMP.optimizer_with_attributes(Cbc.Optimizer, "logLevel"=>0)

@testset "PowerModelsWildfire" begin

    @testset "OPS" begin
        @testset "test case5_risk_sys1 consistency" begin
            # sufficient gen on all load buses -> turn off all branches
            case = PowerModels.parse_file("./networks/case5_risk_sys1.m")
            result = PowerModelsWildfire.run_ops(case, PowerModels.DCPPowerModel, mip_optimizer);
            @test result["termination_status"] == OPTIMAL

            PowerModels.update_data!(case,result["solution"])

            @test isapprox(calc_total_risk(case), 24.6, atol=1e-4)
            @test isapprox(calc_load(case), 10.0, atol=1e-4)

            @test isapprox(case["branch"]["1"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["2"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["3"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["4"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["5"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["6"]["br_status"], 0, atol=1e-4)
        end

        @testset "test case5_risk_sys2 consistency" begin
            # standard case5 network
            case = PowerModels.parse_file("./networks/case5_risk_sys2.m")
            result = PowerModelsWildfire.run_ops(case, PowerModels.DCPPowerModel, mip_optimizer);
            @test result["termination_status"] == OPTIMAL

            PowerModels.update_data!(case,result["solution"])

            @test isapprox(calc_total_risk(case), 24.5, atol=1e-4)
            @test isapprox(calc_load(case), 7.0, atol=1e-4)

            @test isapprox(case["branch"]["1"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["2"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["3"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["4"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["5"]["br_status"], 1, atol=1e-4)
            @test isapprox(case["branch"]["6"]["br_status"], 0, atol=1e-4)
        end

        @testset "test case14 consistency" begin
            # case14
            case = PowerModels.parse_file("./networks/case14_risk.m")
            result = PowerModelsWildfire.run_ops(case, PowerModels.DCPPowerModel, mip_optimizer);
            @test result["termination_status"] == OPTIMAL

            PowerModels.update_data!(case,result["solution"])

            @test isapprox(calc_total_risk(case), 38.16263269, atol=1e-4)
            @test isapprox(calc_load(case), 0.369533, atol=1e-4)

            @test isapprox(case["branch"]["1"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["2"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["3"]["br_status"], 1, atol=1e-4)
            @test isapprox(case["branch"]["4"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["5"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["6"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["7"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["8"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["9"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["10"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["11"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["12"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["13"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["14"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["15"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["16"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["17"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["18"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["19"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["20"]["br_status"], 0, atol=1e-4)
        end
    end

    @testset "Heuristics" begin
        @testset "Voltage Heuristic" begin
            case = PowerModels.parse_file("./networks/case5_risk_sys2.m")
            result = PowerModelsWildfire.run_voltage_shutoff_heuristic(case, PowerModels.DCPPowerModel, mip_optimizer, risk_threshold=2.0);
            # @test result["termination_status"] == OPTIMAL  :: result is data dict, not a solution dict
            PowerModels.update_data!(case,result)

            @test isapprox(calc_total_risk(case), 20.45, atol=1e-4)
            @test isapprox(calc_load(case), 4.0, atol=1e-4)

            @test isapprox(case["branch"]["1"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["2"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["3"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["4"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["5"]["br_status"], 0, atol=1e-4)
            @test isapprox(case["branch"]["6"]["br_status"], 0, atol=1e-4)
        end
        @testset "Area Heuristic" begin
            case = PowerModels.parse_file("./networks/case14_risk.m")
            result = PowerModelsWildfire.run_area_shutoff_heuristic(case, PowerModels.DCPPowerModel, mip_optimizer, risk_threshold=20.0);
            # @test result["termination_status"] == OPTIMAL :: result is data dict, not a solution dict
            PowerModels.update_data!(case,result)

            @test isapprox(calc_load(case), 2.59, atol=1e-4)

            ## Results from the MLD portion of the problem are degenerate, and can have different solutions on different systems.
            # @test isapprox(calc_total_risk(case), 56.2, atol=1e-4)

            # @test isapprox(case["branch"]["1"]["br_status"], 1, atol=1e-4)
            # @test isapprox(case["branch"]["2"]["br_status"], 0, atol=1e-4)
            # @test isapprox(case["branch"]["3"]["br_status"], 1, atol=1e-4)
            # @test isapprox(case["branch"]["4"]["br_status"], 1, atol=1e-4)
            # @test isapprox(case["branch"]["5"]["br_status"], 0, atol=1e-4)
            # @test isapprox(case["branch"]["6"]["br_status"], 1, atol=1e-4)
            # @test isapprox(case["branch"]["7"]["br_status"], 1, atol=1e-4)
            # @test isapprox(case["branch"]["8"]["br_status"], 0, atol=1e-4)
            # @test isapprox(case["branch"]["9"]["br_status"], 0, atol=1e-4)
            # @test isapprox(case["branch"]["10"]["br_status"], 1, atol=1e-4)
            # @test isapprox(case["branch"]["11"]["br_status"], 0, atol=1e-4)
            # @test isapprox(case["branch"]["12"]["br_status"], 1, atol=1e-4)
            # @test isapprox(case["branch"]["13"]["br_status"], 1, atol=1e-4)
            # @test isapprox(case["branch"]["14"]["br_status"], 0, atol=1e-4)
            # @test isapprox(case["branch"]["15"]["br_status"], 1, atol=1e-4)
            # @test isapprox(case["branch"]["16"]["br_status"], 1, atol=1e-4)
            # @test isapprox(case["branch"]["17"]["br_status"], 1, atol=1e-4)
            # @test isapprox(case["branch"]["18"]["br_status"], 1, atol=1e-4)
            # @test isapprox(case["branch"]["19"]["br_status"], 0, atol=1e-4)
            # @test isapprox(case["branch"]["20"]["br_status"], 1, atol=1e-4)
        end
    end
end


text:
using Statistics: mean

struct ReachableFitnessLogger <: Logger
    original_organisms_alive::Vector{UInt64}
    direct_children_alive::Vector{UInt64}
    direct_grandchildren_alive::Vector{UInt64}

    parents::Dict
    children::Vector

    function ReachableFitnessLogger(snapshot)
        new(get_organism_ids(snapshot), [], [], Dict(), [])
    end
end

function get_reachable_fitness(snapshot, rel_tolerance, min_samples, max_samples)
    reachable_fitness = estimate(rel_tolerance, min_samples, max_samples) do
        is_beneficial = []

        while true
            logger = ReachableFitnessLogger(snapshot)
            simulate_snapshot!(should_terminate, snapshot, logger)

            if length(logger.children) == 0 continue end

            is_beneficial = []

            for child in logger.children
                parent = logger.parents[child]

                fitness_child = get_fitness(snapshot, child)
                fitness_parent = get_fitness(snapshot, parent)

                push!(is_beneficial, fitness_parent < fitness_child ? 1 : 0)
            end

            if length(is_beneficial) == 0 continue end

            break
        end

        return rand(is_beneficial, 5)
    end

    return reachable_fitness
end

function should_terminate(logger::ReachableFitnessLogger, snapshot)
    isempty(logger.original_organisms_alive) &&
    isempty(logger.direct_children_alive) &&
    isempty(logger.direct_grandchildren_alive)
end

function log_step(::ReachableFitnessLogger, model) end
function save_log(::ReachableFitnessLogger) end

function log_birth(logger::ReachableFitnessLogger, model, child, parents=nothing)
    if parents === nothing return end

    parent = rand(parents)

    child_genotype = get_genotype(model, child)
    parent_genotype = get_genotype(model, parent)

    if child_genotype == parent_genotype return end

    child_id = get_id(model, child)
    parent_id = get_id(model, parent)

    if parent_id in logger.original_organisms_alive
        push!(logger.direct_children_alive, child_id)
        push!(logger.children, child)
        logger.parents[child] = parent
    elseif parent_id in logger.direct_children_alive
        push!(logger.direct_grandchildren_alive, child_id)
    end
end

function log_death(logger::ReachableFitnessLogger, model, organism)
    id = get_id(model, organism)

    if id in logger.original_organisms_alive
        delete!(logger.original_organisms_alive, id)
    elseif id in logger.direct_children_alive
        delete!(logger.direct_children_alive, id)
    elseif id in logger.direct_grandchildren_alive
        delete!(logger.direct_grandchildren_alive, id)
    end
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
module TestDiffMap
	using ManifoldLearning
	using Base.Test

	d = 2
	t = 1
	X, L = swiss_roll()
	Y = transform(DiffMap, X; d=d, t=t)

	@test outdim(Y) == d
	@test size(projection(Y), 2) == size(X, 2)
end


text:
# To add a new cell, type '# %%'
# To add a new markdown cell, type '# %% [markdown]'
# %%
using HDF5
using Plots
using Plots.PlotMeasures
using GR
using QuadGK
using Printf
using LaTeXStrings
using DelimitedFiles
ENV["GKSwstype"] = "100"
gr()
# %%
H = 0.25 # height from galactic plane in kpc
rlim = 40 # limit of radius
kmsinkpcyr = 3.1536e7/3.085677581e16 # 1 km/sec in kpc/yr
# timestep = np.linspace(0.01,0.8,80)
time = 300 # snapshot
binsize = 0.01
rootdirectory = "/home/oku/SimulationData/isogal"

models = [
    # ["Osaka2019_isogal", "Osaka2019 (Shimizu+19)"],
        #   , "geodome_model/geodome_original"\
    # ["geodome_model/ver_19.11.1", "Geodesic dome model & Cioffi+ 1988"],
    # ["geodome_model/OKU2020","Geodesic dome model & Athena fitting"],
    # ["centroid_model/ver07271_NoMomCeiling","Centroid model & Athena fitting (alpha = 0)"],
    # ["centroid_model/ver07272_nfb1","Centroid model & Athena fitting (nfb = 1)"],
    # ["centroid_model/ver07272_SFE001","Centroid model & Athena fitting (SFE = 0.01)"],
    # ["centroid_model/ver07311","Centroid model & Athena fitting (alpha = -1)"],
    # ["centroid_model/ver07311_fdens-2","Centroid model & Athena fitting (alpha = -2)"],
    # ["centroid_model/ver08041_alpha-1","Centroid model & Athena fitting (new, alpha = -1)"],
    # ["centroid_model/ver08041_alpha-2","Centroid model & Athena fitting (new, alpha = -2)"],
    # ["centroid_model/ver08041_NoThermal","Centroid model (alpha = -1)"],
    # ["centroid_model/ver08191","Centroid model (alpha = -0.24)"],
    # ["centroid_model/ver07272_CHEVALIER1974","Centroid model & Cioffi+ 1988"],
    # ["centroid_model/ver09011_n1f2a0","Centroid model"],
    # ["centroid_model/ver09011_n1f2a024","Centroid model (new, alpha = -0.24)"],
    # ["centroid_model/ver09011_n1f2a050","Centroid model (new, alpha = -0.50)"],
    # ["centroid_model/ver09011_n1f2a075","Centroid model (new, alpha = -0.75)"],
    # ["centroid_model/ver09011_n1f2a1","Centroid model (new, alpha = -1.00)"],
    # ["centroid_model/ver12151","Centroid"],
    # ["ss_model/ver01092", "Spherical superbubble"],
    # ["ss_model/ver01121", "Spherical superbubble"],
    # ["ss_model/ver01193", "Spherical superbubble"],
    # ["ss_model/ver03221", "Spherical superbubble"],
    # ["ss_model/ver12211", "Spherical superbubble"],
    # ["ss_model/ver03311", "with cool-off"],
    ["ss_model/ver03311/fiducial", "fiducial"],
    # ["ss_model/ver03311/med", "high-reso"],
    # ["ss_model/ver03311_NoCoolOff", "without cool-off"],
    # ["NoFB", "No feedback"],
    ]
# %%
function getcenter(sph)
  x = y = z = mtot = 0
  for i in 1:length(sph.mass)
    x += sph.mass[i] * sph.pos[1, i]
    y += sph.mass[i] * sph.pos[2, i]
    z += sph.mass[i] * sph.pos[3, i]
    mtot += sph.mass[i]
  end
  x /= mtot
  y /= mtot
  z /= mtot
  return [x, y, z]
end

function main(galpos)
  RDF = zeros(Int(ceil(rlim/binsize)))
  for i in 1:length(star.rshock)
    xi = star.pos[1,i] - galpos[1]
    yi = star.pos[2,i] - galpos[2]
    zi = star.pos[3,i] - galpos[3]
    for j in i+1:length(star.rshock)
      xj = star.pos[1,j] - galpos[1]
      yj = star.pos[2,j] - galpos[2]
      zj = star.pos[3,j] - galpos[3]

      r2 = (xi - xj)^2 + (yi - yj)^2 + (zi - zj)^2
      r = sqrt(r2)
      index = Int(ceil(r/binsize))
      if index == 0
        @show r, r/binsize, Int(ceil(r/binsize)), i, j
      end
      RDF[index] += 1/(4*pi*r^2)
    end
  end
  norm = sum(RDF)
  RDF = RDF/norm/binsize
  @show RDF
  return RDF
end
struct Star
    pos::Array{Float32, 2}
    rshock::Array{Float32, 1}
    sft::Array{Float32, 1}
    sniis::Array{Float32, 1}
    sniie::Array{Float32, 1}
end
struct Sph
    pos::Array{Float32, 2}
    mass::Array{Float32, 1}
end
# %%
Volume = [[[0.0 for j in 1:time] for k in [1, 2]] for i in 1:length(models)]
for j in 1:time
  for i in 1:length(models)
        path2snapshot = @sprintf("%s/%s/snapshot_%03d/snapshot_%03d.hdf5", rootdirectory, models[i][1], j, j)
        h5open(path2snapshot, "r") do fp
            posstar = read(fp,"PartType4/Coordinates")
            possph = read(fp,"PartType0/Coordinates")
            mass = read(fp,"PartType0/Masses")
            rshock = read(fp, "PartType4/Rshock")
            age = read(fp, "PartType4/StellarFormationTime")
            sniis = read(fp, "PartType4/SNIIStartTime")
            sniie = read(fp, "PartType4/SNIIENDTime")
            global star = Star(posstar, rshock, age, sniis, sniie)
            global sph = Sph(possph, mass)
        end
        galpos = getcenter(sph)

        result = main(galpos)
        Plots.plot(result, xlim=(1e-1, 1e2), ylim=(1e-3, 1e2), xscale=:log10, yscale=:log10)
        Plots.savefig("tmp.png")
  end
  @show j
end

gr()
font = Plots.font(16)
lfont = Plots.font(12)
Plots.plot(legend=:bottomleft,xlim=(0, time), ylim=(1e-8, 1e-1), yscale=:log10, xlabel="Time [Gyr]", ylabel=L"f_{\rm hot}", xtickfont=font, ytickfont=font, xlabelfont=font, ylabelfont=font, legendfont=lfont, guidefont=font, topmargin=3mm, bottommargin=3mm, rightmargin=3mm)
for i in 1:length(models)
Plots.plot!(Volume[i][1], label="hot phase", markersize=5, markercolor=:red)
Plots.plot!(Volume[i][2], label="cold phase", markersize=5, markercolor=:blue)
end
Plots.savefig("Energyloading.png")
# Plots.savefig("results/Energyloading.pdf")


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
#-----------------------------------------------------------------------------
# load packages
@info "load packages"
using ArgParse
using CovidSurvey
using Random, Dates, Turing
using Base.Threads
using StatsPlots
using PrettyTables
using DataFrames
using BSON
using OrderedCollections

using Memoization
using ReverseDiff
Turing.setrdcache(true)
setadbackend(:reversediff)

plotlyjs()
@info "number of threads available: $(nthreads())"
#-----------------------------------------------------------------------------
# read input
argtable = ArgParseSettings(
    description="This script samples from Regional.model using NUTS."
)
@add_arg_table! argtable begin
    "--chunksize"
        help = "chunksize to be used by ForwardDiff.jl"
        arg_type = Int
        default = 40
    "--num-samples", "-n"
        help = "number of samples"
        arg_type = Int
        default = 10
    "--rwstep", "-r"
        help = "random walk step size. Defaults to 7"
        arg_type = Int
        default = 7
    "--num-chains", "-t"
        help = "number of chains sampled with multi-threading"
        arg_type = Int
        default = 1
    "--num-warmup", "-w"
        help = "number of samples to use for warmup/adaptation"
        arg_type = Int
        default = 10
    "--observations-end", "-o"
        help = "end of observation date as string yyyy-mm-dd. Use full data set by defaults"
        arg_type = String
        default = "2021-02-01"
    "--cases-start", "-s"
        help = "start date for case observations as string yyyy-mm-dd. case observations start on 2020-06-01 by default"
        arg_type = String
        default = "2020-06-01"
    "--add-covariates", "-c"
        help = "list of covariate names (default empty string = no covariates). Choose: C: contacts, R:"
        arg_type = String
        default = nothing
    "--model", "-m"
        help = "choose from: 'hospit', 'cases', 'deaths'. Defaul: 'hospit'"
        arg_type = String
        default = "hospit"
    "--standardize", "-z"
        help = "normalize predictors to unit standard deviation"
        arg_type = Bool
        default = false
    "--plot-results", "-p"
        arg_type = Bool
        default = false
    "--foldername", "-f"
        arg_type = String
        default = "$(today())"
    "--name-prefix", "-u"
        arg_type = String
        default = ""
    "--include-early-dynamics", "-i"
        arg_type = Bool
        default = false
    "--semiparametric", "-S"
        arg_type = Bool
        default = true
    "--seed"
        help = "random seed to use"
        arg_type = Int
        default = nothing
end
parsed_args = parse_args(ARGS, argtable)

#----------------------------------------------------------------------------
# load data
@info "load data"
data_params = (
      observationsend   = parsed_args["observations-end"]
    , predictors        = parsed_args["add-covariates"] |> CovidSurvey.parse_predictors
    # , hospitmodel       = Regional.HospitInit(obs_stop="2020-07-01")
    , casemodel         = Regional.CaseInit(obs_start=parsed_args["cases-start"])
    , seromodel         = Regional.SeroInit(delay=0, std=.5)
    , rwstep           = parsed_args["rwstep"]
    , epidemicstart    = 20
    , numimpute        = 6
    , include_early_dynamic = parsed_args["include-early-dynamics"]
    , link              = KLogistic(4.5)
    , invlink           = KLogit(4.5)
    , covariates_kwargs = Dict(
      :semiparametric => parsed_args["semiparametric"],
      :fname          => [projectdir("data/contacts/dk/", "averaged_contact_rates_region=$(r).csv" ) for r in 1:Regional.nregions],#[projectdir("data", "mean_contact_rates_region=$r.csv" ) for r in 1:Regional.nregions],
      :shift          => -1,
      :conditions     => :date => x->x>Date("2020-11-10"), # >= 1000 cases / day
      :datecol        => "dates",
      :aggregation    => parsed_args["rwstep"],
      :standardize    => parsed_args["standardize"],
      :normalize      => true,
      :mobility       => [projectdir("data/mobility/mobility_region=$(Regional.regions[i]).csv") for i in 1:Regional.nregions]
      )
    )

data = Regional.load_data(; data_params... )
turing_data = data.turing_data;
#----------------------------------------------------------------------------
# prepare model
# ps = OrderedDict( zip( keys(data_params), values(data_params) ) )
# ps[:regional] = true
# ps[:model]  = parsed_args["model"]
# ps[:steps]  = parsed_args["num-samples"]
# ps[:warmup] = parsed_args["num-warmup"]
# ps[:chains] = parsed_args["num-chains"]
# ps[:rwstep] = parsed_args["rwstep"]
# ps[:seed]   = parsed_args["seed"]
# ps[:prefix] = isempty(parsed_args["name-prefix"]) ? "" : parsed_args["name-prefix"]*"-"
# ps = NamedTuple(ps)

ps = (
        regional = true,
        model    = parsed_args["model"],
        preds    = parsed_args["add-covariates"],
        steps    = parsed_args["num-samples"],
        warmup   = parsed_args["num-warmup"],
        observ   = parsed_args["observations-end"],
        chains   = parsed_args["num-chains"],
        rwstep   = parsed_args["rwstep"],
        seed     = parsed_args["seed"],
        std      = parsed_args["standardize"],
        semipar  = parsed_args["semiparametric"],
        prefix   = isempty(parsed_args["name-prefix"]) ? "" : parsed_args["name-prefix"]*"-"
)
# cases    = parsed_args["cases-start"],
# predstart= parsed_args["covariates-start"]

const model = Regional.name2model[ps.model]
!isnothing(ps.seed) && Random.seed!(ps.seed);
@info ps
m = model(turing_data, vectorize=true)
m();
## ==========================================================================
# sampling
@time chain = let
    thinning = 10
    if ps.chains > 1
        sample(m, NUTS(ps.warmup, 0.99; max_depth=6), MCMCThreads(), ps.steps + ps.warmup, ps.chains; progress=true, thinning, save_state=true)
    else
        sample(m, NUTS(ps.warmup, 0.99; max_depth=6), ps.steps + ps.warmup; progress=true, thinning, save_state=true) #; max_depth=15
    end
end
## ==========================================================================
@info "save"
fdir = projectdir("reports", parsed_args["foldername"])
fname = Regional.save_results(fdir, ps, data, chain)

## ==========================================================================
@info "post-processing"
using InlineStrings
using WeakRefStrings
Regional.postprocessing(fname)

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
function chaos(index,curr_iter,max_iter,Value)
    x = zeros(max_iter + 1)
    x[1]=0.7

    G = zeros(max_iter)
    if index == 1
        # Chebyshev map
        for i=1:max_iter
            x[i+1]=cos(i*acos(x[i]))
            G[i]=((x[i]+1)*Value)/2
        end
    elseif index == 2
        # Circle map
        a=0.5
        b=0.2
        for i=1:max_iter
            x[i+1]=mod(x[i]+b-(a/(2*pi))*sin(2*pi*x[i]),1)
            G[i]=x[i]*Value
        end
    elseif index == 3
        # Gauss/mouse map
        for i=1:max_iter
            if x[i]==0
                x[i+1]=0
            else
                x[i+1]=mod(1/x[i],1)
            end
            G[i]=x[i]*Value
        end
    elseif index == 4
        # Iterative map
        a=0.7
        for i=1:max_iter
            x[i+1]=sin((a*pi)/x[i])
            G[i]=((x[i]+1)*Value)/2
        end
    elseif index == 5
        # Logistic map
        a=4
        for i=1:max_iter
            x[i+1]=a*x[i]*(1-x[i])
            G[i]=x[i]*Value
        end
    elseif index == 6
        # Piecewise map
        P=0.4
        for i=1:max_iter
            if x[i]>=0 && x[i]<P
                x[i+1]=x[i]/P
            end
            if x[i]>=P && x[i]<0.5
                x[i+1]=(x[i]-P)/(0.5-P)
            end
            if x[i]>=0.5 && x[i]<1-P
                x[i+1]=(1-P-x[i])/(0.5-P)
            end
            if x[i]>=1-P && x[i]<1
                x[i+1]=(1-x[i])/P
            end
            G[i]=x[i]*Value
        end

    elseif index == 7
        # Sine map
        for i=1:max_iter
            x[i+1] = sin(pi*x[i])
            G[i]=(x[i])*Value
        end
    elseif index == 8
        # Singer map
        u=1.07
        for i=1:max_iter
            x[i+1] = u*(7.86*x[i]-23.31*(x[i]^2)+28.75*(x[i]^3)-13.302875*(x[i]^4))
            G[i]=(x[i])*Value
        end
    elseif index == 9
        # Sinusoidal map
        for i=1:max_iter
            x[i+1] = 2.3*x[i]^2*sin(pi*x[i])
            G[i]=(x[i])*Value
        end

    elseif index == 10
        # Tent map
        x[1]=0.6
        for i=1:max_iter
            if x[i]<0.7
                x[i+1]=x[i]/0.7
            end
            if x[i]>=0.7
                x[i+1]=(10/3)*(1-x[i])
            end
            G[i]=(x[i])*Value
        end

    end
    return G[curr_iter]

end


text:
using Documenter
using Tomography

makedocs(
    sitename = "Tomography.jl",
    format = Documenter.HTML(),
    modules = [Tomography],
    pages = Any[
        "Home" => "index.md",
        "Phantom" => "phantom.md",
        "Radon Transform" => "radon.md"
        ]
)

deploydocs(
    repo = "github.com/KNU-MATH-AI/Tomography.jl.git",
    target = "build",
    make = nothing,
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    structmatches, structmatch, isstructmatch,
    substructmatches, substructmatch, issubstructmatch,
    querymatch, isquerymatch,
    fastquerymatches, isqueryidentical, issmartsgroupmatch

"""
    structmatches(mol1::UndirectedGraph, mol2::UndirectedGraph; kwargs...) -> Iterator

Generate molecular graph match mappings between `mol1` and `mol2`. If no match
found, return nothing.
"""
function structmatches(mol1::UndirectedGraph, mol2::UndirectedGraph;
        prefilter=true, kwargs...)
    !prefilter || fastidentityfilter(mol1, mol2) || return []
    afunc = atommatch(mol1, mol2)
    bfunc = bondmatch(mol1, mol2)
    return graphmatches(mol1, mol2,
        nodematcher=afunc, edgematcher=bfunc; kwargs...)
end

function structmatch(mol1, mol2; kwargs...)
    res = iterate(structmatches(mol1, mol2; kwargs...))
    return res === nothing ? nothing : res[1]
end

isstructmatch(mol1, mol2; kwargs...
    ) = structmatch(mol1, mol2; kwargs...) !== nothing

"""
    substructmatches(mol1::UndirectedGraph, mol2::UndirectedGraph; kwargs...) -> Iterator

Generate substructure match mappings between `mol1` and `mol2`. If no match
found, return nothing.

The mapping is based on only edge induced subgraph isomorphism and therefore
it does not care disconnected single atom matches. This function is intended
for use in substructure search. If you need single atom SMARTS match
(ex. [#16;X2;!R]), see [`querymatch`](@ref).
"""
function substructmatches(mol1::UndirectedGraph, mol2::UndirectedGraph;
        prefilter=true, kwargs...)
    !prefilter || fastsubstrfilter(mol1, mol2) || return []
    afunc = atommatch(mol1, mol2)
    bfunc = bondmatch(mol1, mol2)
    return edgesubgraphmatches(mol1, mol2,
        nodematcher=afunc, edgematcher=bfunc; kwargs...)
end

function substructmatch(mol1, mol2; kwargs...)
    res = iterate(substructmatches(mol1, mol2; kwargs...))
    return res === nothing ? nothing : res[1]
end

issubstructmatch(mol1, mol2; kwargs...
    ) = substructmatch(mol1, mol2; kwargs...) !== nothing

"""
    querymatch(mol::UndirectedGraph, query::QueryMol; kwargs...) -> Dict{Int,Int}

Generate substructure match mappings between `mol1` and `mol2`. If no match
found, return nothing.

This accepts also disconnected single atom but returns only the first match.
This function is intended for use in SMARTS query search

"""
function querymatch(mol::UndirectedGraph, query::QueryMol; kwargs...)
    mnodeset = nodeset(mol)
    qnodeset = nodeset(query)
    afunc = atommatch(mol, query)
    bfunc = bondmatch(mol, query)
    if edgecount(query) != 0
        # Edge induced subgraph mapping
        for emap in edgesubgraphmatches(
                mol, query, nodematcher=afunc, edgematcher=bfunc; kwargs...)
            # Isolated node mapping
            msub = edgesubgraph(mol, Set(keys(emap)))
            qsub = edgesubgraph(query, Set(values(emap)))
            miso = setdiff(mnodeset, nodeset(msub))
            qiso = setdiff(qnodeset, nodeset(qsub))
            mrem = nodesubgraph(mol, miso)
            qrem = nodesubgraph(query, qiso)
            nmap = maxcardmap(nodeset(mrem), nodeset(qrem), afunc)
            # TODO: connectivity filter
            if length(nmap) == nodecount(qrem)
                return (emap, nmap)
            end
        end
    elseif nodecount(query) != 0
        # Isolated nodes only
        nmap = maxcardmap(mnodeset, qnodeset, afunc)
        # TODO: connectivity filter
        if length(nmap) == nodecount(query)
            return (Dict(), nmap)
        end
    end
    return
end

"""
    isquerymatch(mol, query; kwargs...)

Return whether mol matches with the query.
"""
isquerymatch(mol, query; kwargs...
    ) = querymatch(mol, query; kwargs...) !== nothing

"""
    fastquerymatches(mol::UndirectedGraph, query::QueryMol; kwargs...
        ) -> Dict{Int,Int}

Generate query match mappings between `mol` and `query`. If no match
found, return nothing.

The `query` should not have any component level expression that means it should
not have any dots (`.`). This is intended for use in functional group detection.

"""
function fastquerymatches(mol::UndirectedGraph, query::QueryMol; kwargs...)
    isempty(query.connectivity) || throw(
        ErrorException("Component level query is disallowed"))
    afunc = atommatch(mol, query)
    bfunc = bondmatch(mol, query)
    nodecount(query) == 0 && return ()
    if nodecount(query) == 1
        # node match
        match = Iterators.filter(nodeset(mol)) do ma
            return afunc(ma, 1)
        end
        return ((Dict(), Dict(ma => 1)) for ma in match)
    elseif edgecount(query) == 1
        # edge match
        (qu, qv) = getedge(query, 1)
        match = Iterators.filter(edgeset(mol)) do mb
            (u, v) = getedge(mol, mb)
            return (
                ((afunc(u, qu) && afunc(v, qv))
                || (afunc(u, qv) && afunc(v, qu)))
                && bfunc(mb, 1)
            )
        end
        return ((Dict(mb => 1), Dict()) for mb in match)
    else
        # subgraph match
        return ((emap, Dict()) for emap in edgesubgraphmatches(
            mol, query, nodematcher=afunc, edgematcher=bfunc; kwargs...))
    end
end

function isqueryidentical(mol::UndirectedGraph, query::QueryMol; kwargs...)
    isempty(query.connectivity) || throw(
        ErrorException("Component level query is disallowed"))
    afunc = atommatch(mol, query)
    bfunc = bondmatch(mol, query)
    nodecount(query) == 0 && return ()
    if nodecount(query) == 1
        for n in nodeset(mol)
            afunc(n, 1) && return true
        end
        return false
    elseif edgecount(query) == 1
        (qu, qv) = getedge(query, 1)
        for e in edgeset(mol)
            (u, v) = getedge(mol, e)
            pred = ((afunc(u, qu) && afunc(v, qv))
                || (afunc(u, qv) && afunc(v, qu))) && bfunc(e, 1)
            pred && return true
        end
        return false
    else
        return isgraphmatch(
            mol, query, nodematcher=afunc, edgematcher=bfunc; kwargs...)
    end
end

function issmartsgroupmatch(
        mol::UndirectedGraph, query::QueryMol, root::Int; kwargs...)
    # For recursive SMARTS match
    isempty(query.connectivity) || throw(
        ErrorException("Component level query is disallowed"))
    afunc = atommatch(mol, query)
    bfunc = bondmatch(mol, query)
    @assert root in nodeset(mol) "node $(root) does not exist"
    if nodecount(query) == 1
        # node match
        return afunc(root, 1)
    elseif edgecount(query) == 1
        # edge match
        for (inc, adj) in neighbors(mol, root)
            if afunc(root, 1) && afunc(adj, 2) && bfunc(inc, 1)
                return true
            end
        end
        return false
    else
        # subgraph match
        for n in incidences(mol, root)
            if isedgesubgraphmatch(
                    mol, query, nodematcher=afunc, edgematcher=bfunc,
                    mandatory=Dict(n=>1); kwargs...)
                return true
            end
        end
        return false
    end
end

function atommatch(mol1::UndirectedGraph, mol2::UndirectedGraph)
    sym1 = atomsymbol(mol1)
    sym2 = atomsymbol(mol2)
    pi1 = pielectron(mol1)
    pi2 = pielectron(mol2)
    return (a1, a2) -> sym1[a1] == sym2[a2] && pi1[a1] == pi2[a2]
end

function atommatch(mol::UndirectedGraph, querymol::QueryMol)
    return (a, qa) -> querymatchtree(nodeattr(querymol, qa).query, mol, a)
end

function bondmatch(mol1::UndirectedGraph, mol2::UndirectedGraph)
    # TODO: need bond attribute matching?
    return (b1, b2) -> true
end

function bondmatch(mol::UndirectedGraph, querymol::QueryMol)
    return (b, qb) -> querymatchtree(edgeattr(querymol, qb).query, mol, b)
end

function querymatchtree(query::Pair, mol::UndirectedGraph, i::Int)
    if query.first == :any
        return true
    elseif query.first == :and
        return all(querymatchtree(q, mol, i) for q in query.second)
    elseif query.first == :or
        return any(querymatchtree(q, mol, i) for q in query.second)
    elseif query.first == :not
        return !querymatchtree(query.second, mol, i)
    elseif query.first == :stereo
        # TODO: stereo not implemented yet
        return true
    elseif query.first == :recursive
        subq = parse(SMARTS, query.second)
        return issmartsgroupmatch(mol, subq, i)
    else
        if query.first == :atom_sssrsizes
            return query.second in atom_sssrsizes(mol)[i]
        else
            return mol[query.first][i] == query.second
        end
    end
end

function fastidentityfilter(mol1::UndirectedGraph, mol2::UndirectedGraph)
    if nodecount(mol1) != nodecount(mol2)
        return false
    elseif edgecount(mol1) != edgecount(mol2)
        return false
    elseif circuitrank(mol1) != circuitrank(mol2)
        return false
    end
    return true
end

function fastsubstrfilter(mol1::UndirectedGraph, mol2::UndirectedGraph)
    if nodecount(mol1) < nodecount(mol2)
        return false
    elseif edgecount(mol1) < edgecount(mol2)
        return false
    elseif circuitrank(mol1) < circuitrank(mol2)
        return false
    end
    return true
end


text:
"""
    statistical_moments(results, n; field_apply=real)::Vector{Matrix}

Calculate moments up to `n` of results at each position and wavenumber/time, after applying `field_apply`.
"""
function statistical_moments(results::AbstractVector{SimRes}, num_moments::Int; field_apply=real) where {T,SimRes<:SimulationResult{T}}

    # Number of positions and wavenumbers/time points sampled
    X, K = size(results[1])

    # Number of realisations
    R = length(results)

    moments = [Matrix{T}(undef,X,K) for i=1:num_moments]

    # For each wavenumber or timestep, calculate each moment up to num_moments
    for i = 1:X
        for j = 1:K
            responses = [field_apply(field(results[r],i,j)) for r=1:R]
            μ = mean(responses)
            moments[1][i,j] = μ
            for m=2:num_moments
                moment = sum((responses .- μ).^m)/(R-1)
                moments[m][i,j] = sign(moment) * abs(moment)^(1.0/m)
            end
        end
    end
    return moments
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "wannier90"
version = v"3.1.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/wannier-developers/wannier90/archive/v3.1.0.tar.gz", "40651a9832eb93dec20a8360dd535262c261c34e13c41b6755fa6915c936b254")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/wannier90-*/

cat > make.inc << EOF
    F90    = $FC
    FCOPTS = $FFLAGS -fPIC -O3
    LDOPTS = $LDFLAGS -fPIC
    LIBS   = -lopenblas -L${libdir}
    DYNLIBRARYEXTENSION = $dlext
EOF

make -j${nproc} wannier dynlib

# Installation in the Makefile is broken
install -d ${bindir} ${libdir}
for exe in wannier90; do
    install -m755 ${exe}.x ${bindir}/${exe}${exeext}
done
for lib in libwannier.*; do
    install -m644 $lib ${libdir}/$lib
done
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_gfortran_versions(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libwannier", :libwannier),
    ExecutableProduct("wannier90", :wannier90)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="OpenBLAS32_jll", uuid="656ef2d0-ae68-5445-9ca0-591084a874a2")),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# Simulate an inhomogeneous Poisson point process on a rectangle
# This is done by simulating a homogeneous Poisson process, which is then
# thinned according to a (spatially *dependent*) p-thinning.
# The intensity function is lambda(x,y)=exp(-(x^2+y^2)/s^2), where s>0.
# Author: H. Paul Keeler, 2019.
# Website: hpaulkeeler.com
# Repository: github.com/hpaulkeeler/posts
# For more details, see the post:
# hpaulkeeler.com/testing-the-julia-language-with-point-process-simulations/

#Note: Need the .+ for adding a scalar to an array
#Also need . for sqrt, exp, cos, sin etc and assinging scalars to arrays
#Best to use vectors instead of 1-D matrices eg x=rand(n),  NOT x=rand(n,1).

using Distributions #for random simulations
using Plots #for plotting
using BlackBoxOptim #for blackbox optimizing

#Simulation window parameters
xMin=-1;xMax=1;
yMin=-1;yMax=1;
xDelta=xMax-xMin;yDelta=yMax-yMin; #rectangle dimensions
areaTotal=xDelta*yDelta;

s=0.5; #scale parameter

#Point process parameters
function fun_lambda(x,y)
    100*exp.(-(x.^2+y.^2)/s^2); #intensity function
end

###START -- find maximum lambda -- START ###
#For an intensity function lambda, given by function fun_lambda,
#finds the maximum of lambda in a rectangular region given by
#[xMin,xMax,yMin,yMax].
#NOTE: Need xMin, xMax, yMin, yMax to be floats eg xMax=1. See boundSearch

function fun_Neg(x)
     -fun_lambda(x[1],x[2]); #negative of lambda
end
xy0=[(xMin+xMax)/2.0,(yMin+yMax)/2.0];#initial value(ie centre)

#Find largest lambda value
boundSearch=[(1.0xMin,1.0xMax), (1.0yMin, 1.0yMax)];
#WARNING: Values of boundSearch cannot be integers!
resultsOpt=bboptimize(fun_Neg;SearchRange = boundSearch);
lambdaNegMin=best_fitness(resultsOpt); #retrieve minimum value found by bboptimize
lambdaMax=-lambdaNegMin;
###END -- find maximum lambda -- END ###

#define thinning probability function
function fun_p(x,y)
    fun_lambda(x,y)/lambdaMax;
end

#Simulate a Poisson point process
numbPoints=rand(Poisson(areaTotal*lambdaMax)); #Poisson number of points
xx=xDelta*rand(numbPoints).+xMin;#x coordinates of Poisson points
yy=yDelta*(rand(numbPoints)).+yMin;#y coordinates of Poisson points

#calculate spatially-dependent thinning probabilities
p=fun_p(xx,yy);
#Generate Bernoulli variables (ie coin flips) for thinning
booleRetained=rand(numbPoints,1).<p; #points to be retained
xxRetained=xx[booleRetained]; yyRetained=yy[booleRetained];

#Plotting
plot1=scatter(xxRetained,yyRetained,xlabel ="x",ylabel ="y", leg=false);
display(plot1);

max_stars_repo_name: hpkeeler/posts

max_issues_repo_name: hpkeeler/posts

max_forks_repo_name: hpkeeler/posts


text:
"""

```
init_velocities(sys::System{T},opt::Options) where T
```

Initialize velocities with random distribution.

"""
function init_velocities(sys::System{T},opt::Options) where T
  @unpack n = sys
  @unpack kT = opt
  # type of random distribution to be created
  if opt.initial_velocities == :normal
    v = randn(T,n)
  elseif opt.initial_velocities == :flat
    v = rand(T,n)
    @. v = random_rotate(v/norm(v))
    @. v = rand()*v
  elseif opt.initial_velocities == :zero
    v = zeros(T,n)
  else
    error("Invalid option for Options().initial_velocities.")
  end
  # Remove possible remaining drift
  remove_drift!(v)
  # Adjust average to desidred temperature  for i in 1:n
  kavg = kinetic(v)/n
  @. v = sqrt((dim(T)*kT/2)/kavg)*v
  return v
end

"""

```
function random_rotate(v)
```

Randomly rotate vector `v`, which might be a `Point2D` or a `Point3D`

"""
function random_rotate(v::Point2D)
  θ = 2π*rand()
  M = @SMatrix[ cos(θ) -sin(θ)
                sin(θ)  cos(θ) ]
  return M*v
end

function random_rotate(v::Point3D)
  θ = 2π*rand()
  ϕ = 2π*rand()
  ψ = 2π*rand()
  M = @SMatrix[ cos(ψ)*cos(θ)*cos(ϕ)-sin(ψ)*sin(ϕ)   -cos(ψ)*cos(θ)*sin(ϕ)-sin(ψ)*cos(ϕ)  cos(ψ)*sin(θ)
                sin(ψ)*cos(θ)*cos(ϕ)+cos(ψ)*sin(ϕ)   -sin(ψ)*cos(θ)*sin(ϕ)+cos(ψ)*cos(ϕ)  sin(ψ)*sin(θ)
                    -sin(θ)*cos(ϕ)                          sin(θ)*sin(ϕ)                    cos(θ)      ]
  return M*v
end


text:
using EmpiricalModeDecomposition
using Test

emd_dir = joinpath(dirname(pathof(EmpiricalModeDecomposition)), "..")

# Conduct Testing
include("$(emd_dir)/test/basic_emd.jl")
include("$(emd_dir)/test/basic_eemd.jl")


text: log_temperature(t::Real) = 1 / log(t)


text:
### MonitorTuner

# MonitorTuner is the most elemental tuner type
# It only monitors the acceptance rate of the current period

immutable MonitorTuner <: AbstractTuner
    period::Int # Tuning period over which acceptance rate is computed
    verbose::Bool

    function MonitorTuner(period::Int,verbose::Bool)
      @assert period > 0 "Tuning period should be positive"
      new(period,verbose)
    end
end

_tuner(::Type{Val{:monitor}},period::Int;verbose::Bool =true) = MonitorTuner(period,verbose)

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

type MonitorTunerState <: AbstractTunerState
    accepted::Vector{Int} # Number of accepted MCMC samples during current tuning period
    proposed::Vector{Int} # Number of proposed MCMC samples during current tuning period
    index::Int

    function MonitorTunerState(accepted::Vector{Int},proposed::Vector{Int},index::Int)
        @assert index > 0 "Total number of proposed MCMC samples should be non-negative"
        new(accepted,proposed,index)
    end
end

function _tunerstate{T<:AbstractFloat}(tuner::MonitorTuner,nburnin::Int,::Type{T})
    nsteps = ceil(Int,nburnin/tuner.period)
    MonitorTunerState(zeros(Int,nsteps),zeros(Int,nsteps),1)
end

show(io::IO,tuner::MonitorTuner) = println(io,"MonitorTuner: period = $(tuner.period), verbose = $(tuner.verbose)")

function show(io::IO,state::MonitorTunerState)
    println(io,"MonitorTunerState: ")
    println(io," accepted = $(accepted(state))")
    println(io," proposed = $(proposed(state))")
    println(io," acceptance rate = $(round(rate(state),3))")
    println(io," total proposed = $(total(state))")
end

function showstep(t::MonitorTuner,state::MonitorTunerState)
    if verbose(t) && state.index <= numtunesteps(state)
        a,p = current(state)
        println("  accepted/proposed = $(a)/$(p)")
        println("  acceptance rate = $(round(a/p,3))")
    end
end


text:
module SGP4

using Dates, Parameters, StaticArrays
using SatelliteToolbox.SatelliteToolboxTLE

################################################################################
#                                  Constants
################################################################################

# Earth Equatorial radius [m].
const R0 = 6378137.0

# Julian Day of J2000.0 epoch.
const JD_J2000 = 2451545.0

################################################################################
#                                    Types
################################################################################

include("types.jl")

################################################################################
#                                   Includes
################################################################################

include("gmst.jl")
include("sgp4_model.jl")
include("helpers.jl")

end


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

@doc raw"""DaemonSetUpdateStrategy is a struct used to control the update strategy for a DaemonSet.

    IoK8sApiAppsV1DaemonSetUpdateStrategy(;
        rollingUpdate=nothing,
        type=nothing,
    )

    - rollingUpdate::IoK8sApiAppsV1RollingUpdateDaemonSet : Rolling update config params. Present only if type = \"RollingUpdate\".
    - type::String : Type of daemon set update. Can be \"RollingUpdate\" or \"OnDelete\". Default is RollingUpdate.
"""
mutable struct IoK8sApiAppsV1DaemonSetUpdateStrategy <: SwaggerModel
    rollingUpdate::Any # spec type: Union{ Nothing, IoK8sApiAppsV1RollingUpdateDaemonSet } # spec name: rollingUpdate
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function IoK8sApiAppsV1DaemonSetUpdateStrategy(;rollingUpdate=nothing, type=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1DaemonSetUpdateStrategy, Symbol("rollingUpdate"), rollingUpdate)
        setfield!(o, Symbol("rollingUpdate"), rollingUpdate)
        validate_property(IoK8sApiAppsV1DaemonSetUpdateStrategy, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type IoK8sApiAppsV1DaemonSetUpdateStrategy

const _property_map_IoK8sApiAppsV1DaemonSetUpdateStrategy = Dict{Symbol,Symbol}(Symbol("rollingUpdate")=>Symbol("rollingUpdate"), Symbol("type")=>Symbol("type"))
const _property_types_IoK8sApiAppsV1DaemonSetUpdateStrategy = Dict{Symbol,String}(Symbol("rollingUpdate")=>"IoK8sApiAppsV1RollingUpdateDaemonSet", Symbol("type")=>"String")
Base.propertynames(::Type{ IoK8sApiAppsV1DaemonSetUpdateStrategy }) = collect(keys(_property_map_IoK8sApiAppsV1DaemonSetUpdateStrategy))
Swagger.property_type(::Type{ IoK8sApiAppsV1DaemonSetUpdateStrategy }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAppsV1DaemonSetUpdateStrategy[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1DaemonSetUpdateStrategy }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1DaemonSetUpdateStrategy[property_name]

function check_required(o::IoK8sApiAppsV1DaemonSetUpdateStrategy)
    true
end

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


text:
# Example of DandelionWebSocket.jl:
# Send some text and binary frames to ws://echo.websocket.org,
# which echoes them back.

using Requests: URI

using DandelionWebSockets

import DandelionWebSockets: on_text, on_binary,
                            state_connecting, state_open, state_closing, state_closed

mutable struct MockHandler <: WebSocketHandler
    client::WSClient
    stop_channel::Channel{Any}
end

# These are called when you get text/binary frames, respectively.
on_text(::MockHandler, s::String)  = println("Received text")
on_binary(::MockHandler, data::Vector{UInt8}) = println("Received data")

function state_closed(::MockHandler)
    println("State: CLOSED")

    # Signal the script that the connection is closed.
    put!(stop_chan, true)
end

stop_chan = Channel{Any}(1)

# Create a WSClient, which we can use to connect and send frames.
client = WSClient()
handler = MockHandler(client, stop_chan)

uri = URI("ws://localhost:$(ARGS[1])")
println("Connecting to $uri... ")

wsconnect(client, uri, handler)
println("Connected.")
# The second message on `stop_chan` indicates that the connection is closed, so we can exit.
take!(stop_chan)