text:
mutable struct Context <: Handle
    handle
    function Context()
        ctx_ptr = null_ptr_ref()
        @check spvc_context_create(ctx_ptr)
        ctx = new(ctx_ptr[])
        finalizer(spvc_context_destroy, ctx)
    end
end

struct IR <: Handle
    handle
end

struct Compiler <: Handle
    handle
end

struct ShaderResources <: Handle
    handle
end

struct ReflectedResource
    id::Int
    base_type::Int
    type::Int
    name::String
end

Base.convert(T::Type{ReflectedResource}, x::spvc_reflected_resource) = T(x.id, x.base_type_id, x.type_id, unsafe_string(x.name))

function IR(ctx::Context, file)
    parsed_ir_ptr = null_ptr_ref()

    filesize = stat(file).size
    code = Array{UInt8}(undef, cld(filesize, 4))
    open(file) do io
        readbytes!(io, code, filesize)
    end
    code_final = reinterpret(UInt32, code)

    @check spvc_context_parse_spirv(ctx, code, length(code_final), parsed_ir_ptr)

    IR(parsed_ir_ptr[])
end

function Compiler(ctx::Context, ir::IR; backend=SPVC_BACKEND_GLSL, capture_mode=SPVC_CAPTURE_MODE_TAKE_OWNERSHIP)
    compiler_ptr = null_ptr_ref()
    @check spvc_context_create_compiler(ctx, backend, ir, capture_mode, compiler_ptr)
    Compiler(compiler_ptr[])
end

function ShaderResources(compiler::Compiler)
    resources_ptr = null_ptr_ref()
    @check spvc_compiler_create_shader_resources(compiler, resources_ptr)
    ShaderResources(resources_ptr[])
end

function get_resource_list(resources::ShaderResources; resource_type=SPVC_RESOURCE_TYPE_UNIFORM_BUFFER)
    list = Ref{Ptr{spvc_reflected_resource}}()
    count = Ref{Int32}()
    @check spvc_resources_get_resource_list_for_type(resources, resource_type, list, count)
    # vec = Vector{spvc_reflected_resource}(undef, count[])
    # @check spvc_resources_get_resource_list_for_type(resources, resource_type, vec, count)
    convert.(ReflectedResource, unsafe_wrap(Array, list[], (count[],)))
end


text:
# Implements trust region code

using LinearAlgebra, SparseArrays

@doc raw"""
    dampingInitialization(B, g, delta, epsilon, a = 0.5)

Finds an appropriate damping parameter to start the search.

Specifically,
```math
λ ∈ [0, ||B|| + (1 + ϵ) \frac{||g||}{Δ}]
```

## Arguments
- `B`: The (approximate or true) Hessian Matrix
- `g`: The (approximate or true) gradient vector
- `epsilon`: Parameter (> 0) weighes the gradient versus the hessian in
determining the damping ratio size
- `a`: Parameter (> 0) How quickly to damp the Hessian

See also: [`findDamping`](@ref)
"""
function dampingInitialization(B, g, delta, epsilon, IPD, a = 0.5,
                                    basedamping = 1e-6)

    dampingMax = (norm(B) + ((1 + epsilon) * norm(g) / delta))

    if isposdef(B)
        damping = basedamping
        Bdamped = B + damping * IPD
        cholBdamp = cholesky(Bdamped, check = false)

        if issuccess(cholBdamp)
            return damping, dampingMax, cholBdamp
        end
    end

    while true
        damping = a * dampingMax

        Bdamped = B + damping * IPD

        # Get the Cholesky decomposition without throwing an exception
        luB = lu(Bdamped, check = false)

        # Checks positive definite and non-singular
        if issuccess(luB)
            return damping, dampingMax, Bdamped
        else
            # If the it fails, increase damping factor. Converge to a -> 1
            a = (a + 1)/ 2
            println("More Damping...")

            # Once "a" gets close to 1, might as well stop
            if a > 0.99
                return dampingMax, dampingMax, Bdamped
            end
        end
    end
end

"""
    findDamping(B, g, delta = 1, gamma = 1.5, epsilon = 0.5,
                      condNumMax = 1e7, a = 1e-6, verbose = false)

Implements a search for the appropriate trust region.

## Arguments:
- `B`: The (approximate or true) Hessian Matrix
- `g`: The (approximate or true) gradient vector
- `delta`: is the initial trust region size
- `condNumMax`: The max condition number on the damped Hessian

for `gamma`, `epsilon`, and `a`, see [`dampingInitialization`](@ref)

returns the damping factor `damping`, the newton step `dk`, and the damped
hessian represented as a Cholesky Decomposition `rCho`

This corresponds to algorithm 2.6 in Nocedal et Yuan (1998)

"""
function findDamping(B, g, dualSize = 0, delta = 1, gamma = 1.5, epsilon = 0.5,
                        condNumMax = 1e7, a = 1e-6, verbose = false)

    # For Primal-Dual the effective Identity is [I 0; 0 -I] (albeit with
    # correct dimensions)
    getIPD(p, d) = Diagonal([ones(p); -ones(d)])
    IPD = getIPD(size(g, 1) - dualSize, dualSize)

    # Initialize the damping factor and find the initial damped hessian
    # represented as a Cholesky Decomposition (rCho)
    # damping, dampingMax, rCho = dampingInitialization(B, g, delta, epsilon,
    #                                                   IPD, a)
    damping, dampingMax, BDamp = dampingInitialization(B, g, delta, epsilon,
                                                      IPD, a)

    # Use the Cholesky decomposition to more easily solve the problem.
    # LCho = sparse(rCho.L)
    # dk = - LCho' \ (LCho \ g)
    dk = - BDamp \ g

    if verbose
        println("Damping Start = $damping and dampingMax = $dampingMax")
        println("dk = $dk")
    end

    counter = 1

    # Have to use Array to convert briefly to a dense matrix in order
    # to determine the condition number.
    # while (norm(dk) > delta) && (cond(Array(LCho * LCho')) > condNumMax)
    while (norm(dk) > delta) && (cond(Array(BDamp)) > condNumMax)
        if verbose
            println("Increasing Damping")
            println("Cholesky Factorization:")
            display(rCho)
            println("dk")
            display(dk)
        end

        LCho = lu(BDamp).L

        qk = LCho \ dk # Problematic

        updateNumerator = norm(dk)^2 * (gamma * norm(dk) - delta)
        updateDenominator = norm(qk)^2 * delta

        damping = min(damping + updateNumerator / updateDenominator, dampingMax)

        if verbose
            println("qk = $qk")
            println("New Damping: $damping")
        end

        # rCho = lu(B + damping * IPD, check = false)
        BDamp = B + damping * IPD
        luB = lu(BDamp, check = false)

        if issuccess(luB)
            dk = - BDamp \ g
            counter += 1

            if damping == dampingMax
                break
            end
        else
            # Need to really increase the damping factor
            damping *= 10
        end

        if counter >= 10
            break
        end
    end

    println("Condition Number: $(cond(Array(BDamp)))")

    return damping, dk, BDamp

end

ext: jl

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
### Some prime group generation algorithms. References:
# + https://crypto.stackexchange.com/questions/820/how-does-one-calculate-a-primitive-root-for-diffie-hellman
# + https://math.stackexchange.com/questions/124408/finding-a-primitive-root-of-a-prime-number

# Eventually will need to revisit this implmentation to comply with CryptoGroups standart

using Primes: isprime, nextprime
using Random: AbstractRNG

function n_bit_random_number(rng::AbstractRNG, len::Integer)
    max_n = ( BigInt(1) << len ) - 1
    if len > 2
        min_n = BigInt(1) << (len - 1)
        return rand(rng, min_n:max_n)
    end
    return rand(rng, 1:max_n)
end

function nbit_prime_of_size(rng::AbstractRNG, n_bits::Integer)
    # generate a random nbit number
    r = n_bit_random_number(rng, n_bits)
    return nextprime(r)
end

rngprime(rng::AbstractRNG, N) = nbit_prime_of_size(rng, N)
rngint(rng::AbstractRNG, N) = n_bit_random_number(rng, N)

### I could manage to get 512 security
# Something wrong with this one
function sophie_germain_group(rng::AbstractRNG, g::Integer, t::Int)
    @assert g!=1
    while true
        q = rngprime(rng, 2*t)
        p = 2*q + 1
        if g!=p-1 && isprime(p)
            return MODP(; p, g, q)
        end
    end
end

function generate_qp(rng::AbstractRNG, t)
    while true
        q = rngprime(rng, t)
        r = rngint(rng, t)

        p = q*r + 1
        if isprime(p)
            return q, p
        end
    end
end

function dsa_standart_group(rng::AbstractRNG, tp::Int, tg::Int)
    q, p = generate_qp(rng, 2*tp)
    while true
        u = mod(rngint(rng, tg), p)
        g = powermod(u, div((p-1), q), p)

        if g!=1 || g!=0
            return MODP(; p, g, q)
        end
    end
end


text:
#NPV TEST MODEL FOR MC HAMMER
#MAY 2019
clearconsole()

using Distributions
using Dates
using Gadfly
using StatsBase
using MCHammer
using DataFrames

#Setup the Date Range for the analysis
dr = collect(Date(2019,1,01):Dates.Year(1):Date(2023,01,01))

#Setup Global Inputs
ForecastYrs = 5
Units = [5000, 10000, 17000, 18000, 20000]
InitialInvestment = 250000
Investment = [100000, 0, 0, 25000,0] #fill(0,ForecastYrs)

#Setup Outputs
Sensitivity_Tbl = []

ProjectNPV = []
USP = []
USC =[]
DR = []
OP =[]
Annual_CashFlows =[]

for i = 1:10000

    UnitSellPrice = GBMM(80, 0.2, 0.1, ForecastYrs)
    UnitCost = GBMM(40, 0.1, 0.05, ForecastYrs)

    #Each period the discount rate is independent. If you use an additive method instead of multiplicative, you can end up with differences. These may or may not impact the decision. For simulation it is best to use the risk free rate.

        #Multiplicative Method
        DiscountRate = cumprod(rand(Normal(0.02,0.0075),ForecastYrs)+fill(1,ForecastYrs))#accumulate(+,rand(Normal(0.02,0.0075),ForecastYrs))+fill(1,ForecastYrs)

        #Additive Method
        #DiscountRate = accumulate(+,rand(Normal(0.02,0.0075),ForecastYrs))+fill(1,ForecastYrs)

        #a static DR
        #DiscountRate = accumulate(+, fill(0.02,ForecastYrs))+fill(1,ForecastYrs)

    #print(DiscountRate)

    #DCF Elements
    Annual_Sales = UnitSellPrice .* Units
    Annual_COGS = UnitCost .* Units
    OPEX = rand(TriangularDist(.2,0.5,0.35),ForecastYrs) .* Annual_Sales

    #Constant Dollar Cashflow
    #CashFlow_C = (Annual_Sales - Annual_COGS - OPEX - Investment)

    #Discounted CashFLow
    CashFlow = (Annual_Sales - Annual_COGS - OPEX - Investment) ./ DiscountRate

    Trial_NPV = sum(CashFlow)-InitialInvestment

#Setup output tables
    push!(ProjectNPV, Trial_NPV)
    push!(USC, mean(UnitCost))
    push!(USP, mean(UnitSellPrice))
    push!(DR,  mean(DiscountRate))
    push!(OP,  mean(OPEX))
    push!(Annual_CashFlows,  CashFlow)

end
Sensitivity_Tbl = DataFrame(hcat(ProjectNPV, USC, USP, DR, OP))
rename!(Sensitivity_Tbl, [:ProjectNPV, :USC, :USP, :DR, :OP])
NPV_Sensitivity = cormat(Sensitivity_Tbl,1)

#Stats

plot(ProjectNPV, x=ProjectNPV, Geom.density)
print("Project Mean: ", mean(ProjectNPV),"\n")
print("Project Std.Dev: ", std(ProjectNPV),"\n")
print("Prob. of Neg. NPV: ", GetCertainty(ProjectNPV,0,0),"\n")
print("NPV p10, p50, p90 : ", quantile(collect(Float64, ProjectNPV),[0.1,0.5,0.9]),"\n")
println("")
println("OUTPUTS: Annual_CashFlows, ProjectNPV, Sensitivity_Tbl")
println("date range = dr")


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    PlanePartition(normal; tol=1e-6)

A method for partitioning spatial data into a family of hyperplanes defined
by a `normal` direction. Two points `x` and `y` belong to the same
hyperplane when `(x - y) ⋅ normal < tol`.
"""
struct PlanePartition{T,N} <: SPredicatePartitionMethod
  normal::SVector{N,T}
  tol::Float64

  function PlanePartition{T,N}(normal, tol) where {N,T}
    new(normalize(normal), tol)
  end
end

PlanePartition(normal::SVector{N,T}; tol=1e-6) where {T,N} =
  PlanePartition{T,N}(normal, tol)

PlanePartition(normal::NTuple{N,T},; tol=1e-6) where {T,N} =
  PlanePartition(SVector(normal), tol=tol)

(p::PlanePartition)(x, y) = abs((x - y) ⋅ p.normal) < p.tol


text:
using GraphSignals

"""
    AbstractTrajectoryState

The `TrajectoryState` is a component of the `AbstractStateRepresentation` specifying how to store
this representation in the `trajectory` of the RL agent for learning.

It might be necessary to implement your own when defining a new state representation or
when another data formatting is necessary for learning. To implement a subtype of
`AbstractTrajectoryState` one must provide:
- a struct subtype of `TabularTrajectoryState` or `NonTabularTrajectoryState` depending on you usage.
- a constructor with signature `YourTrajectoryState(sr::YourStateRepresentation)`.
- a function `Flux.functor(::Type{YourTrajectoryState}, ts)` (see Flux.functor for further information).
- an NN model taking a `YourTrajectoryState` instance as an input.
"""
abstract type AbstractTrajectoryState end
Flux.functor(t::Type{<:AbstractTrajectoryState}, ts) = throw(ErrorException("missing function Flux.functor($(t), ::Any)."))

"""
    TabularTrajectoryState

Abstract subtype of `AbstractTrajectoryState` for array based representations of state.

This type isn't currently used, but it will certainly be a lot clearer if one can easily
distinguish between array based representations, and more exotic ones (e.g. graphs, named tuples...).
"""
abstract type TabularTrajectoryState <: AbstractTrajectoryState end

"""
    NonTabularTrajectoryState

Abstract subtype of `AbstractTrajectoryState` for non array based representations of state.

This type is th only one currently used, but it will certainly be a lot clearer if one can easily
distinguish between array based representations, and more exotic ones (e.g. graphs, named tuples...).
"""
abstract type NonTabularTrajectoryState <: AbstractTrajectoryState end
abstract type GraphTrajectoryState <: NonTabularTrajectoryState end

"""
    AbstractStateRepresentation{TS}

The AbstractStateRepresentation is the abstract type of the structures representing the internal state of the
CPModel and eventually the state of the search in a way that will be as expressive as possible.

It requires a specific `AbstractTrajectoryState` to make the bridge between the internal state representation
and one the RL agent can use to decide of the value to be assigned when branching.

A user can use the `DefaultStateRepresentation` with the `DefaultTrajectoryState` provided by the package but
he has the possibility to define his own.

To define a new one, the user must provide:
- a new structure, subtype of `AbstractStateRepresentation` with the appropriate subtype of `AbstractTrajectoryState` and all its dependencies.
- a constructor from a `CPModel`, with keyword argument `action_space`.
- a function `update_representation!(::AbstractStateRepresentation, ::CPModel, ::AbstractIntVar)` to update the state representation at each step.

Look at the DefaultStateRepresentation to get inspired.
"""
abstract type AbstractStateRepresentation{TS <: AbstractTrajectoryState} end
update_representation!(sr::AbstractStateRepresentation, m::CPModel, x::AbstractIntVar) = throw(ErrorException("missing function update_representation!(::$(typeof(sr)), ::$(typeof(m)), ::$(typeof(x)))."))

"""
    trajectoryState(sr::AbstractStateRepresentation{TS})

Return a TrajectoryState based on the present state represented by `sr`.

The type of the returned object is defined by the `TS` parametric type defined in `sr`.
"""
trajectoryState(sr::AbstractStateRepresentation{TS}) where {TS} = TS(sr)

"""
    AbstractFeaturization

Every subtype of `FeaturizedStateRepresentation{F}` requires an `AbstractFeaturization`.

This type gives the possibility to characterise these feature based representations, and thus gives
the ability to easily define new ones. To implement your own featurization, one must provide:
- a struct subtype of `AbstractFeaturization` (no field is used by the `DefaultStateRepresentation`).
- a function `featurize(::FeaturizedStateRepresentation{YourFeaturization, TS}) where TS` returning a feature Matrix with features stored columnwise.
- a function `feature_length(::Type{<:FeaturizedStateRepresentation{YourFeaturization, TS} where TS}` returning the size of a feature vector.
- _optionally_ a function `update_features!(::FeaturizedStateRepresentation{YourFeaturization, TS}, ::CPModel) where TS` to update your feature vectors based on statistics gathered by the CP Model.
"""
abstract type AbstractFeaturization end

"""
    FeaturizedStateRepresentation{F, TS}

Abstract subtype of `AbstractStateRepresentation` specializing the type of features used in the representation.

When a user wants to try a new featurization with the same organisation of the featurized elements, instead of
having to completely redefine a new type of AbstractStateRepresentation, he can keep the same and
just use a new `AbstractFeaturization`.
"""
abstract type FeaturizedStateRepresentation{F <: AbstractFeaturization, TS} <: AbstractStateRepresentation{TS} end

featurize(sr::FeaturizedStateRepresentation) = throw(ErrorException("missing function featurize(::$(typeof(sr)))."))
function global_featurize(sr::FeaturizedStateRepresentation) end
function update_features!(::FeaturizedStateRepresentation, ::CPModel) end
feature_length(sr::Type{<:FeaturizedStateRepresentation}) = throw(ErrorException("missing function feature_length(::$(sr))."))
global_feature_length(sr::Type{<:FeaturizedStateRepresentation}) = 0

include("default/defaulttrajectorystate.jl")
include("default/defaultstaterepresentation.jl")
include("tsptw/tsptwstaterepresentation.jl")


text:
module GDB

using Distributed

# @everywhere begin
using DistributedArrays
using SharedArrays
using LinearAlgebra
using IterativeSolvers
using SparseArrays
using HDF5
using ForwardDiff
using RecipesBase
using Logging
using ProgressMeter

include("./domain/grid.jl")
include("./domain/barrier.jl")
include("./domain/interpolation.jl")
include("./domain/operators.jl")
include("./domain/ghosts.jl")

include("./linsolve/gmres.jl")
include("./linsolve/relaxation.jl")
include("./linsolve/solvers.jl")

include("./magflux/trace.jl")
include("./magflux/maps.jl")
include("./magflux/annalus.jl")
include("./magflux/solovev.jl")

include("./physics/formulas.jl")
include("./physics/timestep.jl")

# end

end


text:
immutable Body3D <: Body

end


text:
struct SRIConstantCache{VType1,VType2,MType,uType} <: StochasticDiffEqConstantCache
  c₀::VType1
  c₁::VType1
  A₀::MType
  A₁::MType
  B₀::MType
  B₁::MType
  α::VType2
  β₁::VType2
  β₂::VType2
  β₃::VType2
  β₄::VType2
  stages::Int
  H0::Vector{uType}
  H1::Vector{uType}
  error_terms::Int
end

function SRIConstantCache(tableau,rate_prototype,error_terms)
  @unpack c₀,c₁,A₀,A₁,B₀,B₁,α,β₁,β₂,β₃,β₄ = tableau
  stages = length(α)
  H0 = Array{typeof(rate_prototype)}(undef,stages)
  H1 = Array{typeof(rate_prototype)}(undef,stages)
  SRIConstantCache(c₀,c₁,A₀',A₁',B₀',B₁',α,β₁,β₂,β₃,β₄,stages,H0,H1,error_terms)
end

function alg_cache(alg::SRI,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  SRIConstantCache(alg.tableau,rate_prototype,alg.error_terms)
end

@cache struct SRICache{randType,uType,rateType,tabType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  H0::Vector{uType}
  H1::Vector{uType}
  A0temp::rateType
  A1temp::rateType
  B0temp::rateType
  B1temp::rateType
  A0temp2::rateType
  A1temp2::rateType
  B0temp2::rateType
  B1temp2::rateType
  atemp::rateType
  btemp::rateType
  E₁::rateType
  E₂::rateType
  E₁temp::rateType
  ftemp::rateType
  gtemp::rateType
  chi1::randType
  chi2::randType
  chi3::randType
  tmp::uType
  tab::tabType
end

function alg_cache(alg::SRI,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  H0 = Vector{typeof(u)}()
  H1 = Vector{typeof(u)}()
  tab = SRIConstantCache(alg.tableau,rate_prototype,alg.error_terms)
  for i = 1:tab.stages
    push!(H0,zero(u))
    push!(H1,zero(u))
  end
  #TODO Reduce memory
  A0temp = zero(rate_prototype); A1temp = zero(rate_prototype)
  B0temp = zero(rate_prototype); B1temp = zero(rate_prototype)
  A0temp2 = zero(rate_prototype); A1temp2 = zero(rate_prototype)
  B0temp2 = zero(rate_prototype); B1temp2 = zero(rate_prototype)
  atemp = zero(rate_prototype); btemp = zero(rate_prototype)
  E₁ = zero(rate_prototype); E₂ = zero(rate_prototype); E₁temp = zero(rate_prototype)
  ftemp = zero(rate_prototype); gtemp = zero(rate_prototype)

  if typeof(ΔW) <: Union{SArray,Number}
    chi1 = copy(ΔW)
    chi2 = copy(ΔW)
    chi3 = copy(ΔW)
  else
    chi1 = zero(ΔW)
    chi2 = zero(ΔW)
    chi3 = zero(ΔW)
  end
  tmp = zero(u)
  SRICache(u,uprev,H0,H1,A0temp,A1temp,B0temp,B1temp,A0temp2,A1temp2,B0temp2,B1temp2,
    atemp,btemp,E₁,E₂,E₁temp,ftemp,gtemp,chi1,chi2,chi3,tmp,tab)
end

struct SRIW1ConstantCache <: StochasticDiffEqConstantCache end
alg_cache(alg::SRIW1,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}}) = SRIW1ConstantCache()

@cache struct SRIW1Cache{randType,uType,rateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  chi1::randType
  chi2::randType
  chi3::randType
  fH01o4::rateType
  g₁o2::rateType
  H0::uType
  H11::uType
  H12::uType
  H13::uType
  g₂o3::rateType
  Fg₂o3::rateType
  g₃o3::rateType
  Tg₃o3::rateType
  mg₁::rateType
  E₁::rateType
  E₂::rateType
  fH01::rateType
  fH02::rateType
  g₁::rateType
  g₂::rateType
  g₃::rateType
  g₄::rateType
  tmp::uType
end

function alg_cache(alg::SRIW1,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    chi1 = copy(ΔW)
    chi2 = copy(ΔW)
    chi3 = copy(ΔW)
  else
    chi1 = zero(ΔW)
    chi2 = zero(ΔW)
    chi3 = zero(ΔW)
  end
  fH01o4 = zero(rate_prototype)
  g₁o2 = zero(rate_prototype)
  H0 = zero(u)
  H11 = zero(u)
  H12 = zero(u)
  H13 = zero(u)
  g₂o3 = zero(rate_prototype)
  Fg₂o3 = zero(rate_prototype)
  g₃o3 = zero(rate_prototype)
  Tg₃o3 = zero(rate_prototype)
  mg₁ = zero(rate_prototype)
  E₁ = zero(rate_prototype)
  E₂ = zero(rate_prototype)
  fH01 = zero(rate_prototype); fH02 = zero(rate_prototype)
  g₁ = zero(rate_prototype); g₂ = zero(rate_prototype); g₃ = zero(rate_prototype); g₄ = zero(rate_prototype)
  tmp = zero(u)
  SRIW1Cache(u,uprev,chi1,chi2,chi3,fH01o4,g₁o2,H0,H11,H12,H13,g₂o3,Fg₂o3,g₃o3,Tg₃o3,mg₁,E₁,E₂,fH01,fH02,g₁,g₂,g₃,g₄,tmp)
end

struct FourStageSRIConstantCache{T,T2} <: StochasticDiffEqConstantCache
  a021::T
  a031::T
  a032::T
  a041::T
  a042::T
  a043::T
  a121::T
  a131::T
  a132::T
  a141::T
  a142::T
  a143::T
  b021::T
  b031::T
  b032::T
  b041::T
  b042::T
  b043::T
  b121::T
  b131::T
  b132::T
  b141::T
  b142::T
  b143::T
  α1::T
  α2::T
  α3::T
  α4::T
  c02::T2
  c03::T2
  c04::T2
  c11::T2
  c12::T2
  c13::T2
  c14::T2
  beta11::T
  beta12::T
  beta13::T
  beta14::T
  beta21::T
  beta22::T
  beta23::T
  beta24::T
  beta31::T
  beta32::T
  beta33::T
  beta34::T
  beta41::T
  beta42::T
  beta43::T
  beta44::T
end

function SRIW2ConstantCache(T::Type, T2::Type)
  a021 = convert(T, 1)
  a031 = convert(T, 1//4)
  a032 = convert(T, 1//4)
  a041 = convert(T, 0)
  a042 = convert(T, 0)
  a043 = convert(T, 0)
  a121 = convert(T, 1//4)
  a131 = convert(T, 1)
  a132 = convert(T, 0)
  a141 = convert(T, 0)
  a142 = convert(T, 0)
  a143 = convert(T, 1//4)
  b021 = convert(T, 0)
  b031 = convert(T, 1)
  b032 = convert(T, 1//2)
  b041 = convert(T, 0)
  b042 = convert(T, 0)
  b043 = convert(T, 0)
  b121 = convert(T, -1//2)
  b131 = convert(T, 1)
  b132 = convert(T, 0)
  b141 = convert(T, 2)
  b142 = convert(T, -1)
  b143 = convert(T, 1//2)
  α1 = convert(T, 1//6)
  α2 = convert(T, 1//6)
  α3 = convert(T, 2//3)
  α4 = convert(T, 0)
  c02 = convert(T2, 1)
  c03 = convert(T2, 1//2)
  c04 = convert(T2, 0)
  c11 = convert(T2, 0)
  c12 = convert(T2, 1//4)
  c13 = convert(T2, 1)
  c14 = convert(T2, 1//4)
  beta11 = convert(T, -1)
  beta12 = convert(T, 4//3)
  beta13 = convert(T, 2//3)
  beta14 = convert(T, 0)
  beta21 = convert(T, 1)
  beta22 = convert(T, -4//3)
  beta23 = convert(T, 1//3)
  beta24 = convert(T, 0)
  beta31 = convert(T, 2)
  beta32 = convert(T, -4//3)
  beta33 = convert(T, -2//3)
  beta34 = convert(T, 0)
  beta41 = convert(T, -2)
  beta42 = convert(T, 5//3)
  beta43 = convert(T, -2//3)
  beta44 = convert(T, 1)
  FourStageSRIConstantCache(a021,a031,a032,a041,a042,a043,a121,a131,a132,a141,
  a142,a143,b021,b031,b032,b041,b042,b043,b121,b131,b132,b141,b142,b143,α1,α2,
  α3,α4,c02,c03,c04,c11,c12,c13,c14,beta11,beta12,beta13,beta14,beta21,beta22,
  beta23,beta24,beta31,beta32,beta33,beta34,beta41,beta42,beta43,beta44)
end

function alg_cache(alg::SRIW2,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  SRIW2ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function SOSRIConstantCache(T::Type, T2::Type)
  a021 = convert(T, -0.04199224421316468)
  a031 = convert(T, 2.842612915017106)
  a032 = convert(T, -2.0527723684000727)
  a041 = convert(T, 4.338237071435815)
  a042 = convert(T, -2.8895936137439793)
  a043 = convert(T, 2.3017575594644466)
  a121 = convert(T, 0.26204282091330466)
  a131 = convert(T, 0.20903646383505375)
  a132 = convert(T, -0.1502377115150361)
  a141 = convert(T, 0.05836595312746999)
  a142 = convert(T, 0.6149440396332373)
  a143 = convert(T, 0.08535117634046772)
  b021 = convert(T, -0.21641093549612528)
  b031 = convert(T, 1.5336352863679572)
  b032 = convert(T, 0.26066223492647056)
  b041 = convert(T, -1.0536037558179159)
  b042 = convert(T, 1.7015284721089472)
  b043 = convert(T, -0.20725685784180017)
  b121 = convert(T, -0.5119011827621657)
  b131 = convert(T, 2.67767339866713)
  b132 = convert(T, -4.9395031322250995)
  b141 = convert(T, 0.15580956238299215)
  b142 = convert(T, 3.2361551006624674)
  b143 = convert(T, -1.4223118283355949)
  α1 = convert(T, 1.140099274172029)
  α2 = convert(T, -0.6401334255743456)
  α3 = convert(T, 0.4736296532772559)
  α4 = convert(T, 0.026404498125060714)
  c02 = convert(T2, -0.04199224421316468)
  c03 = convert(T2, 0.7898405466170333)
  c04 = convert(T2, 3.7504010171562823)
  c11 = convert(T2, 0.0)
  c12 = convert(T2, 0.26204282091330466)
  c13 = convert(T2, 0.05879875232001766)
  c14 = convert(T2, 0.758661169101175)
  beta11 = convert(T, -1.8453464565104432)
  beta12 = convert(T, 2.688764531100726)
  beta13 = convert(T, -0.2523866501071323)
  beta14 = convert(T, 0.40896857551684956)
  beta21 = convert(T, 0.4969658141589478)
  beta22 = convert(T, -0.5771202869753592)
  beta23 = convert(T, -0.12919702470322217)
  beta24 = convert(T, 0.2093514975196336)
  beta31 = convert(T, 2.8453464565104425)
  beta32 = convert(T, -2.688764531100725)
  beta33 = convert(T, 0.2523866501071322)
  beta34 = convert(T, -0.40896857551684945)
  beta41 = convert(T, 0.11522663875443433)
  beta42 = convert(T, -0.57877086147738)
  beta43 = convert(T, 0.2857851028163886)
  beta44 = convert(T, 0.17775911990655704)
  FourStageSRIConstantCache(a021,a031,a032,a041,a042,a043,a121,a131,a132,a141,
  a142,a143,b021,b031,b032,b041,b042,b043,b121,b131,b132,b141,b142,b143,α1,α2,
  α3,α4,c02,c03,c04,c11,c12,c13,c14,beta11,beta12,beta13,beta14,beta21,beta22,
  beta23,beta24,beta31,beta32,beta33,beta34,beta41,beta42,beta43,beta44)
end

function alg_cache(alg::SOSRI,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  SOSRIConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function SOSRI2ConstantCache(T::Type, T2::Type)
  a021 = convert(T, 0.13804532298278663)
  a031 = convert(T, 0.5818361298250374)
  a032 = convert(T, 0.4181638701749618)
  a041 = convert(T, 0.4670018408674211)
  a042 = convert(T, 0.8046204792187386)
  a043 = convert(T, -0.27162232008616016)
  a121 = convert(T, 0.45605532163856893)
  a131 = convert(T, 0.7555807846451692)
  a132 = convert(T, 0.24441921535482677)
  a141 = convert(T, 0.6981181143266059)
  a142 = convert(T, 0.3453277086024727)
  a143 = convert(T, -0.04344582292908241)
  b021 = convert(T, 0.08852381537667678)
  b031 = convert(T, 1.0317752458971061)
  b032 = convert(T, 0.4563552922077882)
  b041 = convert(T, 1.73078280444124)
  b042 = convert(T, -0.46089678470929774)
  b043 = convert(T, -0.9637509618944188)
  b121 = convert(T, 0.6753186815412179)
  b131 = convert(T, -0.07452812525785148)
  b132 = convert(T, -0.49783736486149366)
  b141 = convert(T, -0.5591906709928903)
  b142 = convert(T, 0.022696571806569924)
  b143 = convert(T, -0.8984927888368557)
  α1 = convert(T, -0.15036858140642623)
  α2 = convert(T, 0.7545275856696072)
  α3 = convert(T, 0.686995463807979)
  α4 = convert(T, -0.2911544680711602)
  c02 = convert(T2, 0.13804532298278663)
  c03 = convert(T2, 0.9999999999999992)
  c04 = convert(T2, 0.9999999999999994)
  c11 = convert(T2, 0.0)
  c12 = convert(T2, 0.45605532163856893)
  c13 = convert(T2, 0.999999999999996)
  c14 = convert(T2, 0.9999999999999962)
  beta11 = convert(T, -0.45315689727309133)
  beta12 = convert(T, 0.8330937231303951)
  beta13 = convert(T, 0.3792843195533544)
  beta14 = convert(T, 0.24077885458934192)
  beta21 = convert(T, -0.4994383733810986)
  beta22 = convert(T, 0.9181786186154077)
  beta23 = convert(T, -0.25613778661003145)
  beta24 = convert(T, -0.16260245862427797)
  beta31 = convert(T, 1.4531568972730915)
  beta32 = convert(T, -0.8330937231303933)
  beta33 = convert(T, -0.3792843195533583)
  beta34 = convert(T, -0.24077885458934023)
  beta41 = convert(T, -0.4976090683622265)
  beta42 = convert(T, 0.9148155835648892)
  beta43 = convert(T, -1.4102107084476505)
  beta44 = convert(T, 0.9930041932449877)
  FourStageSRIConstantCache(a021,a031,a032,a041,a042,a043,a121,a131,a132,a141,
  a142,a143,b021,b031,b032,b041,b042,b043,b121,b131,b132,b141,b142,b143,α1,α2,
  α3,α4,c02,c03,c04,c11,c12,c13,c14,beta11,beta12,beta13,beta14,beta21,beta22,
  beta23,beta24,beta31,beta32,beta33,beta34,beta41,beta42,beta43,beta44)
end

function alg_cache(alg::SOSRI2,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  SOSRI2ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

@cache struct FourStageSRICache{uType,randType,tabType,rateNoiseType,rateType,possibleRateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  chi1::randType
  chi2::randType
  chi3::randType
  tab::tabType
  g1::rateNoiseType
  g2::rateNoiseType
  g3::rateNoiseType
  g4::rateNoiseType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  E₁::rateType
  E₂::rateType
  tmp::rateType
  H02::possibleRateType
  H03::possibleRateType
end

function alg_cache(alg::SRIW2,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    chi1 = copy(ΔW)
    chi2 = copy(ΔW)
    chi3 = copy(ΔW)
  else
    chi1 = zero(ΔW)
    chi2 = zero(ΔW)
    chi3 = zero(ΔW)
  end
  tab = SRIW2ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype); g2 = zero(noise_rate_prototype)
  g3 = zero(noise_rate_prototype); g4 = zero(noise_rate_prototype)
  k1 = zero(rate_prototype); k2 = zero(rate_prototype)
  k3 = zero(rate_prototype); k4 = zero(rate_prototype)
  E₁ = zero(rate_prototype); E₂ = zero(rate_prototype)
  tmp = zero(rate_prototype)
  FourStageSRICache(u,uprev,chi1,chi2,chi3,tab,g1,g2,g3,g4,k1,k2,k3,k4,E₁,E₂,tmp,tmp,tmp)
end

function alg_cache(alg::SOSRI,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    chi1 = copy(ΔW)
    chi2 = copy(ΔW)
    chi3 = copy(ΔW)
  else
    chi1 = zero(ΔW)
    chi2 = zero(ΔW)
    chi3 = zero(ΔW)
  end
  tab = SOSRIConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype); g2 = zero(noise_rate_prototype)
  g3 = zero(noise_rate_prototype); g4 = zero(noise_rate_prototype)
  k1 = zero(rate_prototype); k2 = zero(rate_prototype)
  k3 = zero(rate_prototype); k4 = zero(rate_prototype)
  E₁ = zero(rate_prototype); E₂ = zero(rate_prototype)
  tmp = zero(rate_prototype)
  FourStageSRICache(u,uprev,chi1,chi2,chi3,tab,g1,g2,g3,g4,k1,k2,k3,k4,E₁,E₂,tmp,tmp,tmp)
end

function alg_cache(alg::SOSRI2,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    chi1 = copy(ΔW)
    chi2 = copy(ΔW)
    chi3 = copy(ΔW)
  else
    chi1 = zero(ΔW)
    chi2 = zero(ΔW)
    chi3 = zero(ΔW)
  end
  tab = SOSRI2ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype); g2 = zero(noise_rate_prototype)
  g3 = zero(noise_rate_prototype); g4 = zero(noise_rate_prototype)
  k1 = zero(rate_prototype); k2 = zero(rate_prototype)
  k3 = zero(rate_prototype); k4 = zero(rate_prototype)
  E₁ = zero(rate_prototype); E₂ = zero(rate_prototype)
  tmp = zero(rate_prototype); H02 = zero(rate_prototype)
  H03 = zero(rate_prototype)
  FourStageSRICache(u,uprev,chi1,chi2,chi3,tab,g1,g2,g3,g4,k1,k2,k3,k4,E₁,E₂,tmp,H02,H03)
end


text:
struct Point
	x::Float32
	y::Float32
	z::Float32
end

@enum Plane begin
	undef_plane
    x_plane
    y_plane
	z_plane
end

struct LOR
	P1::Point
	P2::Point
	timeDiff::Float32
end

struct Event
	lor::LOR
	mainPlane::Plane
end

max_stars_repo_name: michielvangendt/PET-Julia

max_issues_repo_name: michielvangendt/PET-Julia

max_forks_repo_name: michielvangendt/PET-Julia


text:
using Documenter, CCBlade

makedocs(
    format = Documenter.HTML(),
    sitename = "CCBlade Documentation",
    pages = [
        "index.md",
        "tutorial.md"
    ],
    modules = [CCBlade]
)

# deploydocs(
#     repo = "github.com/byuflowlab/CCBlade.jl.git"
# )


text:
"""
    NLevelBasis(N)

Basis for a system consisting of N states.
"""
struct NLevelBasis{T} <: Basis
    shape::Vector{T}
    N::T
    function NLevelBasis(N::T) where T<:Integer
        if N < 1
            throw(DimensionMismatch())
        end
        new{T}([N], N)
    end
end

==(b1::NLevelBasis, b2::NLevelBasis) = b1.N == b2.N

"""
    transition(b::NLevelBasis, to::Integer, from::Integer)

Transition operator ``|\\mathrm{to}⟩⟨\\mathrm{from}|``.
"""
function transition(b::NLevelBasis, to::Integer, from::Integer)
    if to < 1 || b.N < to
        throw(BoundsError("'to' index has to be between 1 and b.N"))
    end
    if from < 1 || b.N < from
        throw(BoundsError("'from' index has to be between 1 and b.N"))
    end
    op = SparseOperator(b)
    op.data[to, from] = 1.
    op
end

"""
    nlevelstate(b::NLevelBasis, n::Integer)

State where the system is completely in the n-th level.
"""
function nlevelstate(b::NLevelBasis, n::Integer)
    if n < 1 || b.N < n
        throw(BoundsError("n has to be between 1 and b.N"))
    end
    basisstate(b, n)
end


text:
##### temperature function
@inline function tempFunc(temp, p)
    k = exp(-p.Ea/(8.3145*(temp+273.15)))*(1.0-exp(temp+273.15 - p.T⁺))
    k = max(0.0, k)
    OGT_rate = exp(-p.Ea/(8.3145*(p.T⁺-2)))
    return k/OGT_rate
end

##### calculate photosynthesis rate (mmolC/individual/second)
@inline function calc_PS(par, temp, Chl, Bm, p)
    αI  = par * p.α * p.Φ
    PCm = p.PCmax * tempFunc(temp, p)
    PS  = PCm * (1.0 - exp(-αI / max(1.0e-30, PCm) * Chl / max(1.0e-30, Bm))) * Bm
    return PS
end

@kernel function calc_inorganic_uptake_kernel!(plank, nuts, p)
    i = @index(Global)
    @inbounds plank.PS[i] = calc_PS(nuts.par[i], nuts.T[i], plank.Chl[i], plank.Bm[i], p) * plank.ac[i]

end

function calc_inorganic_uptake!(plank, nuts, p, arch::Architecture)
    kernel! = calc_inorganic_uptake_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, nuts, p)
    wait(device(arch), event)
    return nothing
end

##### calculate respiration (mmolC/individual/second)
@kernel function calc_respir_kernel!(plank, T, p)
    i = @index(Global)
    @inbounds plank.resp[i] = p.respir_a * plank.Bm[i] * tempFunc(T[i], p) * plank.ac[i]
end
function calc_respir!(plank, T, p, arch)
    kernel! = calc_respir_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, T, p)
    wait(device(arch), event)
    return nothing
end

##### update C quotas
@kernel function update_quotas_kernel!(plank, ΔT)
    i = @index(Global)
    @inbounds plank.Bm[i]  += (plank.PS[i] - plank.resp[i]) * ΔT
    @inbounds plank.age[i] += ΔT / 3600.0 * plank.ac[i]
end
function update_quotas!(plank, ΔT, arch)
    kernel! = update_quotas_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, ΔT)
    wait(device(arch), event)
    return nothing
end

##### update cell size
@kernel function update_cellsize_kernel!(plank, p)
    i = @index(Global)
    @inbounds plank.Sz[i]  = plank.Bm[i] / (p.Cquota * p.Nsuper)
    @inbounds plank.Chl[i] = plank.Bm[i] * p.Chl2C
end
function update_cellsize!(plank, p, arch)
    kernel! = update_cellsize_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, p)
    wait(device(arch), event)
    return nothing
end

##### track temperature history
@kernel function calc_thermal_history_kernel!(plank, nuts, p, ΔT)
    i = @index(Global)
    # @inbounds plank.Th[i] += max(0, nuts.T[i] + 273.15 - p.T⁺) * ΔT / 3600
    @inbounds plank.Th[i] += ΔT/3600 * exp(0.22 * max(0, nuts.T[i] + 273.15 - p.T⁺)) * isless(p.T⁺, nuts.T[i]+273.15)
end
function calc_thermal_history!(plank, nuts, p, ΔT, arch)
    kernel! = calc_thermal_history_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, nuts, p, ΔT)
    wait(device(arch), event)
    return nothing
end


text:
#  Copyright 2016, Liye Zhang and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuSO
# Surrogate Optimization package for Julia
# See http://github.com/zhangliye/JuSO.jl
#############################################################################
module JuSOT

export
# Classes
  SoModel,
#function
  test_function1, test_ut1

include("utils.jl")

type SoModel
  name::ASCIIString
end

function test_function1()
  println("i am test_function1()")
  return 1
end

end # module


text:
"""
`ind = nearestbelow(array, nearto)` return the indices for the elements in `array` that are nearest below `nearto`.
The nearest elements are `array[ind]`.
# Example
julia> nearestbelow([1,2,3,4,5,6],[1.9, 3.3])
2-element Vector{Float64}:
 1
 3
"""
function nearestbelow(array::Vector{<:Real}, nearto::Vector{<:Real})
	ind = Vector{Int64}(undef,length(nearto));
	for (i,target) in enumerate(nearto)
		ind[i] = findlast(array .<= target);
	end
	return ind
end


text:
const t0 = TreatmentTerm(:g, TestTreatment(:t, 0), TestParallel(0))
const t1 = treat(:g, TestTreatment(:t, 0), TestParallel(0))
const t2 = treat(:g, TestTreatment(:t, 1), TestParallel(0))
const t3 = treat(:g, TestTreatment(:t, 0), TestParallel(1))

@testset "TermSet" begin
    @test TermSet() == termset()
    @test TermSet([term(:a)]) == termset([term(:a)]) == TermSet(Set{AbstractTerm}([term(:a)]))
    @test TermSet(1, :a, term(:a)) == termset(1, :a, term(:a)) ==
        TermSet(Set{AbstractTerm}([term(1), term(:a)]))

    ts0 = TermSet()
    @test isempty(ts0)
    @test length(ts0) == 0
    ts1 = TermSet(term(:a))
    @test isempty(ts1) == false
    @test length(ts1) == 1
    @test term(:a) in ts1
    @test push!(ts0, term(:a)) == ts1
    @test pop!(ts0, term(:a)) == term(:a)
    @test isempty(ts0)
    @test pop!(ts0, term(:a), 0) === 0
    delete!(ts1, term(:a))
    @test ts1 == ts0
    ts1 = TermSet(term(:a))
    @test empty!(ts1) == ts0

    @test eltype(TermSet) == AbstractTerm
    ts1 = TermSet(term(:a))
    @test iterate(ts1)[1] == term(:a)
    @test Base.emptymutable(ts1) == ts0
end

@testset "==" begin
    @test term(:x) + term(:y) == term(:y) + term(:x)
    @test term(:x) + term(:y) + term(:y) == term(:y) + term(:x) + term(:x)

    @test term(:x) & term(:y) == term(:y) & term(:x)
    @test term(:x) & term(:y) + term(:y) & term(:x) == term(:y) & term(:x)
    @test term(:x) & term(:y) + term(:z) == term(:z) + term(:y) & term(:x)
    @test lag(term(:x),1) & term(:z) == term(:z) & lag(term(:x),1)
    # StatsModels does not enforce uniqueness of terms in interactions
    @test term(:x) & term(:y) & term(:x) != term(:x) & term(:y)

    @test @formula(y ~ x + z) == @formula(y ~ z + x)
    @test @formula(y ~ t & z + x) == @formula(y ~ x + z & t)
    @test @formula(y ~ lag(x, 1)&z + t) == @formula(y ~ t + z&lag(x, 1))
end

@testset "treat" begin
    @test t1 == t0
    @test t2 != t0
    @test t3 != t0
end

@testset "hastreat" begin
    @test hastreat(t0)

    f0 = @formula(y ~ x)
    f1 = @formula(y ~ treat(g, ttreat(t, 0), tpara(0)))
    f2 = @formula(y ~ treat(g, ttreat(t, 0), tpara(0)) & x)

    @test hastreat(f0.rhs) == false
    @test hastreat(f1.rhs)
    @test hastreat(f2.rhs)

    @test hastreat(f0) == false
    @test hastreat(f1)
    @test hastreat(f2)
end

@testset "parse_treat" begin
    @testset "with @formula" begin
        f = @formula(y ~ x)
        @test_throws ArgumentError parse_treat(f)

        f = @formula(y ~ treat(g, ttreat(t, 0), tpara(0)) + x & z)
        t = parse_treat(f)
        @test t == (t1, TermSet(), termset(term(:x)&term(:z)))

        f = @formula(y ~ treat(g, ttreat(t, 0), tpara(0)) & x & z)
        t = parse_treat(f)
        @test t == (t1, TermSet(term(:x), term(:z)), TermSet())

        f = @formula(y ~ treat(g, ttreat(t, 0), tpara(0)) & treat(g, ttreat(t, 0), tpara(0)))
        @test_throws ArgumentError parse_treat(f)

        f = @formula(y ~ treat(g, ttreat(t, 0), tpara(0)) & treat(g, ttreat(t, 1), tpara(0)))
        @test_throws ArgumentError parse_treat(f)

        # + checks uniqueness of terms
        f = @formula(y ~ treat(g, ttreat(t, 0), tpara(0)) + treat(g, ttreat(t, 0), tpara(0)))
        @test f == @formula(y ~ treat(g, ttreat(t, 0), tpara(0)))

        f = @formula(y ~ treat(g, ttreat(t, 0), tpara(0)) + treat(g, ttreat(t, 1), tpara(0)))
        @test_throws ArgumentError parse_treat(f)

        f = @formula(y ~ treat(g, ttreat(t, 0), tpara(0)) & x + treat(g, ttreat(t, 0), tpara(0)) & z)
        @test_throws ArgumentError parse_treat(f)
    end

    @testset "without @formula" begin
        f = term(:y) ~ term(:x)
        @test_throws ArgumentError parse_treat(f)

        f = term(:y) ~ treat(:g, TR, PR) + term(:x) & term(:z) + lag(term(:x),1)
        t = parse_treat(f)
        @test t == (t1, TermSet(), TermSet((term(:x) & term(:z) + lag(term(:x),1))))

        f = term(:y) ~ treat(:g, TR, PR) & term(:x) & term(:z)
        t = parse_treat(f)
        @test t == (t1, TermSet((term(:x), term(:z))), TermSet())

        f = term(:y) ~ treat(:g, TR, PR) & treat(:g, TR, PR)
        @test_throws ArgumentError parse_treat(f)

        f = term(:y) ~ treat(:g, TR, PR) & treat(:g, TestTreatment(:t, 1), PR)
        @test_throws ArgumentError parse_treat(f)

        # + checks uniqueness of terms
        f = term(:y) ~ treat(:g, TR, PR) + treat(:g, TR, PR)
        @test f == (term(:y) ~ treat(:g, TR, PR))

        f = term(:y) ~ treat(:g, TR, PR) + treat(:g, TestTreatment(:t, 1), PR)
        @test_throws ArgumentError parse_treat(f)

        f = term(:y) ~ treat(:g, TR, PR) & term(:x) +
            treat(:g, TR, PR) & term(:z)
        @test_throws ArgumentError parse_treat(f)
    end
end

@testset "parse_intercept!" begin
    @test isintercept(term(1))
    @test isomitsintercept(term(0))
    @test isomitsintercept(term(-1))

    @test parse_intercept!(TermSet()) == (false, false)
    @test parse_intercept!(TermSet(term(:x))) == (false, false)
    ts = TermSet(term(1))
    @test parse_intercept!(ts) == (true, false)
    @test isempty(ts)
    ts = TermSet((term(0), term(-1)))
    @test parse_intercept!(ts) == (false, true)
    @test isempty(ts)
    ts = TermSet((term(1), term(0), term(:x)))
    @test parse_intercept!(ts) == (true, true)
    @test collect(ts) == [term(:x)]
end

@testset "schema" begin
    hrs = exampledata("hrs")
    sc = schema(termset(:wave, :oop_spend), hrs)
    @test sc[term(:wave)] isa ContinuousTerm
    @test sc[term(:oop_spend)] isa ContinuousTerm
    sc = schema(termset(:wave, :oop_spend), hrs, Dict(:wave=>CategoricalTerm))
    @test sc[term(:wave)] isa CategoricalTerm
    @test sc[term(:oop_spend)] isa ContinuousTerm

    cols = VecColumnTable(hrs)
    @test concrete_term(term(:wave), cols, nothing) isa ContinuousTerm
    @test concrete_term(term(:wave), cols, Dict(:wave=>CategoricalTerm)) isa CategoricalTerm
    @test concrete_term(term(:wave), cols, term(:wave)) isa Term

    @test termvars(termset()) == Symbol[]
    @test termvars(termset(:wave, :oop_spend)) in ([:wave, :oop_spend], [:oop_spend, :wave])
end


text:
## This file defines basic mathematical functions and operators for
## Timenum and Timemtr types. These operators and functions make sense
## only for objects that are restricted to numeric values. For type
## Timedata, they are not defined, similarly to DataFrames.
##
## Operators usually try to build on the same operator for DataArrays.
## Only where it makes sense, Timematr operators will be build on
## operators that are defined on numeric Arrays.

#######################################################
## type preserving symmetric operators and functions ##
#######################################################

## mathematical operators are defined elementwise
##
## if two Timematr / Timenum objects are involved, indices and column
## names must coincide.
##
## elementwise operators involving two different TimeData types are
## not allowed.
##
## function headers
## f(inst::NewType, inst2::NewType)
## f(inst::NewType, b::Number)
## f(b::Number, inst::NewType)

## type preserving: applied to Timematr, it will again return Timematr
## object

## macro is parametrized with respect to type, because this way
## Timematr implementations involving core() can easily be re-defined
## and related to DataArray operators as well
##
## if time indices do not coincide, time indices will be inferred from
## the first TimeData instance involved in the operation!
macro pres_symmetric_timenum(f, myType)
    esc(quote
        function $(f)(inst::$(myType), inst2::$(myType))
            if !similarMeta(inst, inst2)
                error("Timematr variables not similar")
            end

            dfResult = DataFrame()
            for col in names(inst)
                dfResult[col] = $(f)(inst.vals[col], inst2.vals[col])
            end
            return $(myType)(dfResult, idx(inst))
        end

        function $(f)(inst::$(myType), val::Number)

            dfResult = DataFrame()
            for col in names(inst)
                dfResult[col] = $(f)(inst.vals[col], val)
            end
            return $(myType)(dfResult, idx(inst))
        end

        function $(f)(val::Number, inst::$(myType))

            dfResult = DataFrame()
            for col in names(inst)
                dfResult[col] = $(f)(val, inst.vals[col])
            end
            return $(myType)(dfResult, idx(inst))
        end

    end)
end

macro pres_symmetric_timematr(f)
    ## can be related to Array operators through core function
    esc(quote
        ## tm1 .+ tm2
        function $(f)(inst::Timematr, inst2::Timematr)
            if !similarMeta(inst, inst2)
                error("Timematr variables not similar")
            end
            return Timematr($(f)(core(inst), core(inst2)),
                            names(inst), idx(inst))
        end

        ## tm .+ 3
        function $(f)(inst::Timematr,b::Number)
            return Timematr($(f)(core(inst), b), names(inst),
                            idx(inst))
        end

        ## 3 .+ tm
        function $(f)(b::Number,inst::Timematr)
            return Timematr($(f)(b,core(inst)), names(inst),
                            idx(inst))
        end
    end)
end

importall Base
t = :Timenum
for f in [:(.+), :(.-), :(.*), :(./), :(.^)]
    eval(macroexpand(:(@pres_symmetric_timematr($f))))
    eval(macroexpand(:(@pres_symmetric_timenum($f, $t))))
end

## other function with two numbers as input
for f in [:(div), :(mod), :(fld), :(rem)]
    eval(macroexpand(:(@pres_symmetric_timematr($f))))
    eval(macroexpand(:(@pres_symmetric_timenum($f, $t))))
end

## only implemented for Timematr
const matrix_mult = [:(*), :(/)]
for f in matrix_mult
    eval(macroexpand(:(@pres_symmetric_timematr($f))))
end

#################################################
## type preserving functions without arguments ##
#################################################

## f(inst::NewType)
const unary_operators = [:(+), :(-)]    # [:(*), :(/)]
const mathematical_functions = [:abs, :sign,
                                ## :acos, :acosh, :asin, :asinh,
                                ## :atan, :atanh, :sin, :sinh, :cos,
                                ## :cosh, :tan, :tanh,
                                ## :exp2, :expm1, :log10, :log1p,
                                ## :log2, :exponent,
                                ## :lgamma, :digamma, :erf, :erfc,
                                :exp, :log, :sqrt, :gamma
                                ]

pres_msUnitary_functions = [unary_operators, mathematical_functions]

macro pres_msUnitary_timenum(f, myType)
    esc(quote
        function $(f)(inst::$(myType))

            dfResult = DataFrame()
            for col in names(inst)
                dfResult[col] = $(f)(inst.vals[col])
            end

            $(myType)(dfResult, idx(inst))
        end
    end)
end

macro pres_msUnitary_timematr(f)
    esc(quote
        $(f)(inst::Timematr) =
            Timematr($(f)(core(inst)), names(inst), idx(inst))
    end)
end

t = :Timenum
for f in pres_msUnitary_functions
    eval(macroexpand(:(@pres_msUnitary_timematr($f))))
    eval(macroexpand(:(@pres_msUnitary_timenum($f, $t))))
end

function !(td::Timedata)
    ## negate boolean values
    df = DataFrame()
    for (nam, col) in eachcol(td.vals)
        df[nam] = !col
    end
    return Timedata(df, td.idx)
end

######################################################
## non-preserving symmetric functions and operators ##
######################################################

## f(inst::NewType, inst2::NewType)
## f(inst::NewType, b::Number)
## f(b::Number, inst::NewType)

## for timedata and timenum objects apply function to DataArrays
macro nonpres_msSymmetric_timedata(f)
    esc(quote
        function $(f)(inst::AbstractTimedata, inst2::AbstractTimedata)
            if !equMeta(inst, inst2)
                error("TimeData objects not similar")
            end

            df2 = DataFrame()
            for nam in names(inst)
                df2[nam] = $(f)(inst.vals[nam], inst2.vals[nam])
            end
            return Timedata(df2, idx(inst))
        end

        function $(f)(inst::AbstractTimedata,b::Union(String,Number))
            df2 = DataFrame()
            for col in eachcol(inst.vals)
                df2[col[1]] = $(f)(col[2], b)
            end
            return Timedata(df2, idx(inst))
        end

        function $(f)(b::Union(String,Number),inst::AbstractTimedata)
            df2 = DataFrame()
            for col in eachcol(inst.vals)
                df2[col[1]] = $(f)(col[2], b)
            end
            return Timedata(df2, idx(inst))
        end
    end)
end

## for timematr objects apply function to Array
## for comparison, only Numbers need to be considered
macro nonpres_msSymmetric_timematr(f)
    esc(quote
        function $(f)(inst::Timematr, inst2::Timematr)
            if !equMeta(inst, inst2)
                error("Timematr objects not similar")
            end

            return Timedata(convert(Array{Bool, 2}, $(f)(core(inst), core(inst2))),
                            names(inst), idx(inst))
        end

        function $(f)(inst::Timematr, b::Number)
            return Timedata(convert(Array{Bool, 2}, $(f)(core(inst),b)),
                            names(inst), idx(inst))
        end

        function $(f)(b::Union(String,Number),inst::Timematr)
            return Timedata(convert(Array{Bool, 2}, $(f)(b,core(inst))),
                            names(inst), idx(inst))
        end
    end)
end

for f in [:(.==), :(.!=), :(.>), :(.>=), :(.<), :(.<=)]
    eval(macroexpand(:(@nonpres_msSymmetric_timematr($f))))
    eval(macroexpand(:(@nonpres_msSymmetric_timedata($f))))
end

#########################################################
## type preserving functions, single or extra argument ##
#########################################################

## f(td::NewType)
## f(td::NewType, i::Integer)

macro pres_msSingle_or_extra_timematr(f)
    esc(quote
        $(f)(inst::Timematr) =
            Timematr($(f)(core(inst)), names(inst), idx(inst))
        $(f)(inst::Timematr, i::Integer) =
            Timematr($(f)(core(inst), i), names(inst), idx(inst))
    end)
end

macro pres_msSingle_or_extra_timenum(f, myType)
    esc(quote
        function $(f)(inst::$(myType))
            dfResult = DataFrame()
            for col in names(inst)
                dfResult[col] = $(f)(inst.vals[col])
            end
            return $(myType)(dfResult, idx(inst))
        end

        function $(f)(inst::Timenum, i::Integer)
            dfResult = DataFrame()
            for col in names(inst)
                dfResult[col] = $(f)(inst.vals[col], i)
            end
            return $(myType)(dfResult, idx(inst))
        end
    end)
end

t = :Timenum
for f in [:round, :ceil, :floor, :trunc]
    eval(macroexpand(:(@pres_msSingle_or_extra_timematr($f))))
    eval(macroexpand(:(@pres_msSingle_or_extra_timenum($f, $t))))
end

######################################################
## non-preserving comparison operators for timedata ##
######################################################

macro time_data_comparison(f)
    esc(quote
        function $(f)(td1::Timedata, td2::Timedata)
            if !equMeta(td1, td2)
                error("Timedata objects not similar")
            end

            dfResult = DataFrame()
            for col in names(td1)
                dfResult[col] = $(f)(td1.vals[col], td2.vals[col])
            end
            return Timedata(dfResult, idx(td1))
        end
    end)
end

for f in [:(&), :(|), :($)]
    eval(macroexpand(:(@time_data_comparison($f))))
end


text:
# Julia wrapper for header: /usr/include/libavutil/tea.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_tea_alloc,
    av_tea_init,
    av_tea_crypt

function av_tea_alloc()
    ccall((:av_tea_alloc, libavutil), Ptr{AVTEA}, ())
end

function av_tea_init(ctx, key::NTuple{16, UInt8}, rounds::Integer)
    ccall((:av_tea_init, libavutil), Cvoid, (Ptr{AVTEA}, NTuple{16, UInt8}, Cint), ctx, key, rounds)
end

function av_tea_crypt(ctx, dst, src, count::Integer, iv, decrypt::Integer)
    ccall((:av_tea_crypt, libavutil), Cvoid, (Ptr{AVTEA}, Ptr{UInt8}, Ptr{UInt8}, Cint, Ptr{UInt8}, Cint), ctx, dst, src, count, iv, decrypt)
end


text:
# csci_textload.jl
#
#	Load Campbell Scientific TOA5 files
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 1.6.0
# Created: 04.11.2013
# Last Edit: 28.07.2021

"""# csci_textload

Load TOA5 data files generated by Campbell Scientific dataloggers

---

### Examples

`data = csci_textload(\"K:\\\\Data\\\\AirTemp.dat\")` Load the data from the file into a dataframe, data, converting time columns to DateTimes

`data = csci_textload([\"K:\\\\Data\\\\AirTemp_20141219-000000.dat\",\"K:\\\\Data\\\\AirTemp_20141220-000000.dat\"])` Load a list of similarly formatted files

`data = csci_textload(\"K:\\\\Data\",\"MetData_AirTemp\")` Load a directory of data which start with \"MetData_AirTemp\"

`data = csci_textload(\"K:\\\\Data\",\"MetData_AirTemp\",DateTime(2014,12,18))` Load a directory of data starting at 18.12.2014

`data = csci_textload(\"K:\\\\Data\",\"MetData_AirTemp\",DateTime(2014,12,18),DateTime(2014,12,19,12))` Load a directory of data from 18.12.2014 to 12:00 on 19.12.2014

`(data,logger_Str,columns_Str,units_Str,type_Str) = csci_textload(\"K:\\\\Data\\\\AirTemp.dat\",headerout = true)`

---

`data = csci_textload(F;headerlines=4,headeroutput=false,verbose=true,strlist=[],intlist=[],timelist=[])`\n
* **data**::DataFrame = Loaded data, including a fully parsed time column
* **F**::String or Vector{String} = File or an array of files to load
* **headerlines**::Int (optional) = Number of header lines, 4 is default
* **headeroutput**::Bool (optional) = Add four header lines to output, FALSE is default
* **verbose**::Bool (optional) = Display information about what is happening when set to TRUE, TRUE is default
* **strlist**::Vector{String} (optional) = List of column names that should be strings, [] is default
* **intlist**::Vector{String} (optional) = List of column names that should be integers, [] is default
* **timelist**::Vector{String} (optional) = List of column names that should be DateTime values, [] is default
* **select**::Vector{Symbol} (optional) = List of column names as Symbols to be loaded, names not included will not be loaded. [] is default
\n\n
`data = csci_textload(F,`rootname`;headerlines=4,headeroutput=false,verbose=true,recur=2,strlist=[],intlist=[],timelist=[])`\n
`data = csci_textload(F,rootname,`mindate`;headerlines=4,headeroutput=false,verbose=true,recur=2,strlist=[],intlist=[],timelist=[])`\n
`data = csci_textload(F,rootname,mindate,`maxdate`;headerlines=4,headeroutput=false,verbose=true,recur=2,strlist=[],intlist=[],timelist=[])`\n
* **F**::String = Directory of files to load
* **rootname**::String = Root of the file names to process, \"Dedelow_CR3000\" is the root of \"Dedelow_CR3000_Soil.dat\", \"Dedelow_CR3000_Temperature\", etc.
* **mindate**::DateTime = Load files including and after this date and time
* **maxdate**::DateTime = Load files up to this date and time"""
function csci_textload(F::String,rootname::String,mindate::DateTime=DateTime(0),maxdate::DateTime=DateTime(9999);headerlines::Int=4,headeroutput::Bool=false,verbose::Bool=true,recur::Int=2,strlist::Vector{String}=String[],intlist::Vector{String}=String[],timelist::Vector{String}=String[],select::Vector{Symbol}=Symbol[])

	# Check
	if !isfile(F) & !isdir(F)
		error("First input must be a valid file or directory. For a single file use csci_textread.")
	end

	# Load a directory of data
	(files,folder) = dirlist(F,regex=r"\.dat$",recur=recur)

	# Find files only with the root name
	filebase = ""
	keepers = Array{Bool}(undef,length(files))
	for i=1:1:length(files)
		keepers[i] = occursin(Regex(rootname),basename(files[i]))
	end
	files = files[keepers]

	# Files Timestamp Conversion
	(mintimes,maxtimes) = csci_times(files)
	f = sortperm(mintimes)
	files = files[f]
	mintimes = mintimes[f]
	maxtimes = maxtimes[f]

	# Remove Files Out of Range
	f1 = Array{Bool}(undef,length(files))
	for i=1:1:length(files)
		f1[i] = (mintimes[i] <= mindate < maxtimes[i]) | (mintimes[i] <= maxdate < maxtimes[i])
	end
	f = (mindate .<= mintimes .< maxdate) .| (mindate .<= maxtimes .< maxdate) .| f1
	files = files[f]
	mintimes = mintimes[f]
	maxtimes = maxtimes[f]

	# Convert to Vector{String,1}
	temp = Vector{String}(undef,length(files))
	for i=1:1:length(files)
		temp[i] = files[i]
	end
	files = temp

	if headeroutput
		(D,loggerStr,colsStr,unitsStr,processingStr) = csci_textload(files,headerlines=headerlines,headeroutput=headeroutput,verbose=verbose,strlist=strlist,intlist=intlist,timelist=timelist,select=select)

		# Remove Values Still Out of Range
		f = find(mindate .<= D[1] .< maxdate)
		D = D[f,:]

		return D, loggerStr, colsStr, unitsStr, processingStr
	else
		D = csci_textload(files,headerlines=headerlines,headeroutput=headeroutput,verbose=verbose,strlist=strlist,intlist=intlist,timelist=timelist,select=select)

		# Remove Values Still Out of Range
		f = findall(mindate .<= D[!,1] .< maxdate)
		D = D[f,:]

		return D
	end
end

function csci_textload(F::Vector{T};headerlines::Int=4,headeroutput::Bool=false,verbose::Bool=true,strlist::Vector{String}=String[],intlist::Vector{String}=String[],timelist::Vector{String}=String[],select::Vector{Symbol}=Symbol[]) where T <: String
	D = DataFrame[]
	loggerStr = ""
	colsStr = ""
	unitsStr = ""
	processingStr = ""

	# Load the list of files
	tempD = DataFrame[]
	for i=1:1:length(F)
		if headeroutput
			(tempD,loggerStr,colsStr,unitsStr,processingStr) = csci_textread(F[i],headerlines=headerlines,headeroutput=headeroutput,verbose=verbose,strlist=strlist,intlist=intlist,timelist=timelist,select=select)
		else
			(tempD) = csci_textread(F[i],headerlines=headerlines,headeroutput=headeroutput,verbose=verbose,strlist=strlist,intlist=intlist,timelist=timelist,select=select)
		end

		if isempty(D)
			D = tempD
		else
			D = [D;tempD]
		end
	end

	# Output Results
	if headeroutput
		return D, loggerStr, colsStr, unitsStr, processingStr
	else
		return D
	end
end

max_stars_repo_name: JeremyRueffer/ClimateDataIO

max_issues_repo_name: JeremyRueffer/ClimateDataIO


text:
# Workaround for JuliaLang/julia/pull/28625
if Base.HOME_PROJECT[] !== nothing
    Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

using Documenter, Parametron

makedocs(
    modules = [Parametron, Parametron.Functions],
    checkdocs = :exports,
    format = :html,
    root = @__DIR__,
    sitename ="Parametron.jl",
    authors = "Twan Koolen and contributors.",
    pages = [
        "Home" => "index.md",
        "API" => [
            "api/functions.md",
            "api/parameter.md",
            "api/lazyexpression.md",
            "api/model.md",
            "api/debug.md",
        ]
    ],
    html_prettyurls = parse(Bool, get(ENV, "CI", "false")),
)

deploydocs(
    repo = "github.com/tkoolen/Parametron.jl.git"
)


text:
using .Turing
export pointwise_log_likelihoods, psis_loo, psis, loo_from_psis

const TURING_MODEL_ARG = """
`model`: A Turing model with data in the form of `model(data)`.
"""

"""
    pointwise_log_likelihoods(model::DynamicPPL.Model, chains::Chains) -> Array

Compute pointwise log-likelihoods from a Turing model.

# Arguments
  - $TURING_MODEL_ARG
  - $CHAINS_ARG

# Returns
  - `Array`: A three dimensional array of pointwise log likelihoods. This array should be
    indexed using `array[data, sample, chains]`.
"""
function pointwise_log_likelihoods(model::DynamicPPL.Model, chains::Chains)
    # subset of chains for mcmc samples
    chain_params = MCMCChains.get_sections(chains, :parameters)
    # compute the pointwise log likelihoods
    log_like_dict = DynamicPPL.pointwise_loglikelihoods(model, chain_params)
    # Size of array (n_steps, n_chains) using first parameter
    dims = size(last(first(log_like_dict)))
    # parse "var[i]" -> i
    ind_from_string(x) = parse(Int, split(split(x, "[")[2], "]")[1])
    # collect variable names
    sorted_keys = sort(collect(keys(log_like_dict)); by=ind_from_string)
    # Convert from dictionary to 3d array
    array = [reshape(log_like_dict[i], 1, dims...) for i in sorted_keys]
    return reduce(vcat, array)
end

"""
    psis_loo(model::DynamicPPL.Model, chains::Chains, args...; kwargs...) -> PsisLoo

Use Pareto-Smoothed Importance Sampling to calculate the leave-one-out cross validation
score from a `chains` object and a Turing model.

# Arguments

  - $CHAINS_ARG
  - $TURING_MODEL_ARG
  - $ARGS [`psis`](@ref).
  - $KWARGS [`psis`](@ref).

See also: [`psis`](@ref), [`loo`](@ref), [`PsisLoo`](@ref).
"""
function psis_loo(model::DynamicPPL.Model, chains::Chains, args...; kwargs...)
    pointwise_log_likes = pointwise_log_likelihoods(model, chains)
    return psis_loo(pointwise_log_likes, args...; kwargs...)
end

"""
    loo_from_psis(model::DynamicPPL.Model, chains::Chains, args...; kwargs...) -> PsisLoo

Use Pareto-Smoothed Importance Sampling to calculate the leave-one-out cross validation
score from a `Chains` object, a Turing model, and a precalculated `Psis` object.

# Arguments

  - $CHAINS_ARG
  - $TURING_MODEL_ARG
  - `psis`: A `Psis` object containing the results of Pareto smoothed importance sampling.

See also: [`psis`](@ref), [`psis_loo`](@ref), [`PsisLoo`](@ref).
"""
function loo_from_psis(model::DynamicPPL.Model, chains::Chains, psis::Psis)
    pointwise_log_likes = pointwise_log_likelihoods(model, chains)
    return loo_from_psis(pointwise_log_likes, psis)
end

"""
    psis(model::DynamicPPL.Model, chains::Chains, args...; kwargs...) -> Psis

Generate samples using Pareto smoothed importance sampling (PSIS).

# Arguments
  - $TURING_MODEL_ARG
  - $CHAINS_ARG

See also: [`psis`](@ref), [`loo`](@ref), [`PsisLoo`](@ref).
"""
function psis(model::DynamicPPL.Model, chains::Chains, args...; kwargs...)
    log_ratios = pointwise_log_likelihoods(model, chains)
    return psis(-log_ratios, args...; kwargs...)
end

function naive_lpd(model::DynamicPPL.Model, chains::Chains, args...; kwargs...)
    log_ratios = pointwise_log_likelihoods(model, chains)
    return naive_lpd(log_ratios, args...; kwargs...)
end

max_stars_repo_name: ParadaCarleton/JuLOOa

max_issues_repo_name: ParadaCarleton/JuLOOa

max_forks_repo_name: ParadaCarleton/JuLOOa


text:
using CUDA
using AMDGPU
using Oceananigans.Architectures: CPU, CUDAGPU, AMDGPU, AbstractMultiArchitecture

import Base: zeros

zeros(FT, ::CPU, N...) = zeros(FT, N...)
zeros(FT, ::CUDAGPU, N...) = CUDA.zeros(FT, N...)
zeros(FT, ::ROCMGPU, N...) = AMDGPU.ROCArray(zeros(FT, N...))

zeros(arch::AbstractArchitecture, grid, N...) = zeros(eltype(grid), arch, N...)
zeros(grid::AbstractGrid, N...) = zeros(eltype(grid), architecture(grid), N...)

@inline Base.zero(grid::AbstractGrid) = zero(eltype(grid))


text:
using LinearAlgebra, Test, Random
using Knet.KnetArrays: A_mul_Bt, At_mul_B, At_mul_Bt, KnetArray
using Knet.Ops20: mat
using CUDA: CUDA, functional
using AutoGrad: gradcheck
include("combinatorics.jl")

#Random.seed!(42)
nsample(a,n)=collect(a)[randperm(length(a))[1:n]]

@testset "linalg" begin
    for t in (Float32,Float64)
        a = rand(t,3,5)
        b = rand(t,5,2)
        mmul(w)=w[1]*w[2]
        @test gradcheck(mmul, (a,b))
        if CUDA.functional()
            c = a * b
            ka = KnetArray(a)
            kb = KnetArray(b)
            kc = ka * kb
            @test isapprox(c, Array(kc))
            @test gradcheck(mmul, (ka,kb))
            # Issue #456: 2nd derivative for MLP
            mmulABt(w)=A_mul_Bt(w[1],w[2])
            mmulAtB(w)=At_mul_B(w[1],w[2])
            mmulAtBt(w)=At_mul_Bt(w[1],w[2])
            @test gradcheck(mmulABt, (ka,kb'))
            @test gradcheck(mmulAtB, (ka',kb))
            @test gradcheck(mmulAtBt, (ka',kb'))
        end

        if CUDA.functional()
            # cannot gradcheck axpy!, scal! overwriting
            d = rand(t,3,5)
            kd = KnetArray(d)
            r = rand()
            @test isapprox(axpy!(r,a,d), axpy!(r,ka,kd))
            @test isapprox(lmul!(r,d), lmul!(r,kd))
            @test isapprox(rmul!(d,r), rmul!(kd,r))
        end

        @test gradcheck(transpose, a)
        at = a'
        bt = b'
        s = a[1]
        mmul1(w)=w[1]'*w[2]
        mmul2(w)=w[1]*w[2]'
        mmul3(w)=w[1]'*w[2]'
        @test gradcheck(mmul1, Any[at,b])
        @test gradcheck(mmul2, Any[a,bt])
        @test gradcheck(mmul3, Any[at,bt])
        @test gradcheck(mmul1, Any[a,s])
        @test gradcheck(mmul1, Any[s,b])
        @test gradcheck(mmul2, Any[a,s])
        @test gradcheck(mmul2, Any[s,b])
        @test gradcheck(mmul3, Any[a,s])
        @test gradcheck(mmul3, Any[s,b])
        if CUDA.functional()
            kat = ka'
            kbt = kb'
            @test isapprox(at, Array(kat))
            @test gradcheck(transpose, ka)
            @test isapprox(kat'*kb, at'*b)
            @test isapprox(ka*kbt', a*bt')
            @test isapprox(kat'*kbt', at'*bt')
            @test gradcheck(mmul1, Any[kat,kb])
            @test gradcheck(mmul2, Any[ka,kbt])
            @test gradcheck(mmul3, Any[kat,kbt])
            @test gradcheck(mmul1, Any[ka,s])
            @test gradcheck(mmul1, Any[s,kb])
            @test gradcheck(mmul2, Any[ka,s])
            @test gradcheck(mmul2, Any[s,kb])
            @test gradcheck(mmul3, Any[ka,s])
            @test gradcheck(mmul3, Any[s,kb])
        end

        @test gradcheck(mat, a)
        if CUDA.functional()
            @test isapprox(mat(a), Array(mat(ka)))
            @test gradcheck(mat, ka)
        end

        for p in collect(permutas(1:2))
            p2(x) = permutedims(x,p)
            @test gradcheck(p2, a)
            if CUDA.functional()
                @test isapprox(p2(a), Array(p2(ka)))
                @test gradcheck(p2, ka)
            end
        end

        a3 = rand(2,3,4)
        if CUDA.functional(); k3 = KnetArray(a3); end
        for p in collect(permutas(1:3))
            p3(x) = permutedims(x,p)
            @test gradcheck(p3, a3)
            if CUDA.functional()
                @test isapprox(p3(a3), Array(p3(k3)))
                @test gradcheck(p3, k3)
            end
        end

        a4 = rand(2,3,4,5)
        if CUDA.functional(); k4 = KnetArray(a4); end
        for p in nsample(permutas(1:4),6)
            p4(x) = permutedims(x,p)
            @test gradcheck(p4, a4)
            if CUDA.functional()
                @test isapprox(p4(a4), Array(p4(k4)))
                @test gradcheck(p4, k4)
            end
        end

        a5 = rand(2,3,4,5,6)
        if CUDA.functional(); k5 = KnetArray(a5); end
        for p in nsample(permutas(1:5),6)
            p5(x) = permutedims(x,p)
            @test gradcheck(p5, a5)
            if CUDA.functional()
                @test isapprox(p5(a5), Array(p5(k5)))
                @test gradcheck(p5, k5)
            end
        end

        # transpose and matmul should work with vectors
        if CUDA.functional()
            a6 = rand(t,3)
            b6 = rand(t,3,3)
            c6 = KnetArray(a6)
            d6 = KnetArray(b6)
            @test c6' ≈ a6'
            @test c6' * c6 ≈ a6' * a6 # (a6'*a6)::Number, (Array(a6')*a6)::Vector (not Matrix), (c6'*c6)::Number.
            @test c6 * c6' ≈ a6 * a6'
            @test b6 * a6 ≈ d6 * c6 # Matrix * Vector => Vector
        end
    end
end

nothing


text:
"""
This module provides `lsqfit` methods for fitting an equation of state
with(out) units.
"""
module NonlinearFitting

using ConstructionBase: constructorof
using LsqFit: curve_fit
using Unitful: AbstractQuantity, upreferred, ustrip, unit

import ..EquationForm
using ..Collections

export lsqfit

"""
    lsqfit(form, eos, xdata, ydata; debug = false, kwargs...)

Fit an equation of state using least-squares fitting method (with the Levenberg-Marquardt algorithm).

# Arguments
- `form::EquationForm`: an `EquationForm` instance. If `EnergyForm`, fit ``E(V)``; if `PressureForm`, fit ``P(V)``; if `BulkModulusForm`, fit ``B(V)``.
- `eos::EquationOfState`: a trial equation of state. If it has units, `xdata` and `ydata` must also have.
- `xdata::AbstractVector`: a vector of volumes (``V``), with(out) units.
- `ydata::AbstractVector`: a vector of energies (``E``), pressures (``P``), or bulk moduli (``B``), with(out) units. It must be consistent with `form`.
- `debug::Bool=false`: if `true`, then an `LsqFit.LsqFitResult` is returned, containing estimated Jacobian, residuals, etc.; if `false`, a fitted `EquationOfState` is returned. The default value is `false`.
- `kwargs`: the rest keyword arguments are the same as that of `LsqFit.curve_fit`. See its [documentation](https://github.com/JuliaNLSolvers/LsqFit.jl/blob/master/README.md)
    and [tutorial](https://julianlsolvers.github.io/LsqFit.jl/latest/tutorial/).
"""
function lsqfit(
    form::EquationForm,
    eos::EquationOfState{<:Real},
    xdata::AbstractVector{<:Real},
    ydata::AbstractVector{<:Real};
    debug = false,
    kwargs...,
)
    E = constructorof(typeof(eos))  # Get the `UnionAll` type
    model = (x, p) -> map(apply(form, E(p...)), x)
    fitted = curve_fit(
        model,
        float(xdata),  # Convert `xdata` elements to floats
        float(ydata),  # Convert `ydata` elements to floats
        float(Collections.fieldvalues(eos));  # TODO: What if these floats are different types?
        kwargs...,
    )
    return debug ? fitted : E(fitted.param...)
end  # function lsqfit
function lsqfit(
    form::EquationForm,
    eos::EquationOfState{<:AbstractQuantity},
    xdata::AbstractVector{<:AbstractQuantity},
    ydata::AbstractVector{<:AbstractQuantity};
    kwargs...,
)
    E = constructorof(typeof(eos))  # Get the `UnionAll` type
    values = Collections.fieldvalues(eos)
    original_units = unit.(values)  # Keep a record of `eos`'s units
    f = x -> map(ustrip ∘ upreferred, x)  # Convert to preferred units and strip the unit
    trial_params = f.(values)
    result = lsqfit(form, E(trial_params...), f.(xdata), f.(ydata); kwargs...)
    if result isa EquationOfState  # i.e., if `debug = false` and no error is thrown
        data = Collections.fieldvalues(result)
        # Convert back to original `eos`'s units
        return E([data[i] * upreferred(u) |> u for (i, u) in enumerate(original_units)]...)
    end
    return result
end  # function lsqfit

end


text:
using LinearAlgebra, Zygote, Optim, IterativeSolvers, Arpack
using CSV, DataFrames

include("model.jl")

function data()
    abalone = CSV.File(open("warpedgp/abalone.csv"))
    abalone = DataFrame(abalone)
    X = Matrix(Matrix(abalone[1:200,2:8])')
    y = Vector{Float64}(abalone[1:200, 9])
    f = Linear()
    g = Identity()
    k = Gaussian()
    d = data(f, X, y)
    return f, g, k, d
end

function testfun()
    f, g, k, d = data()
    unpack = x -> (x[1],)
    fg! = make_fg(unpack) do ℓ
       ks = kernelsystem(f, k, g, d, UniformScaling(ℓ), (1.0,), 1e-6, ())
       -logprob(ks, (), 0)
    end
    return Optim.optimize(Optim.only_fg!(fg!), rand(1), ConjugateGradient())
end

function datavanilla()
    abalone = CSV.File(open("warpedgp/abalone.csv"))
    abalone = DataFrame(abalone)
    X = Matrix(Matrix(abalone[1:200,2:6])')
    y = Vector(abalone[1:200, 7])
    f = Linear()
    g = Identity()
    k = Gaussian()
    d = data(f, X, y)
    return f, g, k, d
end

function testfunvanilla()
    f, g, k, d = datavanilla()
    unpack = x -> (x[1:5], x[6], x[7])
    fg! = make_fg(unpack) do ℓ, ϵ, amp
       ks = kernelsystem(f, k, g, d, Diagonal(ℓ), (amp,), ϵ, ())
       -logprob(ks, (), 0)
    end
    return Optim.optimize(Optim.only_fg!(fg!), rand(7), ConjugateGradient())
end

function datavanillanoτ()
    abalone = CSV.File(open("abalone.csv"))
    abalone = DataFrame(abalone)
    X = Matrix(Matrix(abalone[1:200,2:6])')
    y = Vector(abalone[1:200, 7])
    f = Linear()
    g = Identity()
    k = Gaussian()
    d = data(f, X, y)
    return f, g, k, d
end

function testfunvanillanoτ()
    f, g, k, d = datavanillanoτ()
    unpack = x -> (x[1:5], x[6])
    fg! = make_fg(unpack) do ℓ, ϵ
       ks = kernelsystem(f, k, g, d, Diagonal(ℓ), (1 / (1 + ϵ),), ϵ / (1 + ϵ), ())
       -logprob(ks, (), 0)
    end
    return Optim.optimize(Optim.only_fg!(fg!), rand(6), ConjugateGradient())
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using ArgParse
using DrWatson
@quickactivate
using GenerativeAD

s = ArgParseSettings()
@add_arg_table! s begin
    "exp_prefix"
		arg_type = String
		default = "experiments/tabular"
		help = "Data prefix of experiment files."
	"eval_prefix"
		arg_type = String
		default = "evaluation/tabular"
		help = "Data prefix of evaluation files."
	"-d", "--dry_run"
    	action = :store_true
		help = "Compute diff but do not delete."
end

"""
	synchronize_stats(exp_prefix::String, eval_prefix::String; dry_run=false)

Collects filesnames from evaluation datadir in `eval_prefix` and compares it to
contents of experiments in `exp_prefix` datadir. If `dry_run` is true the function does
not delete the evaluation files for which there are no longer experiment files.
"""
function synchronize_stats(exp_prefix::String, eval_prefix::String; dry_run=false)
    exp_dir = datadir(exp_prefix)
    eval_dir = datadir(eval_prefix)

    @info "Collecting files from $exp_dir and $eval_dir folders."
    exp_files = GenerativeAD.Evaluation.collect_files(exp_dir)
    filter!(x -> !startswith(basename(x), "model"), exp_files)
    eval_files = GenerativeAD.Evaluation.collect_files(eval_dir)

    @info "Collected $(length(exp_files)) files from experiment folder."
    @info "Collected $(length(eval_files)) files from evaluation folder."

    # reverse naming of exp files to evaluation
    exp_eval_files = joinpath.(dirname.(exp_files), "eval_" .* basename.(exp_files))
    exp_eval_files = replace.(exp_eval_files, exp_prefix => eval_prefix)

    dif = setdiff(eval_files, exp_eval_files)
    @info "$(length(dif)) evaluation files lack underlying experiment files."

    if !dry_run
        @info "Removing files."
        rm.(dif)
    end
end

function main(args)
    @unpack exp_prefix, eval_prefix, dry_run = args
	synchronize_stats(exp_prefix, eval_prefix; dry_run=dry_run)
	@info "---------------- DONE -----------------"
end

main(parse_args(ARGS, s))


text:
function dispersion(series)
    S¹ = eltype(series)(0)
    S² = eltype(series)(0)
    D, M = 0, 0
    for (n,a) in enumerate(series)
        S¹ += a
        S² += a^2
        M = S¹/n
        D = S²/n-M^2
    end
    return D, M
end

function currentstd(series, n::Int)
    a = series[1:n]
    D, l = dispersion(a)
    return D
end

function test_currentstd()
    series = randn(Float64, 50)
    D = currentstd(series, 50)
    plot(series)
    return D
end


text:
"""
    Context

Stores common summary information for all Imputor types.

# Fields
* `num::Int`: number of observations
* `count::Int`: number of missing values found
* `limit::Float64`: allowable limit for missing values to impute
* `missing::Function`: returns a Bool if the value counts as missing.
"""
mutable struct Context
    num::Int
    count::Int
    limit::Float64
    missing::Function
end

Context(limit::Float64, missing::Function=ismissing) = Context(0, 0, limit, missing)

Base.copy(x::Context) = Context(x.num, x.count, x.limit, x.missing)

"""
    ismissing(ctx::Context, x) -> Bool

Uses `ctx.missing` to determine if x is missing. If x is a data row or an abstract array
then `ismissing` will return true if `ctx.missing` returns true for any element.
The ctx.count is increased whenever whenever we return true and if `ctx.count / ctx.num`
exceeds our `ctx.limit` we throw an `ImputeError`

# Arguments
* `ctx::Context`: the contextual information about missing information.
* `x`: the value to check (may be an single values, abstract array or row)
"""
function Base.ismissing(ctx::Context, x)
    missing = if isa(x, DataFrameRow)
        any(entry -> ctx.missing(entry[2]), x)
    elseif isa(x, AbstractArray)
        any(ctx.missing, x)
    else
        ctx.missing(x)
    end

    if missing
        ctx.count += 1

        if (ctx.count / ctx.num) > ctx.limit
            throw(ImputeError(
                "More than $(ctx.limit * 100)% of values were missing ()."
            ))
        end

        return true
    else
        return false
    end
end

"""
    findfirst(ctx::Context, data::AbstractVector) -> Int

Returns the first not missing index in `data`.

# Arguments
* `ctx::Context`: the context to pass into `ismissing`
* `data::AbstractVector`: the data array to search

# Returns
* `Int`: the first index in `data` that isn't missing
"""
function Base.findfirst(ctx::Context, data::AbstractVector)
    return findfirst(x -> !ismissing(ctx, x), data)
end

"""
    findlast(ctx::Context, data::AbstractVector) -> Int

Returns the last not missing index in `data`.

# Arguments
* `ctx::Context`: the context to pass into `ismissing`
* `data::AbstractVector`: the data array to search

# Returns
* `Int`: the last index in `data` that isn't missing
"""
function Base.findlast(ctx::Context, data::AbstractVector)
    return findlast(x -> !ismissing(ctx, x), data)
end

"""
    findnext(ctx::Context, data::AbstractVector) -> Int

Returns the next not missing index in `data`.

# Arguments
* `ctx::Context`: the context to pass into `ismissing`
* `data::AbstractVector`: the data array to search

# Returns
* `Int`: the next index in `data` that isn't missing
"""
function Base.findnext(ctx::Context, data::AbstractVector, idx::Int)
    return findnext(x -> !ismissing(ctx, x), data, idx)
end


text:
using Base.Test
using QuantumOptics

@testset "semiclassical" begin

# Test equality
b = GenericBasis(3)
x1 = semiclassical.State(basisstate(b, 1), [0.5im])
x2 = semiclassical.State(basisstate(b, 1), [0.6im])
x3 = semiclassical.State(basisstate(b, 2), [0.5im])

@test x1 == x1
@test x1 != x2
@test x1 != x3

x1 = semiclassical.State(dm(basisstate(b, 1)), [0.5im])
x2 = semiclassical.State(dm(basisstate(b, 1)), [0.6im])
x3 = semiclassical.State(dm(basisstate(b, 2)), [0.5im])

@test x1 == x1
@test x1 != x2
@test x1 != x3

# Test expect and variance
b = FockBasis(10)
a = destroy(b)
n = number(b)

alpha = complex(0.4, 0.3)
psi = coherentstate(b, alpha)
rho = dm(psi)

x_ket = semiclassical.State(psi, [complex(1., 0.)])
x_rho = semiclassical.State(rho, [complex(1., 0.)])

@test expect(a, x_ket) ≈ alpha
@test expect(a, x_rho) ≈ alpha
@test expect(a, [x_ket, x_rho]) ≈ [alpha, alpha]
@test variance(n, x_ket) ≈ abs2(alpha)
@test variance(n, x_rho) ≈ abs2(alpha)
@test variance(n, [x_ket, x_rho]) ≈ [abs2(alpha), abs2(alpha)]

# Test partial trace
b1 = GenericBasis(3)
b2 = GenericBasis(5)
b = b1 ⊗ b2
psi = randstate(b)
x = semiclassical.State(psi, [0.4, -0.3im])
@test ptrace(x, 1).quantum == ptrace(psi, 1)
@test ptrace(x, [2]).quantum == ptrace(psi, 2)

rho = randoperator(b)
x = semiclassical.State(rho, [0.4, -0.3im])
@test ptrace(x, 1).quantum == ptrace(rho, 1)
@test ptrace(x, [2]).quantum == ptrace(rho, 2)

# Test dm function
b = GenericBasis(4)
psi = randstate(b)
u = Complex128[complex(2., 3.)]
state = semiclassical.State(psi, u)
@test dm(state) == semiclassical.State(dm(psi), u)

# Test casting between and semiclassical states and complex vectors
b = GenericBasis(4)
rho = randoperator(b)
u = rand(Complex128, 3)
state = semiclassical.State(rho, u)
state_ = semiclassical.State(randoperator(b), rand(Complex128, 3))
x = Vector{Complex128}(16+3)
semiclassical.recast!(state, x)
semiclassical.recast!(x, state_)
@test state_ == state

# Test master
basis = SpinBasis(1//2)

# Random 2 level Hamiltonian
a1 = 0.5
a2 = 1.9
c = 1.3
d = -4.7

data = [a1 c-1im*d; c+1im*d a2]
H = DenseOperator(basis, data)

a = (a1 + a2)/2
b = (a1 - a2)/2
r = [c d b]

sigma_r = c*sigmax(basis) + d*sigmay(basis) + b*sigmaz(basis)

U(t) = exp(-1im*a*t)*(cos(norm(r)*t)*one(basis) - 1im*sin(norm(r)*t)*sigma_r/norm(r))

# Random initial state
psi0 = randstate(basis)

T = [0:0.5:1;]

fquantum_schroedinger(t, rho, u) = H
fquantum_master(t, rho, u) = H, [], []
function fclassical(t, quantumstate, u, du)
    du[1] = -1*u[1]
end

state0 = semiclassical.State(psi0, Complex128[complex(2., 3.)])
function f(t, state::semiclassical.State{Ket})
    @test 1e-5 > norm(state.quantum - U(t)*psi0)
    @test 1e-5 > abs(state.classical[1] - state0.classical[1]*exp(-t))
end
semiclassical.schroedinger_dynamic(T, state0, fquantum_schroedinger, fclassical; fout=f)
tout, state_t = semiclassical.schroedinger_dynamic(T, state0, fquantum_schroedinger, fclassical)
f(T[end], state_t[end])

function f(t, state::semiclassical.State{DenseOperator})
    @test 1e-5 > tracedistance(state.quantum, dm(U(t)*psi0))
    @test 1e-5 > abs(state.classical[1] - state0.classical[1]*exp(-t))
end
semiclassical.master_dynamic(T, state0, fquantum_master, fclassical; fout=f)
tout, state_t = semiclassical.master_dynamic(T, state0, fquantum_master, fclassical)
f(T[end], state_t[end])

end # testset


text:
using TypeEmulator
using Base.Test

@testset "metadata"     begin include("metadata.jl") end
@testset "types"        begin include("types.jl") end
@testset "linear"       begin include("linear.jl") end
@testset "conversions"  begin include("typeconversions.jl") end
@testset "relations"    begin include("relations.jl") end
@testset "subtype"      begin include("subtype.jl") end


text:
function invariants_check(x)
prim=zeros(10,1);
pv=zeros(21,1);

v=zeros(144,1);
prim[1]=x[1] + x[2] + x[3] + x[4]

prim[2]=x[5] + x[6] + x[7] + x[8] + x[9] + x[10]

prim[3]=x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2

prim[4]=x[5]^2 + x[6]^2 + x[7]^2 + x[8]^2 + x[9]^2 + x[10]^2

prim[5]=x[1]*x[5] + x[1]*x[6] + x[1]*x[7] + x[2]*x[5] + x[2]*x[8] + x[2]*x[9] + x[3]*x[6] + x[3]*x[8] + x[3]*x[10] + x[4]*x[7] +
    x[4]*x[9] + x[4]*x[10]

prim[6]=x[1]^3 + x[2]^3 + x[3]^3 + x[4]^3

prim[7]=x[5]*x[6]*x[7] + x[5]*x[8]*x[9] + x[6]*x[8]*x[10] + x[7]*x[9]*x[10]

prim[8]=x[5]^3 + x[6]^3 + x[7]^3 + x[8]^3 + x[9]^3 + x[10]^3

prim[9]=x[5]^4 + x[6]^4 + x[7]^4 + x[8]^4 + x[9]^4 + x[10]^4

prim[10]=x[1]^4 + x[2]^4 + x[3]^4 + x[4]^4 + x[5]^3*x[6] + x[5]^3*x[7] + x[5]^3*x[8] + x[5]^3*x[9] + x[5]*x[6]^3 +
    x[5]*x[7]^3 + x[5]*x[8]^3 + x[5]*x[9]^3 + x[6]^3*x[7] + x[6]^3*x[8] + x[6]^3*x[10] + x[6]*x[7]^3 +
    x[6]*x[8]^3 + x[6]*x[10]^3 + x[7]^3*x[9] + x[7]^3*x[10] + x[7]*x[9]^3 + x[7]*x[10]^3 + x[8]^3*x[9] +
    x[8]^3*x[10] + x[8]*x[9]^3 + x[8]*x[10]^3 + x[9]^3*x[10] + x[9]*x[10]^3

pv[1]=x[5]*x[6] + x[5]*x[7] + x[5]*x[8] + x[5]*x[9] + x[6]*x[7] + x[6]*x[8] + x[6]*x[10] + x[7]*x[9] + x[7]*x[10] + x[8]*x[9]  + x[8]*x[10] + x[9]*x[10]

pv[2]=x[1]*x[2]*x[10] + x[1]*x[3]*x[9] + x[1]*x[4]*x[8] + x[2]*x[3]*x[7] + x[2]*x[4]*x[6] + x[3]*x[4]*x[5]

pv[3]=x[1]*x[2]*x[5] + x[1]*x[3]*x[6] + x[1]*x[4]*x[7] + x[2]*x[3]*x[8] + x[2]*x[4]*x[9] + x[3]*x[4]*x[10]

pv[4]=x[1]*x[8]^2 + x[1]*x[9]^2 + x[1]*x[10]^2 + x[2]*x[6]^2 + x[2]*x[7]^2 + x[2]*x[10]^2 + x[3]*x[5]^2 + x[3]*x[7]^2  + x[3]*x[9]^2 + x[4]*x[5]^2 + x[4]*x[6]^2 + x[4]*x[8]^2

pv[5]=x[1]*x[5]*x[8] + x[1]*x[5]*x[9] + x[1]*x[6]*x[8] + x[1]*x[6]*x[10] + x[1]*x[7]*x[9] + x[1]*x[7]*x[10] + x[2]*x[5]*x[6] +
    x[2]*x[5]*x[7] + x[2]*x[6]*x[8] + x[2]*x[7]*x[9] + x[2]*x[8]*x[10] + x[2]*x[9]*x[10] + x[3]*x[5]*x[6] + x[3]*x[5]*x[8] + x[3]*x[6]*x[7] + x[3]*x[7]*x[10] + x[3]*x[8]*x[9] + x[3]*x[9]*x[10] + x[4]*x[5]*x[7] + x[4]*x[5]*x[9] +
    x[4]*x[6]*x[7] + x[4]*x[6]*x[10] + x[4]*x[8]*x[9] + x[4]*x[8]*x[10]

pv[6]=x[5]^2*x[6] + x[5]^2*x[7] + x[5]^2*x[8] + x[5]^2*x[9] + x[5]*x[6]^2 + x[5]*x[7]^2 + x[5]*x[8]^2 + x[5]*x[9]^2 +
    x[6]^2*x[7] + x[6]^2*x[8] + x[6]^2*x[10] + x[6]*x[7]^2 + x[6]*x[8]^2 + x[6]*x[10]^2 + x[7]^2*x[9] +
    x[7]^2*x[10] + x[7]*x[9]^2 + x[7]*x[10]^2 + x[8]^2*x[9] + x[8]^2*x[10] + x[8]*x[9]^2 + x[8]*x[10]^2 +
    x[9]^2*x[10] + x[9]*x[10]^2

pv[7]=x[1]^2*x[2]*x[3] + x[1]^2*x[2]*x[4] + x[1]^2*x[3]*x[4] + x[1]*x[2]^2*x[3] + x[1]*x[2]^2*x[4] + x[1]*x[2]*x[3]^2 +
    x[1]*x[2]*x[4]^2 + x[1]*x[3]^2*x[4] + x[1]*x[3]*x[4]^2 + x[2]^2*x[3]*x[4] + x[2]*x[3]^2*x[4] + x[2]*x[3]*x[4]^2

pv[8]=x[1]^2*x[2]*x[10] + x[1]^2*x[3]*x[9] + x[1]^2*x[4]*x[8] + x[1]*x[2]^2*x[10] + x[1]*x[3]^2*x[9] + x[1]*x[4]^2*x[8] +
    x[2]^2*x[3]*x[7] + x[2]^2*x[4]*x[6] + x[2]*x[3]^2*x[7] + x[2]*x[4]^2*x[6] + x[3]^2*x[4]*x[5] + x[3]*x[4]^2*x[5]

pv[9]=x[1]*x[2]*x[5]*x[10] + x[1]*x[3]*x[6]*x[9] + x[1]*x[4]*x[7]*x[8] + x[2]*x[3]*x[7]*x[8] + x[2]*x[4]*x[6]*x[9] +
    x[3]*x[4]*x[5]*x[10]

pv[10]=x[1]*x[2]*x[6]*x[10] + x[1]*x[2]*x[7]*x[10] + x[1]*x[2]*x[8]*x[10] + x[1]*x[2]*x[9]*x[10] + x[1]*x[3]*x[5]*x[9] +
    x[1]*x[3]*x[7]*x[9] + x[1]*x[3]*x[8]*x[9] + x[1]*x[3]*x[9]*x[10] + x[1]*x[4]*x[5]*x[8] + x[1]*x[4]*x[6]*x[8] +
    x[1]*x[4]*x[8]*x[9] + x[1]*x[4]*x[8]*x[10] + x[2]*x[3]*x[5]*x[7] + x[2]*x[3]*x[6]*x[7] + x[2]*x[3]*x[7]*x[9] +
    x[2]*x[3]*x[7]*x[10] + x[2]*x[4]*x[5]*x[6] + x[2]*x[4]*x[6]*x[7] + x[2]*x[4]*x[6]*x[8] + x[2]*x[4]*x[6]*x[10] +
    x[3]*x[4]*x[5]*x[6] + x[3]*x[4]*x[5]*x[7] + x[3]*x[4]*x[5]*x[8] + x[3]*x[4]*x[5]*x[9]

pv[11]=x[1]*x[2]*x[5]^2 + x[1]*x[3]*x[6]^2 + x[1]*x[4]*x[7]^2 + x[2]*x[3]*x[8]^2 + x[2]*x[4]*x[9]^2 + x[3]*x[4]*x[10]^2

pv[12]=x[1]*x[5]*x[10]^2 + x[1]*x[6]*x[9]^2 + x[1]*x[7]*x[8]^2 + x[2]*x[5]*x[10]^2 + x[2]*x[6]^2*x[9] + x[2]*x[7]^2*x[8] +
    x[3]*x[5]^2*x[10] + x[3]*x[6]*x[9]^2 + x[3]*x[7]^2*x[8] + x[4]*x[5]^2*x[10] + x[4]*x[6]^2*x[9] +
    x[4]*x[7]*x[8]^2

pv[13]=x[1]*x[5]*x[8]^2 + x[1]*x[5]*x[9]^2 + x[1]*x[6]*x[8]^2 + x[1]*x[6]*x[10]^2 + x[1]*x[7]*x[9]^2 + x[1]*x[7]*x[10]^2 +
    x[2]*x[5]*x[6]^2 + x[2]*x[5]*x[7]^2 + x[2]*x[6]^2*x[8] + x[2]*x[7]^2*x[9] + x[2]*x[8]*x[10]^2 +
    x[2]*x[9]*x[10]^2 + x[3]*x[5]^2*x[6] + x[3]*x[5]^2*x[8] + x[3]*x[6]*x[7]^2 + x[3]*x[7]^2*x[10] +
    x[3]*x[8]*x[9]^2 + x[3]*x[9]^2*x[10] + x[4]*x[5]^2*x[7] + x[4]*x[5]^2*x[9] + x[4]*x[6]^2*x[7] +
    x[4]*x[6]^2*x[10] + x[4]*x[8]^2*x[9] + x[4]*x[8]^2*x[10]

pv[14]=x[1]*x[8]^2*x[9] + x[1]*x[8]^2*x[10] + x[1]*x[8]*x[9]^2 + x[1]*x[8]*x[10]^2 + x[1]*x[9]^2*x[10] + x[1]*x[9]*x[10]^2  + x[2]*x[6]^2*x[7] + x[2]*x[6]^2*x[10] + x[2]*x[6]*x[7]^2 + x[2]*x[6]*x[10]^2 + x[2]*x[7]^2*x[10] +
    x[2]*x[7]*x[10]^2 + x[3]*x[5]^2*x[7] + x[3]*x[5]^2*x[9] + x[3]*x[5]*x[7]^2 + x[3]*x[5]*x[9]^2 + x[3]*x[7]^2*x[9] + x[3]*x[7]*x[9]^2 + x[4]*x[5]^2*x[6] + x[4]*x[5]^2*x[8] + x[4]*x[5]*x[6]^2 + x[4]*x[5]*x[8]^2 +
    x[4]*x[6]^2*x[8] + x[4]*x[6]*x[8]^2

pv[15]=x[1]^2*x[2]*x[3]*x[5] + x[1]^2*x[2]*x[3]*x[6] + x[1]^2*x[2]*x[4]*x[5] + x[1]^2*x[2]*x[4]*x[7] + x[1]^2*x[3]*x[4]*x[6] +
    x[1]^2*x[3]*x[4]*x[7] + x[1]*x[2]^2*x[3]*x[5] + x[1]*x[2]^2*x[3]*x[8] + x[1]*x[2]^2*x[4]*x[5] +
    x[1]*x[2]^2*x[4]*x[9] + x[1]*x[2]*x[3]^2*x[6] + x[1]*x[2]*x[3]^2*x[8] + x[1]*x[2]*x[4]^2*x[7] +
    x[1]*x[2]*x[4]^2*x[9] + x[1]*x[3]^2*x[4]*x[6] + x[1]*x[3]^2*x[4]*x[10] + x[1]*x[3]*x[4]^2*x[7] +
    x[1]*x[3]*x[4]^2*x[10] + x[2]^2*x[3]*x[4]*x[8] + x[2]^2*x[3]*x[4]*x[9] + x[2]*x[3]^2*x[4]*x[8] +
    x[2]*x[3]^2*x[4]*x[10] + x[2]*x[3]*x[4]^2*x[9] + x[2]*x[3]*x[4]^2*x[10]

pv[16]=x[1]^2*x[2]*x[8]*x[10] + x[1]^2*x[2]*x[9]*x[10] + x[1]^2*x[3]*x[8]*x[9] + x[1]^2*x[3]*x[9]*x[10] + x[1]^2*x[4]*x[8]*x[9] + x[1]^2*x[4]*x[8]*x[10] + x[1]*x[2]^2*x[6]*x[10] + x[1]*x[2]^2*x[7]*x[10] + x[1]*x[3]^2*x[5]*x[9] +
    x[1]*x[3]^2*x[7]*x[9] + x[1]*x[4]^2*x[5]*x[8] + x[1]*x[4]^2*x[6]*x[8] + x[2]^2*x[3]*x[6]*x[7] +
    x[2]^2*x[3]*x[7]*x[10] + x[2]^2*x[4]*x[6]*x[7] + x[2]^2*x[4]*x[6]*x[10] + x[2]*x[3]^2*x[5]*x[7] +
    x[2]*x[3]^2*x[7]*x[9] + x[2]*x[4]^2*x[5]*x[6] + x[2]*x[4]^2*x[6]*x[8] + x[3]^2*x[4]*x[5]*x[7] +
    x[3]^2*x[4]*x[5]*x[9] + x[3]*x[4]^2*x[5]*x[6] + x[3]*x[4]^2*x[5]*x[8]

pv[17]=x[1]*x[2]*x[5]^2*x[10] + x[1]*x[3]*x[6]^2*x[9] + x[1]*x[4]*x[7]^2*x[8] + x[2]*x[3]*x[7]*x[8]^2 + x[2]*x[4]*x[6]*x[9]^2 +
    x[3]*x[4]*x[5]*x[10]^2

pv[18]=x[1]*x[2]*x[5]*x[6]*x[10] + x[1]*x[2]*x[5]*x[7]*x[10] + x[1]*x[2]*x[5]*x[8]*x[10] + x[1]*x[2]*x[5]*x[9]*x[10] +
    x[1]*x[3]*x[5]*x[6]*x[9] + x[1]*x[3]*x[6]*x[7]*x[9] + x[1]*x[3]*x[6]*x[8]*x[9] + x[1]*x[3]*x[6]*x[9]*x[10] +
    x[1]*x[4]*x[5]*x[7]*x[8] + x[1]*x[4]*x[6]*x[7]*x[8] + x[1]*x[4]*x[7]*x[8]*x[9] + x[1]*x[4]*x[7]*x[8]*x[10] +
    x[2]*x[3]*x[5]*x[7]*x[8] + x[2]*x[3]*x[6]*x[7]*x[8] + x[2]*x[3]*x[7]*x[8]*x[9] + x[2]*x[3]*x[7]*x[8]*x[10] +
    x[2]*x[4]*x[5]*x[6]*x[9] + x[2]*x[4]*x[6]*x[7]*x[9] + x[2]*x[4]*x[6]*x[8]*x[9] + x[2]*x[4]*x[6]*x[9]*x[10] +
    x[3]*x[4]*x[5]*x[6]*x[10] + x[3]*x[4]*x[5]*x[7]*x[10] + x[3]*x[4]*x[5]*x[8]*x[10] + x[3]*x[4]*x[5]*x[9]*x[10]

pv[19]=x[1]*x[5]*x[6]*x[9]^2 + x[1]*x[5]*x[6]*x[10]^2 + x[1]*x[5]*x[7]*x[8]^2 + x[1]*x[5]*x[7]*x[10]^2 + x[1]*x[6]*x[7]*x[8]^2  + x[1]*x[6]*x[7]*x[9]^2 + x[2]*x[5]*x[6]^2*x[9] + x[2]*x[5]*x[7]^2*x[8] + x[2]*x[5]*x[8]*x[10]^2 +
    x[2]*x[5]*x[9]*x[10]^2 + x[2]*x[6]^2*x[8]*x[9] + x[2]*x[7]^2*x[8]*x[9] + x[3]*x[5]^2*x[6]*x[10] +
    x[3]*x[5]^2*x[8]*x[10] + x[3]*x[6]*x[7]^2*x[8] + x[3]*x[6]*x[8]*x[9]^2 + x[3]*x[6]*x[9]^2*x[10] +
    x[3]*x[7]^2*x[8]*x[10] + x[4]*x[5]^2*x[7]*x[10] + x[4]*x[5]^2*x[9]*x[10] + x[4]*x[6]^2*x[7]*x[9] +
    x[4]*x[6]^2*x[9]*x[10] + x[4]*x[7]*x[8]^2*x[9] + x[4]*x[7]*x[8]^2*x[10]

pv[20]=x[1]*x[5]*x[8]*x[10]^2 + x[1]*x[5]*x[9]*x[10]^2 + x[1]*x[6]*x[8]*x[9]^2 + x[1]*x[6]*x[9]^2*x[10] + x[1]*x[7]*x[8]^2*x[9] + x[1]*x[7]*x[8]^2*x[10] + x[2]*x[5]*x[6]*x[10]^2 + x[2]*x[5]*x[7]*x[10]^2 + x[2]*x[6]^2*x[7]*x[9] +
    x[2]*x[6]^2*x[9]*x[10] + x[2]*x[6]*x[7]^2*x[8] + x[2]*x[7]^2*x[8]*x[10] + x[3]*x[5]^2*x[7]*x[10] +
    x[3]*x[5]^2*x[9]*x[10] + x[3]*x[5]*x[6]*x[9]^2 + x[3]*x[5]*x[7]^2*x[8] + x[3]*x[6]*x[7]*x[9]^2 +
    x[3]*x[7]^2*x[8]*x[9] + x[4]*x[5]^2*x[6]*x[10] + x[4]*x[5]^2*x[8]*x[10] + x[4]*x[5]*x[6]^2*x[9] +
    x[4]*x[5]*x[7]*x[8]^2 + x[4]*x[6]^2*x[8]*x[9] + x[4]*x[6]*x[7]*x[8]^2

pv[21]=x[5]^2*x[6]*x[10]^2 + x[5]^2*x[7]*x[10]^2 + x[5]^2*x[8]*x[10]^2 + x[5]^2*x[9]*x[10]^2 + x[5]*x[6]^2*x[9]^2 +
    x[5]*x[7]^2*x[8]^2 + x[6]^2*x[7]*x[9]^2 + x[6]^2*x[8]*x[9]^2 + x[6]^2*x[9]^2*x[10] + x[6]*x[7]^2*x[8]^2 +
    x[7]^2*x[8]^2*x[9] + x[7]^2*x[8]^2*x[10]

v[1] = 1
v[2] = pv[1]
v[3] = pv[2]
v[4] = pv[3]
v[5] = pv[4]
v[6] = pv[5]
v[7] = pv[6]
v[8] = pv[1]^2
v[9] = pv[7]
v[10] = pv[8]
v[11] = pv[9]
v[12] = pv[10]
v[13] = pv[11]
v[14] = pv[12]
v[15] = pv[13]
v[16] = pv[14]
v[17] = pv[1]*pv[5]
v[18] = pv[1]*pv[4]
v[19] = pv[1]*pv[3]
v[20] = pv[1]*pv[6]
v[21] = pv[1]*pv[2]
v[22] = pv[15]
v[23] = pv[16]
v[24] = pv[17]
v[25] = pv[18]
v[26] = pv[19]
v[27] = pv[20]
v[28] = pv[21]
v[29] = pv[4]^2
v[30] = pv[3]*pv[6]
v[31] = pv[5]*pv[6]
v[32] = pv[1]*pv[7]
v[33] = pv[4]*pv[6]
v[34] = pv[2]*pv[6]
v[35] = pv[1]*pv[9]
v[36] = pv[6]^2
v[37] = pv[1]*pv[8]
v[38] = pv[1]*pv[11]
v[39] = pv[1]^3
v[40] = pv[3]^2
v[41] = pv[3]*pv[5]
v[42] = pv[1]*pv[12]
v[43] = pv[1]*pv[14]
v[44] = pv[3]*pv[4]
v[45] = pv[2]*pv[3]
v[46] = pv[1]*pv[13]
v[47] = pv[5]*pv[7]
v[48] = pv[6]*pv[13]
v[49] = pv[4]*pv[13]
v[50] = pv[5]*pv[13]
v[51] = pv[5]*pv[10]
v[52] = pv[3]*pv[11]
v[53] = pv[3]*pv[8]
v[54] = pv[3]*pv[9]
v[55] = pv[1]^2*pv[3]
v[56] = pv[6]*pv[11]
v[57] = pv[6]*pv[8]
v[58] = pv[2]*pv[11]
v[59] = pv[1]*pv[18]
v[60] = pv[4]*pv[11]
v[61] = pv[1]*pv[16]
v[62] = pv[6]*pv[9]
v[63] = pv[3]*pv[12]
v[64] = pv[1]*pv[21]
v[65] = pv[6]*pv[12]
v[66] = pv[6]*pv[14]
v[67] = pv[1]*pv[20]
v[68] = pv[12]*pv[13]
v[69] = pv[10]*pv[14]
v[70] = pv[6]*pv[20]
v[71] = pv[10]^2
v[72] = pv[12]*pv[14]
v[73] = pv[4]*pv[17]
v[74] = pv[3]*pv[21]
v[75] = pv[2]*pv[17]
v[76] = pv[7]*pv[10]
v[77] = pv[9]*pv[14]
v[78] = pv[1]*pv[3]*pv[5]
v[79] = pv[7]*pv[12]
v[80] = pv[11]*pv[13]
v[81] = pv[1]^2*pv[13]
v[82] = pv[6]*pv[17]
v[83] = pv[7]*pv[9]
v[84] = pv[5]*pv[15]
v[85] = pv[1]*pv[3]*pv[4]
v[86] = pv[9]*pv[10]
v[87] = pv[5]*pv[21]
v[88] = pv[11]^2
v[89] = pv[9]^2
v[90] = pv[4]*pv[21]
v[91] = pv[6]*pv[21]
v[92] = pv[8]*pv[17]
v[93] = pv[14]*pv[18]
v[94] = pv[1]*pv[4]*pv[7]
v[95] = pv[1]*pv[3]*pv[11]
v[96] = pv[1]*pv[6]*pv[14]
v[97] = pv[1]^3*pv[3]
v[98] = pv[7]*pv[19]
v[99] = pv[3]*pv[6]^2
v[100] = pv[7]*pv[18]
v[101] = pv[10]*pv[21]
v[102] = pv[1]*pv[3]*pv[8]
v[103] = pv[11]*pv[15]
v[104] = pv[1]*pv[6]*pv[7]
v[105] = pv[1]*pv[5]*pv[10]
v[106] = pv[1]^2*pv[15]
v[107] = pv[12]*pv[21]
v[108] = pv[1]*pv[5]*pv[12]
v[109] = pv[13]*pv[20]
v[110] = pv[8]*pv[15]
v[111] = pv[9]*pv[21]
v[112] = pv[4]^2*pv[5]
v[113] = pv[1]*pv[4]*pv[10]
v[114] = pv[10]*pv[16]
v[115] = pv[1]^2*pv[21]
v[116] = pv[1]*pv[6]*pv[12]
v[117] = pv[1]*pv[5]*pv[11]
v[118] = pv[1]*pv[9]*pv[13]
v[119] = pv[1]*pv[10]^2
v[120] = pv[16]*pv[18]
v[121] = pv[4]*pv[6]*pv[10]
v[122] = pv[4]*pv[6]*pv[11]
v[123] = pv[1]^3*pv[8]
v[124] = pv[1]*pv[6]*pv[20]
v[125] = pv[15]*pv[20]
v[126] = pv[1]*pv[3]*pv[21]
v[127] = pv[5]*pv[6]*pv[12]
v[128] = pv[1]*pv[7]*pv[9]
v[129] = pv[17]*pv[18]
v[130] = pv[2]*pv[4]*pv[11]
v[131] = pv[1]*pv[2]*pv[20]
v[132] = pv[4]^2*pv[7]
v[133] = pv[2]*pv[3]*pv[17]
v[134] = pv[1]^2*pv[6]*pv[12]
v[135] = pv[5]*pv[8]*pv[13]
v[136] = pv[4]^2*pv[21]
v[137] = pv[1]*pv[9]*pv[19]
v[138] = pv[3]*pv[4]*pv[17]
v[139] = pv[1]^3*pv[17]
v[140] = pv[1]^3*pv[21]
v[141] = pv[5]*pv[11]*pv[21]
v[142] = pv[1]*pv[4]*pv[6]*pv[7]
v[143] = pv[2]*pv[11]*pv[18]
v[144] = pv[5]*pv[13]*pv[20]
return prim, v, pv

end

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs