text:
module TPT

export TPTSystem,

       # Basic information
       ncomp,
       composition,
       numberdensity,
       totalnumberdensity,
       temperature,

       # Structural properties
       structurefactor,
       paircorrelation,
       cavityfunction,
       nndistance,

       # Interatomic interaction
       pairpotential,
       pairpotential_minimum,
       pairpotential_minimizer,
       hsdiameter_estimate,

       # Thermodynamic properties
       kinetic,
       entropy,
       entropy_gas,
       entropy_conf,
       internal,
       helmholtz,

       # Reference systems
       AHS,
       WCA,

       # Hard-sphere system
       hsdiameter,
       packingfraction,
       totalpackingfraction,
       contactvalue,
       contactgradient,

       # Perturbation
       LennardJones,
       NFE,
       WHTB,
       NFETB,
       BOTB,

       # WCA
       blipfunction,

       # NFE
       Ashcroft,
       BretonnetSilbert,
       fermiwavenumber,
       formfactor,
       screenedformfactor,
       dielectric,
       localfiled,
       wnechar,

       # TB
       bandwidth,

       # utils
       spline,

       # constants
       kB

import Optim
import NLopt
using Dierckx
using Polynomials

include("utils.jl")
include("types.jl")
include("tptsystem.jl")
include("constants.jl")
include(joinpath("reference", "reference.jl"))
include(joinpath("reference", "abstractwca.jl"))
include(joinpath("reference", "ahs.jl"))
include(joinpath("reference", "wca.jl"))
include(joinpath("reference", "lwca.jl"))
include(joinpath("perturbation", "perturbation.jl"))
include(joinpath("perturbation", "nullpert.jl"))
include(joinpath("perturbation", "lennardjones.jl"))
include(joinpath("perturbation", "nfetb.jl"))
include(joinpath("perturbation", "nfe.jl"))
include(joinpath("perturbation", "ashcroft.jl"))
include(joinpath("perturbation", "bretonnet_silbert.jl"))
include(joinpath("perturbation", "harrison.jl"))
include(joinpath("perturbation", "hausleitner.jl"))

end # module


text:
"""
Discretize the acoustic wave equation in 1D using staggered SBP operators.
Weakly impose a homogeneous boundary condition on each boundary.
"""

using SparseArrays
using LinearAlgebra
using Test

import sbp

# Number of grid points
nx = 100

"""
Load staggered grid operators and grids
xp denotes the regular grid with equidistantly spaced grid points, and xm
denotes the cell-centered grid, including boundary points.
"""
ops = sbp.OP2019
desc = ops.load_description()
xp, xm, h, Dp, Dm, Hp, Hm, Pp, Pm = ops.build_operators(desc, nx)

"""

Write the discretization as dq/dt = A * q, q = [v, p], where

A = [ 0, Dp,
      Dm, 0],

v is stored on the xp grid, and p is stored on the xm grid.

"""
rows = [length(xp), length(xm)]
cols = [length(xp), length(xm)]

A = sbp.Sparse.block_matrix(rows, cols)
A12 = Dp
A21 = Dm
A = sbp.Sparse.block_matrix_insert(A, rows, rows, 1, 2, A12)
A = sbp.Sparse.block_matrix_insert(A, rows, rows, 2, 1, A21)

"""
Impose the boundary condition p = 0 on each boundary weakly in an energy
conserving manner.
"""
Bp = sbp.Staggered.boundary_matrix_p(nx)
Hpi = spdiagm(0 => 1 ./diag(Hp))
SAT = -Hpi*Bp

A = sbp.Sparse.block_matrix_insert(A, rows, rows, 1, 2, A12 + SAT)

"""
Check that the discretization is energy conservative to machine precision by
eigenvalue analysis

"""
lam = real(eigvals(Matrix(A)))
@test minimum(lam) > -1e-13
@test maximum(lam) < 1e-13

ext: jl

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
# this user supplied file creates a dictionary of arguments
# if a key is repeated, the last use of the key is used
# it is a little bit dangerous letting the user run arbitrary code
# as part of the solver
# now that this file is read inside a function, it is better encapsulated

arg_dict = Dict{Any, Any}(
# specify physcs, SBP operator, and solve type (steady Newton's method)
"physics" => "Euler",
"operator_type" => "SBPOmega",
"dimensions" => 2,
"run_type" => 5,
"jac_method" => 2,
"jac_type" => 2,
"order" => 1,

# setup discretization
"use_DG" => true,
"volume_integral_type" => 2,
"Volume_flux_name" => "IRFlux",
"face_integral_type" => 2,
"FaceElementIntegral_name" => "ESLFFaceIntegral",
"Flux_name" => "IRFlux",

# use implicit Euler globalization for Newton's method
"newton_globalize_euler" => true,
#"euler_tau" => 1e-1,  # the pseudo-time step can be specified here, but
#                      # the default is fine for this case
"real_time" => false,

# set up initial and boundary conditions
"IC_name" => "ICFreeStream",
"numBC" => 2,
"BC1" => [8],
"BC1_name" => "FreeStreamBC",
"BC2" => [5],
"BC2_name" => "noPenetrationBC",

# problem-specific parameters
"aoa" => 2.0,
"Ma" => 0.5,

# mesh and geometry
"smb_name" => "SRCMESHES/2D_airfoil.smb",
"dmg_name" => "SRCMESHES/2D_airfoil.dmg",

# solve options
"res_abstol" => 1e-10,
"res_reltol" => 1e-9,
"itermax" => 50,
"output_freq" => 1,
"do_postproc" => false,
)


text:
function testvisualizedata()
  vd = FastHenryHelper.VisualizeData(groupfortests())
  verified = [1000.0, 1.0, 0.01, 0.001, 1.0e-6, 0.0254, 2.54e-5]
  for i in eachindex(vd.nodedataarray)
    @test vd.nodedataarray[i].xyz[3] == verified[i]
  end
  @test vd.segmentdataarray[1].n1xyz[3] == verified[3]
  @test vd.segmentdataarray[1].n2xyz[3] == verified[4]
  @test vd.planedataarray[1].c1 ≈ [0.0,0.0,0.0]
  @test vd.planedataarray[1].c2 ≈ [0.000254,0.0,0.0]
  @test vd.planedataarray[1].c3 ≈ [0.000254,0.000254,0.0]
  @test vd.planedataarray[1].node_xyz ≈
    [2.54e-5 0.0001016;5.08e-5 0.000127;7.62e-5 0.0001524]
  FastHenryHelper.todisplayunit!(vd)
  for i in eachindex(vd.nodedataarray)
    @test vd.nodedataarray[i].xyz[3] == 1e-3*verified[i]
  end
  @test vd.segmentdataarray[1].n1xyz[3] == 1e-3*verified[3]
  @test vd.segmentdataarray[1].n2xyz[3] == 1e-3*verified[4]
  @test vd.planedataarray[1].c1 ≈ 1e-3 * [0.0,0.0,0.0]
  @test vd.planedataarray[1].c2 ≈ 1e-3 * [0.000254,0.0,0.0]
  @test vd.planedataarray[1].c3 ≈ 1e-3 * [0.000254,0.000254,0.0]
  @test vd.planedataarray[1].node_xyz ≈
    1e-3 * [2.54e-5 0.0001016;5.08e-5 0.000127;7.62e-5 0.0001524]
end
testvisualizedata()

function testmesh()
  inductor = coilcraft1010vsgroup("1010VS-111ME")
  vd = FastHenryHelper.VisualizeData(inductor)
  meshtest1 = mesh(inductor)
  meshtest2 = mesh(Node(1,1,1))
end
testmesh()


text:
push!(LOAD_PATH, "@stdlib")
import Pkg; Pkg.add("Conda"); using Conda

CC = joinpath(Conda.BINDIR, "gcc")
CXX = joinpath(Conda.BINDIR, "g++")
CMAKE = joinpath(Conda.BINDIR, "cmake")
MAKE = joinpath(Conda.BINDIR, "make")

@info "Install CONDA dependencies..."
pkgs = Conda._installed_packages()
if !("boost" in pkgs)
    Conda.add("boost", channel="anaconda")
end

SRC_DIR = "$(@__DIR__)/../src"
CUR_DIR = @__DIR__
# install dependencies
if !isdir("$SRC_DIR/amgcl")
    download("https://github.com/ddemidov/amgcl/archive/master.zip", "$SRC_DIR/amgcl.zip")
    run(`unzip -o $SRC_DIR/amgcl.zip -d $SRC_DIR`)
    run(`mv $SRC_DIR/amgcl-master $SRC_DIR/amgcl`)
    run(`rm $SRC_DIR/amgcl.zip`)
end

function compile_op(DIR)
    DIR = abspath(DIR)
    cd(DIR)
    if !isdir("build")
        mkdir("build")
    end
    flag = false
    for file in readdir("build")
        if Sys.islinux() && endswith(file, ".so")
            flag = true
            break
        end
        if Sys.isapple() && endswith(file, ".dylib")
            flag = true
            break
        end
    end
    if flag
        @info "Library exists"
        return
    end
    cd("build")
    run(`$CMAKE ..`)
    run(`$MAKE -j`)
    cd(CUR_DIR)
end

try
    run(`nvcc --version`)
    global NVCC = true
catch
    @warn("`nvcc` is not found. The FWI module of `FwiFlow` only has a GPU kernel.")
    global NVCC = false
end

if NVCC
    try
        DIR = joinpath(SRC_DIR, "CppOps/FWI/Src")
        compile_op(DIR)
    catch
        @warn("CppOps/FWI/Src Failed.")
        global NVCC = false
    end
end

if NVCC
    try
        DIR = joinpath(SRC_DIR, "CppOps/FWI")
        compile_op(DIR)
    catch
        @warn("CppOps/FWI Failed.")
        global NVCC = false
    end
end

for name in ["Poisson", "Laplacian", "Upwlap", "Upwps", "Saturation"]
    DIR = joinpath(SRC_DIR, "CppOps/$name")
    compile_op(DIR)
end


text:
using Plots

using simulator.materials.spesifics
using simulator.fresnelltools
using simulator.matrixcore

λ = 1200e-9

n1 = 1.5
nm = Pd_Werner(λ)
n2 = 1.2

θs = 0.7:1e-4:1

function point(θ, n2, which)
    i_1_2_te, i_1_2_tm, θm = FresnellBoundrary(n1, nm, θ)
    bulk = FresnellSlab(nm, 2*π/λ, 50e-9, θm)
    i_2_3_te, i_2_3_tm, _ = FresnellBoundrary(nm, n2, θm)

    system_te = CascadeScattering([i_1_2_te, bulk, i_2_3_te])
    system_tm = CascadeScattering([i_1_2_tm, bulk, i_2_3_tm])

    te_r = abs(system_te[2, 1])^2
    tm_r = abs(system_tm[2, 1])^2

    [te_r, tm_r][which]
end

anim = @animate for n ∈ 1:0.001:1.2
    plot(θs, point.(θs, n, 1), ylim = [0.5, 1], title="Plasmon reflectance in varying index environment", legend=:bottomright, label="TE")
    plot!(θs, point.(θs, n, 2), label="TM")
end

mp4(anim, "~/Downloads/temp.mp4", fps=30)


text:
export
    default_milp_optimizer,
    set_default_milp_optimizer!,
    default_optimizer_attributes,
    set_default_optimizer_attributes!,
    clear_default_optimizer_attributes!

global MILP_OPTIMIZER = nothing
global DEFAULT_OPTIMIZER_ATTRIBUTES = Dict{Union{String,MOI.AbstractOptimizerAttribute},Any}()
"""
    default_milp_optimizer()

Returns the black box optimizer to be use when formulating JuMP models.
"""
default_milp_optimizer() = MILP_OPTIMIZER
"""
    set_default_milp_optimizer!(optimizer)

Set the black box optimizer to be use when formulating JuMP models.
"""
function set_default_milp_optimizer!(optimizer)
    global MILP_OPTIMIZER = optimizer
end

"""
    default_optimizer_attributes()

Return a dictionary of default optimizer attributes.
"""
default_optimizer_attributes() = DEFAULT_OPTIMIZER_ATTRIBUTES

"""
    set_default_optimizer_attributes!(vals)

Set default optimizer attributes.
e.g. `set_default_optimizer_attributes!(Dict("PreSolve"=>-1))`
"""
function set_default_optimizer_attributes!(pair::Pair,pairs...)
    push!(DEFAULT_OPTIMIZER_ATTRIBUTES,pair)
    set_default_optimizer_attributes!(pairs...)
end
set_default_optimizer_attributes!(d::Dict) = set_default_optimizer_attributes!(d...)
set_default_optimizer_attributes!() = nothing

"""
    clear_default_optimizer_attributes!()

Clear the default optimizer attributes.
"""
function clear_default_optimizer_attributes!()
    empty!(DEFAULT_OPTIMIZER_ATTRIBUTES)
end


text:
using Revise
using LinearAlgebra,Random, BenchmarkTools, Zygote
Random.seed!(1111)
#TODO- adapt to the current problem under considration

#value of the function-used for testing purposes

f_Xa2(X) = ((a) -> (x->(x-a)^2).(X))
D_Xa2(X) = ((a)->2a .- 2X)

#projection matrix that is a function of a parameter
function f_PXay(f_Xa, X, y)
  return (
  function PXay(a)
    Xa = f_Xa(a)
    return Xa*(cholesky(Xa'*Xa)\(Xa'*y))
  end)
end

function D_PXayold(f_Xa, D_Xa, X, y)
  return (
  function DPXay(a)
    Xa = f_Xa(a)
    DXa = D_Xa(a)
    XatXa = Xa' * Xa
    XatXainv = cholesky(XatXa)\I
    return (DXa*XatXainv*(Xa'*y) -
      Xa*XatXainv*(DXa'*Xa+Xa'*DXa)*XatXainv * Xa'*y +
      Xa*XatXainv*(DXa'*y))
  end)
end

function D_PXay(f_Xa, D_Xa, X, y)
  return (
  function DPXay(a, X::TM = f_Xa(a), DaX::TM = D_Xa(a)) where TM

    K= size(X,2)
    cuI = TM(I, K, K) #need an identity matrix
    XtXinv = cholesky!(X'*X)\cuI
    XXtXinv = X*XtXinv
    XtXinvXtv = XtXinv * (X'*y)

    return DaX*XtXinvXtv -
      XXtXinv*(DaX'*X+X'*DaX)*XtXinvXtv +
      XXtXinv*(DaX'*y)
  end)
end

function testD(f,D,args...; testname="", erroronfail::Bool = true)
  testval = D(args...) |> sum

  ∇(args) = gradient((args)->sum(f(args...)), args)
  verval = ∇(args)[1][1]

  local ε::Float64
  try
    ε =sum(abs.(testval .- verval))
  catch err
    msg = "Test $testname: failed probably due to dimensional mismatch: $err
      #args: $args
      testval: $testval
      verval: $verval"
    error(msg)
  end
  if args==verval
    @assert εrms==0
    @info "Test $testname: passed with a difference of 0"
  elseif testval≈verval
    @info "Test $testname: passed with absolute error of $(ε)"
  else
    msg = "Test $testname: failed with absolute error of $(ε)
      #args: $args
      testval: $testval
      verval: $verval"

    erroronfail && error(msg)
    @warn msg
  end

  return nothing
end

function testderivs(N=3,K=2)
  X= rand(N,K)
  y=rand(N) + X*rand(K)
  h = 10^-4
  a = 10.0

  Xa2 = f_Xa2(X)
  DXa2 = D_Xa2(X)
  Pa2 = f_PXay(Xa2, X, y)
  DPa2 = D_PXay(Xa2, DXa2, X, y)

  sPa2(a) = sum(Pa2(a))
  fd = (sPa2(a+h) - sPa2(a))/h
  @info "fd: $fd"

  testD(Pa2, DPa2, a, testname="Pay")

end

@time testderivs(10^4,10)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module custom_module
  using LinearAlgebra

  function squareMeBaby(A)
    #  A = unsafe_wrap(Array, data, len)
    ## Square array and return the result
    println(typeof(A))

    #  A[:]=[i*i for i in A]
    #  return A
    B = A.*A
    return B
    #  return map(x -> x*x, A)
  end

  function mutateMeByTen!(A)
    ## Multiple array in place by ten
    lmul!(10, A)
  end

  export dummy
  export squareMeBaby
  export mutateMeByTen!

  function dummy()
    println("Julia says... Hello, world ! Function dummy() from module custom_module has been executed !")
  end

end
## TODO Fix this
const julia_dummy = @cfunction(custom_module.dummy, (Cvoid), ())

function addMeBaby(x, y)
  return x+y
end
const julia_addMeBabyInt = @cfunction(addMeBaby, Cint, (Cint, Cint,))
#  const julia_addMeBabyFloat = @cfunction(addMeBaby, Cfloat, (Cfloat, Cfloat,))

max_stars_repo_name: Clonkk/nimjl

max_issues_repo_name: Clonkk/nimjl

max_forks_repo_name: Clonkk/nimjl


text:
using ParameterInference
using LowRankModels

function rprincipalcomponents(A::AbstractArray, k::Int, scale::Float64=1.0; abs_tol=1e-4^3, rel_tol=1e-4^3, max_iter=10000, kwargs...)
	loss = L1Loss()
	r = ZeroReg()
    # # Start with PCA
    m = GLRM(A, QuadLoss(), ZeroReg(), ZeroReg(), k; kwargs...)
    fit!(m)
    # Then do the other thing
    m = GLRM(A, loss, r, r, k; kwargs...)
    fit!(m, ProxGradParams(;abs_tol, rel_tol, max_iter))
	return m
end

function nanprincipalcomponents(A::AbstractArray, k::Int, scale::Float64=1.0; abs_tol=1e-4^3, rel_tol=1e-4^3, max_iter=10000, kwargs...)
    A = SparseArrays.fkeep!(A, (i,j,x) -> abs(x) < 5) # ! better
	loss = QuadLoss()
	r = ZeroReg()
    m = GLRM(A, loss, r, r, k; kwargs...)
    fit!(m, ProxGradParams(;abs_tol, rel_tol, max_iter))
	return m
end


text:
module NNHelferlein

import Base: iterate, length, summary
using LinearAlgebra
using Unicode
import ZipFile
import HDF5
import JLD2
using Statistics: mean, std
using ProgressMeter, Dates
using IterTools: ncycle, takenth
import CSV
# import OdsIO    # removed because of PyCall incompatibilities!
import DataFrames
import Random
using Printf
import CUDA
using Knet #: KnetArray, Param, @diff
import Images, Colors
import Augmentor
import MLDataUtils
using TensorBoardLogger, Logging

include("types.jl")
include("util.jl")
include("nets.jl")
include("layers.jl")
include("attn.jl")
include("funs.jl")
include("transformers.jl")
include("images.jl")
include("dataframes.jl")
include("texts.jl")
include("train.jl")
include("imagenet.jl")

DATA_DIR = normpath(joinpath(dirname(pathof(@__MODULE__)), "..", "data"))

export DNN, Classifier, Regressor, Chain,          # chains
       ImageLoader, DataLoader, preproc_imagenet,
       get_class_labels,
       iterate, length,
       Layer, Dense, Conv, Pool, Flat, PyFlat,         # layers
       DeConv, UnPool,
       Embed,
       RSeqClassifier, RSeqTagger,
       hidden_states, cell_states,
       Softmax, Dropout, BatchNorm, LayerNorm,
       Linear,
       AttentionMechanism, AttnBahdanau,
       AttnLuong, AttnDot, AttnLocation,
       AttnInFeed,
       leaky_sigm, leaky_relu, leaky_tanh,
       PositionalEncoding, positional_encoding_sincos,  # transformers
       mk_padding_mask, mk_peek_ahead_mask,
       dot_prod_attn,
       dataframe_read, dataframe_split,         # import data
       dataframe_minibatches, mk_class_ids,
       mk_image_minibatch,
       tb_train!,
       predict_top5, predict_imagenet,
       predict, hamming_dist,
       get_imagenet_classes,                    # images
       image2array, array2image, array2RGB,
       WordTokenizer,                           # texts
       get_tatoeba_corpus,
       seq_minibatch, seq2seq_minibatch,
       crop_array, init0, convert2KnetArray,             # utils
       blowup_array, recycle_array,
       de_embed,
       print_network,
       DATA_DIR

end # module


text:
"""
    UniversalFunctions

Universal stability and stability correction
functions for `SurfaceFluxes` module.
"""
module UniversalFunctions

using DocStringExtensions
using CLIMAParameters
using CLIMAParameters.SurfaceFluxes.UniversalFunctions

# Transport Types
export MomentumTransport, HeatTransport

# Types
export AbstractUniversalFunction, Gryanik, Grachev, Businger

# Methods
export phi, psi, Psi

const FTypes = Union{Real, AbstractArray}

abstract type AbstractUniversalFunction{FT <: FTypes} end
const AUF = AbstractUniversalFunction
const APS = AbstractParameterSet

Base.eltype(uf::AbstractUniversalFunction{FT}) where {FT} = FT

#####
##### Interface
#####

abstract type AbstractTransportType end
struct MomentumTransport <: AbstractTransportType end
struct HeatTransport <: AbstractTransportType end

Base.broadcastable(tt::AbstractUniversalFunction) = Ref(tt)
Base.broadcastable(tt::AbstractTransportType) = Ref(tt)

"""
    phi

Universal stability function for wind shear
(ϕ_m) and temperature gradient (ϕ_h)
"""
function phi end

"""
    psi

Universal stability correction function for
momentum (ψ_m) and heat (ψ_h)
"""
function psi end

"""
    Psi

Integral of universal stability correction
function for momentum (ψ_m) and heat (ψ_h)
"""
function Psi end

#####
##### Forwarding methods for free parameters
#####

Pr_0(uf::AUF) = Pr_0(uf.param_set, uf)
a_m(uf::AUF) = a_m(uf.param_set, uf)
a_h(uf::AUF) = a_h(uf.param_set, uf)
b_m(uf::AUF) = b_m(uf.param_set, uf)
b_h(uf::AUF) = b_h(uf.param_set, uf)
c_h(uf::AUF) = c_h(uf.param_set, uf)

π_group(uf::AUF, ::HeatTransport) = Pr_0(uf)
π_group(::AUF, ::MomentumTransport) = 1

#####
##### Businger
#####

"""
    Businger

# Reference

 - [Nishizawa2018](@cite)

# Original research

 - [Businger1971](@cite)

# Equations in reference:

    ϕ_m: Eq. A1
    ϕ_h: Eq. A2
    ψ_m: Eq. A3
    ψ_h: Eq. A4

# Fields

$(DocStringExtensions.FIELDS)
"""
struct Businger{FT, PS} <: AbstractUniversalFunction{FT}
    "Parameter set"
    param_set::PS
    "Monin-Obhukov Length"
    L::FT
    Businger(param_set::PS, L::FT) where {FT, PS} = new{FT, PS}(param_set, L)
end

# CLIMAParameters wrapper
Pr_0(param_set::APS, ::Businger) = Pr_0_Businger(param_set)
a_m(param_set::APS, ::Businger) = a_m_Businger(param_set)
a_h(param_set::APS, ::Businger) = a_h_Businger(param_set)

f_momentum(uf::Businger, ζ) = sqrt(sqrt(1 - 15 * ζ))

f_heat(uf::Businger, ζ) = sqrt(1 - 9 * ζ)

function phi(uf::Businger, ζ, ::MomentumTransport)
    if ζ < 0
        f_m = f_momentum(uf, ζ)
        return 1 / f_m
    else
        FT = eltype(uf)
        _a_m = FT(a_m(uf))
        return _a_m * ζ + 1
    end
end

function phi(uf::Businger, ζ, tt::HeatTransport)
    if ζ < 0
        f_h = f_heat(uf, ζ)
        return 1 / f_h
    else
        FT = eltype(uf)
        _a_h = FT(a_h(uf))
        _π_group = FT(π_group(uf, tt))
        return _a_h * ζ / _π_group + 1
    end
end

function psi(uf::Businger, ζ, ::MomentumTransport)
    FT = eltype(uf.L)
    if ζ < 0
        f_m = f_momentum(uf, ζ)
        log_term = log((1 + f_m)^2 * (1 + f_m^2) / 8)
        return log_term - 2 * atan(f_m) + FT(π) / 2
    else
        _a_m = FT(a_m(uf))
        return -_a_m * ζ
    end
end

function psi(uf::Businger, ζ, tt::HeatTransport)
    if ζ < 0
        f_h = f_heat(uf, ζ)
        return 2 * log((1 + f_h) / 2)
    else
        FT = eltype(uf)
        _a_h = FT(a_h(uf))
        _π_group = FT(π_group(uf, tt))
        return -_a_h * ζ / _π_group
    end
end

function Psi(uf::Businger, ζ, tt::MomentumTransport)
    FT = eltype(uf)
    if abs(ζ) < eps(FT)
        # Psi_m in Eq. A13
        if ζ >= 0
            _a_m = FT(a_m(uf))
            return -_a_m * ζ / 2
        else
            return -FT(15) * ζ / FT(8)
        end
    else
        # Note that "1-f^3" in is a typo, it is
        # supposed to be "1-f_m^3". This was
        # confirmed by communication with the
        # author.
        if ζ < 0
            f_m = f_momentum(uf, ζ)
            log_term = log((1 + f_m)^2 * (1 + f_m^2) / 8)
            π_term = FT(π) / 2
            tan_term = 2 * atan(f_m)
            cubic_term = (1 - f_m^3) / (12 * ζ)
            return log_term - tan_term + π_term - 1 + cubic_term
        else
            _a_m = FT(a_m(uf))
            return -_a_m * ζ / 2
        end
    end
end

function Psi(uf::Businger, ζ, tt::HeatTransport)
    FT = eltype(uf)
    _a_h = FT(a_h(uf))
    if abs(ζ) < eps(typeof(uf.L))
        # Psi_h in Eq. A14
        if ζ >= 0
            _π_group = FT(π_group(uf, tt))
            return -_a_h * ζ / (2 * _π_group)
        else
            return -9 * ζ / 4
        end
    else
        f_h = f_heat(uf, ζ)
        if uf.L >= 0
            _π_group = FT(π_group(uf, tt))
            return -_a_h * ζ / (2 * _π_group)
        else
            log_term = 2 * log((1 + f_h) / 2)
            return log_term + 2 * (1 - f_h) / (9 * ζ) - 1
        end
    end
end

#####
##### Gryanik
#####

"""
    Gryanik <: AbstractUniversalFunction{FT}

# References
 - [Gryanik2020](@cite)

# Equations in reference:

    ϕ_m: Eq. 13
    ϕ_h: Eq. 13
    ψ_m: Eq. 14
    ψ_h: Eq. 14

# Fields

$(DocStringExtensions.FIELDS)
"""
struct Gryanik{FT, PS} <: AbstractUniversalFunction{FT}
    "Parameter set"
    param_set::PS
    "Monin-Obhukov Length"
    L::FT
    Gryanik(param_set::PS, L::FT) where {FT, PS} = new{FT, PS}(param_set, L)
end

# CLIMAParameters wrapper
Pr_0(param_set::APS, ::Gryanik) = Pr_0_Gryanik(param_set)
a_m(param_set::APS, ::Gryanik) = a_m_Gryanik(param_set)
a_h(param_set::APS, ::Gryanik) = a_h_Gryanik(param_set)
b_m(param_set::APS, ::Gryanik) = b_m_Gryanik(param_set)
b_h(param_set::APS, ::Gryanik) = b_h_Gryanik(param_set)

function phi(uf::Gryanik, ζ, tt::MomentumTransport)
    if 0 < ζ
        FT = eltype(uf)
        _a_m = FT(a_m(uf))
        _b_m = FT(b_m(uf))
        return 1 + (_a_m * ζ) / (1 + _b_m * ζ)^(FT(2 / 3))
    else
        return phi(Businger(uf), ζ, tt)
    end
end

function phi(uf::Gryanik, ζ, tt::HeatTransport)
    if 0 < ζ
        FT = eltype(uf)
        _Pr_0 = FT(Pr_0(uf))
        _a_h = FT(a_h(uf))
        _b_h = FT(b_h(uf))
        return _Pr_0 * (1 + (_a_h * ζ) / (1 + _b_h * ζ))
    else
        return phi(Businger(uf), ζ, tt)
    end
end

function psi(uf::Gryanik, ζ, tt::MomentumTransport)
    if 0 < ζ
        FT = eltype(uf)
        _a_m = FT(a_m(uf))
        _b_m = FT(b_m(uf))
        return -3 * (_a_m / _b_m) * ((1 + _b_m * ζ)^(FT(1 / 3)) - 1)
    else
        return psi(Businger(uf), ζ, tt)
    end
end

function psi(uf::Gryanik, ζ, tt::HeatTransport)
    if 0 < ζ
        FT = eltype(uf)
        _Pr_0 = FT(Pr_0(uf))
        _a_h = FT(a_h(uf))
        _b_h = FT(b_h(uf))
        return -_Pr_0 * (_a_h / _b_h) * log1p(_b_h * ζ)
    else
        return psi(Businger(uf), ζ, tt)
    end
end

#####
##### Grachev
#####

"""
    Grachev <: AbstractUniversalFunction{FT}

# References
 - [Grachev2007](@cite)

Equations in reference:

    ϕ_m: Eq. 13
    ϕ_h: Eq. 13
    ψ_m: Eq. 14
    ψ_h: Eq. 14

# Fields

$(DocStringExtensions.FIELDS)
"""
struct Grachev{FT, PS} <: AbstractUniversalFunction{FT}
    "Parameter set"
    param_set::PS
    "Monin-Obhukov Length"
    L::FT
    Grachev(param_set::PS, L::FT) where {FT, PS} = new{FT, PS}(param_set, L)
end

# CLIMAParameters wrapper
Pr_0(param_set::APS, ::Grachev) = Pr_0_Grachev(param_set)
a_m(param_set::APS, ::Grachev) = a_m_Grachev(param_set)
a_h(param_set::APS, ::Grachev) = a_h_Grachev(param_set)
b_m(param_set::APS, ::Grachev) = b_m_Grachev(param_set)
b_h(param_set::APS, ::Grachev) = b_h_Grachev(param_set)
c_h(param_set::APS, ::Grachev) = c_h_Grachev(param_set)

function phi(uf::Grachev, ζ, tt::MomentumTransport)
    if 0 < ζ
        FT = eltype(uf)
        _a_m = FT(a_m(uf))
        _b_m = FT(b_m(uf))
        return 1 + _a_m * ζ * (1 + ζ)^FT(1 / 3) / (1 + _b_m * ζ)
    else
        return phi(Businger(uf), ζ, tt)
    end
end

function phi(uf::Grachev, ζ, tt::HeatTransport)
    if 0 < ζ
        FT = eltype(uf)
        _a_h = FT(a_h(uf))
        _b_h = FT(b_h(uf))
        _c_h = FT(c_h(uf))
        return 1 + (_a_h * ζ + _b_h * ζ^2) / (1 + _c_h * ζ + ζ^2)
    else
        return phi(Businger(uf), ζ, tt)
    end
end

function psi(uf::Grachev, ζ, tt::MomentumTransport)
    if 0 < ζ
        FT = eltype(uf)
        _a_m = FT(a_m(uf))
        _b_m = FT(b_m(uf))
        B_m = cbrt(1 / _b_m - 1)
        x = cbrt(1 + ζ)
        sqrt3 = FT(sqrt(3))
        # Note: there is a mismatch between
        # Gryanik, Eq. 26, and Grachev Eq. 12.
        # We use the Grachev Eq. 12.
        linear_term = -3 * (_a_m / _b_m) * (x - 1)
        log_term_1 = 2 * log((x + B_m) / (1 + B_m))
        log_term_2 = log((x^2 - x * B_m + B_m^2) / (1 - B_m + B_m^2))
        atan_term_1 = atan((2 * x - B_m) / (sqrt3 * B_m))
        atan_term_2 = atan((2 - B_m) / (sqrt3 * B_m))
        atan_terms = atan_term_1 - atan_term_2
        bracket_term = log_term_1 - log_term_2 + 2 * sqrt3 * atan_terms
        return linear_term + _a_m * B_m / (2 * _b_m) * bracket_term
    else
        return psi(Businger(uf), ζ, tt)
    end
end

function psi(uf::Grachev, ζ, tt::HeatTransport)
    if 0 < ζ
        FT = eltype(uf)
        _Pr_0 = FT(Pr_0(uf))
        _a_h = FT(a_h(uf))
        _b_h = FT(b_h(uf))
        _c_h = FT(c_h(uf))
        B_h = sqrt(_c_h^2 - 4)
        coeff = _a_h / B_h - _b_h * _c_h / (2 * B_h)
        log_term_1 = log((2 * ζ + _c_h - B_h) / (2 * ζ + _c_h + B_h))
        log_term_2 = log((_c_h - B_h) / (_c_h + B_h))
        log_terms = log_term_1 - log_term_2
        term_2 = _b_h / 2 * log1p(_c_h * ζ + ζ^2)
        return -coeff * log_terms - term_2
    else
        return psi(Businger(uf), ζ, tt)
    end
end

#####
##### Conversions
#####

Businger(uf::Grachev) = Businger(uf.param_set, uf.L)
Businger(uf::Gryanik) = Businger(uf.param_set, uf.L)

Grachev(uf::Businger) = Grachev(uf.param_set, uf.L)
Grachev(uf::Gryanik) = Grachev(uf.param_set, uf.L)

Gryanik(uf::Businger) = Gryanik(uf.param_set, uf.L)
Gryanik(uf::Grachev) = Gryanik(uf.param_set, uf.L)

end # module

max_forks_repo_name: jm-c/CLIMA


text:
function initial_update!(integ, cache, init::ClassicSolverInit)
    @unpack u, f, p, t = integ
    @unpack d, x, Proj = cache
    q = integ.alg.order

    if q > 5
        @warn "ClassicSolverInit might be unstable for high orders"
    end

    @unpack ddu, du, x_tmp, x_tmp2, m_tmp, K1 = cache

    # Initialize on u0
    condition_on!(x, Proj(0), view(u, :), m_tmp, K1, x_tmp.Σ, x_tmp2.Σ.mat)
    if isinplace(f)
        f(du, u, p, t)
    else
        du .= f(u, p, t)
    end
    integ.destats.nf += 1
    condition_on!(x, Proj(1), view(du, :), m_tmp, K1, x_tmp.Σ, x_tmp2.Σ.mat)

    if q < 2
        return
    end

    # Use a jac or autodiff to initialize on ddu0
    if integ.alg.initialization.init_on_du
        if isinplace(f)
            dfdt = copy(u)
            ForwardDiff.derivative!(dfdt, (du, t) -> f(du, u, p, t), du, t)

            if !isnothing(f.jac)
                f.jac(ddu, u, p, t)
            else
                ForwardDiff.jacobian!(ddu, (du, u) -> f(du, u, p, t), du, u)
            end
        else
            dfdt = ForwardDiff.derivative((t) -> f(u, p, t), t)
            if !isnothing(f.jac)
                ddu .= f.jac(du, u, p, t)
            else
                ddu .= ForwardDiff.jacobian(u -> f(u, p, t), u)
            end
        end
        ddfddu = ddu * view(du, :) + view(dfdt, :)
        condition_on!(x, Proj(2), ddfddu, m_tmp, K1, x_tmp.Σ, x_tmp2.Σ.mat)
        if q < 3
            return
        end
    end

    # Compute the other parts with classic solvers
    t0 = integ.sol.prob.tspan[1]
    dt =
        10 * OrdinaryDiffEq.ode_determine_initdt(
            u,
            t,
            1,
            integ.opts.dtmax,
            integ.opts.abstol,
            integ.opts.reltol,
            integ.opts.internalnorm,
            integ.sol.prob,
            integ,
        )
    integ.destats.nf += 2

    nsteps = q + 2
    tmax = t0 + nsteps * dt
    tstops = t0:dt:tmax
    alg = integ.alg.initialization.alg
    sol = solve(
        remake(integ.sol.prob, tspan=(t0, tmax)),
        alg,
        dense=false,
        save_start=false,
        abstol=integ.opts.abstol / 100,
        reltol=integ.opts.reltol / 100,
        saveat=tstops,
    )
    # This is necessary in order to fairly account for the cost of initialization!
    integ.destats.nf += sol.destats.nf
    integ.destats.njacs += sol.destats.njacs
    integ.destats.nsolve += sol.destats.nsolve
    integ.destats.nw += sol.destats.nw
    integ.destats.nnonliniter += sol.destats.nnonliniter
    integ.destats.nnonlinconvfail += sol.destats.nnonlinconvfail
    integ.destats.ncondition += sol.destats.ncondition

    # Filter & smooth to fit these values!
    us = [view(u, :) for u in sol.u]
    return rk_init_improve(integ, cache, sol.t, us, dt)
end

function rk_init_improve(integ, cache::GaussianODEFilterCache, ts, us, dt)
    @unpack A, Q = cache
    @unpack x, x_pred, x_filt, measurement = cache

    # Predict forward:
    make_preconditioners!(cache, dt)
    @unpack P, PI = cache

    mul!(x, P, x)

    preds = []
    filts = [copy(x)]

    # Filter through the data forwards
    for (i, (t, u)) in enumerate(zip(ts, us))
        predict_mean!(x_pred, x, A)
        predict_cov!(x_pred, x, A, Q, cache.C1, cache.default_diffusion)
        push!(preds, copy(x_pred))

        H = cache.E0 * PI
        measurement.μ .= H * x_pred.μ .- u
        X_A_Xt!(measurement.Σ, x_pred.Σ, H)

        update!(x_filt, x_pred, measurement, H, cache.K1, cache.x_tmp2.Σ.mat, cache.m_tmp)
        push!(filts, copy(x_filt))

        x = x_filt
    end

    # Smooth backwards
    for i in length(filts):-1:2
        xf = filts[i-1]
        xs = filts[i]
        xp = preds[i-1] # Since `preds` is one shorter
        smooth!(xf, xs, A, Q, integ, 1)
    end

    mul!(cache.x, PI, filts[1])

    return nothing
end


text:
module Twiddle

"""
    repeatpattern{T<:Unsigned}(::Type{T}, pattern::Unsigned)

Repeats the bitpattern of the `pattern` integer throughout the an
Unsigned integer of type `T`. Note that this assumed the size of T (in bits)
is larger than that of the `pattern` integer.

This is useful for bit-twiddling code that is to work on any word size.

Let's use a very simple example to show why this is useful.
Let's say you have a function in which a byte is masked with a bitmask
0x33.

```julia
f(x::UInt8) = x & 0x33
```

Suppose you wanted this to work for any word size, so 8 bit, 16 bit, 32 bit, 64
bit, and 128 bit words? One way is to write `f` as a set of seperate methods:

```julia
f(x::UInt8) = x & 0x33
f(x::UInt16) = x & 0x3333
f(x::UInt32) = x & 0x33333333
f(x::UInt64) = x & 0x3333333333333333
f(x::UInt128) = x & 0x33333333333333333333333333333333
```

This is clearly not ideal, you would be writing out the same method multiplet
times, but with a differently sized literal for the mask each time.

Alternatively, you could write one parametric function with repeatbyte:

```julia
f2(x::T) where {T<:Unsigned} = x & Twiddle.repeatpattern(T, 0x33)
# Or a non-parametric version using `typeof`.
f3(x::Unsigned) = x & Twiddle.repeatpattern(typeof(x), 0x33)
```

You might expect this to be less efficient - `repeatpattern` uses several operations
to generate the values 0x33, 0x3333 and so on, whereas in the seperate methods,
those literal values are hard coded. However, thanks to constant folding during
compilation, those operations are done once at compilation time and so the
native instructions generated are identical. You can check this with `@code_llvm`
or `@code_native`.
"""
@inline function repeatpattern(::Type{T}, pattern::Unsigned) where {T<:Unsigned}
    return div(typemax(T), typemax(pattern)) * pattern
end

@inline function repeatbyte(::Type{T}, byte::UInt8) where {T<:Unsigned}
    Base.depwarn("repeatbyte is deprecated, use repeatpattern instead", :repeatbyte)
    return repeatpattern(T, byte)
end

@inline mergebits(a, b, mask) = a ⊻ ((a ⊻ b) & mask)

include("bitswapping.jl")
include("masking.jl")
include("nibbles.jl")
include("counting.jl")

end


text:
using Robotics

using Test

include("checkpuma.jl")


text:
include("attributes.jl")

shift_and_rotate(X, o, M) = (X .- o[1,1:size(X, 2)]')*M'

function rastrigin(X)
    X = 0.0512*X;
    return sum(X.^2 - 10cos.(2π*X) .+ 10, dims = 2);
end

elliptic(X) = sum((1e6).^((0:size(X,2)-1)'/(size(X,2)-1+1e-6)) .* X.^2, dims = 2)

function schaffer(X)
    X = X .^ 2
    mask = vcat(2:size(X, 2),1)
    return sum(0.5.+(sin.(sqrt.(X.+X[:,mask])).^2 .-0.5)./(1 .+ 0.001*(X+X[:,mask])).^2, dims = 2)
end

function HGBat(X)
    X = 0.05X .- 1;
    sq = sum(X.^2, dims = 2)
    s = sum(X, dims = 2)
    return sqrt.(abs.(sq.^2 - s.^2)) + (0.5*sq + s)/size(X,2) .+ 0.5;
end

function rosenbrock(X)
    X = 0.02048*X .+ 1;
    sum(100*(X[:,1:end-1].^2-X[:,2:end]).^2+(X[:,1:end-1] .- 1).^2, dims = 2);
end

function griewank(X)
    X = 6X;
    return sum(X.^2,dims=2)/4000 - prod(cos.(X) ./ sqrt.(1:size(X,2))', dims=2) .+ 1;
end

ackley(X) = -20exp.(-0.2*sqrt.(mean(X.^2,dims=2))) - exp.(mean(cos.(2*pi*X),dims=2)) .+(20 + exp(1))

function happycat(X)
    X = 0.05*X;
    sq = sum(X.^2, dims=2)
    D = size(X,2)
    return abs.(sq.-D).^0.25 + (0.5*sq+sum(X,dims=2))/D .+ 0.5;
end

discus(X) = 1e6*X[:,1].^2 + sum(X[:,2:end].^2, dims=2);

########################3

# base function here
bent_cigar(X) = 100 .+ X[:,1].^2 + 1e6*sum(X[:,2:end].^2, dims = 2)

function schwefel(X)
    D = size(X,2)
    Z = 10X .+ 4.2097e2
    g = Z .* sin.(sqrt.(abs.(Z)))
    mask = Z .> 500
    temp = 500 .- mod.(Z[mask],500)
    g[mask] = temp.*sin.(sqrt.(abs.(temp))) - (Z[mask].-500).^2/10000/D

    mask = Z .< -500
    temp = mod.(abs.(Z[mask]),500) .- 500
    g[mask] = temp.*sin.(sqrt.(abs.(temp))) - (Z[mask] .- 500)/10000/D
    return 418.9829*D .- sum(g, dims = 2)
end

function bi_rastrigin(X, O, M)
    D = size(X,2)
    O = O[1,1:D]'
    s   = 1 - 1/(2*sqrt(D+20)-8.2)
    μ0 = 2.5
    μ1 = -sqrt((μ0^2-1)/s)
    Y   = (X .- O)/10
    tmp = 2sign.(O).*Y .+ μ0
    Z   = (tmp .- μ0) * M'
    F = 700 .+ min.(sum((tmp .- μ0).^2, dims=2), D .+ s*sum((tmp .- μ1).^2, dims = 2))
    F += 10*(D .- sum(cos.(2π*Z), dims = 2))
    return F
end

function rosenbrock_plus_griewangk(X)
    Y = 0.05X
    Z = Y .+ 1;
    f5   = 100*(Z.^2-Z[:, vcat(2:end,1)]).^2 + (Z .- 1).^2
    return 1900 .+ sum(f5.^2/4000-cos.(f5) .+ 1, dims = 2)
end

function hybrid1(X, S)
    1700 .+ schwefel( X[:, S[1]]) .+ rastrigin(X[:, S[2]]) .+ elliptic(X[:, S[3]])
end
function hybrid2(X, S)
    return 1600 .+
    schaffer(X[:,S[1]]) .+
    HGBat(X[:,S[2]]) .+
    rosenbrock(X[:,S[3]]) .+
    schwefel(X[:,S[4]])
end

function hybrid3(Z, S)
    return 2100 .+
    schaffer(Z[:,S[1]]) .+
    HGBat(Z[:,S[2]]) .+
    rosenbrock(Z[:,S[3]]) .+
    schwefel(Z[:,S[4]]) .+
    elliptic(Z[:,S[5]])
end

function composition1(X, O, M)
    λ  = [1 10 1]
    δ  = [10 20 30]
    bias   = [0 100 200]
    W      = zeros(size(X,1),3);
    F      = zeros(size(W));
    D = size(X,2)

    for (i, f) in enumerate([rastrigin, griewank, schwefel])
        o = O[i,1:D]'
        tmp    = sum((X .- o).^2, dims = 2)
        W[:,i] = 1 ./ (sqrt.(tmp) .+ 1e-10) .* exp.(-tmp/2/D/δ[i]^2)
        F[:,i] = f((X .- o)*M[(i-1)*D+1:i*D,:]')
    end

    W = W ./ sum(W, dims=2);
    return 2200 .+ sum(W.*(λ.*F .+ bias), dims=2);
end

function composition2(X, O, M)
    λ = [10 1e-6 10 1]
    δ  = [10 20 30 40]
    bias   = [0 100 200 300]

    W      = zeros(size(X,1),4)
    F      = zeros(size(W))
    D = size(X,2)

    for (i, f) in enumerate([ackley,elliptic,griewank,rastrigin])
        o = O[i,1:D]'
        Z = (X .- o)
        M_ = M[(i-1)*D+1:i*D,:]'
        tmp    = sum(Z.^2, dims = 2)
        W[:,i] = 1 ./ (sqrt.(tmp) .+ 1e-10) .* exp.(-tmp/2/D/δ[i]^2)
        F[:,i] = f(Z*M_)
    end

    W = W ./ sum(W, dims=2);
    return 2400 .+ sum(W.*(λ.*F .+ bias), dims=2);
end

function composition3(X, O, M)
    λ = [10 1 10 1e-6 1]
    δ  = [10 20 30 40 50]
    bias   = [0 100 200 300 400]
    W      = zeros(size(X,1),5)
    F      = zeros(size(W))
    D = size(X,2)

    for (i, f) in enumerate([rastrigin,happycat,ackley,discus,rosenbrock])
        o = O[i,1:D]'
        tmp    = sum((X .- o).^2, dims = 2)
        W[:,i] = 1 ./ (sqrt.(tmp) .+ 1e-10) .* exp.(-tmp/2/D/δ[i]^2)
        F[:,i] = f((X .- o)*M[(i-1)*D+1:i*D,:]')
    end

    W = W ./ sum(W, dims=2);
    return 2500 .+ sum(W.*(λ.*F .+ bias), dims=2);
end

function cec2020_f1(x)
    D = (size(x,2))
    o, M, _ = load_data_from_dict(1, D)
    bent_cigar(shift_and_rotate(x, o, M))
end

function cec2020_f2(x)
    D = (size(x,2))
    o, M, _ = load_data_from_dict(2, D)
    1100 .+ schwefel(shift_and_rotate(x, o, M))
end

function cec2020_f3(x)
    D = (size(x,2))
    o, M, _ = load_data_from_dict(3, D)
    bi_rastrigin(x, o, M)
end

function cec2020_f4(x)
    D = (size(x,2))
    o, M, _ = load_data_from_dict(4, D)
    rosenbrock_plus_griewangk(shift_and_rotate(x, o, M))
end

function cec2020_f5(x)
    D = (size(x,2))
    o, M, S_ = load_data_from_dict(5, D)
    p     = ceil.(Int, D*[0.3, 0.3, 0.4])
    p[1]  = D - sum(p[2:end])
    p     = vcat(0,cumsum(p))
    S = [S_[p[i]+1:p[i+1]] for i in 1:length(p)-1]
    hybrid1(shift_and_rotate(x, o, M), S)
end

function cec2020_f6(x)
    D = size(x,2)
    o, M, S_ = load_data_from_dict(6, D)

    if D == 5
        p = [1, 1, 1, 2]
    else
        p = ceil.(Int, D*[0.2, 0.2, 0.3, 0.3])
        p[1]  = D - sum(p[2:end])
    end

    p     = vcat(0,cumsum(p))
    S = [S_[p[i]+1:p[i+1]] for i in 1:length(p)-1]

    hybrid2(shift_and_rotate(x, o, M), S)
end

function cec2020_f7(x)
    D = (size(x,2))
    o, M, S_ = load_data_from_dict(7, D)

    if D == 5
        p = ones(Int, 5)
    else
        p = ceil.(Int, D*[0.1, 0.2, 0.2, 0.2, 0.3])
        p[1]  = D - sum(p[2:end])
    end

    p     = vcat(0,cumsum(p))
    S = [S_[p[i]+1:p[i+1]] for i in 1:length(p)-1]
    hybrid3(shift_and_rotate(x, o, M), S)
end

function cec2020_f8(x)
    D = (size(x,2))
    O, M, _ = load_data_from_dict(8, D)
    composition1(x, O, M)
end

function cec2020_f9(x)
    D = (size(x,2))
    O, M, _ = load_data_from_dict(9, D)
    composition2(x, O, M)
end

function cec2020_f10(x)
    D = (size(x,2))
    O, M, _ = load_data_from_dict(10, D)
    composition3(x, O, M)
end

cec2020_f1(x::AbstractVector{T}) where T <: Real = cec2020_f1(x')[1]
cec2020_f2(x::AbstractVector{T}) where T <: Real = cec2020_f2(x')[1]
cec2020_f3(x::AbstractVector{T}) where T <: Real = cec2020_f3(x')[1]
cec2020_f4(x::AbstractVector{T}) where T <: Real = cec2020_f4(x')[1]
cec2020_f5(x::AbstractVector{T}) where T <: Real = cec2020_f5(x')[1]
cec2020_f6(x::AbstractVector{T}) where T <: Real = cec2020_f6(x')[1]
cec2020_f7(x::AbstractVector{T}) where T <: Real = cec2020_f7(x')[1]
cec2020_f8(x::AbstractVector{T}) where T <: Real = cec2020_f8(x')[1]
cec2020_f9(x::AbstractVector{T}) where T <: Real = cec2020_f9(x')[1]
cec2020_f10(x::AbstractVector{T}) where T <: Real = cec2020_f10(x')[1]

export cec2020_f1, cec2020_f2, cec2020_f3, cec2020_f4, cec2020_f5
export cec2020_f6, cec2020_f7, cec2020_f8, cec2020_f9, cec2020_f10
export  get_cec2020_minimum, get_cec2020_problem


text:
module URIs2

import URIs

export URI, uri2filepath, filepath2uri, @uri_str

struct URI
    scheme::Union{String,Nothing}
    authority::Union{String,Nothing}
    path::String
    query::Union{String,Nothing}
    fragment::Union{String,Nothing}
end

function percent_decode(str::AbstractString)
    return URIs.unescapeuri(str)
end

function URI(value::AbstractString)
    m = match(r"^(([^:/?#]+?):)?(\/\/([^/?#]*))?([^?#]*)(\?([^#]*))?(#(.*))?", value)

    m===nothing && error("Invalid argument.")

    return URI(
        m.captures[2],
        m.captures[4]===nothing ? nothing : percent_decode(m.captures[4]),
        m.captures[5]===nothing ? nothing : percent_decode(m.captures[5]),
        m.captures[7]===nothing ? nothing : percent_decode(m.captures[7]),
        m.captures[9]===nothing ? nothing : percent_decode(m.captures[9])
    )
end

URI(uri::URI) = uri

function URI(;
    scheme::Union{AbstractString,Nothing}=nothing,
    authority::Union{AbstractString,Nothing}=nothing,
    path::AbstractString="",
    query::Union{AbstractString,Nothing}=nothing,
    fragment::Union{AbstractString,Nothing}=nothing
    )
    return URI(scheme, authority, path, query, fragment)
end

@inline function is_rfc3986_unreserved(c::Char)
    return 'A' <= c <= 'Z' ||
        'a' <= c <= 'z' ||
        '0' <= c <= '9' ||
        c == '-' ||
        c == '.' ||
        c == '_' ||
        c == '~'
end

@inline function is_rfc3986_sub_delim(c::Char)
    return c == '!' ||
        c == '$' ||
        c == '&' ||
        c == '\'' ||
        c == '(' ||
        c == ')' ||
        c == '*' ||
        c == '+' ||
        c == ',' ||
        c == ';' ||
        c == '='
end

@inline function is_rfc3986_pchar(c::Char)
    return is_rfc3986_unreserved(c) ||
        is_rfc3986_sub_delim(c) ||
        c == ':' ||
        c == '@'
end

@inline function is_rfc3986_query(c::Char)
    return is_rfc3986_pchar(c) || c=='/' || c=='?'
end

@inline function is_rfc3986_fragment(c::Char)
    return is_rfc3986_pchar(c) || c=='/' || c=='?'
end

@inline function is_rfc3986_userinfo(c::Char)
    return is_rfc3986_unreserved(c) ||
        is_rfc3986_sub_delim(c) ||
        c == ':'
end

@inline function is_rfc3986_reg_name(c::Char)
    return is_rfc3986_unreserved(c) ||
        is_rfc3986_sub_delim(c)
end

function encode(io::IO, s::AbstractString, issafe::Function)
    for c in s
        if issafe(c)
            print(io, c)
        else
            print(io, '%')
            print(io, uppercase(string(Int(c), base=16, pad=2)))
        end
    end
end

@inline function is_ipv4address(s::AbstractString)
    if length(s)==1
        return '0' <= s[1] <= '9'
    elseif length(s)==2
        return '1' <= s[1] <= '9' && '0' <= s[2] <= '9'
    elseif length(s)==3
        return (s[1]=='1' && '0' <= s[2] <= '9' && '0' <= s[3] <= '9') ||
            (s[1]=='2' && '0' <= s[2] <= '4' && '0' <= s[3] <= '9') ||
            (s[1]=='2' && s[2] == '5' && '0' <= s[3] <= '5')
    else
        return false
    end
end

@inline function is_ipliteral(s::AbstractString)
    # TODO Implement this
    return false
end

function encode_host(io::IO, s::AbstractString)
    if is_ipv4address(s) || is_ipliteral(s)
        print(io, s)
    else
        # The host must be a reg-name
        encode(io, s, is_rfc3986_reg_name)
    end
end

function encode_path(io::IO, s::AbstractString)
    # TODO Write our own version
    print(io, URIs.escapepath(s))
end

function Base.print(io::IO, uri::URI)
    scheme = uri.scheme
    authority = uri.authority
    path = uri.path
    query = uri.query
    fragment = uri.fragment

 	if scheme!==nothing
        print(io, scheme)
        print(io, ':')
 	end

 	if authority!==nothing
        print(io, "//")

		idx = findfirst("@", authority)
		if idx !== nothing
			# <user>@<auth>
			userinfo = SubString(authority, 1:idx.start-1)
			host_and_port = SubString(authority, idx.start + 1)
			encode(io, userinfo, is_rfc3986_userinfo)
            print(io, '@')
        else
            host_and_port = SubString(authority, 1)
		end

		idx3 = findfirst(":", host_and_port)
		if idx3 === nothing
            encode_host(io, host_and_port)
		else
			# <auth>:<port>
            encode_host(io, SubString(host_and_port, 1:idx3.start-1))
			print(io, SubString(host_and_port, idx3.start))
        end
     end

     # Append path
     encode_path(io, path)

    if query!==nothing
        print(io, '?')
        encode(io, query, is_rfc3986_query)
    end

 	if fragment!==nothing
        print(io, '#')
        encode(io, fragment, is_rfc3986_fragment)
    end

    return nothing
end

function Base.string(uri::URI)
    io = IOBuffer()

    print(io, uri)

    return String(take!(io))
end

macro uri_str(ex)
    return URI(ex)
end

include("uri_helpers.jl")

end


text:
using Test
import ParallelNeighbors
using Random: MersenneTwister, shuffle!
using NearestNeighbors: BruteTree, knn, Euclidean, SqEuclidean
using LinearAlgebra: LowerTriangular, Diagonal
using MultivariateStats: fit, predict, MDS
using Suppressor: @suppress_err

# TODO: change to `as` import once we drop Julia 1.5 support
const knn_parallel = ParallelNeighbors.knn

function create_unique_data(dim, size, rng = MersenneTwister(0))
    # Use MDS to calculate samples that are of pairwise-unique distance to each other
    # https://www.stat.pitt.edu/sungkyu/course/2221Fall13/lec8_mds_combined.pdf
    base_vector = Float64.(collect(1:size*size))
    unique_dists = shuffle!(rng, base_vector) # unique distances
    Dfull = reshape(unique_dists, size, size)
    D = LowerTriangular(Dfull)
    D = D + D' - Diagonal(D)
    # we don't care if some of the eigenvalues are degenerate in testing
    mds = @suppress_err fit(MDS, D; distances = true, maxoutdim = dim)
    return predict(mds)
end

function create_unique_split(dim, n_train, n_test)
    X = create_unique_data(dim, n_train + n_test)
    return X[:, 1:n_train], X[:, n_train+1:end]
end

function reference_knn(Xtrain, Xtest, k; metric = Euclidean())
    knn(BruteTree(Xtrain, metric), Xtest, k, true)
end

@testset "ParallelNeighbors.jl" begin
    ks = 1:5
    dims = 1:3
    batch_sizes = 1:5
    train_sizes = 1:8
    test_sizes = 1:8

    # we cannot sample less than k training points (otherwise how should we determine the k nearest neighbors)
    for k in ks, dim in dims,
        batch_size in batch_sizes,
        test_size in test_sizes,
        train_size in train_sizes[k:end]

        # test distances with discrete data and indices with continuous data
        Xtrain, Xtest = create_unique_split(dim, train_size, test_size)

        pointwise_batch_size = min(batch_size, test_size)
        batched_batch_size = max(min(batch_size, train_size), k)

        # evaluate with euclidean distance
        results_euclidean = (reference_knn(Xtrain, Xtest, k; metric = Euclidean()),
            knn_parallel(Xtrain, Xtest, k; metric = Euclidean(), algorithm = :full),
            knn_parallel(Xtrain, Xtest, k, pointwise_batch_size; metric = Euclidean(), algorithm = :hybrid_batch_test),
            knn_parallel(Xtrain, Xtest, k, batched_batch_size; metric = Euclidean(), algorithm = :hybrid_batch_all))

        # TODO: evaluate with squared euclidean distance once the following pull request is merged
        # https://github.com/KristofferC/NearestNeighbors.jl/pull/137
        # results_sqeuclidean = (reference_knn(Xtrain, Xtest, k; metric = SqEuclidean()),
        #     knn_full(Xtrain, Xtest, k; metric = SqEuclidean()),
        #     knn_pointwise(Xtrain, Xtest, k, pointwise_batch_size; metric = SqEuclidean()),
        #     knn_batch(Xtrain, Xtest, k, batched_batch_size; metric = SqEuclidean()))

        ref_idx, full_idx, point_idx, batch_idx = map(first, results_euclidean)
        @test all(ref_idx .== point_idx)
        @test all(ref_idx .== batch_idx)
        @test all(ref_idx .== full_idx)

        ref_dist, full_dist, point_dist, batch_dist = map(last, results_euclidean)
        @test all(ref_dist .≈ point_dist)
        @test all(ref_dist .≈ batch_dist)
        @test all(ref_dist .≈ full_dist)
    end
end


text:
function plot_result(plot_dictionary, input_data, spatial_data, spatial_data_matrices, optimal_path)
    A = zeros(3, size(plot_dictionary["overlay_image"],1), size(plot_dictionary["overlay_image"],2) )
    look_up_table = [0 255 0; 176 224 230;255 255 255;255 255 0; 218 165 32;255 160 122;255 0 0] ./ 255

    for i in 1:size(plot_dictionary["overlay_image"], 1)
        for j in 1:size(plot_dictionary["overlay_image"], 2)
            A[:, i, j] = look_up_table[Int.(plot_dictionary["overlay_image"][i, j]), :]'
        end
    end
    p = Images.colorview(ColorTypes.RGB, A)
    Plots.plot(p)
    for idx = 1 : size(optimal_path, 1) - 1
        x1 = spatial_data["nodes"][optimal_path[idx], 2]
        y1 = spatial_data["nodes"][optimal_path[idx], 3]
        x2 = spatial_data["nodes"][optimal_path[idx + 1], 2]
        y2 = spatial_data["nodes"][optimal_path[idx + 1], 3]
        ohl_cable = ohl_or_cable(spatial_data, spatial_data_matrices, optimal_path[idx], optimal_path[idx + 1])
        if input_data["technology"] == "dc"
            if ohl_cable == "ugc"
                p = Plots.plot!([x1, x2], [y1, y2], color = :brown, linewidth = 2, label = "")
            else
                p = Plots.plot!([x1, x2], [y1, y2], color = :brown, linestyle = :dot, linewidth = 2, label = "")
            end
        else
            if ohl_cable == "ugc"
                p = Plots.plot!([x1, x2], [y1, y2], color = :black, linewidth = 2, label = "")
            else
                p = Plots.plot!([x1, x2], [y1, y2], color = :black, linestyle = :dot, linewidth = 2, label = "")
            end
        end
    end
    Plots.plot!(xlims=(input_data["boundaries"]["xmin"],input_data["boundaries"]["xmax"]), ylims = (input_data["boundaries"]["ymin"], input_data["boundaries"]["ymax"]))
    plot_file = join(["plot_test_",input_data["strategy"],"_",input_data["technology"],".pdf"])
    Plots.savefig(p, plot_file)
 end

 function ohl_or_cable(spatial_data, spatial_data_matrices, node_id1, node_id2)
    index_s1 = findall(node_id1 .== spatial_data["segments"][:, 2])
    index_s2 = findall(node_id2 .== spatial_data["segments"][:, 3])
    index_s = intersect(index_s1, index_s2)
    if isempty(index_s)
        index_s1 = findall(node_id2 .== spatial_data["segments"][:, 2])
        index_s2 = findall(node_id1 .== spatial_data["segments"][:, 3])
        index_s = intersect(index_s1, index_s2)
    end
    if any(index_s[1] .== spatial_data_matrices["ohl"]["segments"][:, 1])
        ohl_cable = "ohl"
    else
        ohl_cable = "ugc"
    end
    return ohl_cable
 end


text:
using Test, Random
using NiLang, NiLang.AD

function _sce(x::AbstractArray{T,N}, p) where {T,N}
    x = x .- maximum(x; dims=N)  # avoid data overflow
    rho = exp.(x)
    Z = sum(rho; dims=N)
    return dropdims(sum((log.(Z) .- x) .* p; dims=N), dims=N)
end

@testset "softmax_crossentropy" begin
    Random.seed!(2)
    x = randn(10)
    x0 = copy(x)
    p = randn(10); p=p./maximum(p)
    res = _sce(x, p)
    imax = 0
    Z = 0.0
    out = 0.0
    xmax = 0.0
    x_ = x
    p_ = p
    @instr i_softmax_crossentropy(x_, p_, imax, xmax, Z, out)
    @show Z
    @test isapprox(imax, argmax(x0), atol=1e-8)
    @test isapprox(out, res[], atol=1e-8)
    @instr (~i_softmax_crossentropy)(x_, p_, imax, xmax, Z, out)
    args = x_, p_, imax, xmax, Z, out
    @test check_inv(i_softmax_crossentropy, args)
    args = x_, p_, imax, xmax, Z, out
    @test check_grad(i_softmax_crossentropy, args; iloss=6, verbose=true)
end

@testset "logsumexp" begin
	function logsumexp2(x)
	  	mx = maximum(x)
	  	log.(sum(exp.(x .- mx))) .+ mx
	end

	x = randn(100)
	@test i_ascending!(Float64[], Int[], x)[1][end] == maximum(x)
	@test i_logsumexp(0.0, 0.0, Float64[], Int[], x)[1] ≈ logsumexp2(x)
	@test check_inv(i_logsumexp, (0.0, 0.0, Float64[], Int[], x))
end


text:
module rings
using LinearAlgebra, DSP
#include("..\\cheb.jl")
include("cheb.jl")
using .cheb
export ring

"""

UNITS: cgs

Σ is an even function of x
Φ is an odd function of x
U is an odd function of x
"""

# Consts
G  = 6.67430e-8 # Gravitational constant, cgs
M♄ = 5.6834e29  # mass of Saturn, in grams
A₁ = 6.7187     # Goldreich-Tremaine constant
Mm = 1.0e8      # Megameter. Typical orbit is 100Mm.

Pan     = (4.95e19, 1.33584e10) # (mass, orbital radius) of Pan (Encke gap)
Daphnis = (8.40e16, 1.36504e10) # (mass, orbital radius) of Daphnis (Keeler gap)

mutable struct ring # an abstract, idealized 1D (in r) ring, with periodic BC's and a single forcing moon
    X₀::Float64  # THE key nondimensional parameter.
    ξ₀::Float64 # the other key non-dimensional parameter: sets the size of the domain (in dim'less units)
    N::Integer  # number of nodes for a full ring
    β::Float64 # exponent
    nn::Integer # half-ring
    wr::Any
    wc::Any
    ξ::Any
    Ξ::Any
    σ::Any
    Ar::Any
    a::Any
    Br::Any
    b::Any
    γ::Float64
end
function ring(X₀,ξ₀,N,β=0.0,γ=0.0;h=0.0) # h is the Hill radius in units of Lcrit
    @assert X₀ > 2^24 # 2^24 is critical value for X; below this, the method does not apply
    @assert ξ₀ > 3   # needs to be greater than a few; this is a bit of a fuzzy boundary, but certainly >3, I believe.
    ξvisc = X₀^(1/24)
    #
    # Find nn, which is half of N, since we treat only half the domain.
    # (Why not just specify nn to start? So we can test/try both odd and even N for the full domain.)
    nn = N÷2 + N%2 # nn: "non-negative" (positions)
    #
    # Get collocation points and differentiation matrices
    (y,DM,DeM,DoM) = chebdif(N,1)
    D1 = DM[:,:,1]   # full derivative matrix
    De = DeM[:,:,1]  # derivative of even function
    Do = DoM[:,:,1]  # derivative of odd function
    #
    # Get (dim'less) velocity and reduced (to half domain) velocity
    # Velocity with zero Hill radius:
    w0  = (ξ₀-1)^-4 .* ( (y .+ ξ₀/(ξ₀-1)).^-4 - (y .- ξ₀/(ξ₀-1)).^-4 )
    # Velocity including Hill radius effect:
        aa = 0.711557; bb = -7.58607
        function g(xoh)
            max(2.5,abs(xoh))
        end
        function f(xoh)
            (1+ aa/g(xoh) + bb/g(xoh)/g(xoh))^(-1)
        end
        function ww(x)
            sign(x)/x^4 * f(x/h)
        end
        xL = y * (ξ₀-1) .+ ξ₀
        xR = y * (ξ₀-1) .- ξ₀
        wL = ww.(xL)
        wR = ww.(xR)
        ww  = wL + wR
    wr = w0[1:nn]      # w (right)
    wc = 0.0 * (ww-w0)[1:nn] # w "correction"
    #
    # Form full A matrix. Could use sparse matrix for diagm() but not really needed honestly; we are just constructing
    # the matrix, and memory is not an issue.
    A = diagm((1+β) .-wr.^2) * De
    #
    # Get reduced A and column a.
    Ar = A[:,2:end-1]
    a  = A[:,end]
    @assert size(Ar) == (nn,nn-2)
    @assert size(a) == (nn,)
    #
    # Repeat for B1.
    τ = X₀^(-1/8)
    print("Dimless relaxation time: $τ\n")
    B1 = diagm(wr) ./ (ξ₀-1) ./ τ
    B1r = B1[:,2:end-1]
    b1  = B1[:,end]
    @assert size(B1r) == (nn,nn-2)
    @assert size(b1) == (nn,)
    #
    # And B2.
    B2 = 2 * diagm(wr) * diagm(Do*wr)
    B2r = B2[:,2:end-1]
    b2  = B2[:,end]
    @assert size(B2r) == (nn,nn-2)
    @assert size(b2) == (nn,)
    #
    # Now solve.
    Br = B1r + B2r
    b  = b1  + b2
    σ = (Ar - Br) \ (b - a)
    #
    # And return.
    yr = y[1:nn]
    σ = [0, σ..., 1]
    σσ = [0, σ...] # for plotting purposes
    ξ = ξ₀ .- yr * (ξ₀-1)
    ξξ = [0, ξ...] # for plotting purposes
    ΞΞ = ξξ / ξvisc
    return ring(X₀,ξ₀,N,β,nn,wr,wc,ξξ,ΞΞ,σσ,Ar,a,Br,b,γ)
end

function getA!(r::ring;hrel=0.0)
    N = r.N; σ = r.σ[2:end]; β = r.β; wr = r.wr; wc = r.wc; γ=r.γ
    (y,DM,DeM,DoM) = chebdif(N,1)
    D1 = DM[:,:,1]   # full derivative matrix
    De = DeM[:,:,1]  # derivative of even function
    Do = DoM[:,:,1]  # derivative of odd function
    #
    @assert 0.0 ≤ hrel ≤ 1.0
    w = wr #+ hrel*wc
    A = diagm((1+β)*σ.^(β+γ) .- w.^2) * De
    #
    r.Ar = A[:,2:end-1]
    r.a  = A[:,end]
    return 0
end

function getB!(r::ring;hrel=0.0)
    if r.γ == 0
        return 0
    else
        N = r.N; σ = r.σ[2:end]; β = r.β; wr = r.wr; γ=r.γ; ξ₀=r.ξ₀
        (y,DM,DeM,DoM) = chebdif(N,1)
        D1 = DM[:,:,1]   # full derivative matrix
        De = DeM[:,:,1]  # derivative of even function
        Do = DoM[:,:,1]  # derivative of odd function
        #
        @assert 0.0 ≤ hrel ≤ 1.0
        w = wr #+ hrel*r.wc
        τ = r.X₀^(-1/8)
        B1 = diagm(w.*σ.^γ) ./ (ξ₀-1) ./ τ
        B1r = B1[:,2:end-1]
        b1  = B1[:,end]
        #@assert size(B1r) == (nn,nn-2)
        #@assert size(b1) == (nn,)
        #
        # And B2.
        B2 = 2 * diagm(wr) * diagm(Do*wr)
        B2r = B2[:,2:end-1]
        b2  = B2[:,end]
        #@assert size(B2r) == (nn,nn-2)
        #@assert size(b2) == (nn,)
        #
        # Now solve.
        Br = B1r + B2r
        b  = b1  + b2
        r.Br=Br; r.b=b
    end
    return 0
end

function iterate!(r::ring)
    getA!(r)
    getB!(r)
    Br = r.Br;  Ar = r.Ar;  b = r.b;  a = r.a;
    σ = (Ar - Br) \ (b - a)
    r.σ = [0, 0, σ..., 1]
    return 0
end

#(X₀,ξ₀,N,β=1.0)

struct rring
    r::ring
    moon::Tuple # mass, a0
    disk::Tuple # nu, srel
    X₀::Float64 # fundamental parameter
    β::Float64
    γ::Float64
    Lset::Tuple # lengthscales
    gap::Tuple # width, Δwidth
    h::Float64  # Hill radius
end
function rring(m,β,γ;N=1000,a0=1.0e10,ν=1.0,srel=30,ξ₀=10)#,h=0.0)
#    function gimmearing(m,ν=1.0,srel=30,β=2,a₀=100Mm,N=1000,ξ₀=10;γ=0)
        function Ω(a₀)
            return √(G*M♄/a₀^3)
        end
        a₀=a0
        Ω = Ω(a₀)
        α = A₁^2 / 18π * Ω * (m/M♄)^2 * a₀^5
        s₀ = srel / Ω
        κ₀ = 3ν
        X₀ = α^2/κ₀^5/s₀^3
#        return
    h = a0 * (m/3/M♄)^(1/3)

    Lvisc = (α/κ₀)^(1/3)
    w₀ = √(κ₀/s₀)
    Lcrit = (α/w₀)^(1/4)
    Lαs = (α*s₀)^0.2          # another lengthscale constructed from α, ν, s
    Lκs = √(κ₀*s₀)
    Lset = (Lvisc,Lcrit,Lαs,Lκs)
        r = ring(X₀,ξ₀,N,β,γ;h=h/Lcrit)
#    end
#    r = gimmering(m,ν,srel,β,a₀,N,ξ₀=10)
    for i in 1:5
        iterate!(r)
    end
    # Smooth:
    r.σ = conv(r.σ,[1,4,6,4,1]./16)[3:end-2]
    #

    i10 = findfirst(x->x>0.1,r.σ)
    ξ10 = r.ξ[i10]

    i50 = findfirst(x->x>0.5,r.σ)
    ξ50 = r.ξ[i50]

    i90 = findfirst(x->x>0.9,r.σ)
    ξ90 = r.ξ[i90]

    width = 2*(ξ50 * Lcrit)
    Δwidth = 2*(ξ90-ξ10)*Lcrit
    #width=(0.0,0.0)
    #Δwidth=(0.0,0.0)
    h*=0
    return rring(r,(m,a0),(ν,srel),X₀,β,γ,Lset,(width,Δwidth,2*ξ90*Lcrit,2*ξ10*Lcrit),h)
    # Get 10%, 50%, 90% points
end

end#module

max_stars_repo_name: diakosmos/razor

max_issues_repo_name: diakosmos/razor

max_forks_repo_name: diakosmos/razor


text:
using Pkg
using Pollen

using GraphPlot
using FluxTraining
const PACKAGE = FluxTraining

# Create target folder
DIR = abspath(mkpath(ARGS[1]))

# Create Project
m = PACKAGE
ms = [m, m.Events, m.Phases, m.Loggables]

@info "Creating project..."
project = Project(
    Pollen.Rewriter[
        Pollen.DocumentFolder(pkgdir(m), prefix = "documents"),
        Pollen.ParseCode(),
        Pollen.ExecuteCode(),
        Pollen.PackageDocumentation(ms),
        Pollen.DocumentGraph(),
        Pollen.SearchIndex(),
        Pollen.SaveAttributes((:title,)),
        Pollen.LoadFrontendConfig(pkgdir(m))
    ],
)

@info "Rewriting documents..."
Pollen.rewritesources!(project)

@info "Writing to disk at \"$DIR\"..."
builder = Pollen.FileBuilder(
    Pollen.JSON(),
    DIR,
)
Pollen.build(
    builder,
    project,
)


text:
abstract type VTPRRuleModel <: MixingRule end

struct VTPRRule{γ} <: VTPRRuleModel
    components::Array{String,1}
    activity::γ
    references::Array{String,1}
end

@registermodel VTPRRule
export VTPRRule
function VTPRRule(components::Vector{String}; activity = UNIFAC, userlocations::Vector{String}=String[],activity_userlocations::Vector{String}=String[], verbose::Bool=false)
    init_activity = activity(components;userlocations = activity_userlocations,verbose)

    references = ["10.1016/S0378-3812(01)00626-4"]
    model = VTPRRule(components, init_activity,references)
    return model
end

function ab_premixing(::Type{PR},mixing::VTPRRule,Tc,pc,kij)
    Ωa, Ωb = ab_consts(PR)
    _Tc = Tc.values
    _pc = pc.values
    a = epsilon_LorentzBerthelot(SingleParam(pc, @. Ωa*R̄^2*_Tc^2/_pc),kij)
    bi = @. Ωb*R̄*_Tc/_pc
    bij = ((bi.^(3/4).+bi'.^(3/4))/2).^(4/3)
    b = PairParam("b",Tc.components,bij)
    return a,b
end

function mixing_rule(model::PRModel,V,T,z,mixing_model::VTPRRuleModel,α,a,b,c)
    n = sum(z)
    invn = (one(n)/n)
    invn2 = invn^2
    lnγ_res_ = lnγ_res(mixing_model.activity,1e5,T,z)
    g_E_res = sum(z[i]*R̄*T*lnγ_res_[i] for i ∈ @comps)*invn
    b̄ = dot(z,Symmetric(b),z) * invn2
    c̄ = dot(z,c)/n
    Σab = invn*sum(z[i]*a[i,i]*α[i]/b[i,i]/(R̄*T) for i ∈ @comps)
    ā = b̄*R̄*T*(Σab-1/0.53087*(g_E_res/(R̄*T)))
    return ā,b̄,c̄
end


text:
function bubble_sort(seq::Vector{Int64})::Vector{Int64}
    l = length(seq)
    for _ = 0:l-1
        for n = 1:l-1
            if seq[n+1] < seq[n]
                seq[n], seq[n+1] = (seq[n+1], seq[n])
            end
        end
    end
    return seq
end

if abspath(PROGRAM_FILE) == @__FILE__
    unsorted = [14, 11, 19, 5, 16, 10, 19, 12, 5, 12]
    expected = [5, 5, 10, 11, 12, 12, 14, 16, 19, 19]
    @assert(bubble_sort(unsorted) == expected)
    println("OK")
end


text:
export rule

# NOTE: We removed this rule in favor of Cateforical distribution, Gaussian Mixture node always returns Categorical for switch distribution
# In case of Bernoulli prior later on multiplication of Categorical and Bernoulli results in Bernoulli posterior marginal
# @rule NormalMixture{2}(:switch, Marginalisation) (q_out::Any, q_m::NTuple{2, NormalMeanVariance}, q_p::NTuple{2, Gamma}) = begin
#     U1 = -score(AverageEnergy(), NormalMeanPrecision, Val{ (:out, :μ, :τ) }, map(as_marginal, (q_out, q_m[1], q_p[1])), nothing)
#     U2 = -score(AverageEnergy(), NormalMeanPrecision, Val{ (:out, :μ, :τ) }, map(as_marginal, (q_out, q_m[2], q_p[2])), nothing)
#     return Bernoulli(clamp(softmax((U1, U2))[1], tiny, 1.0 - tiny))
# end

@rule NormalMixture{N}(:switch, Marginalisation) (q_out::Any, q_m::NTuple{N, UnivariateNormalDistributionsFamily}, q_p::NTuple{N, GammaDistributionsFamily}) where { N } = begin
    U = map(zip(q_m, q_p)) do (m, p)
        return -score(AverageEnergy(), NormalMeanPrecision, Val{ (:out, :μ, :τ) }, map((q) -> Marginal(q, false , false), (q_out, m, p)), nothing)
    end
    return Categorical(clamp.(softmax(U), tiny, 1.0 - tiny))
end

@rule NormalMixture{N}(:switch, Marginalisation) (q_out::Any, q_m::NTuple{N, MultivariateNormalDistributionsFamily}, q_p::NTuple{N, Wishart}) where { N } = begin
    U = map(zip(q_m, q_p)) do (m, p)
        return -score(AverageEnergy(), MvNormalMeanPrecision, Val{ (:out, :μ, :Λ) }, map((q) -> Marginal(q, false, false), (q_out, m, p)), nothing)
    end
    return Categorical(clamp.(softmax(U), tiny, 1.0 - tiny))
end


text:
############################################################################################
"""
$(TYPEDEF)

Default arguments for SMC constructor.

# Fields
$(TYPEDFIELDS)
"""
struct SMCDefault{
    F<:Function,
    T<:BaytesCore.ResamplingMethod,
    B<:BaytesCore.UpdateBool,
    U<:BaytesCore.UpdateBool
}
    "Number of tuning steps used when constructing sampler."
    Ntuning::Int64
    "Threshold for resampling chains."
    resamplingmethod::T
    "Threshold for resampling chains. Set to 1.0 if resampling should always be applied."
    resamplingthreshold::Float64
    "Function that is applied to correlation vector of parameter to determine if jittering can be stopped."
    jitterfun::F
    "Boolean if fixed number of jittersteps are applied or jittering is based on parameter correlation."
    jitteradaption::B
    "Stopping threshold against mean correlation of jittered model parameter."
    jitterthreshold::Float64
    "Minimum number of jittering steps."
    jittermin::Int64
    "Maximum number of jittering steps."
    jittermax::Int64
    "Boolean if generate(_rng, objective) for corresponding model is stored in PF Diagnostics."
    generated::U
    function SMCDefault(;
        Ntuning=10,
        resamplingmethod::T = BaytesFilters.Systematic(),
        resamplingthreshold=0.9,
        jitterfun::F = maximum,
        jitteradaption::B = BaytesCore.UpdateTrue(),
        jitterthreshold=0.9,
        jittermin=1,
        jittermax=10,
        generated=BaytesCore.UpdateFalse(),
    ) where {
    F<:Function,
    T<:BaytesCore.ResamplingMethod,
    B<:BaytesCore.UpdateBool
    }
        ArgCheck.@argcheck 0 < Ntuning
        ArgCheck.@argcheck 0.0 <= resamplingthreshold <= 1.0
        ArgCheck.@argcheck 0.0 <= jitterthreshold <= 1.0
        ArgCheck.@argcheck 0 < jittermin <= jittermax
        return new{F, T, B, typeof(generated)}(
            Ntuning,
            resamplingmethod,
            resamplingthreshold,
            jitterfun,
            jitteradaption,
            jitterthreshold,
            jittermin,
            jittermax,
            generated
            )
    end
end

################################################################################
function update!(particles::SMCParticles, model::ModelWrapper, tagged::Tagged, upd::BaytesCore.UpdateFalse)
    return nothing
end
function update!(particles::SMCParticles, model::ModelWrapper, tagged::Tagged, upd::BaytesCore.UpdateTrue)
    # Update parameter particles that are not tagged
    for iter in eachindex(particles.model)
        ModelWrappers.fill!(particles.model[iter],
            merge(model.val, ModelWrappers.subset(particles.model[iter].val, tagged.parameter))
        )
    end
    return nothing
end

############################################################################################
"""
$(TYPEDEF)

SMC Algorithm.

# Fields
$(TYPEDFIELDS)
"""
struct SMC{A<:SMCParticles,B<:SMCTune} <: AbstractAlgorithm
    particles::A
    tune::B
    function SMC(particles::A, tune::B) where {A<:SMCParticles,B<:SMCTune}
        return new{A,B}(particles, tune)
    end
end

function SMC(
    _rng::Random.AbstractRNG,
    Jitterkernel::J,
    objective::Objective,
    default::SMCDefault=SMCDefault(),
    info::BaytesCore.SampleDefault = BaytesCore.SampleDefault()
) where {
    J<:AbstractConstructor
}
    ## Print assumptions to user
    @unpack chains = info
    @unpack Ntuning, resamplingmethod, resamplingthreshold, jitterfun, jitteradaption, jitterthreshold, jittermin, jittermax, generated = default
    Ndata = maximum(size(objective.data))
    #!NOTE: Check if steps after first jitter call can be captured
    capture = Jitterkernel isa MCMCConstructor ? BaytesCore.UpdateFalse() : BaytesCore.UpdateTrue()
    ## Initialize tune
    tagged = ModelWrappers.Tagged(objective.model, BaytesCore.get_sym(Jitterkernel))
    tune = SMCTune(
        tagged,
        chains,
        Ndata,
        Ntuning,
        resamplingthreshold,
        resamplingmethod,
        capture,
        jitterfun,
        jitteradaption,
        jitterthreshold,
        jittermin,
        jittermax,
        generated
    )
    ## Initialize SMCParticles
    particles = SMCParticles(_rng, Jitterkernel, objective, info, tune)
    ## Return SMC container
    return SMC(particles, tune)
end

############################################################################################
"""
$(SIGNATURES)
Propose new parameter with smc sampler. If update=true, objective function will be updated with input model and data.

# Examples
```julia
```

"""
function propose!(
    _rng::Random.AbstractRNG,
    smc::SMC,
    model::ModelWrapper,
    data::D,
    temperature::F = model.info.reconstruct.default.output(1.0),
    update::U=BaytesCore.UpdateTrue()
) where {D,F<:AbstractFloat,U<:BaytesCore.UpdateBool}
    ## Update kernel parameter values with non-tagged parameter from other sampler
    update!(smc.particles, model, smc.tune.tagged, update)
    ## Set back tune.jitter and update buffer to store correct information for diagnostics
    update!(smc.tune)
    update!(smc.particles.buffer)
    ## If latent θ trajectory increasing - propagate forward
    propagate!(_rng, smc.particles, smc.tune, data, temperature)
    ## Predict new data given current particles
    predict!(_rng, smc.particles, smc.tune, data, temperature)
    ## Adjust weights with log likelihood INCREMENT at time t
    weight!(_rng, smc.particles, smc.tune, data, temperature)
    ## Resample and jitter particles if criterion fulfilled
    #!NOTE: Cannot use resample before propagate! (as in pf) as previous temperature otherwise unknown
    ESS, resampled = resample!(_rng, smc.particles, smc.tune, data, temperature)
    ## Choose proposal model parameter to update model.val
    path = BaytesFilters.draw!(_rng, smc.particles.weights)
    ModelWrappers.fill!(model,
        merge(model.val, ModelWrappers.subset(smc.particles.model[path].val, smc.tune.tagged.parameter))
    )
    ## Pack and return Models and particles
    return model.val,
    SMCDiagnostics(
        smc.particles,
        temperature,
        ESS,
        resampled,
        smc.tune.jitter.Nsteps.current,
        smc.tune.iter.current,
        ModelWrappers.generate(_rng, Objective(model, data, smc.tune.tagged, temperature), smc.tune.generated)
    )
end

############################################################################################
#export
export SMC, SMCDefault, propose!


text:
"""
    MultivariateStdNormal{T, d} <: MultivariateStandardizedDistribution{T, d}
The multivariate standard normal distribution.
"""
struct MultivariateStdNormal{T, d} <: MultivariateStandardizedDistribution{T, d}
    coefs::Vector{T}
end
MultivariateStdNormal{T, d}() where {T, d} = MultivariateStdNormal{T, d}(T[])
MultivariateStdNormal(T::Type, d::Int) = MultivariateStdNormal{T, d}()
MultivariateStdNormal(v::Vector{T}, d::Int) where {T} = MultivariateStdNormal{T, d}()
MultivariateStdNormal{T}(d::Int) where {T} = MultivariateStdNormal{T, d}()
MultivariateStdNormal(d::Int) = MultivariateStdNormal{Float64, d}(Float64[])
rand(rng::AbstractRNG, ::MultivariateStdNormal{T, d}) where {T, d} = randn(rng, T, d)
nparams(::Type{<:MultivariateStdNormal}) = 0
coefnames(::Type{<:MultivariateStdNormal}) = String[]
distname(::Type{<:MultivariateStdNormal}) = "Multivariate Normal"


text:
function SphereGenerator()
    return SphereGenerator(())
end

function ball_on_support(obj::SphereGenerator, arg0::List)
    return jcall(obj, "ballOnSupport", EnclosingBall, (List,), arg0)
end


text:
module Math

include("multiindex.jl")

end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -12.503038214935698 <= q <= 11.138318577023822)

                     add_NL_constraint(m, :(swish(-0.2825986234777189 + -0.06477303738755547*swish(0.91111659897616 + 0.32194568310871174*$(x[1]) + -0.04079743189253193*$(x[2]) + -0.7834407613653229*$(x[3]) + -0.5105294854228863*$(x[4])) + 0.23435836886254213*swish(-0.27698098400696525 + 0.765733032903245*$(x[1]) + 0.6374561618281618*$(x[2]) + -0.006430667300674209*$(x[3]) + -0.6595935184751247*$(x[4])) + 0.7516401137686608*swish(0.7789283911146407 + -0.8853253179157212*$(x[1]) + -0.6613182381252614*$(x[2]) + 0.5007525575333553*$(x[3]) + 0.7960733427996214*$(x[4]))) + swish(0.12692347233027856 + 0.060984281016326136*swish(0.91111659897616 + 0.32194568310871174*$(x[1]) + -0.04079743189253193*$(x[2]) + -0.7834407613653229*$(x[3]) + -0.5105294854228863*$(x[4])) + 0.4767325000240725*swish(-0.27698098400696525 + 0.765733032903245*$(x[1]) + 0.6374561618281618*$(x[2]) + -0.006430667300674209*$(x[3]) + -0.6595935184751247*$(x[4])) + -0.9515049437955176*swish(0.7789283911146407 + -0.8853253179157212*$(x[1]) + -0.6613182381252614*$(x[2]) + 0.5007525575333553*$(x[3]) + 0.7960733427996214*$(x[4]))) + swish(-0.022625917828582143 + -0.9740734403326825*swish(0.91111659897616 + 0.32194568310871174*$(x[1]) + -0.04079743189253193*$(x[2]) + -0.7834407613653229*$(x[3]) + -0.5105294854228863*$(x[4])) + -0.9306328050206538*swish(-0.27698098400696525 + 0.765733032903245*$(x[1]) + 0.6374561618281618*$(x[2]) + -0.006430667300674209*$(x[3]) + -0.6595935184751247*$(x[4])) + 0.61848968844355*swish(0.7789283911146407 + -0.8853253179157212*$(x[1]) + -0.6613182381252614*$(x[2]) + 0.5007525575333553*$(x[3]) + 0.7960733427996214*$(x[4]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
"""
Include a file to your lintContext.
If file does not exist or has already been included, raise a lint warning.
"""
function lintinclude(ctx::LintContext, f::AbstractString)
    if ispath(f)
        # TODO: use paths instead of files
        if f in ctx.included
            msg(ctx, :W356, f, "file included more than once")
        else
            push!(ctx.included, f)
        end

        oldloc = location(ctx)
        str = open(readstring, f)
        location!(ctx, Location(f, 1))

        # TODO: make sure to perform the include at top level
        _lintstr(str, ctx)

        location!(ctx, oldloc)
    else
        msg(ctx, :W357, f, "included file doesn't exist")
    end
end


text:
# This module interface the data retrived from Rath, Alexander.
# “Characterisation of Cell Growth, Metabolism and
# Recombinant Protein Production during Transient and
# Steady State Conditions for the Human Cell Line AGE1.HN-AAT,”
# 2017. https://pure.mpg.de/pubman/item/item_2508673_4.

# The AGE1.HN cells express only marker genes characteristic
# for neuronal cells and lack markers characteristic for glial cells.

module RathData

    import CSV

    using DataFrames
    using Serialization
    using ProjAssistant
    @gen_sub_proj

    include("const.jl")
    include("raw_dat.jl")
    include("proc_dat.jl")
    include("a1at_aa_rel_ab.jl")
    include("data_interface.jl")

    function __init__()
        @create_proj_dirs

        _load_rath_bundle()
        _define_interface()
    end

end


text:
#-------------split-op.jl------------------------------------------------------#
#
# Purpose: A simple split-operator method for solving the schrodinger equation
#          The split operator method involves continuous transformations from
#              real space to phase space, before each tranformation, we multiply
#              by operators that exist either completely in real space or in
#              phase space. This allows us to solve Psi = Psi_0*e^(iHt)
#
#------------------------------------------------------------------------------#

using PyPlot

# Module to hold all parameters for simulation
module Par
    const xmax = 10.0
    const res = 512
    const dt = 0.05
    const timesteps = 1000
    const dx = 2*xmax / res
    const x = -xmax+dx:dx:xmax
    const dk = pi / xmax
    const k = vcat(0:res/2-1, -res/2:-1)*dk
end

# Type to hold all operators
type Operators
    V::Vector{Complex{Float64}}
    PE::Vector{Complex{Float64}}
    KE::Vector{Complex{Float64}}
    wfc::Vector{Complex{Float64}}
end

# Function to initialize the wfc and potential
function init(voffset::Float64, wfcoffset::Float64)
    V = 0.5 * (Par.x - voffset).^2
    wfc = 3* exp(-(Par.x-wfcoffset).^2/2)
    PE = exp(-0.5*im*V*Par.dt)
    KE = exp(-0.5*im*Par.k.^2*Par.dt)

    opr = Operators(V, PE, KE, wfc)
end

# Function for the split-operator loop
function split_op(opr::Operators)

    for i = 1:Par.timesteps
        # Half-step in real space
        opr.wfc = opr.wfc.*opr.PE

        # fft to phase space
        opr.wfc = fft(opr.wfc)

        # Full step in phase space
        opr.wfc = opr.wfc.*opr.KE

        # ifft back
        opr.wfc = ifft(opr.wfc)

        # final half-step in real space
        opr.wfc = opr.wfc.*opr.PE

        # plotting density and potential
        density = abs2(opr.wfc)

        clf()
        plot(density)
        plot(opr.V)
        sleep(0.01)
    end
end

# main function
function main()
    opr = init(2.0,-1.0)
    split_op(opr)
end

main()

max_stars_repo_name: oist/skillpill-julia

max_issues_repo_name: oist/skillpill-julia

max_forks_repo_name: oist/skillpill-julia


text:
struct CallAtChoiceMap{K,T} <: ChoiceMap
    key::K
    submap::T
end

Base.isempty(choices::CallAtChoiceMap) = isempty(choices.submap)

function get_address_schema(::Type{T}) where {T<:CallAtChoiceMap}
    SingleDynamicKeyAddressSchema()
end

function get_submap(choices::CallAtChoiceMap{K,T}, addr::K) where {K,T}
    choices.key == addr ? choices.submap : EmptyChoiceMap()
end

get_submap(choices::CallAtChoiceMap, addr::Pair) = _get_submap(choices, addr)
get_value(choices::CallAtChoiceMap, addr::Pair) = _get_value(choices, addr)
has_value(choices::CallAtChoiceMap, addr::Pair) = _has_value(choices, addr)
get_submaps_shallow(choices::CallAtChoiceMap) = ((choices.key, choices.submap),)
get_values_shallow(::CallAtChoiceMap) = ()

# TODO optimize CallAtTrace using type parameters

struct CallAtTrace <: Trace
    gen_fn::GenerativeFunction # the ChoiceAtCombinator (not the kernel)
    subtrace::Any
    key::Any
end

function Base.getindex(trace::CallAtTrace, addr::Pair)
    first, rest = addr
    if first == trace.key
        return trace.subtrace[rest]
    else
        error("Address prefix $addr not found.")
    end
end
function Base.getindex(trace::CallAtTrace, addr)
    if addr == trace.key
        return trace.subtrace[]
    else
        error("Address $addr not found.")
    end
end

get_args(trace::CallAtTrace) = (get_args(trace.subtrace)..., trace.key)
get_retval(trace::CallAtTrace) = get_retval(trace.subtrace)
get_score(trace::CallAtTrace) = get_score(trace.subtrace)
get_gen_fn(trace::CallAtTrace) = trace.gen_fn

function get_choices(trace::CallAtTrace)
    CallAtChoiceMap(trace.key, get_choices(trace.subtrace))
end

struct CallAtCombinator{T,U,K} <: GenerativeFunction{T, CallAtTrace}
    kernel::GenerativeFunction{T,U}
end

function call_at(kernel::GenerativeFunction{T,U}, ::Type{K}) where {T,U,K}
    CallAtCombinator{T,U,K}(kernel)
end

function accepts_output_grad(gen_fn::CallAtCombinator)
    accepts_output_grad(gen_fn.kernel)
end

unpack_call_at_args(args) = (args[end], args[1:end-1])

function assess(gen_fn::CallAtCombinator, args::Tuple, choices::ChoiceMap)
    (key, kernel_args) = unpack_call_at_args(args)
    if length(collect(get_submaps_shallow(choices))) > 1 || !isempty(get_values_shallow(choices))
        error("Not all constraints were consumed")
    end
    submap = get_submap(choices, key)
    assess(gen_fn.kernel, kernel_args, submap)
end

function propose(gen_fn::CallAtCombinator, args::Tuple)
    (key, kernel_args) = unpack_call_at_args(args)
    (submap, weight, retval) = propose(gen_fn.kernel, kernel_args)
    choices = CallAtChoiceMap(key, submap)
    (choices, weight, retval)
end

function simulate(gen_fn::CallAtCombinator, args::Tuple)
    (key, kernel_args) = unpack_call_at_args(args)
    subtrace = simulate(gen_fn.kernel, kernel_args)
    CallAtTrace(gen_fn, subtrace, key)
end

function generate(gen_fn::CallAtCombinator{T,U,K}, args::Tuple,
                  choices::ChoiceMap) where {T,U,K}
    (key, kernel_args) = unpack_call_at_args(args)
    submap = get_submap(choices, key)
    (subtrace, weight) = generate(gen_fn.kernel, kernel_args, submap)
    trace = CallAtTrace(gen_fn, subtrace, key)
    (trace, weight)
end

function project(trace::CallAtTrace, selection::Selection)
    subselection = selection[trace.key]
    project(trace.subtrace, subselection)
end

function update(trace::CallAtTrace, args::Tuple, argdiffs::Tuple,
                choices::ChoiceMap)
    (key, kernel_args) = unpack_call_at_args(args)
    key_changed = (key != trace.key)
    submap = get_submap(choices, key)
    if key_changed
        (subtrace, weight) = generate(trace.gen_fn.kernel, kernel_args, submap)
        weight -= get_score(trace.subtrace)
        discard = get_choices(trace)
        retdiff = UnknownChange()
    else
        (subtrace, weight, retdiff, subdiscard) = update(
            trace.subtrace, kernel_args, argdiffs[1:end-1], submap)
        discard = CallAtChoiceMap(key, subdiscard)
    end
    new_trace = CallAtTrace(trace.gen_fn, subtrace, key)
    (new_trace, weight, retdiff, discard)
end

function regenerate(trace::CallAtTrace, args::Tuple, argdiffs::Tuple,
                    selection::Selection)
    (key, kernel_args) = unpack_call_at_args(args)
    key_changed = (key != trace.key)
    subselection = selection[key]
    if key_changed
        if !isempty(subselection)
            error("Cannot select addresses under new key $key in regenerate")
        end
        (subtrace, weight) = generate(trace.gen_fn.kernel, kernel_args, EmptyChoiceMap())
        weight -= project(trace.subtrace, EmptySelection())
        retdiff = UnknownChange()
    else
        (subtrace, weight, retdiff) = regenerate(
            trace.subtrace, kernel_args, argdiffs[1:end-1], subselection)
    end
    new_trace = CallAtTrace(trace.gen_fn, subtrace, key)
    (new_trace, weight, retdiff)
end

function choice_gradients(trace::CallAtTrace, selection::Selection, retval_grad)
    subselection = selection[trace.key]
    (kernel_input_grads, value_submap, gradient_submap) = choice_gradients(
        trace.subtrace, subselection, retval_grad)
    input_grads = (kernel_input_grads..., nothing)
    value_choices = CallAtChoiceMap(trace.key, value_submap)
    gradient_choices = CallAtChoiceMap(trace.key, gradient_submap)
    (input_grads, value_choices, gradient_choices)
end

function accumulate_param_gradients!(trace::CallAtTrace, retval_grad)
    kernel_input_grads = accumulate_param_gradients!(trace.subtrace, retval_grad)
    (kernel_input_grads..., nothing)
end

export call_at


text:
using Decs
using DecFP
using Decimals
using BenchmarkTools

# Test performance of comparing < and == of numbers with the same or different scales,
# same or different signs, conversion to/from string, +, -, *, /

for T in (Dec, Decimal, Dec64, Dec128, Float64, BigFloat)

    for (sgn, num, scl) in ((0, 0, 0), (1, 0, 0), (0, 123456789, -8),
                            (0, 123400000, -5), (0, 123400000, 0),
                            (0, 1234, 0), (0, 1234, 5))

        s = string(sgn == 0 ? "" : "-", num) * (scl == 0 ? "" : string('e', scl))
        println(T, ":\t", s)
        n = parse(T, s)
        print("parse: ")
        @btime parse($T, $s)
        print("string: ")
        @btime string($n)

        for vs in (-5, 0, 5)
            sv = scl + vs
            vv = string(sgn == 0 ? "" : "-", num) * (sv == 0 ? "" : string('e', sv))
            println("value: ", vv)
            vn = parse(T, vv)
            print("==:     ")
            @btime $n == $vn
            print("<:      ")
            @btime $n < $vn
            print("+:      ")
            @btime $n + $vn
            print("*:      ")
            @btime $n * $vn
            if !iszero(vn)
                print("/:      ")
                @btime $n / $vn
            end
        end
    end
end


text:
function test_sgd_solver(backend)
    println("-- Testing simple SGD solver call")
    registry_reset(backend)
    srand(12345678)
    ############################################################
    # Prepare Random Data
    ############################################################
    N = 5    # works with arbitrary minibatch size as long as
             # N == batch_size in MemoryDataLayer so it cycles through
             # and gets the same data during forward()
    M = 10
    P = 4

    X = rand(M, N)
    W = rand(M, P)
    B = rand(P, 1)

    Y = (W'*X .+ B)
    Y = Y + 0.01*randn(size(Y))

    ############################################################
    # Define network
    ############################################################
    data_layer = MemoryDataLayer(batch_size=N, data=Array[X,Y])

    w1 = InnerProductLayer(neuron=Neurons.Sigmoid(), name="ip1",output_dim=20, tops=[:a], bottoms=[:data])
    w2 = InnerProductLayer(neuron=Neurons.Identity(), name="ip2",output_dim=4, tops=[:b], bottoms=[:a])
    loss_layer = SquareLossLayer(name="loss", bottoms=[:b, :label] )

    net = Net("TEST", backend, [w1, w2, loss_layer, data_layer])

    # Make a Solver with max iterations 2
    method = Adadelta()
    params = make_solver_parameters(method,
                                    max_iter=2,
                                    beta1=0.9,
                                    beta2=0.999,
                                    epsilon=1e-8,
                                    load_from="")
    solver = Solver(method, params)
    solve(solver, net)
    # TODO check gradient updates
    # TODO check snapshot loading
    # TODO check statistic saving
    # TODO check other lr policies
    # TODO check other mom policies
    destroy(net)
end

if test_cpu
    test_sgd_solver(backend_cpu)
end

if test_gpu
    test_sgd_solver(backend_gpu)
end


text:
input = "9,12,1,4,17,0,18"
invals = parse.(Int,split(input, ","))

using BenchmarkTools

starting = invals
function part1(starting, finish)
    turntracker = Dict{Int,Array{Int}}()
    numbertracker = zeros(Int, 1,finish)
    for (i,v) in enumerate(starting)
        turntracker[v] = [i]
        numbertracker[i] = v
    end
    prevval = starting[end]
    prevturn = length(starting)

    for turn in (length(starting) + 1):finish
        if !haskey(turntracker, prevval)
            error("Could not find $prevval spoken on turn $turn")
        elseif length(turntracker[prevval]) == 1
            # number only spoken on previous turn, and not before
            val = 0
        else
            # number spoken before, need diff between previous two turns where it was spoken
            pt = turntracker[prevval][(end-1):end]
            val = diff(pt)[1]
        end
        if haskey(turntracker, val)
            push!(turntracker[val], turn)
        else
            turntracker[val] = [turn]
        end
        numbertracker[turn] = val
        prevval = val
    end
    return numbertracker, turntracker
end

n, t = part1(starting, 2020)
n[2020]
using Profile
@profview part1(starting, 2020)

@profview part2n, part2t = part1(starting, 30000000)


text:
using Documenter, Projections

ENV["DOCUMENTER_DEBUG"] = "true"

makedocs(
    format = :html,
    sitename = "Projections",
    modules = [Projections]
)

deploydocs(
    repo   = "github.com/Hua-Zhou/Projections.jl.git",
    target = "build"
)


text:
module KymatioWrapper

using DataFrames
using Interpolations
using PyCall

export filter_bank, scattering_filter_factory, Scattering, computeJ, numpyify

function __init__()
    py"""
    import torch
    from kymatio.scattering1d.filter_bank import scattering_filter_factory
    from kymatio import Scattering1D
    from kymatio import Scattering2D
    from kymatio.scattering2d.filter_bank import filter_bank

    torch.cuda.is_available()
    """
end

struct Scattering{D}
    scatter::PyObject
    useGpu::Bool
end

######################### various utils #########################
next2(N) = 2^ceil.(Int, log2(N))
function padTo(x, T;dims=ndims(x))
    szx = size(x)
    toFill = zeros(eltype(x), szx[1:dims - 1]..., T - szx[dims], szx[dims + 1:end]...)
    cat(x, toFill, dims=dims) # assumes the *last* dimension is the one you're using
end
function padTo(x, N, M; dims=ndims(x))
    w = padTo(x, N, dims=ndims(x) - 1)
    padTo(w, M, dims=ndims(x))
end
computeJ(N) = max(floor(Int, min(log2.(N)...) / 2), 3)
"""
Convert any torch pyobjects to numpy (and thus julia) arrays
"""
function numpyify(a::AbstractArray)
    A = copy(a)
    for d in A
        for x in d
            if length(x) == 2 && typeof(x[2]) <: PyObject
                if x[2].is_cuda
                    xp = x[2].clone().detach().cpu()
                else
                    xp = x
                end
                d[x[1]] = (xp.numpy())[:,1]
            end
        end
    end
    return A
end
function numpyify(d::Dict)
    D = copy(d)
    for x in D
        if length(x) == 2 && typeof(x[2]) <: PyObject
            if x[2].is_cuda
                xp = x[2].clone().detach().cpu()
            else
                xp = x[2]
            end
            D[x[1]] = (xp.numpy())[:,1]
        end
    end
    return D
end
################################################################

"""
    ψ₁, ψ₂, ϕ₃, σξ = filter_bank(N::Real, J=floor(log2(N)/2), Q=8)
    ψ, params, ϕₙ = filter_bank(N::Real, J=floor(log2(N)/2), Q=8)

Slightly different format than the native one. Note that it rounds your length
    up to the next power of 2. ψ₁ is (length×nframes), and starts with the
    averaging (phi_f[0] in the original), and increases in frequency, as in
    Wavelets.jl, which is backwards from the native version. Same for ψ₂ with
    phi_f[1]. The final averaging stands alone in ϕ₃, and the parameters are in
    the array of dataframes σξ.
"""

function filter_bank(N::Real, J::Real=computeJ(N), Q::Real=8)
    T = ceil(Int, log2(N))
    phi_f, psi1_f, psi2_f, _ = py"scattering_filter_factory($T,$J,$Q)";
    ψ₁, params1 = extractDictionaries(phi_f, psi1_f, 0)
    ψ₂, params2 = extractDictionaries(phi_f, psi2_f, 1)
    if 2 in keys(phi_f)
        ϕ₃ = upInterp(phi_f[2], 2^2)
    else
        ϕ₃ = zeros(T)
    end
    σξ = [DataFrame(params1'), DataFrame(params2')]
    [rename!(x, [:σ,:j,:ξ]) for x in σξ]
    return (ψ₁, ψ₂, ϕ₃, σξ)
end

function filter_bank(ky::Scattering{1})
    sc = ky.scatter
    phi_f, psi1_f, psi2_f = (numpyify(sc.phi_f), numpyify(sc.psi1_f),
                             numpyify(sc.psi2_f))
    ψ₁, params1 = extractDictionaries(phi_f, psi1_f, 0)
    ψ₂, params2 = extractDictionaries(phi_f, psi2_f, 1)
    if 2 in keys(phi_f)
        ϕ₃ = upInterp(phi_f[2], 2^2)
    else
        ϕ₃ = zeros(2^(sc.J_pad))
    end
    σξ = [DataFrame(params1'), DataFrame(params2')]
    [rename!(x, [:σ,:j,:ξ]) for x in σξ]
    return (ψ₁, ψ₂, ϕ₃, σξ)
end
function filter_bank(ky::Scattering{2})
    N = 641
    filter_bank(N, ky.scatter.J, ky.scatter.L)
end

function filter_bank(N, J=computeJ(N), L=8)
    filters_set = py"filter_bank($(N[1]),$(N[2]),$J,L=$L)"
    nfilters = length(filters_set["psi"]) + 1
    ψ = zeros(N..., nfilters)
    params = -1 .* ones(nfilters, 2)
    ψ[:,:,1] = py"($(filters_set)['phi'][0][..., 0])" .* (py"($(filters_set)['phi'][0][..., 0])")'
    params[1,1] = filters_set["psi"][1]["j"]
    params[1,2] = NaN
    for ii = 1:(nfilters - 1)
        ψ[:,:,nfilters - ii + 1] = py"($(filters_set)['psi'][$(ii-1)][0][..., 0])" .* (py"($(filters_set)['psi'][$(ii-1)][0][..., 0])")'
        params[nfilters - ii + 1, 1] = filters_set["psi"][ii]["j"]
        params[nfilters - ii + 1, 2] = filters_set["psi"][ii]["theta"]
    end
    params = DataFrame(params)
    rename!(params, [:j,:θ])
    ϕₙ = Dict()
    for ii = 1:J
        ϕₙ[ii] = py"($(filters_set)['phi'][$ii-1][..., 0])"
    end
    return ψ, params, ϕₙ
end

scattering_filter_factory(N,J,Q) = filter_bank(N, J, Q)

function extractDictionaries(phi_f, psi_f, m)
    nfreq = length(psi_f) + 1
    ψ = zeros(length(psi_f[1][0]), nfreq)
    params = zeros(3, nfreq)
    ψ[:,1] = upInterp(phi_f[m], 2^m)
    params[1,1] = phi_f["sigma"]
    params[2,1] = phi_f["j"]
    params[3,1] = phi_f["xi"]
    for (ii, ψi) in enumerate(psi_f)
        ψ[:, nfreq - ii + 1] .= ψi[0]
        params[1,nfreq - ii + 1] = ψi["sigma"]
        params[2,nfreq - ii + 1] = ψi["j"]
        params[3,nfreq - ii + 1] = ψi["xi"]
    end
    return (ψ, params)
end

"""
    upInterp(ϕ,k)
given a vector of values, do a periodic cubic interpolation with k times as
    many entries
"""
function upInterp(ϕ::AbstractArray{T,1}, k) where T
    itr = interpolate(ϕ, BSpline(Cubic(Periodic(OnGrid()))))
     x = range(1, length(ϕ), length=k * length(ϕ))
    return itr(x)
end
upInterp(ϕ::AbstractArray{T,2}, k) where T = upInterp(ϕ[:,1], k + 1)

# 1D version
function Scattering(N::Int64; J=computeJ(N), Q=8,
                    max_order=2, average=true, oversampling=0,
                    vectorize=true, useGpu=true)
    if J < 2.13
        @warn "silently does the wrong thing if J<2.13."
    end
    T = N
    scatterer = py"Scattering1D($J, $T, $Q, $(max_order), $(average), $(oversampling), $(vectorize),backend='torch')"
    if useGpu
        cuVersion = scatterer.cuda()
    end
    Scattering{1}(scatterer, useGpu)
end

"""
    compute_meta_scattering(J::Int,Q,max_order=2)
    compute_meta_scattering(s::Scattering{1})
what's on the tin.
"""
function compute_meta_scattering(J::Int, Q, max_order=2)
    dictMeta = py"Scattering1D.compute_meta_scattering($(J), $(Q), $(max_order))"
    dictMeta
    meta = DataFrame(j1=Float64[], j2=Float64[], key=Tuple{Vararg{Int64}}[],
                    σ1=Float64[], σ2=Float64[], order=Int[], ξ1=Float64[],
                    ξ2=Float64[], n1=Int[], n2=Int[])
    σ = dictMeta["sigma"].numpy()
    σ = DataFrame(σ₁=σ[:,1], σ₂=σ[:,2])
    order = DataFrame(order=dictMeta["order"].numpy())
    n = dictMeta["n"].numpy()
    n = DataFrame(n1=n[:,1], n2=n[:,2])
    j = dictMeta["j"].numpy()
    j = DataFrame(j1=j[:,1], j2=j[:,2])
    key = DataFrame(key=dictMeta["key"])
    return hcat(σ, n, j, key, order)
end
function compute_meta_scattering(s::Scattering{1})
    compute_meta_scattering(s.scatter.J, s.scatter.Q, s.scatter.max_order)
end

function (s::Scattering{1})(x)
    x = Float32.(permutedims(x, (2:ndims(x)..., 1)))
    flipped = true
    if s.useGpu
        res = s.scatter.forward(py"torch.from_numpy($x).cuda().contiguous()")
    else
        res = s.scatter.forward(py"torch.from_numpy($x).contiguous()")
    end
    res = res.cpu().numpy()
    n = ndims(res)
    if flipped
        n = ndims(res)
        return permutedims(res, (n, n - 1, 1:(n - 2)...))
    else
        return res
    end
end
################################################################
# 2D version
################################################################

function Scattering(N::Tuple{Int64,Int64}; J=computeJ(N),
                    L=8, max_order=2, useGpu=true)
    # if you're using Julia size conventions, switch that around
    T = next2.(N)
    scatterer = py"Scattering2D($J, shape=$T, L=$L, max_order=$max_order)"
    if useGpu
        scatterer.cuda()
    end
    Scattering{2}(scatterer, useGpu)
end

function (s::Scattering{2})(x)
    @assert maximum(s.scatter.shape .% 2 .== 0) "s has shape = $(s.scatter.shape)"
    flipped = false
    if next2(size(x, 1)) == s.scatter.shape[1]
        x = Float32.(permutedims(x, ((3:ndims(x))..., 1, 2)))
        flipped = true
    end
    x = padTo(x, s.scatter.shape...)
    torchX =
    if s.useGpu
        res = s.scatter.forward(py"torch.from_numpy($x).cuda()")
    else
        res = s.scatter.forward(py"torch.from_numpy($x).contiguous()")
    end

    res = res.cpu().numpy()
    if flipped
        n = ndims(res)
        return permutedims(res, (n, n - 1, n - 2, 1:(n - 3)...))
    else
        return res
    end
end

end


text:
module ASHBenchmark

using AverageShiftedHistograms, BenchmarkTools
versioninfo()

# compile function
ash(randn(10), linspace(-4, 4, 100), m = 3)

x = randn(1_000_000)
@benchmark ash(x, linspace(-6, 6, 1000), m = 5)
end #module


text:
using IcanProj
using PyPlot

opt = get_options()

wsh_info = get_wsh_info()

# for i in 1:length(wsh_info)

i = 1

fig = figure("pyplot_histogram",figsize=(7,5))

h = plt[:hist](wsh_info[i].elev[:], 20, alpha = 0.3)

xlabel("Elevation (m)")
ylabel("Number of gridcells")

# end