text:
"""
    QuantisedDiatomic

This module exports two user facing functions:

- `generate_configurations`
    Creates a set of velocities and positions for diatomic molecule with specified
    vibrational `ν` and rotational `J` quantum numbers.

- `quantise_diatomic`
    Obtains vibrational `ν` and rotational `J` quantum numbers for a diatomic molecule
    with a given set of velocities and positions.

The central concept of this module is the EBK procedure which is nicely detailed here:
[Larkoski2006](@cite)

Inspired by VENUS96: [Hase1996](@cite)
"""
module QuantisedDiatomic

using QuadGK: QuadGK
using Rotations: Rotations
using LinearAlgebra: norm, normalize!, nullspace, cross, I
using UnitfulAtomic: austrip
using Distributions: Uniform
using UnicodePlots: lineplot, lineplot!, DotCanvas, histogram
using Optim: Optim
using Roots: Roots
using TimerOutputs: TimerOutputs, @timeit
using Interpolations: interpolate, BSpline, Cubic, Line, OnGrid, scale, hessian

using NQCDynamics: Simulation, Calculators, DynamicsUtils, masses
using NQCBase: Atoms, PeriodicCell, InfiniteCell
using NQCModels: NQCModels
using NQCModels.AdiabaticModels: AdiabaticModel

export generate_configurations
export quantise_diatomic

const TIMER = TimerOutputs.TimerOutput()

include("energy_evaluation.jl")
include("binding_curve.jl")
include("random_configuration.jl")

struct EffectivePotential{T,B,F}
    μ::T
    J::Int
    binding_curve::BindingCurve{T,B,F}
end

function (effective_potential::EffectivePotential)(r)
    (;μ, J, binding_curve) = effective_potential
    rotational = J*(J+1) / (2μ*r^2)

    potential = binding_curve.fit(r)
    return rotational + potential
end

struct RadialMomentum{T,B,F}
    total_energy::T
    V::EffectivePotential{T,B,F}
end

function (radial_momentum::RadialMomentum)(r)
    (;total_energy, V) = radial_momentum
    return sqrt_avoid_negatives(2V.μ * (total_energy - V(r)))
end

"""
    generate_configurations(sim, ν, J; samples=1000, height=10, normal_vector=[0, 0, 1],
        translational_energy=0, direction=[0, 0, -1], position=[0, 0, height])

Generate positions and momenta for given quantum numbers

`translational_energy`, `direction` and `position` specify the kinetic energy in
a specific direction with the molecule placed with centre of mass at `position`.

Keyword arguments `height` and `normal_vector` become relevant if the potential
requires specific placement of the molecule.
These allow the molecule to be placed at a distance `height` in the direction
`normal_vector` when performing potential evaluations.
"""
function generate_configurations(sim, ν, J;
    samples=1000,
    height=10.0,
    surface_normal=[0, 0, 1.0],
    translational_energy=0.0,
    direction=[0, 0, -1.0],
    atom_indices=[1, 2],
    r=zeros(size(sim)),
    bond_lengths=0.5:0.01:5.0
)

    μ = reduced_mass(masses(sim)[atom_indices])
    _, slab = separate_slab_and_molecule(atom_indices, r)
    generation = GenerationParameters(direction, austrip(translational_energy), samples)
    environment = EvaluationEnvironment(atom_indices, size(sim), slab, austrip(height), surface_normal)

    binding_curve = calculate_binding_curve(bond_lengths, sim.calculator.model, environment)
    plot_binding_curve(binding_curve.bond_lengths, binding_curve.potential, binding_curve.fit)

    V = EffectivePotential(μ, J, binding_curve)

    total_energy, bounds = find_total_energy(V, ν)

    radial_momentum = RadialMomentum(total_energy, V)
    bonds, momenta = select_random_bond_lengths(bounds, radial_momentum, samples)
    velocities = momenta ./ μ

    plot_distributions(bonds, velocities)

    @info "Generating the requested configurations..."
    configure_diatomic.(sim, bonds, velocities, J, environment, generation, μ)
end

function generate_1D_vibrations(model::AdiabaticModel, μ::Real, ν::Integer;
    samples=1000, bond_lengths=0.5:0.01:5.0
)

    J = 0
    environment = EvaluationEnvironment([1], (1,1), zeros(1,0), 0.0, [1.0])
    binding_curve = calculate_binding_curve(bond_lengths, model, environment)
    plot_binding_curve(binding_curve.bond_lengths, binding_curve.potential, binding_curve.fit)

    V = EffectivePotential(μ, J, binding_curve)

    total_energy, bounds = find_total_energy(V, ν)

    radial_momentum = RadialMomentum(total_energy, V)
    bonds, momenta = select_random_bond_lengths(bounds, radial_momentum, samples)
    velocities = momenta ./ μ

    plot_distributions(bonds, velocities)

    return bonds, velocities
end

"""
    sqrt_avoid_negatives(x)

Same as `sqrt` but returns `zero(x)` if `x` is negative.
Used here just in case the endpoints are a little off.
"""
function sqrt_avoid_negatives(x)
    if x < zero(x)
        @warn """
        `sqrt(x)` of negative number attempted: x = $x, returning zero instead.\
        This is acceptable if x is very close to zero. Otherwise something has gone wrong.
        """
        return zero(x)
    else
        return sqrt(x)
    end
end

"""
    find_total_energy(V::EffectivePotential, ν)

Returns the energy associated with the specified quantum numbers
"""
function find_total_energy(V::EffectivePotential, ν)

    μ = V.μ
    k = calculate_force_constant(V.binding_curve)
    E = guess_initial_energy(k, μ, ν, V.binding_curve.potential_minimum)

    "Evaluate vibrational quantum number."
    function nᵣ(E, r₁, r₂)
        kernel(r) = sqrt_avoid_negatives(E - V(r))
        integral, _ = QuadGK.quadgk(kernel, r₁, r₂; maxevals=100)
        return sqrt(2μ)/π * integral - 1/2
    end

    "Derivative of the above"
    function ∂nᵣ∂E(E, r₁, r₂)
        kernel(r) = 1 / sqrt_avoid_negatives(E - V(r))
        integral, _ = QuadGK.quadgk(kernel, r₁, r₂; maxevals=100)
        return sqrt(2μ)/2π * integral
    end

    "Second derivative of the above"
    function ∂²nᵣ∂E²(E, r₁, r₂)
        kernel(r) = 1 / sqrt_avoid_negatives(E - V(r))^3
        integral, _ = QuadGK.quadgk(kernel, r₁, r₂; maxevals=100)
        return -sqrt(2μ)/4π * integral
    end

    """
    Evaluates the function, gradient and hessian as described in the Optim documentation.
    """
    function fgh!(F,G,H,x)
        E = x[1]
        r₁, r₂ = find_integral_bounds(E, V)

        νtmp = nᵣ(E, r₁, r₂)
        difference = νtmp - ν

        if G !== nothing
            G[1] = ∂nᵣ∂E(E, r₁, r₂) * sign(difference)
        end
        if H !== nothing
            H[1] = ∂²nᵣ∂E²(E, r₁, r₂) * sign(difference)
        end
        if F !== nothing
            return abs(difference)
        end
    end

    @info "Converging the total energy to match the chosen quantum numbers..."
    optim = Optim.optimize(Optim.only_fgh!(fgh!), [E], Optim.Newton(),
        Optim.Options(show_trace=true, x_tol=1e-3)
    )
    E = Optim.minimizer(optim)[1]
    bounds = find_integral_bounds(E, V)

    E, bounds
end

function find_integral_bounds(total_energy::Real, V::EffectivePotential)
    energy_difference(r) = total_energy - V(r)

    r₀ = V.binding_curve.equilibrium_bond_length
    minimum_bond_length = first(V.binding_curve.bond_lengths)
    maximum_bond_length = last(V.binding_curve.bond_lengths)
    r₁ = @timeit TIMER "Lower bound" Roots.find_zero(energy_difference, (minimum_bond_length, r₀))
    r₂ = @timeit TIMER "Upper bound" Roots.find_zero(energy_difference, (r₀, maximum_bond_length))

    return r₁, r₂
end

"""
    quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix;
        height=10, normal_vector=[0, 0, 1])

Quantise the vibrational and rotational degrees of freedom for the specified
positions and velocities

When evaluating the potential, the molecule is moved to `height` in direction `normal_vector`.
If the potential is independent of centre of mass position, this has no effect.
Otherwise, be sure to modify these parameters to give the intended behaviour.
"""
function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix;
    height=10.0, surface_normal=[0, 0, 1.0], atom_indices=[1, 2],
    show_timer=false, reset_timer=false,
    bond_lengths=0.5:0.01:5.0
)

    reset_timer && TimerOutputs.reset_timer!(TIMER)

    r, slab = separate_slab_and_molecule(atom_indices, r)
    environment = EvaluationEnvironment(atom_indices, size(sim), slab, austrip(height), surface_normal)

    r_com = subtract_centre_of_mass(r, masses(sim)[atom_indices])
    v_com = subtract_centre_of_mass(v, masses(sim)[atom_indices])
    p_com = v_com .* masses(sim)[atom_indices]'

    k = DynamicsUtils.classical_kinetic_energy(sim, v_com)
    p = calculate_diatomic_energy(bond_length(r_com), sim.calculator.model, environment)
    E = k + p

    L = total_angular_momentum(r_com, p_com)
    J = round(Int, (sqrt(1+4*L^2) - 1) / 2) # L^2 = J(J+1)ħ^2

    μ = reduced_mass(masses(sim)[atom_indices])

    binding_curve = calculate_binding_curve(bond_lengths, sim.calculator.model, environment)

    V = EffectivePotential(μ, J, binding_curve)

    r₁, r₂ = @timeit TIMER "Finding bounds" find_integral_bounds(E, V)

    function nᵣ(E, r₁, r₂)
        kernel(r) = sqrt_avoid_negatives(E - V(r))
        integral, _ = QuadGK.quadgk(kernel, r₁, r₂; maxevals=100)
        return sqrt(2μ)/π * integral - 1/2
    end

    ν = @timeit TIMER "Calculating ν" nᵣ(E, r₁, r₂)

    TimerOutputs.complement!(TIMER)
    show_timer && show(TIMER)

    return round(Int, ν), J
end

subtract_centre_of_mass(x, m) = x .- centre_of_mass(x, m)
@views centre_of_mass(x, m) = (x[:,1].*m[1] .+ x[:,2].*m[2]) ./ (m[1]+m[2])
reduced_mass(atoms::Atoms) = reduced_mass(atoms.masses)

function reduced_mass(m::AbstractVector)
    if length(m) == 1
        return m[1]
    elseif length(m) == 2
        m[1]*m[2]/(m[1]+m[2])
    else
        throw(error("Mass vector of incorrect length."))
    end
end

@views bond_length(r) = norm(r[:,1] .- r[:,2])
@views total_angular_momentum(r, p) = norm(cross(r[:,1], p[:,1]) + cross(r[:,2], p[:,2]))

end # module


text:
## Do inference in parallel

using Distributed

## add number of processors desired

addprocs(5)

## Load relevant code into scope for each worker

@everywhere begin
    using Pkg;Pkg.activate(".")
end
@everywhere using JLD2,LinearAlgebra,TransformVariables,LogDensityProblems,MCMCChains,Dates
@everywhere import KenyaCoVSD

## Load data

@load("data/linelist_data_with_pos_neg_20feb_to_27apr_c_age.jld2")
@load("data/serological_data_with_20feb_to_31dec_age.jld2")
@load("data/deaths_data20210521.jld2")
@load("data/N_kenya.jld2")

## Parameter transformations and initial KE guess for HMC

trans = as((ct_min1 = as(Real, 0.4, 0.999),
        R₀ = as(Real, 0.0, 7.0),
        ϵ = as(Real, 0.7, 0.999),
        χ₁=as(Real, 0.1, 30.0),χ₂=as(Real, 0.1, 30.0),
        p_test₁=as(Real, 5e-6, 5e-3),p_test₂ = as(Real, 5e-6, 5e-3),
        p_choose1=as(Real, 0.1, 0.999),
        P_eff=as(Real, 0.5, 0.999),
        schooleffect = as(Real, 0.5, 0.999),
        extra_transmissibility = as(Real, 0.5, 2.5),
        influx_exposed_new_variant = as(Real, 0.0, 10e3),
        p_test_boost = as(Real, 1.0, 3.0),
        χ_boost= as(Real, 1.0, 3.0),
        E₀=as(Real, 0.0, 10e3)))

D = Diagonal(0.1*ones(15))

## Name for each county with informed priors for detection rates

semiurbancounties_lowserodata = ["Kajiado", "Kiambu", "Machakos","Isiolo" , "Kirinyaga", "Murang'a",
                                "Trans Nzoia","Nyandarua","Meru","Bungoma","Nyamira","Kakamega",
                                "Laikipia", "Taita Taveta", "Makueni","Elgeyo-Marakwet",
                                "Tharaka-Nithi", "Kericho","Bomet", "Homa Bay", "Busia","Migori",
                                "Nandi","Vihiga"]

ruralcounties_lowserodata = [ "Garissa", "Kitui",
                            "Lamu", "Marsabit","Narok", "Tana River", "Turkana", "Wajir",
                            "Baringo", "Mandera","Samburu", "West Pokot"]

counties_with_informed_priors = [semiurbancounties_lowserodata;ruralcounties_lowserodata]

priors =[fill(KenyaCoVSD.priors_twogroup_newvariant_semi_urbanrural_fitted_detection,length(semiurbancounties_lowserodata));
        fill(KenyaCoVSD.priors_twogroup_newvariant_rural_fitted_detection,length(ruralcounties_lowserodata))]

## parallelised loop for fitting each county
@distributed for i = 1:36
    name = counties_with_informed_priors[i]
    prior = priors[i]
    model = KenyaCoVSD.CoVAreaModel(name,KenyaCoVSD.ll_twogroup_newvariant_infboost,prior;
                case_data = linelist_data_with_pos_neg,
                sero_data = serological_data,
                death_data = deaths_data,
                pop_data = N_kenya)

    KenyaCoVSD.inferparameters!(model,2000,trans,0.05,D;serowaningrate = 0.0)

    @save("modelfits/$(name)_model.jld2",model)
end

ext: jl


text:
export bode

"""
### function bode(sys::StateSpace, w::Array=[]; plot=true)

Compute the bode diagram.

##### Args

* sys: System in state space representation.
* w: Array of frequencies in which the bode diagram will be computed. If empty,
  the frequencies will be chosen according to the system poles.

##### Parameters

* plot: Plot the bode diagram (default = `true`).

##### Returns

Gain (dB), Phase (deg), w (freq. array) for each input/output pair.

"""

function bode(sys::StateSpace, w::Array=[]; plot=true)
    A, B, C, D = sys.A, sys.B, sys.C, sys.D

    # Number of inputs.
    num_u = size(B,2)

    # Number of outputs.
    num_y = size(C,1)

    # Number of states.
    num_states = size(A,1)

    # Identity matrix.
    eye_s = eye(num_states)

    # Check if the vector of frequencies was provided.
    if (isempty(w))
        # If not, compute the bode diagram from a frequency two decades lower
        # than the minimum frequency and two decades higher than the maximum
        # frequency.
        poles = pole(sys)

        # Get the minimum and maximum frequency of the system removing the poles
        # at origin.
        highest_freq  = maximum(abs.(real(poles[poles .!= 0])))
        smallest_freq = minimum(abs.(real(poles[poles .!= 0])))

        w_initial = (smallest_freq != 0) ? smallest_freq/100.0 : 1e-2
        w_final   = (highest_freq != 0)  ?  highest_freq*100.0 : 1e+2

        # Compute 500 logarithmically spaced points.
        w = logspace(log10(w_initial), log10(w_final), 500)
    end

    # Frequency response for each pair input/output.
    gain  = Vector{Vector{Float64}}(length(w))
    phase = Vector{Vector{Float64}}(length(w))

    for i = 1:length(w)
        s = Complex(0,w[i])

        # `fr` is always a vector. Hence, we can squeeze one dimension here.
        fr = squeeze(C*inv(s*eye_s-A)*B+D,1)

        gain[i]  = 20*log10.([ norm(f) for f in fr ])
        phase[i] = [ angle(f)*180.0/pi for f in fr ]

        if i > 1
            # Make sure that the phase plot is continuous. This can be done by
            # analysing the difference between the current phase and the latest
            # one and checking if it is higher than +180.0 or lower than -180.0.
            Δ = phase[i] - phase[i-1]
            phase[i] += (Δ .> 180.0) * (-360.0) + (Δ .< -180.0) * (+360.0)
        end
    end

    ( plot ) && plot_bode(gain, phase, w, num_u, num_y)

    gain, phase, w
end

"""
### function bode(G::TransferFunction, w::Array=[]; plot=true)

Compute the bode diagram.

##### Args

* G: Transfer function of the dynamic system.
* w: Array of frequencies in which the bode diagram will be computed. If empty,
  the frequencies will be chosen according to the system poles.

##### Parameters

* plot: Plot the bode diagram (default = `true`).

##### Returns

Gain (dB), Phase (deg), w (freq. array) for each input/output pair.

"""

function bode(G::TransferFunction, w::Array=[]; plot=true)
    # Convert the transfer function to state space and compute the Bode diagram.
    bode(tf2ss(G), w; plot=plot)
end


text:
function read_isdf( filename )  #=

reads an industrial system described in the Industrial System Description Format (iSDF) from the file whose file name is designated by the argument "filename".

This function requires that the Julia package EzXML has been loaded prior to the execution.

Author: Taku Noda
Started on Oct. 25, 2020
Update on June 29, 2021

Copyright 2021 Central Research Institute of Electric Power Industry

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

=#

# reading a system description from file:
sys = readxml( filename )

# reading the root element:
sys_root = sys.root
path = nodepath( sys_root )
tag = sys_root.name
if tag != "system"
  isdf_error( "The tag of the root iSDF element must be <system>. The tag <$tag>, which is not valid, has been found.", path )
end
if !haskey( sys_root, "name" )
  isdf_error( "The <system> ... </system> section does not have a \"name\" attribute. A <system> ... </system> section must have a \"name\" attribute so as to identify the system by its name.", path )
end
name = sys_root[ "name" ]
if isempty( name )
  isdf_error( "The \"name\" attribute of the <system> ... </system> section is empty. It cannot be empty.", path )
end

# reading the component and subsystem elements:
num_subsyss = 0
for elm in eachelement( sys_root )
  tag = elm.name
  if tag == "component"
    check_component( elm )
  elseif tag == "subsystem"
    num_subsyss = num_subsyss + 1
    if num_subsyss > 1
      isdf_error( "Currently, only one subsystem is allowed in the definition of an industrial system. More than one subsystem is found.", nodepath( elm ) )
    end
    check_component( elm, true )
  end
end

# returning the entire system as an EzXML document:
return sys

end

max_stars_repo_name: takunoda/omegasim

max_issues_repo_name: takunoda/omegasim

max_forks_repo_name: takunoda/omegasim


text:
structeq(a::T, b::T) where {T} = all(f->getfield(a, f) == getfield(b, f), fieldnames(T))

@testset "parsefunction" begin
    @test structeq(MOIU.parsefunction(:x),
                   MOIU.ParsedSingleVariable(:x))
    @test structeq(MOIU.parsefunction(:([x,y,z])),
                   MOIU.ParsedVectorOfVariables([:x,:y,:z]))

    @test structeq(MOIU.parsefunction(:(x + y + 2.0)),
                   MOIU.ParsedScalarAffineFunction(MOIU.ParsedScalarAffineTerm.([1.0,1.0],[:x,:y]), 2.0))
    @test structeq(MOIU.parsefunction(:(x + -3y + 2.0)),
                   MOIU.ParsedScalarAffineFunction(MOIU.ParsedScalarAffineTerm.([1.0,-3.0],[:x,:y]), 2.0))
    @test structeq(MOIU.parsefunction(:(2*x*y + y + 1.0)),
                   MOIU.ParsedScalarQuadraticFunction(MOIU.ParsedScalarAffineTerm.([1.0],[:y]), MOIU.ParsedScalarQuadraticTerm.([2.0],[:x],[:y]), 1.0))

    @test_throws AssertionError MOIU.parsefunction(:(x - y))

    @test structeq(MOIU.parsefunction(:([x, 2x+y+5.0])),
                   MOIU.ParsedVectorAffineFunction(MOIU.ParsedVectorAffineTerm.([1,2,2],MOIU.ParsedScalarAffineTerm.([1.0,2.0,1.0],[:x,:x,:y])), [0.0,5.0]))
    @test structeq(MOIU.parsefunction(:([x, 2x+y+5.0, 1*x*x])),
                   MOIU.ParsedVectorQuadraticFunction(MOIU.ParsedVectorAffineTerm.([1,2,2],MOIU.ParsedScalarAffineTerm.([1.0,2.0,1.0],[:x,:x,:y])), MOIU.ParsedVectorQuadraticTerm.([3],MOIU.ParsedScalarQuadraticTerm.([2.0],[:x],[:x])), [0.0,5.0,0.0]))

end

@testset "separatelabel" begin
    @test MOIU.separatelabel(:(variables: x)) == (:variables, :x)
    @test MOIU.separatelabel(:(variables: x, y)) == (:variables, :((x,y)))
    @test MOIU.separatelabel(:(minobjective: x + y)) == (:minobjective, :(x+y))
    @test MOIU.separatelabel(:(con1: 2x <= 1)) == (:con1, :(2x <= 1))
    @test MOIU.separatelabel(:(con1: [x,y] in S)) == (:con1, :([x,y] in S))
end

@MOIU.model GeneralModel (MOI.ZeroOne, MOI.Integer) (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval) (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.SecondOrderCone, MOI.RotatedSecondOrderCone, MOI.PositiveSemidefiniteConeTriangle) () (MOI.SingleVariable,) (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction) (MOI.VectorOfVariables,) (MOI.VectorAffineFunction,)

@testset "loadfromstring" begin
    @testset "one variable" begin
        s = """
        variables: x
        bound: x >= 1.0
        """
        model = GeneralModel{Float64}()
        x = MOI.add_variable(model)
        MOI.set(model, MOI.VariableName(), x, "x")
        bound = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(1.0))
        MOI.set(model, MOI.ConstraintName(), bound, "bound")

        model2 = GeneralModel{Float64}()
        MOIU.loadfromstring!(model2, s)
        MOIU.test_models_equal(model, model2, ["x"], ["bound"])
    end

    @testset "linear constraints" begin
        s = """
        variables: x, y
        linear1: x + y >= 1.0
        linear2: x + y <= 1.0
        linear3: x + y == 1.0
        """
        model = GeneralModel{Float64}()
        x = MOI.add_variable(model)
        y = MOI.add_variable(model)
        MOI.set(model, MOI.VariableName(), x, "x")
        MOI.set(model, MOI.VariableName(), y, "y")
        linear1 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0), MOI.GreaterThan(1.0))
        linear2 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0), MOI.LessThan(1.0))
        linear3 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0), MOI.EqualTo(1.0))
        MOI.set(model, MOI.ConstraintName(), linear1, "linear1")
        MOI.set(model, MOI.ConstraintName(), linear2, "linear2")
        MOI.set(model, MOI.ConstraintName(), linear3, "linear3")

        model2 = GeneralModel{Float64}()
        MOIU.loadfromstring!(model2, s)
        MOIU.test_models_equal(model, model2, ["x", "y"], ["linear1", "linear2", "linear3"])
    end

    @testset "minimization: linear objective" begin
        s = """
        variables: x, y
        minobjective: x + -2y + 1.0
        """
        model = GeneralModel{Float64}()
        x = MOI.add_variable(model)
        y = MOI.add_variable(model)
        MOI.set(model, MOI.VariableName(), x, "x")
        MOI.set(model, MOI.VariableName(), y, "y")
        MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -2.0], [x, y]), 1.0))
        MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

        model2 = GeneralModel{Float64}()
        MOIU.loadfromstring!(model2, s)
        MOIU.test_models_equal(model, model2, ["x", "y"], String[])
    end

    @testset "maximization: linear objective" begin
        s = """
        variables: x, y
        maxobjective: x + -2y + 1.0
        """
        model = GeneralModel{Float64}()
        x = MOI.add_variable(model)
        y = MOI.add_variable(model)
        MOI.set(model, MOI.VariableName(), x, "x")
        MOI.set(model, MOI.VariableName(), y, "y")
        MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -2.0], [x, y]), 1.0))
        MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

        model2 = GeneralModel{Float64}()
        MOIU.loadfromstring!(model2, s)
        MOIU.test_models_equal(model, model2, ["x", "y"], String[])
    end

    @testset "SOC constraints" begin
        s = """
        variables: x, y, z
        varsoc: [x,y,z] in SecondOrderCone(3)
        affsoc: [2x,y+1,-1*z] in SecondOrderCone(3)
        affsoc2: [1.0,2.0,3.0] in SecondOrderCone(3)
        """
        model = GeneralModel{Float64}()
        x = MOI.add_variable(model)
        y = MOI.add_variable(model)
        z = MOI.add_variable(model)
        MOI.set(model, MOI.VariableName(), x, "x")
        MOI.set(model, MOI.VariableName(), y, "y")
        MOI.set(model, MOI.VariableName(), z, "z")
        varsoc = MOI.add_constraint(model, MOI.VectorOfVariables([x,y,z]), MOI.SecondOrderCone(3))
        affsoc = MOI.add_constraint(model, MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1,2,3], MOI.ScalarAffineTerm.([2.0,1.0,-1.0], [x,y,z])), [0.0,1.0,0.0]), MOI.SecondOrderCone(3))
        affsoc2 = MOI.add_constraint(model, MOI.VectorAffineFunction(MOI.VectorAffineTerm{Float64}[],[1.0,2.0,3.0]), MOI.SecondOrderCone(3))
        MOI.set(model, MOI.ConstraintName(), varsoc, "varsoc")
        MOI.set(model, MOI.ConstraintName(), affsoc, "affsoc")
        MOI.set(model, MOI.ConstraintName(), affsoc2, "affsoc2")

        model2 = GeneralModel{Float64}()
        MOIU.loadfromstring!(model2, s)
        MOIU.test_models_equal(model, model2, ["x", "y", "z"], ["varsoc", "affsoc", "affsoc2"])
    end

end


text:
import Dates
# Conway's Game of Life
using Distributions
row = 10
col = 10
generations = 15
pop_percent = 0.619

function generate_input()
    grid = zeros(row,col)
    for i in 1:row
        for j in 1:col
            val = rand(Uniform(0.00,1.00))
            if val < pop_percent
                grid[i,j]=0
            else
                grid[i,j]=1
            end
        end
    end
    println(grid)
    println("####################################")
    return grid
end

# Display the new generation
function output_test(row,col,future)
    for i in 1:row
        y=""
        for j in 1:col
            if future[i,j]==0.0
                y=string(y, ". ")
            else future[i,j]==1.0
                y=string(y, "@ ")
            end
        end
        println(y)
    end
end

# Simulate Game of Life for given number of Generations
function game_of_life_test(grid)
    future = deepcopy(grid)
    temp = deepcopy(grid)
    r,c = size(future)

    for x in 1:generations
        println("Generation:",x)

        # Loop through each cell
        for i in 1:r
            for j in 1:c
                alive_neighbors=0

                # Finding number of neighbours for current cell
                for a in -1:1
                    for b in -1:1

                        # Variables for counting Neighbors
                        p = 0
                        q = 0

                        # Current number of rows and columns
                        np = r
                        nq = c

                        # Corner Cases including cyclic neighbors
                        # Top Left Corner
                        if i == 1 && j == 1
                            if i+a < 1 && j+b < 1
                                p = np
                                q = nq
                            elseif i+a < 1
                                p = np
                                q = j + b
                            elseif j+b < 1
                                p = i + a
                                q = nq
                            else
                                p = i + a
                                q = j + b
                            end

                        # Top Right Corner
                        elseif i == 1 && j == nq
                            if i+a < 1 && j+b > nq
                                p = np
                                q = 1
                            elseif i+a < 1
                                p = np
                                q = j + b
                            elseif j+b > nq
                                p = i + a
                                q = 1
                            else
                                p = i + a
                                q = j + b
                            end

                        # Bottom Right Corner
                        elseif i == np && j == nq
                            if i+a > np && j+b > nq
                                p = 1
                                q = 1
                            elseif i+a > np
                                p = 1
                                q = j + b
                            elseif j+b > nq
                                p = i + a
                                q = 1
                            else
                                p = i + a
                                q = j + b
                            end

                        # Bottom Left Corner
                        elseif i == np && j == 1
                            if i+a > np && j+b < 1
                                p = 1
                                q = nq
                            elseif i+a > np
                                p = 1
                                q = j + b
                            elseif j+b < 1
                                p = i + a
                                q = nq
                            else
                                p = i + a
                                q = j + b
                            end

                        # Edge Cases
                        # Top Edge
                        elseif i == 1 && i+a < 1
                            p = np
                            q = j + b

                        # Bottom Edge
                        elseif i == np && i+a > np
                            p = 1
                            q = j + b

                        # Left Edge
                        elseif j == 1 && j+b < 1
                            p = i + a
                            q = nq

                        # Right Edge
                        elseif j == nq && j+b > nq
                            p = i + a
                            q = 1

                        # Internal cells
                        else
                            p = i + a
                            q = j + b
                        end

                        alive_neighbors+=temp[p,q]
                    end
                end

                # Current cell to be subtracted to avoid duplicates
                alive_neighbors = alive_neighbors - temp[i,j] - 1

                # Lonely Cell dies
                if temp[i,j]==1 && alive_neighbors<2
                    future[i,j]=1

                # Cell dies due to overcrowding
                elseif temp[i,j]==1 && alive_neighbors>3
                    future[i,j]=0

                # A new cell is born if there are 3 neighbors
                elseif temp[i,j]==0 && alive_neighbors==3
                    future[i,j]=0

                end
            end
        end
        output_test(r,c,future)
    end
    return future
end

function runtest(inp,out,num)
    test_out = game_of_life_test(inp)
    #orig_out = game_of_life(inp)
    println("*************************")
    println(test_out)
    println("*************************")
    #println(orig_out)
    #println("*************************")

    if test_out == out
        open("log.txt","a") do io
            println(io,Dates.now(), "  - Time the Program was Executed - Test Case ",num," | Result - Pass")
	    println(io,"Output=",test_out)
        end
    else
        open("log.txt","a") do io
            println(io,Dates.now(), "  - Time the Program was Executed - Test Case ",num," | Result - Fail")
            println(io,"Output=",test_out)
        end
    end
end

#runtest(generate_input())
#runtest(generate_input())
#runtest(generate_input())

Input1 = [0 0 0 0 0; 0 0 1 0 0; 0 0 1 0 0; 0 0 1 0 0; 0 0 0 0 0]
Input2 = [0 0 0 0 0 0; 0 0 0 1 0 0; 0 1 0 1 0 0; 0 0 1 1 0 0; 0 0 0 0 0 0]
Input3 = [0 0 0 0 0 0; 0 1 1 0 0 0; 0 1 0 0 0 0; 0 0 0 0 1 0; 0 0 0 1 1 0; 0 0 0 0 0 0]

Output1 = [0 0 0 0 0; 0 0 0 0 0; 0 1 1 1 0; 0 0 0 0 0; 0 0 0 0 0]
Output2 = [0 1 0 0 0 1; 1 1 0 0 0 0; 1 0 0 0 0 0; 0 0 0 0 0 0; 0 0 0 0 0 0]
Output3 = [0 0 0 0 0 0; 0 1 1 0 0 0; 0 1 1 0 0 0; 0 0 0 1 1 0; 0 0 0 1 1 0; 0 0 0 0 0 0]

runtest(Input1, Output1,1)
runtest(Input2, Output2,2)
runtest(Input3, Output3,3)


text:
include("DataDictionary.jl")
include("FluxDriver.jl")
include("FVA.jl")
include("Utility.jl")
using GLPK
using DelimitedFiles


text:
# Run tests on functions in src/by_reference_doc/albums.jl

using Spotify.Albums

@testset verbose = true "GET-request endpoints for albums" begin

    # Input argument from composite types defined in src/types.jl
    album_id = SpAlbumId()

    # Cycle through different input keywords for testing
    markets = ["US", "NL", "DE"]
    offsets = [5, 18]

    @testset "For market = $(market_id)" for market_id in markets

        @test ~isempty(album_get_single(album_id, market = market_id )[1])
        @test ~isempty(album_get_tracks(album_id, market = market_id)[1])
        @test ~isempty(album_get_multiple("$(album_id),5XgEM5g3xWEwL4Zr6UjoLo,2rpT0freJsmUmmPluVWqg5",
                                          market = market_id))

    end

    @testset "For offset = $(num_offset)" for num_offset in offsets

        @test ~isempty(album_get_tracks(album_id, offset = num_offset)[1])

        @test ~isempty(album_get_saved(offset = num_offset)[1])

    end

    @test ~isempty(album_get_contains("$(album_id),5XgEM5g3xWEwL4Zr6UjoLo"))

end


text:
function factoredFastMarching(pEik::EikonalParam, mem::EikonalTempMemory)
# time1 = 0.0;
# time2 = 0.0;
# time3 = 0.0;

kappaSquared 	= pEik.kappaSquared;
h 				= pEik.Mesh.h;
src 			= pEik.src;
n 				= pEik.Mesh.n+1;
HO 				= pEik.HO;

if isempty(pEik.T1)
	pEik.T1 = zeros(Float64,tuple(n...));
end
T 		 = pEik.T1;
if isempty(pEik.ordering)
	pEik.ordering = zeros(EikIdxType,prod(n));
end
ordering = pEik.ordering;
if isempty(pEik.OP)
	pEik.OP = zeros(Int8,tuple(n...));
end
OP       = pEik.OP;
OP[:] 	 = 0;

T[:] = Inf;
Done = mem.Done;
Done[:] = false;

heapArr = mem.V;
heapIdx = mem.J;
frontHeap = initHeap(heapArr,heapIdx);

monoThresh = max(h[1],h[2])^2;

# Done = zeros(Bool,n[1]+4,n[2]+4); # the +4 is actually only for HO. For FO 2 is enough but need to change access to Done in helpFuncs

src_cs = loc2cs(src,n);

kappaSrc = sqrt(kappaSquared[src_cs]);
T[src_cs] = kappaSrc;
insertToHeap(frontHeap,0.0,src_cs);

k_order = 1;
Neighbors2D = zeros(Int64,4,2);
Boundary2D = zeros(Bool,4);
curr_loc = zeros(Int64,2);
GradLoc = zeros(Float64,2);
A = zeros(Float64,2);
B = zeros(Float64,2);
OpLoc = zeros(Int8,2);
curr_val = zeros(Float64,1);

hinv = 1./h;

while frontHeap.size > 0
	# tic();
	curr = getMin(frontHeap,curr_val);
	cs2loc(curr_loc,curr,n);
	while getPaddedArr2D(Done,curr_loc[1],curr_loc[2]) && frontHeap.size>0
		curr = getMin(frontHeap,curr_val);
		cs2loc(curr_loc,curr,n);
	end
	# time1 += toq();
	if getPaddedArr2D(Done,curr_loc[1],curr_loc[2]) && frontHeap.size==0
        break;
    end
	ordering[k_order] = curr;
    k_order += 1;
	setPaddedArr2D(Done,curr_loc[1],curr_loc[2],true);
	get2DNeighbors(curr_loc,n,Neighbors2D,Boundary2D);
	for k = 1:4
		if Boundary2D[k]
			@inbounds nei_loc1 = Neighbors2D[k,1];
			@inbounds nei_loc2 = Neighbors2D[k,2];

			if !getPaddedArr2D(Done,nei_loc1,nei_loc2)
				# tic();
				T0loc = analyticLocal(h,src,nei_loc1,nei_loc2,GradLoc);

				T_k = getUpdatedFactoredVal_2D(HO,kappaSquared,T,Done,nei_loc1,nei_loc2,h,hinv,src,T0loc,GradLoc,OpLoc,A,B);

				# @inbounds if T_k*T0loc < curr_val[1]
					# @inbounds if abs(T_k*T0loc-curr_val[1]) > monoThresh
						 # println("Monotonicity not fulfilled by ",(T_k*T0loc-curr_val)," - not corrected. HO = ",HO);
					# end
				# end
				# time2 += toq();
				# tic();

				if T_k < T[nei_loc1,nei_loc2]
					@inbounds T[nei_loc1,nei_loc2] = T_k;
					@inbounds OP[nei_loc1,nei_loc2] = getOpCode(OpLoc[1],OpLoc[2],0) ;
					@inbounds insertToHeap(frontHeap,T0loc*T_k,loc2cs(nei_loc1,nei_loc2,n));
				end
				# time3 += toq();
			end

		end
	end
end

# println("time heap: ",time1, " time quad: ",time2, " time heapInsert: ",time3);
return T,ordering;
end;


text:
# using Flux: conv
using DSP
using LinearAlgebra

prod0__(x, y) = cat(
    [x .* y for y in eachslice(y, dims = length(size(y)))]...,
    dims = length(size(x)),
)
prod_0_(x, y) = cat(
    [x .* y for x in eachslice(x, dims = length(size(x)))]...,
    dims = length(size(x)),
)
prod__0(x, y) = cat(
    sum([
        x .* y
        for
        (x, y) in zip(
            eachslice(x, dims = length(size(x))),
            eachslice(y, dims = length(size(y))),
        )
    ]),
    dims = length(size(x)),
)
# prod111(x, y) = join([
#     x[:, :, :, 2] .* y[:, :, :, 3] - x[:, :, :, 3] .* y[:, :, :, 2],
#     x[:, :, :, 3] .* y[:, :, :, 1] - x[:, :, :, 1] .* y[:, :, :, 3],
#     x[:, :, :, 1] .* y[:, :, :, 2] - x[:, :, :, 2] .* y[:, :, :, 1],
# ])

function myconv(x, y)
    DSP.conv(x, y)
    # Flux.conv(cat(x,dims=5),cat(y,dims=5);pad=(size(y).-1).÷2)[:,:,:,1,1]
end

function conv0__(x, y)
    dims = length(size(y))
    cat([myconv(dropdims(x; dims), y) for y in eachslice(y; dims)]...; dims)
end

function conv_0_(x, y)
    dims = length(size(x))
    cat([myconv(x, dropdims(y; dims)) for x in eachslice(x; dims)]; dims)
end

conv__0(x, y) = cat(
    sum([
        myconv(x, y)
        for
        (x, y) in zip(
            eachslice(x, dims = length(size(x))),
            eachslice(y, dims = length(size(y))),
        )
    ]),
    dims = length(size(x)),
)
# conv111 =
#     (x, y) -> join([
#         myconv(x[:, :, :, 2], y[:, :, :, 3]) -
#         myconv(x[:, :, :, 3], y[:, :, :, 2]),
#         myconv(x[:, :, :, 3], y[:, :, :, 1]) -
#         myconv(x[:, :, :, 1], y[:, :, :, 3]),
#         myconv(x[:, :, :, 1], y[:, :, :, 2]) -
#         myconv(x[:, :, :, 2], y[:, :, :, 1]),
#     ])
# conv121 =
#     (x, y) -> join([
#         myconv(x[:, :, :, 1], -y[:, :, :, 3] / sqrt(3) .- y[:, :, :, 5]) +
#         myconv(x[:, :, :, 2], y[:, :, :, 2]) +
#         myconv(x[:, :, :, 3], y[:, :, :, 1]),
#         myconv(x[:, :, :, 1], y[:, :, :, 2]) +
#         myconv(x[:, :, :, 2], 2 / sqrt(3) * y[:, :, :, 3]) +
#         myconv(x[:, :, :, 3], y[:, :, :, 4]),
#         myconv(x[:, :, :, 1], y[:, :, :, 1]) +
#         myconv(x[:, :, :, 2], y[:, :, :, 4]) +
#         myconv(x[:, :, :, 3], -y[:, :, :, 3] / sqrt(3) .+ y[:, :, :, 5]),
#     ])

function FieldProd(l1, l2, l)
    # l1, l2, l = ranks
    if l1 == 0
        return prod0__
    elseif l2 == 0
        return prod_0_
    elseif l1 == l2 && l == 0
        return prod__0
    elseif l1 == l2 == l == 1
        return prod111
    end

    # (x, y) -> yield(f(x, y), x.grid)
    # (x,y)->yield
end

function FieldConv(l1, l2, l)
    if l1 == 0
        f = conv0__
    elseif l2 == 0
        f = conv_0_
    elseif l1 == l2 && l == 0
        f = conv__0
    elseif l1 == l2 == l == 1
        f = conv111
    elseif (l1, l2, l) == (1, 2, 1)
        f = conv121
    end
    f
    # (x, y) -> tfconv(f,x,y)
end

"""
    ⊗(x,y)
    fieldprod(x::AbstractArray, y::AbstractArray)
    fieldprod(x::AbstractArray, y::AbstractArray,l)

    Alias of fieldprod(u, v), `u ⊗ v` computes the appropriate pointwise product between 2 scalar or vector fields (inferred as scalar-scalar, scalar-vector, dot, cross). For greater clarity one may also write `u ⋅ v` for dot and `u ⨉ v` for cross
"""
function fieldprod(x::AbstractArray, y::AbstractArray, l)
    lx, ly = getl.([x, y])
    FieldProd(lx, ly, l)(x, y)
end
function fieldprod(x::AbstractArray, y::AbstractArray)
    lx, ly = getl.([x, y])
    l = abs(lx - ly)
    FieldProd(lx, ly, l)(x, y)
end
function getl(x)
    n = size(x)[end]
    if n == 1
        return 0
    end

    dims = ndims(x) - 1
    if dims == 2
        return 1
    elseif dims == 3
        return (n - 1) ÷ 2
    end

end
function ⊗(x,y)
    fieldprod(x,y)
end
#
# function LinearAlgebra.:⋅(x::Array{T,3}, y::Array{T,3}) where T<:AbstractFloat
#     fieldprod(x,y,0)
# end
# function LinearAlgebra.:⋅(x::Array{T,4}, y::Array{T,4}) where T<:AbstractFloat
#     fieldprod(x,y,0)
# end
# function ⨉(x::Array{T,4}, y::Array{T,4}) where T<:AbstractFloat
#     fieldprod(x,y,1)
# end
#
# function LinearAlgebra.adjoint(x::AbstractArray)
#     x
# end
"""
"""
function field_conv(x::AbstractArray, y::AbstractArray; rank = nothing)
    lx, ly = l.([x, y])
    FieldConv(lx, ly, rank === nothing ? lx + ly : rank)(x, y)
end

"""
    function field_norm(field::AbstractArray)
"""
function fieldnorm(x::AbstractArray)
    # n=length(size(x))-1
    if size(x)[end] == 1
        return abs.(x)
    end
    # r=sum([x[:,:,:,i:i] .^ 2 for i in 1:size(x,4)])
    # sqrt.(r.+.01).-.1
    # sqrt.(sum([x .^ 2 for x in eachslice(x, dims = length(size(x)))]))
    cat(
        sqrt.(sum([x .^ 2 for x in eachslice(x, dims = ndims(x))])),
        dims = ndims(x),
    )
end

function Δ(x, y)
    sum(abs.(x .- y)) / sum(abs.(y))
end

function nae(yhat, y; sumy = sum(abs.(y)))
    if sumy == 0
        error()
    end
    sum(abs.(yhat .- y)) / sumy
end


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

using BinaryBuilder, Dates, LibGit2

name = "xdmf"
repo = LibGit2.clone("https://gitlab.kitware.com/xdmf/xdmf.git", tempname())
hash_long = string(LibGit2.revparseid(repo, "master"))
build = Dates.format(Dates.now(), "yyyymmdd")
version = VersionNumber("3+$build")

# Collection of sources required to build xdmf
sources = [
    "https://gitlab.kitware.com/xdmf/xdmf.git" =>
    hash_long,

    "https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.10/hdf5-1.10.3/src/hdf5-1.10.3.tar.gz" =>
    "b600d7c914cfa80ae127cd1a1539981213fee9994ac22ebec9e3845e951d9b39",

    "https://dl.bintray.com/boostorg/release/1.71.0/source/boost_1_71_0.tar.gz" =>
    "96b34f7468f26a141f6020efb813f1a2f3dfb9797ecf76a7d7cbd843cc95f5bd",
]

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

cd hdf5-1.10.3/
./configure --prefix=$prefix --host=$target
make
make install
cd ..

export BOOST_ROOT=$WORKSPACE/srcdir/boost_1_71_0
cd xdmf
mkdir build
cd build
cmake .. -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DBUILD_SHARED_LIBS=1
make
make install
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libXdmf", :libXdmf),
    LibraryProduct(prefix, "libXdmfCore", :libXdmfCore)
]

# Dependencies that must be installed before this package can be built
dependencies = [
	"https://raw.githubusercontent.com/JuliaIO/LightXML.jl/master/deps/build_XML2Builder.v2.9.9.jl"
]

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


text:
using Base.Test
using Elemental

for T in (Float32, Float64, Complex{Float32}, Complex{Float64})
    A = Elemental.Matrix(T)
    Elemental.gaussian!(A, 10)
    @test isa(similar(A), Elemental.Matrix{T})
    @test isa(similar(A, Elemental.ElInt), Elemental.Matrix{Elemental.ElInt})
    @test isa(similar(A, Float32), Elemental.Matrix{Float32})
    @test isa(similar(A, Elemental.ElInt, (5,5)), Elemental.Matrix{Elemental.ElInt})
    @test isa(similar(A, Float32, (5,5)), Elemental.Matrix{Float32})
end


text:
type EvalError
  err
  bt
end

EvalError(err) = EvalError(err, [])

# @render Inline e::EvalError Text(sprint(Base.showerror, e.err, e.bt))

rendererr(err) = strong(".error-description", err)

function splitlink(path)
  m = match(r"^(.?:?[^:]*)(?::(\d+))$", path)
  m == nothing && return path, 0
  m.captures[1], parse(Int, m.captures[2])
end

function btlines(bt)
  isempty(bt) && return []
  length(bt) > 100 && (bt = bt[1:100])
  trace = Base.stacktrace(bt)
  t_trace = map(trace) do frame
    (split(sprint(Base.StackTraces.show_spec_linfo, frame), '(')[1],
    string(frame.file),
    frame.line,
    frame.inlined ? " [inlined]" : "")
  end
  # clip stack traces
  internal_ind = findlast(t_trace) do t_frame
    ismatch(r"loading\.jl", t_frame[2]) &&
    ismatch(r"include_string", t_frame[1])
  end
  internal_ind == 0 ? t_trace : t_trace[1:internal_ind-1]
end

highlights(lines) =
  @>> lines filter(x->x[3]>0) map(x->d(:file=>fullpath(x[2]),:line=>x[3]))

highlights(e::EvalError) = highlights(btlines(e.bt))

function renderbt(ls)
  span(".error-trace",
       [div(".trace-entry",
        c(fade("in "), f, fade(" at "),
          render(Inline(), Copyable(baselink(loc, li))), fade(inlined)))
        for (f, loc, li, inlined) in ls])
end

function render(::Editor, e::EvalError)
  header = sprint(showerror, e.err)
  header = split(header, '\n')
  trace = btlines(e.bt)
  view = if isempty(trace)
    if length(header) == 1
      rendererr(header[1])
    else
      Tree(rendererr(header[1]), [rendererr(join(header[2:end], '\n'))])
    end
  else
    Tree(rendererr(header[1]), [rendererr(join(header[2:end], '\n')); renderbt(trace)])
  end

  d(:type => :error,
    :view => render(Inline(),
                    Copyable(view, sprint(showerror, e.err, e.bt))),
    :highlights => highlights(e))
end

render(::Console, e::EvalError) =
  @msg result(render(Editor(), e))


text:
#AbstractMatrix

using CairoMakie

data = cumsum(randn(4, 101), dims = 2)

fig, ax, sp = series(data, labels=["label $i" for i in 1:4])
axislegend(ax)

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
# rotatez.jl

export imrotate!, imrotate_adj!
export imrotate, imrotate_adj

using LinearInterpolators: LinearSpline
using LinearInterpolators: SparseInterpolator, AffineTransform2D, rotate
using LinearInterpolators: TwoDimensionalTransformInterpolator

"""
    linearinterp!(A, s, e, len)
Assign key values in `SparseInterpolator` (linear) `A` that are calculated from `s`, `e` and `len`.
`s` means start (x[1])
`e` means end (x[end])
`len` means length (length(x))
"""
function linearinterp!(
    A::SparseInterpolator{<:AbstractFloat},
    s::RealU,
    e::RealU,
    len::Int,
)

    dec = ceil(Int, s) - s
    ncoeff = length(A.C)
    ncol = len
    for i = 1:ncoeff
        if isodd(i)
            A.C[i] = dec
        else
            A.C[i] = 1 - dec
        end
    end

    if e <= ncol
        A.J[end] = ceil(Int, e)
        for i = ncoeff-1:-1:1
            A.J[i] = max(1, A.J[end] - ceil(Int, (ncoeff - i) / 2))
        end
    else
        A.J[1] = floor(Int, s)
        for i = 2:ncoeff
            A.J[i] = min(ncol, A.J[1] + ceil(Int, (i - 1) / 2))
        end
    end
    return A
end

"""
    rotate_x!(output, img, tan_θ, interp)
"""
function rotate_x!(
    output::AbstractMatrix{<:RealU},
    img::AbstractMatrix{<:RealU},
    tan_θ::RealU,
    interp::SparseInterpolator,
)

    len = size(img, 1) # size of a square image
    c_y = (len + 1)/2 # center of yi
    idx = 1:len

    for i in idx
        s = (idx[i] - c_y) * tan_θ + 1
        e = (idx[i] - c_y) * tan_θ + len
        linearinterp!(interp, s, e, len)
        mul!((@view output[:, i]), interp, (@view img[:, i])) # need mul! to avoid allocating
    end
    return output
end

"""
    rotate_x_adj!(output, img, tan_θ, interp)
"""
function rotate_x_adj!(
    output::AbstractMatrix{<:RealU},
    img::AbstractMatrix{<:RealU},
    tan_θ::RealU,
    interp::SparseInterpolator,
)

    len = size(img, 1)
    c_y = (len + 1)/2 # center of yi
    idx = 1:len

    for i in idx
        s = (idx[i] - c_y) * tan_θ + 1
        e = (idx[i] - c_y) * tan_θ + len
        linearinterp!(interp, s, e, len)
        mul!((@view output[:, i]), interp', (@view img[:, i])) # need mul! to avoid allocating
    end
    return output
end

"""
    rotate_y!(output, img, sin_θ, interp)
"""
function rotate_y!(
    output::AbstractMatrix{<:RealU},
    img::AbstractMatrix{<:RealU},
    sin_θ::RealU,
    interp::SparseInterpolator,
)

    len = size(img, 1)
    c_x = (len + 1)/2 # center of xi
    idx = 1:len

    for i in idx
        s = (idx[i] - c_x) * sin_θ + 1
        e = (idx[i] - c_x) * sin_θ + len
        linearinterp!(interp, s, e, len)
        mul!((@view output[i, :]), interp, (@view img[i, :]))
    end
    return output
end

"""
    rotate_y_adj!(output, img, sin_θ, interp)
"""
function rotate_y_adj!(
    output::AbstractMatrix{<:RealU},
    img::AbstractMatrix{<:RealU},
    sin_θ::RealU,
    interp::SparseInterpolator,
)

    len = size(img, 1)
    c_x = (len + 1)/2 # center of xi
    idx = 1:len

    for i in idx
        s = (idx[i] - c_x) * sin_θ + 1
        e = (idx[i] - c_x) * sin_θ + len
        linearinterp!(interp, s, e, len)
        mul!((@view output[i, :]), interp', (@view img[i, :]))
    end
    return output
end

"""
    rotl90!(B::AbstractMatrix, A::AbstractMatrix)
In place version of `rotl90`, returning rotation of `A` in `B`.
"""
function rotl90!(B::AbstractMatrix, A::AbstractMatrix)
    ind1, ind2 = axes(A)
    n = first(ind2) + last(ind2)
    for i = axes(A, 1), j = ind2
        B[n - j, i] = A[i, j]
    end
    return B
end

"""
    rotr90!(B::AbstractMatrix, A::AbstractMatrix)
In place version of `rotr90`, returning rotation of `A` in `B`.
"""
function rotr90!(B::AbstractMatrix, A::AbstractMatrix)
    ind1, ind2 = axes(A)
    m = first(ind1) + last(ind1)
    for i = ind1, j = axes(A, 2)
        B[j, m - i] = A[i, j]
    end
    return B
end

"""
    rot180!(B::AbstractMatrix, A::AbstractMatrix)
In place version of `rot180`, returning rotation of `A` in `B`.
"""
function rot180!(B::AbstractMatrix, A::AbstractMatrix)
    ind1, ind2 = axes(A,1), axes(A,2)
    m, n = first(ind1)+last(ind1), first(ind2)+last(ind2)
    for j=ind2, i=ind1
        B[m-i,n-j] = A[i, j]
    end
    return B
end

"""
    rot_f90!(output, img, m)
In-place version of rotating an image by 90/180/270 degrees
"""
function rot_f90!(output::AbstractMatrix, img::AbstractMatrix, m::Int)
    if m == 0
        output .= img
    elseif m == 1
        rotl90!(output, img)
    elseif m == 2
        rot180!(output, img)
    elseif m == 3
        rotr90!(output, img)
    else
        throw("invalid m!")
    end
    return output
end

"""
    rot_f90_adj!(output, img, m)
    The adjoint of rotating an image by 90/180/270 degrees
"""
function rot_f90_adj!(output::AbstractMatrix, img::AbstractMatrix, m::Int)
    if m == 0
        output .= img
    elseif m == 1
        rotr90!(output, img)
    elseif m == 2
        rot180!(output, img)
    elseif m == 3
        rotl90!(output, img)
    else
        throw("invalid m!")
    end
    return output
end

"""
    imrotate!(output, img, θ, plan)
Rotate a 2D image `img` by angle `θ ∈ [0,2π]`
in counter-clockwise direction (opposite to `imrotate` in Julia)
using 3-pass 1d linear interpolations.
"""
function imrotate!(
    output::AbstractMatrix{<:RealU},
    img::AbstractMatrix{<:RealU},
    θ::RealU,
    plan::PlanRotate{<:Number,RotateMode{:one}},
)
    @boundscheck size(output) == size(img) || throw("size")
    @boundscheck size(img, 1) == size(img, 2) || throw("row != col")

    if mod(θ, 2π) ≈ 0
        output .= img
        return output # todo: check!
    end

    m = mod(floor(Int, 0.5 + θ/(π/2)), 4)

    N = size(img, 1) # M = N!

    if θ ≈ m * (π/2)
        padzero!(plan.workmat2, img, (plan.padsize, plan.padsize, plan.padsize, plan.padsize))
        rot_f90!(plan.workmat1, plan.workmat2, m)
    else
        mod_theta = θ - m * (π/2) # make sure it is between -45 and 45 degree
        tan_mod_theta = tan(mod_theta / 2)
        sin_mod_theta = - sin(mod_theta)
        padzero!(plan.workmat1, img, (plan.padsize, plan.padsize, plan.padsize, plan.padsize))
        rot_f90!(plan.workmat2, plan.workmat1, m)
        rotate_x!(plan.workmat1, plan.workmat2, tan_mod_theta, plan.interp)
        rotate_y!(plan.workmat2, plan.workmat1, sin_mod_theta, plan.interp)
        rotate_x!(plan.workmat1, plan.workmat2, tan_mod_theta, plan.interp)
    end

    output .= (@view plan.workmat1[plan.padsize + 1 : plan.padsize + N,
                                   plan.padsize + 1 : plan.padsize + N])

    return output
end

"""
    imrotate_adj!(output, img, θ, plan)
Adjoint of `imrotate!` for a 2D image
using 3-pass 1d linear interpolations.
"""
function imrotate_adj!(
    output::AbstractMatrix{<:RealU},
    img::AbstractMatrix{<:RealU},
    θ::RealU,
    plan::PlanRotate{<:Number,RotateMode{:one}},
)

    @boundscheck size(output) == size(img) || throw("size")
    @boundscheck size(img, 1) == size(img, 2) || throw("row != col")

    if mod(θ, 2π) ≈ 0
        output .= img
        return output
    end

    m = mod(floor(Int, 0.5 + θ/(π/2)), 4)
    N = size(img, 1) # M = N!

    if θ ≈ m * (π/2)
        padzero!(plan.workmat2, img, (plan.padsize, plan.padsize, plan.padsize, plan.padsize))
        rot_f90_adj!(plan.workmat1, plan.workmat2, m)
    else
        mod_theta = θ - m * (π/2) # make sure it is between -45 and 45 degree
        tan_mod_theta = tan(mod_theta / 2)
        sin_mod_theta = - sin(mod_theta)
        padzero!(plan.workmat1, img, (plan.padsize, plan.padsize, plan.padsize, plan.padsize))
        rotate_x_adj!(plan.workmat2, plan.workmat1, tan_mod_theta, plan.interp)
        rotate_y_adj!(plan.workmat1, plan.workmat2, sin_mod_theta, plan.interp)
        rotate_x_adj!(plan.workmat2, plan.workmat1, tan_mod_theta, plan.interp)
        rot_f90_adj!(plan.workmat1, plan.workmat2, m) # must be two different arguments
    end
    output .= (@view plan.workmat1[plan.padsize + 1 : plan.padsize + N,
                                   plan.padsize + 1 : plan.padsize + N])
    return output
end

"""
    imrotate!(output, img, θ, plan)
Rotate a 2D image `img` by angle `θ ∈ [0,2π]`
in counter-clockwise direction (opposite to `imrotate` in Julia)
using 2d bilinear interpolation.
"""
function imrotate!(
    output::AbstractMatrix{<:RealU},
    img::AbstractMatrix{<:RealU},
    θ::RealU,
    plan::PlanRotate{T,RotateMode{:two}},
) where {T}

    @boundscheck size(output) == size(img) || throw("size")
    @boundscheck size(img, 1) == size(img, 2) || throw("row != col")

    if mod(θ, 2π) ≈ 0
        output .= img
        return output
    end

    N = size(img, 1) # M = N!

    ker = LinearSpline(T)
    rows = size(plan.workmat2)
    cols = size(plan.workmat1)
    c = ((1 + rows[1]) / 2, (1 + rows[2]) / 2)
    R = c + rotate(2π - θ, AffineTransform2D{T}() - c)
    A = TwoDimensionalTransformInterpolator(rows, cols, ker, R)

    padzero!(plan.workmat1, img, (plan.padsize, plan.padsize, plan.padsize, plan.padsize))
    mul!(plan.workmat2, A, plan.workmat1)

    output .= (@view plan.workmat2[plan.padsize + 1 : plan.padsize + N,
                                   plan.padsize + 1 : plan.padsize + N])
    return output
end

"""
    imrotate(img, θ; method::Symbol=:two)
Rotate a 2D image `img` by angle `θ ∈ [0,2π]`
in counter-clockwise direction (opposite to `imrotate` in Julia)
using either 2d linear interpolation (for `:two`)
or 3-pass 1D interpolation (for `:one`)
"""
function imrotate(
    img::AbstractMatrix{I},
    θ::RealU;
    method::Symbol=:two,
    T::DataType = promote_type(I, Float32),
) where {I <: Number}
    output = similar(Matrix{T}, size(img))
    plan = plan_rotate(size(img, 1); T, nthread = 1, method)[1]
    imrotate!(output, img, θ, plan)
    return output
end

"""
    imrotate_adj!(output, img, θ, plan)
Adjoint of `imrotate!` for a 2D image
using 2D bilinear interpolation.
"""
function imrotate_adj!(
    output::AbstractMatrix{<:RealU},
    img::AbstractMatrix{<:RealU},
    θ::RealU,
    plan::PlanRotate{T,RotateMode{:two}},
) where {T}

    @boundscheck size(output) == size(img) || throw("size")
    @boundscheck size(img, 1) == size(img, 2) || throw("row != col")

    N = size(img, 1) # M = N!

    ker = LinearSpline(T)
    rows = size(plan.workmat2)
    cols = size(plan.workmat1)
    c = ((1 + rows[1]) / 2, (1 + rows[2]) / 2)
    R = c + rotate(2π - θ, AffineTransform2D{T}() - c)
    A = TwoDimensionalTransformInterpolator(rows, cols, ker, R)

    padzero!(plan.workmat1, img, (plan.padsize, plan.padsize, plan.padsize, plan.padsize))
    mul!(plan.workmat2, A', plan.workmat1) # todo: internals of A' ?

    output .= (@view plan.workmat2[plan.padsize + 1 : plan.padsize + N,
                                   plan.padsize + 1 : plan.padsize + N])
    return output
end

"""
    imrotate_adj(img, θ; method::Symbol=:two)
Adjoint of rotating a 2D image `img` by angle `θ ∈ [0,2π]`
in counter-clockwise direction (opposite to `imrotate` in Julia)
using either 2d linear interpolations or 3-pass 1D interpolation.
"""
function imrotate_adj(
    img::AbstractMatrix{I},
    θ::RealU;
    method::Symbol=:two,
    T::DataType = promote_type(I, Float32),
) where {I <: Number}
    output = similar(Matrix{T}, size(img))
    plan = plan_rotate(size(img, 1); T, nthread = 1, method)[1]
    imrotate_adj!(output, img, θ, plan)
    return output
end

function _task3(z, fun::Function, output, image3, θ, plans)
    id = Threads.threadid()
#   1 ≤ id ≤ length(plans) || throw("bug: id=$id nplan=$(length(plans))")
    return fun(
        (@view output[:,:,z]),
        (@view image3[:,:,z]),
        θ,
        plans[id],
    )
end

# rotate image3 using foreach
function _imrotate!(
    fun::Function, # imrotate! or imrotate_adj!
    output::AbstractArray{<:RealU,3},
    image3::AbstractArray{<:RealU,3},
    θ::RealU,
    plans::Vector{<:PlanRotate}, # must be ≥ Threads.nthreads()
    ntasks::Int = Threads.nthreads(), # ≥ 1
)

    size(output) == size(image3) || throw(DimensionMismatch())
    size(image3,1) == plans[1].nx || throw("nx")
    length(plans) == Threads.nthreads() || throw("#threads")

    task = z -> _task3(z, fun, output, image3, θ, plans)
    Threads.foreach(task, foreach_setup(1:size(image3,3)); ntasks)

    return output
end

"""
    imrotate!(output, image3, θ, plans, ntasks=nthreads)
In-place version of rotating a 3D `image3` by `θ`
in counter-clockwise direction (opposite to `imrotate` in Julia)
using `foreach` with `ntasks`.
"""
function imrotate!(
    output::AbstractArray{<:RealU,3},
    image3::AbstractArray{<:RealU,3},
    θ::RealU,
    plans::Vector{<:PlanRotate},
    ntasks::Int = Threads.nthreads(),
)
    return _imrotate!(imrotate!, output, image3, θ, plans, ntasks)
end

"""
    imrotate_adj!(output, image3, θ, plans, ntasks=nthreads)
Adjoint of `imrotate!` using `foreach`.
"""
function imrotate_adj!(
    output::AbstractArray{<:RealU,3},
    image3::AbstractArray{<:RealU,3},
    θ::RealU,
    plans::Vector{<:PlanRotate},
    ntasks::Int = Threads.nthreads(),
)
    return _imrotate!(imrotate_adj!, output, image3, θ, plans, ntasks)
end

# rotate image3 using `Threads.@threads`
function _imrotate!(
    fun::Function, # imrotate! or imrotate_adj!
    output::AbstractArray{<:RealU,3},
    image3::AbstractArray{<:RealU,3},
    θ::RealU,
    plans::Vector{<:PlanRotate}, # must be ≥ Threads.nthreads()
    ::Symbol, # :thread
)

    size(output) == size(image3) || throw(DimensionMismatch())
    size(image3,1) == plans[1].nx || throw("nx")
    length(plans) == Threads.nthreads() || throw("#threads")

    Threads.@threads for z = 1:size(image3,3) # 1:nz
        id = Threads.threadid() # thread id
        fun(
            (@view output[:, :, z]),
            (@view image3[:, :, z]),
            θ,
            plans[id],
        )
    end

    return output
end

"""
    imrotate!(output, image3, θ, plans, :thread)
Mutating version of rotating a 3D `image3` by `θ`
in counter-clockwise direction (opposite of `imrotate` in Julia)
writing result into `output`,
using `Threads.@threads`.
"""
function imrotate!(
    output::AbstractArray{<:RealU,3},
    image3::AbstractArray{<:RealU,3},
    θ::RealU,
    plans::Vector{<:PlanRotate},
    ::Symbol, # :thread
)
    return _imrotate!(imrotate!, output, image3, θ, plans, :thread)
end

"""
    imrotate_adj!(output, image3, θ, plans, :thread)
Adjoint of `imrotate!` using `@threads`.
"""
function imrotate_adj!(
    output::AbstractArray{<:RealU,3},
    image3::AbstractArray{<:RealU,3},
    θ::RealU,
    plans::Vector{<:PlanRotate},
    ::Symbol, # :thread
)
    return _imrotate!(imrotate_adj!, output, image3, θ, plans, :thread)
end


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

export
    av_opencl_get_device_list,
    av_opencl_free_device_list,
    av_opencl_set_option,
    av_opencl_get_option,
    av_opencl_free_option,
    av_opencl_alloc_external_env,
    av_opencl_free_external_env,
    av_opencl_errstr,
    av_opencl_register_kernel_code,
    av_opencl_init,
    av_opencl_compile,
    av_opencl_get_command_queue,
    av_opencl_buffer_create,
    av_opencl_buffer_write,
    av_opencl_buffer_read,
    av_opencl_buffer_write_image,
    av_opencl_buffer_read_image,
    av_opencl_buffer_release,
    av_opencl_uninit,
    av_opencl_benchmark

function av_opencl_get_device_list(device_list)
    ccall((:av_opencl_get_device_list,libavutil),Cint,(Ptr{Ptr{AVOpenCLDeviceList}},),device_list)
end

function av_opencl_free_device_list(device_list)
    ccall((:av_opencl_free_device_list,libavutil),Void,(Ptr{Ptr{AVOpenCLDeviceList}},),device_list)
end

function av_opencl_set_option(key,val)
    ccall((:av_opencl_set_option,libavutil),Cint,(Cstring,Cstring),key,val)
end

function av_opencl_get_option(key,out_val)
    ccall((:av_opencl_get_option,libavutil),Cint,(Cstring,Ptr{Ptr{UInt8}}),key,out_val)
end

function av_opencl_free_option()
    ccall((:av_opencl_free_option,libavutil),Void,())
end

function av_opencl_alloc_external_env()
    ccall((:av_opencl_alloc_external_env,libavutil),Ptr{AVOpenCLExternalEnv},())
end

function av_opencl_free_external_env(ext_opencl_env)
    ccall((:av_opencl_free_external_env,libavutil),Void,(Ptr{Ptr{AVOpenCLExternalEnv}},),ext_opencl_env)
end

function av_opencl_errstr(status::cl_int)
    ccall((:av_opencl_errstr,libavutil),Cstring,(cl_int,),status)
end

function av_opencl_register_kernel_code(kernel_code)
    ccall((:av_opencl_register_kernel_code,libavutil),Cint,(Cstring,),kernel_code)
end

function av_opencl_init(ext_opencl_env)
    ccall((:av_opencl_init,libavutil),Cint,(Ptr{AVOpenCLExternalEnv},),ext_opencl_env)
end

function av_opencl_compile(program_name,build_opts)
    ccall((:av_opencl_compile,libavutil),cl_program,(Cstring,Cstring),program_name,build_opts)
end

function av_opencl_get_command_queue()
    ccall((:av_opencl_get_command_queue,libavutil),cl_command_queue,())
end

function av_opencl_buffer_create(cl_buf,cl_buf_size::Csize_t,flags::Integer,host_ptr)
    ccall((:av_opencl_buffer_create,libavutil),Cint,(Ptr{cl_mem},Csize_t,Cint,Ptr{Void}),cl_buf,cl_buf_size,flags,host_ptr)
end

function av_opencl_buffer_write(dst_cl_buf::cl_mem,src_buf,buf_size::Csize_t)
    ccall((:av_opencl_buffer_write,libavutil),Cint,(cl_mem,Ptr{UInt8},Csize_t),dst_cl_buf,src_buf,buf_size)
end

function av_opencl_buffer_read(dst_buf,src_cl_buf::cl_mem,buf_size::Csize_t)
    ccall((:av_opencl_buffer_read,libavutil),Cint,(Ptr{UInt8},cl_mem,Csize_t),dst_buf,src_cl_buf,buf_size)
end

function av_opencl_buffer_write_image(dst_cl_buf::cl_mem,cl_buffer_size::Csize_t,dst_cl_offset::Integer,src_data,plane_size,plane_num::Integer)
    ccall((:av_opencl_buffer_write_image,libavutil),Cint,(cl_mem,Csize_t,Cint,Ptr{Ptr{UInt8}},Ptr{Cint},Cint),dst_cl_buf,cl_buffer_size,dst_cl_offset,src_data,plane_size,plane_num)
end

function av_opencl_buffer_read_image(dst_data,plane_size,plane_num::Integer,src_cl_buf::cl_mem,cl_buffer_size::Csize_t)
    ccall((:av_opencl_buffer_read_image,libavutil),Cint,(Ptr{Ptr{UInt8}},Ptr{Cint},Cint,cl_mem,Csize_t),dst_data,plane_size,plane_num,src_cl_buf,cl_buffer_size)
end

function av_opencl_buffer_release(cl_buf)
    ccall((:av_opencl_buffer_release,libavutil),Void,(Ptr{cl_mem},),cl_buf)
end

function av_opencl_uninit()
    ccall((:av_opencl_uninit,libavutil),Void,())
end

function av_opencl_benchmark(device,platform::cl_platform_id,benchmark)
    ccall((:av_opencl_benchmark,libavutil),Int64,(Ptr{AVOpenCLDeviceNode},cl_platform_id,Ptr{Void}),device,platform,benchmark)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
using Mimi

# when using @defcomp, regions index must be specified
# in the run_timestep function, regions must be specified
# set_dimension! must be used to specify your regions
# update_param! each row corresponds to a time step / each column corresponds to a region
# with regionalized models, easing to save each component as separate file

include("MyModel.jl")
using .MyModel

m = construct_MyModel()
run(m)

getdataframe(m, :emissions, :E_Global)
explore(m)

max_stars_repo_name: tlorans/SustainableFinance

max_issues_repo_name: tlorans/SustainableFinance

max_forks_repo_name: tlorans/SustainableFinance


text:
# Autogenerated wrapper script for SCALAPACK32_jll for powerpc64le-linux-gnu-libgfortran3
export libscalapack32

using OpenBLAS32_jll
using MPICH_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("SCALAPACK32")
JLLWrappers.@declare_library_product(libscalapack32, "libscalapack32.so")
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS32_jll, MPICH_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libscalapack32,
        "lib/libscalapack32.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export ImaginaryFiber2D

"""
    ImaginaryFiber2D(f, y=(0., 0.))

Construct the fiber `ImaginaryFiber2D` of the bivariate polynomial `f` at `y`.
"""
function ImaginaryFiber2D(p::MP.AbstractPolynomialLike{S}, y=(0.0, 0.0)) where S
    vars = MP.variables(p)
    @assert (length(vars) == 2) "Expected a bivariate polynomial, but got $(p)."

    T = realified_coefficient_type(S)
    f_re, f_im = SP.Polynomial.( one(T) .* realify(p))

    ImaginaryFiber2D(f_re, f_im, y, zeros(T, 2, 4))
end

function update_fiber!(F::ImaginaryFiber2D{T}, y) where T
    F.y = (y[1], y[2])
    F
end

function evaluate(fiber::ImaginaryFiber2D{T}, x::SVector{2,T}) where T
    xy = SVector(x[1], x[2], fiber.y[1], fiber.y[2])
    out = SVector(SP.evaluate(fiber.f_re, xy), SP.evaluate(fiber.f_im, xy))
    out
end

function jacobian(fiber::ImaginaryFiber2D, x::AbstractVector, precomputed=false)
    U = fiber.U
    xy = SVector(x[1], x[2], fiber.y[1], fiber.y[2])
    @inbounds begin
        ∇f_re = SP.gradient(fiber.f_re, xy)
        for i=1:4
            U[1, i] = ∇f_re[i]
        end
        U[2, 3] = U[1, 1]
        U[2, 4] = U[1, 2]
        U[2, 1] = -U[1, 3]
        U[2, 2] = -U[1, 4]
    end

    # SMatrix has column-major storage
    SMatrix{2,2}(U[1,1], U[2,1], U[1, 2], U[2, 2])
end


text:
using Test
using Downloads

using OpenTelemetryUber

OpenTelemetryUber.init(
    meter_provider = MeterProvider(),
    tracer_provider = TracerProvider(;
        span_processor = CompositSpanProcessor(SimpleSpanProcessor(ConsoleExporter())),
    ),
)

console_metric_reader = MetricReader(ConsoleExporter())
console_metric_reader()

in_memory_metric_reader = MetricReader(InMemoryExporter())
in_memory_metric_reader()

prometheus_metric_reader = MetricReader(PrometheusExporter())

@test_throws Downloads.RequestError with_span("test download") do
    Downloads.download("https://julialang.org")
    Downloads.download("http://localhost:9966")  # Prometheus
    Downloads.download("http://localhost:1234")  # throw exception
end


text:
using Printf
using TimerOutputs
using Oceananigans
using Oceananigans.Advection
using Oceananigans.Utils

include("benchmark_utils.jl")

#####
##### Benchmark setup and parameters
#####

const timer = TimerOutput()

Nt = 10  # Number of iterations to use for benchmarking time stepping.

   float_types = [Float64]       # Float types to benchmark.
         archs = [CPU()]         # Architectures to benchmark on.
@hascuda archs = [CPU(), GPU()]  # Benchmark GPU on systems with CUDA-enabled GPUs.

time_steppers = (:QuasiAdamsBashforth2, :RungeKutta3)

#####
##### Run benchmarks
#####

for arch in archs, FT in float_types, ts in time_steppers
    N = arch isa CPU ? 64 : 256

    grid = RegularCartesianGrid(FT, size=(N, N, N), extent=(1, 1, 1))
    model = IncompressibleModel(architecture=arch, float_type=FT, grid=grid, timestepper=ts)

    time_step!(model, 1)  # precompile

    bn =  benchmark_name((N, N, N), string(ts), arch, FT)
    @printf("Running benchmark: %s...\n", bn)
    for i in 1:Nt
        @timeit timer bn time_step!(model, 1)
    end
end

#####
##### Print benchmark results
#####

println()
println(oceananigans_versioninfo())
println(versioninfo_with_gpu())
print_timer(timer, title="Time stepper benchmarks", sortby=:name)
println()


text:
# include(pwd() * "/src/vizinanigans.jl")
using QuickVizExample
states = [randn(10,10,10)]
statenames = ["Random Data"]
volumeslice(states, statenames = statenames)


text:
using PoissonFE
using Documenter

makedocs(;
    modules=[PoissonFE],
    authors="Evan Wright <evan.btbucket@gmail.com> and contributors",
    repo="https://github.com/ew-git/PoissonFE.jl/blob/{commit}{path}#L{line}",
    sitename="PoissonFE.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://ew-git.github.io/PoissonFE.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/ew-git/PoissonFE.jl",
)


text:
using EstiMate
using Test

@test true != false


text: Q, R = qr([12 -51 4; 6 167 -68; -4 24 -41])

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
import ..BifurcationsBase

BifurcationsBase.aliasof(::Codim1Sweep) = "Codim1Sweep"
BifurcationsBase.aliasof(::Codim1Solution) = "Codim1Solution"
BifurcationsBase.aliasof(::Codim1Solver) = "Codim1Solver"


text:
using Test
using Indexes

@testset "Indexes" begin
    # TODO
    # @test GenomicFeatures.Indexes.Tabix === GenomicFeatures.Indexes.Tabix
end


text:
########################################################################
# Expected latent info -  for 1-D latent and 1PL pars                  #
########################################################################
"""
```julia
_latent_information(latent::Latent1D, parameters::Parameters1PL)
```

# Description

It computes the information (second derivative of the likelihood) with respect to the 1-dimensional latent variable under the 1PL model.
It follows the parametrization ``a(θ - b)``.

# Arguments

- **`latent::Latent1D`** : Required. A 1-dimensional `Latent1D` latent variable.
- **`parameters::Parameters1PL`** : Required. A 1-parameter logistic parameters object.

# Output

A `Float64` scalar.
"""
function _latent_information(latent::Latent1D, parameters::Parameters1PL)
    p = _probability(latent, parameters)
    return p * (1 - p)::Float64
end

########################################################################
# Expected latent info -  for 1-D latent and 2PL pars                  #
########################################################################
"""
```julia
_latent_information(latent::Latent1D, parameters::Parameters2PL)
```

# Description

It computes the information (-second derivative of the likelihood) with respect to the 1-dimensional latent variable under the 2PL model.
It follows the parametrization ``a(θ - b)``.

# Arguments

- **`latent::Latent1D`** : Required. A 1-dimensional `Latent1D` latent variable.
- **`parameters::Parameters2PL`** : Required. A 2-parameter logistic parameters object.

# Output

A `Float64` scalar.
"""
function _latent_information(latent::Latent1D, parameters::Parameters2PL)
    p = _probability(latent, parameters)
    return Float64(p * (1 - p) * parameters.a^2)::Float64
end

########################################################################
# Expected latent info -  for 1-D latent and 3PL pars                  #
########################################################################
"""
```julia
_latent_information(latent::Latent1D, parameters::Parameters3PL)
```

# Description

It computes the information (-second derivative of the likelihood) with respect to the 1-dimensional latent variable under the 3PL model.
It follows the parametrization ``a(θ - b)``.

# Arguments

- **`latent::Latent1D`** : Required. A 1-dimensional `Latent1D` latent variable.
- **`parameters::Parameters3PL`** : Required. A 3-parameter logistic parameters object.

# Output

A `Float64` scalar.
"""
function _latent_information(latent::Latent1D, parameters::Parameters3PL)
    p = _probability(latent, parameters)
    return (p - parameters.c)^2 * (1 - p) * parameters.a^2 / (1 - parameters.c)^2 / p::Float64
end


text:
@testset "Testing SOAPlite Descriptor with dc Si" begin
@info("Testing SOAPlite Descriptor for dc Si.")
using DescriptorZoo, JuLIP, Test

at = bulk(:Si, cubic=true)
desc = soaplite(at, 6.5, n_max=4, l_max=4)
soaplite_ref = [ 0.1971183445703821,     1.4490657343533957,
                 0.5107929461890574,     1.8172037012517137,
                10.652440832199687,      3.7549653599475037,
                13.358714134715191,      1.3236182276446533,
                 4.70892160955927,      16.75252143092065,
                 0.017463358091459907,  -0.02626119328542633,
                 0.11471972827481619,   -0.4090557633886685,
                 0.03949127476872734,   -0.1725141832457655,
                 0.6151332642099767,     0.7536131359456912,
                -2.6871559169449846,     9.581583145988054,
                 0.004708428029764531,  -0.011058683554778901,
                 0.023747935881985295,  -0.16566775992816207,
                 0.02597352687386267,   -0.055776770152986335,
                 0.38910381993589627,    0.11977765298519145,
                -0.8355797976767374,     5.829080641376547,
                 0.0006486216581383904, -0.0025158961247938256,
                 0.005301754501410711,  -0.009553657612107339,
                 0.010103897726600844,  -0.023776685376270507,
                 0.011598805156323224,   0.07322789608095953,
                 0.15882882691545966,    2.018502009027488,
                 0.0031142829543606703, -0.014809827760354713,
                 0.03513043211547648,   -0.009589132006098629,
                 0.07045157270748875,   -0.16756672316768023,
                 0.04181981593527481,    0.40688151258454447,
                -0.0289359019056896,     0.6220462845081933]
soaplite_now = vcat(desc[1,:]...)
println(@test soaplite_now  ≈  soaplite_ref)
end


text:
### WARNING This is an old ODE solver which is not based on an iterator implementation. We keep it until LSODA has an iterator implementation

include("chvdiffeq.jl")

function solve(problem::PDMPProblem, algo::CHV{Tode}; verbose::Bool = false, ind_save_d = -1:1, ind_save_c = -1:1, n_jumps = Inf64, reltol = 1e-7, abstol = 1e-9, save_positions = (false, true), save_rate = false, finalizer = finalize_dummy) where {Tode <: Symbol}
	verbose && println("#"^30)
	ode = algo.ode
	@assert ode in [:cvode, :lsoda, :adams, :bdf]
	verbose && printstyled(color=:red,"--> Start CHV method\n")

	# initialise the problem. If I call twice this solve function, it should give the same result...
	init!(problem)

	# we declare the characteristics for convenience
	caract = problem.caract
	ratecache = caract.ratecache

	ti, tf = problem.tspan
	n_jumps  += 1 # to hold initial vector
	nsteps  = 1 # index for the current jump number

	xc0 = caract.xc0
	xd0 = caract.xd0

	# Set up initial simulation time
	t = ti

	X_extended = similar(xc0, length(xc0) + 1)
	for ii in eachindex(xc0)
		X_extended[ii] = xc0[ii]
	end
	X_extended[end] = ti

	#useful to use the same array, as it can be used in CHV(ode)
	Xd = caract.xd
	if ind_save_c[1] == -1
		ind_save_c = 1:length(xc0)
	end

	if ind_save_d[1] == -1
		ind_save_d = 1:length(xd0)
	end
	xc_hist = VectorOfArray([copy(xc0)[ind_save_c]])
	xd_hist = VectorOfArray([copy(xd0)[ind_save_d]])

	res_ode = zeros(2, length(X_extended))

	nsteps += 1

	# define the ODE flow, this leads to big memory saving
	if ode == :cvode || ode == :bdf
		Flow = (X0_,Xd_,Δt,r_) -> Sundials.cvode((tt,x,xdot) -> algo(xdot, x, caract, tt), X0_, [0., Δt], abstol = abstol, reltol = reltol, integrator = :BDF)
	elseif	ode==:adams
		Flow = (X0_,Xd_,Δt,r_) -> Sundials.cvode((tt,x,xdot) -> algo(xdot, x, caract, tt), X0_, [0., Δt], abstol = abstol, reltol = reltol, integrator = :Adams)
	elseif ode==:lsoda
		Flow = (X0_,Xd_,Δt,r_) -> LSODA.lsoda((tt,x,xdot,data) -> algo(xdot, x, caract, tt), X0_, [0., Δt], abstol = abstol, reltol = reltol)
	end

	# we use the first time interval from the one generated by the constructor PDMPProblem
	δt = problem.simjptimes.tstop_extended

	# Main loop
	while (t < tf) && (nsteps < n_jumps)

		verbose && println("--> t = ", t," - δt = ", δt, ",nstep =  ", nsteps)

		res_ode .= Flow(X_extended, Xd, δt, ratecache.rate)

		verbose && println("--> ode solve is done!")

		# this holds the new state of the continuous component
		@inbounds for ii in eachindex(X_extended)
			X_extended[ii] = res_ode[end, ii]
		end

		# this is the next jump time
		t = res_ode[end, end]

		caract.R(ratecache.rate, X_extended, Xd, caract.parms, t, false)

		# jump time:
		if (t < tf) && nsteps < n_jumps
			# Update event
			ev = pfsample(ratecache.rate)

			# we perform the jump
			affect!(caract.pdmpjump, ev, X_extended, Xd, caract.parms, t)

			verbose && println("--> Which reaction? => ", ev)
			verbose && println("--> xd = ", Xd)

			# save state, post-jump
			pushTime!(problem, t)
			push!(xc_hist, X_extended[ind_save_c])
			push!(xd_hist, Xd[ind_save_d])

			save_rate && push!(problem.rate_hist, caract.R(ratecache.rate, X_extended, Xd, caract.parms, t, true)[1])

			finalizer(ratecache.rate, caract.xc, caract.xd, caract.parms, t)

			δt = - log(rand())

		else
			if ode in [:cvode, :bdf, :adams]
				res_ode_last = Sundials.cvode((tt, x, xdot)->caract.F(xdot,x,Xd,caract.parms,tt), xc_hist[end], [problem.time[end], tf], abstol = 1e-9, reltol = 1e-7)
			else#if ode==:lsoda
				res_ode_last = LSODA.lsoda((tt, x, xdot, data)->caract.F(xdot,x,Xd,caract.parms,tt), xc_hist[end], [problem.time[end], tf], abstol = 1e-9, reltol = 1e-7)
			end
			t = tf

			# save state
			pushTime!(problem, tf)
			push!(xc_hist, res_ode_last[end,ind_save_c])
			push!(xd_hist, Xd[ind_save_d])
		end
		nsteps += 1
	end
	verbose && println("--> Done")
	verbose && println("--> xc = ", xd_hist[:,1:nsteps-1])
	return PDMPResult(problem.time, xc_hist, xd_hist, problem.rate_hist, save_positions, length(problem.time), 0)
end


text:
"""
Provides interface for loading parameteres of components and mixtures from
*separated values* files.

See `ComponentDatabase`, `BinaryMixtureDatabase`, `getentry`, `getmatrix`.

Out-of-box database stored in `CubicEoSDatabase.Data` module.
"""
module CubicEoSDatabase
    using DelimitedFiles

    export Data
    export ComponentDatabase, MixtureDatabase
    export data, header, reference, source
    export getentry, getmatrix

    include("types.jl")
    include("errors.jl")
    include("requests.jl")
    include("Data.jl")
end


text:
@with_kw struct DeepPolyHeuristic <: Solver
    max_vars::Int64 = 10
    var_frac::Float64 = 0.5
    # switch to hard coded function later, variable function is slow performance wise
    get_fresh_var_idxs = get_fresh_var_idxs_largest_range
end

function forward_linear(solver::DeepPolyHeuristic, L::LayerNegPosIdx, input::SymbolicIntervalHeur)
    output_Low, output_Up = interval_map(L.W_neg, L.W_pos, input.sym.Low, input.sym.Up)
    output_Up[:, end] += L.bias
    output_Low[:, end] += L.bias
    sym = SymbolicInterval(output_Low, output_Up, domain(input))
    return SymbolicIntervalHeur(sym, input.lbs, input.ubs, input.var_los, input.var_his, input.max_vars, input.importance)
end

function forward_act(solver::DeepPolyHeuristic, L::LayerNegPosIdx{ReLU}, input::SymbolicIntervalHeur)
    n_node = n_nodes(L)
    n_sym = size(input.sym.Low, 2) - 1 # last column is constant term
    n_in = dim(domain(input))
    current_n_vars = n_sym - n_in

    subs_sym_lo, subs_sym_hi = substitute_variables(input.sym.Low, input.sym.Up,
                                                    input.var_los, input.var_his,
                                                    n_in, current_n_vars)
    # TODO: write better function for bounds instead of this monstrosity!!!
    tmp_sym = SymbolicInterval(subs_sym_lo, subs_sym_hi, domain(input))
    los, his = bounds(tmp_sym, input.lbs[L.index], input.ubs[L.index])

    ##### Calculate importance score ######
    crossing = is_crossing.(los, his)
    layer_importance = sum(abs.(subs_sym_lo[crossing, :]), dims=1) .+ sum(abs.(subs_sym_hi[crossing, :]), dims=1)
    importance = input.importance .+ layer_importance[1:end-1]  # constant term doesn't need importance

    fv_idxs = solver.get_fresh_var_idxs(input.max_vars, current_n_vars, los, his, solver.var_frac)
    #n_vars = min(input.max_vars - current_n_vars, floor(Int, solver.var_frac * n_node))
    #fv_idxs = fresh_var_largest_range(los, his, n_vars)
    n_vars = length(fv_idxs)

    out_Low, out_Up = zeros(n_node, n_sym + 1 + n_vars), zeros(n_node, n_sym + 1 + n_vars)

    slopes = relaxed_relu_gradient.(los, his)
    out_Up[:, (1:n_sym) ∪ [end]] .= input.sym.Up .* slopes
    out_Up[:, end] .+= slopes .* max.(-los, 0)

    # also apply relaxation to substituted variables
    subs_sym_hi .= subs_sym_hi .* slopes
    subs_sym_hi[:, end] .+= slopes .* max.(-los, 0)

    out_Low[:, (1:n_sym) ∪ [end]] .= input.sym.Low .* relaxed_relu_gradient_lower.(los, his)

    # also apply relaxation to substituted variables
    subs_sym_lo .= subs_sym_lo .* relaxed_relu_gradient_lower.(los, his)

    for (i, v) in enumerate(fv_idxs)
        # store symbolic bounds on fresh variables
        input.var_los[current_n_vars + i, :] .= subs_sym_lo[v, :]
        input.var_his[current_n_vars + i, :] .= subs_sym_hi[v, :]

        # set corresponding entry to unit-vec
        out_Low[v,:] .= unit_vec(n_sym + i, n_sym + 1 + n_vars)
        out_Up[v,:]  .= unit_vec(n_sym + i, n_sym + 1 + n_vars)
    end

    sym = SymbolicInterval(out_Low, out_Up, domain(input))
    output = SymbolicIntervalHeur(sym, input.lbs, input.ubs, input.var_los, input.var_his, input.max_vars, importance)
    output.lbs[L.index] .= los
    output.ubs[L.index] .= his
    return output
end

function forward_act(solver::DeepPolyHeuristic, L::LayerNegPosIdx{Id}, input::SymbolicIntervalHeur)
    n_node = n_nodes(L)
    n_sym = size(input.sym.Low, 2) - 1 # last column is constant term
    n_in = dim(domain(input))
    current_n_vars = n_sym - n_in

    subs_sym_lo, subs_sym_hi = substitute_variables(input.sym.Low, input.sym.Up,
                                                    input.var_los, input.var_his,
                                                    n_in, current_n_vars)
    # TODO: write better function for bounds instead of this monstrosity!!!
    tmp_sym = SymbolicInterval(subs_sym_lo, subs_sym_hi, domain(input))
    los, his = bounds(tmp_sym, input.lbs[L.index], input.ubs[L.index])

    sym = SymbolicInterval(input.sym.Low, input.sym.Up, domain(input))
    output = SymbolicIntervalHeur(sym, input.lbs, input.ubs, input.var_los, input.var_his, input.max_vars, input.importance)
    output.lbs[end] .= los
    output.ubs[end] .= his
    return output

    return input
end

function get_fresh_var_idxs_largest_range(max_vars, current_n_vars, los, his, var_frac)
    n_node = length(los)
    n_vars = min(max_vars - current_n_vars, floor(Int, var_frac * n_node))
    return fresh_var_largest_range(los, his, n_vars)
end


text:
# define the api for users

"""
    $SIGNATURES

Approximate the pointwise the product of functionals on manifolds using KernelDensityEstimate.

Example
-------
```julia
using ApproxManifoldProducts

# two densities on a cylinder
p = manikde!(randn(2,100), (:Euclid, :Circular) )

pts2a = 3.0*randn(1,100).+5.0
pts2b = TransformUtils.wrapRad.(0.5*randn(1,100).+pi)
q = manikde!([pts2a;pts2b], (:Euclid, :Circular) )

# approximate the product between hybrid manifold densities
pq = manifoldProduct([p;q], (:Euclid, :Circular))

# convenient plotting (work in progress...)
# TODO update docs

# direct histogram plot
using Gadfly
plot( x=getPoints(pq)[1,:], y=getPoints(pq)[2,:], Geom.histogram2d )
```

"""
function manifoldProduct(ff::Vector{BallTreeDensity},
                         manif::T;
                         Niter=1  )::BallTreeDensity where {T <: Tuple}
  #

  ndims = Ndim(ff[1])
  N = Npts(ff[1])

  addopT, diffopT, getManiMu, getManiLam = buildHybridManifoldCallbacks(manif)

  bws = ones(ndims)

  dummy = kde!(rand(ndims,N), bws, addopT, diffopT );

  pGM, = prodAppxMSGibbsS(dummy, ff,
                          nothing, nothing, Niter=1,
                          addop=addopT,
                          diffop=diffopT,
                          getMu=getManiMu  );
  #

  bws[:] = getKDEManifoldBandwidths(pGM, manif)
  kde!(pGM, bws, addopT, diffopT)
end

#


text:
@doc raw"""
`SPSDE`: Stratonovich Split Partitioned Stochastic Differential Equation

Defines a partitioned stochastic differential initial value problem
```math
\begin{aligned}
dq (t) &=   v(t, q(t)) \, dt + B(t, q(t)) \circ dW , & q(t_{0}) &= q_{0} , \\
dp (t) &= [ f_1(t, q(t)) + f_2(t, q(t)) ] \, dt + [ G_1(t, q(t)) + G_2(t, q(t)) ] \circ dW , & p(t_{0}) &= p_{0} ,
\end{aligned}
```
with the drift vector fields ``v`` and ``f_i``, diffusion matrices ``B`` and ``G_i``,
initial conditions ``q_{0}`` and ``p_{0}``, the dynamical variables ``(q,p)`` taking
values in ``\mathbb{R}^{d} \times \mathbb{R}^{d}``, and the m-dimensional Wiener process W

### Parameters

* `DT <: Number`: data type
* `TT <: Real`: time step type
* `AT <: AbstractArray{DT}`: array type
* `vType <: Function`: type of `v`
* `f1Type <: Function`: type of `f1`
* `f2Type <: Function`: type of `f2`
* `BType <: Function`: type of `B`
* `G1Type <: Function`: type of `G1`
* `G2Type <: Function`: type of `G2`
* `pType <: Union{NamedTuple,Nothing}`: parameters type

### Fields

* `d`:  dimension of dynamical variable ``q`` and the vector fields ``vi``
* `m`:  dimension of the Wiener process
* `ni`: number of initial conditions
* `ns`: number of sample paths
* `v` :  function computing the drift vector field for the position variable ``q``
* `f1`:  function computing the drift vector field for the momentum variable ``p``
* `f2`:  function computing the drift vector field for the momentum variable ``p``
* `B` :  function computing the d x m diffusion matrix for the position variable ``q``
* `G1`:  function computing the d x m diffusion matrix for the momentum variable ``p``
* `G2`:  function computing the d x m diffusion matrix for the momentum variable ``p``
* `t₀`: initial time
* `q₀`: initial condition for dynamical variable ``q`` (may be a random variable itself)
* `p₀`: initial condition for dynamical variable ``p`` (may be a random variable itself)
* `parameters`: either a `NamedTuple` containing the equations parameters or `nothing`
* `periodicity`: determines the periodicity of the state vector `q` for cutting periodic solutions

The functions `v`, `f`, `B` and `G`, providing the drift vector fields and diffusion matrices, take four arguments,
`v(t, q, p, v)`, `f(t, q, p, f)`, `B(t, q, p,  B)` and `G(t, q, p, G)`, where `t` is the current time, `(q, p)` is the
current solution vector, and `v`, `f`, `B` and `G` are the variables which hold the result
of evaluating the vector fields ``v``, ``f`` and the matrices ``B``, ``G`` on `t` and `(q,p)`.

### Constructors

```julia
SPSDE(m, ns, v, f1, f2, B, G1, G2, t₀, q₀, p₀; parameters=NullParameters(), periodicity=zero(q₀[begin]))
SPSDE(m, ns, v, f1, f2, B, G1, G2, q₀::StateVector, p₀::StateVector; kwargs...) = SPSDE(m, ns, v, f1, f2, B, G1, G2, 0.0, q₀, p₀; kwargs...)
SPSDE(m, ns, v, f1, f2, B, G1, G2, t₀, q₀::State, p₀::State; kwargs...) = SPSDE(m, ns, v, f1, f2, B, G1, G2, t₀, [q₀], [p₀]; kwargs...)
SPSDE(m, ns, v, f1, f2, B, G1, G2, q₀::State, p₀::State; kwargs...) = SPSDE(m, ns, v, f1, f2, B, G1, G2, 0.0, q₀, p₀; kwargs...)
```

### Example

```julia
    function v(λ, t, q, v)
        v[1] = λ*q[1]
        v[2] = λ*q[2]
    end

    function B(μ, t, q, B)
        B[1] = μ*q[1]
        B[2] = μ*q[2]
    end

    t₀ = 0.
    q₀ = [1., 1.]
    λ  = 2.
    μ  = 1.

    v_sde = (t, q, v) -> v(λ, t, q, v)
    B_sde = (t, q, B) -> B(μ, t, q, B)

    sde = SDE(v_sde, B_sde, t₀, q₀)
```
"""
struct SPSDE{dType <: Number, tType <: Real, arrayType <: AbstractArray{dType},
             vType <: Function, f1Type <: Function, f2Type <: Function,
             BType <: Function, G1Type <: Function, G2Type <: Function,
             pType <: OptionalParameters} <: AbstractEquationPSDE{Nothing,Nothing,Nothing}
    d::Int
    m::Int
    ns::Int
    v ::vType
    f1::f1Type
    f2::f2Type
    B ::BType
    G1::G1Type
    G2::G2Type
    t₀::tType
    q₀::Vector{arrayType}
    p₀::Vector{arrayType}
    parameters::pType
    periodicity::arrayType

    function SPSDE(m, ns,
                    v::vType, f1::f1Type, f2::f2Type,
                    B::BType, G1::G1Type, G2::G2Type,
                    t₀::tType, q₀::Vector{arrayType}, p₀::Vector{arrayType};
                    parameters=NullParameters(), periodicity=zero(q₀[begin])) where {
                        dType <: Number, tType <: Real, arrayType <: AbstractArray{dType},
                        vType <: Function, f1Type <: Function, f2Type <: Function,
                        BType <: Function, G1Type <: Function, G2Type <: Function}

        d  = length(q₀[begin])
        ni = length(q₀)

        @assert length(q₀) == length(p₀)
        @assert all(length(q) == d for q in q₀)
        @assert all(length(p) == d for p in p₀)

        @assert ns ≥ 1
        @assert ni == 1 || ns == 1
        # either multiple deterministic initial conditions and one sample path
        # or one deterministic initial condition and multiple sample paths

        new{dType,tType,arrayType,vType,f1Type,f2Type,BType,G1Type,G2Type,typeof(parameters)}(d, m, ns, v, f1, f2, B, G1, G2, t₀, q₀, p₀, parameters, periodicity)
    end
end

SPSDE(m, ns, v, f1, f2, B, G1, G2, q₀::StateVector, p₀::StateVector; kwargs...) = SPSDE(m, ns, v, f1, f2, B, G1, G2, 0.0, q₀, p₀; kwargs...)
SPSDE(m, ns, v, f1, f2, B, G1, G2, t₀, q₀::State, p₀::State; kwargs...) = SPSDE(m, ns, v, f1, f2, B, G1, G2, t₀, [q₀], [p₀]; kwargs...)
SPSDE(m, ns, v, f1, f2, B, G1, G2, q₀::State, p₀::State; kwargs...) = SPSDE(m, ns, v, f1, f2, B, G1, G2, 0.0, q₀, p₀; kwargs...)

const SPSDEPT{PT,DT,TT,AT,VT,F1T,F2T,BT,G1T,G2T} = SPSDE{DT,TT,AT,VT,F1T,F2T,BT,G1T,G2T,PT} # type alias for dispatch on parameters type parameter

Base.hash(sde::SPSDE, h::UInt) = hash(sde.d, hash(sde.m, hash(sde.ns,
                                 hash(sde.v, hash(sde.f1, hash(sde.f2,
                                 hash(sde.B, hash(sde.G1, hash(sde.G2,
                                 hash(sde.t₀, hash(sde.q₀, hash(sde.p₀,
                                 hash(sde.parameters, hash(sde.periodicity, h))))))))))))))

Base.:(==)(sde1::SPSDE, sde2::SPSDE) = (
                                sde1.d == sde2.d
                             && sde1.m == sde2.m
                             && sde1.ns == sde2.ns
                             && sde1.v  == sde2.v
                             && sde1.f1 == sde2.f1
                             && sde1.f2 == sde2.f2
                             && sde1.B  == sde2.B
                             && sde1.G1 == sde2.G1
                             && sde1.G2 == sde2.G2
                             && sde1.t₀ == sde2.t₀
                             && sde1.q₀ == sde2.q₀
                             && sde1.p₀ == sde2.p₀
                             && sde1.parameters  == sde2.parameters
                             && sde1.periodicity == sde2.periodicity)

function Base.similar(equ::SPSDE, t₀::Real, q₀::StateVector, p₀::StateVector, ns::Int=equ.ns)
    @assert all([length(q) == ndims(equ) for q in q₀])
    @assert all([length(p) == ndims(equ) for p in p₀])
    SPSDE(equ.m, ns, equ.v, equ.f1, equ.f2, equ.B, equ.G1, equ.G2, t₀, q₀, p₀; parameters=equ.parameters, periodicity=equ.periodicity)
end

Base.similar(equ::SPSDE, t₀::Real, q₀::State, p₀::State, ns::Int=equ.ns; kwargs...) = similar(equ, t₀, [q₀], [p₀], ns; kwargs...)
Base.similar(equ::SPSDE, q₀::Union{State,StateVector}, p₀::Union{State,StateVector}, ns::Int=equ.ns; kwargs...) = similar(equ, equ.t₀, q₀, p₀, ns; kwargs...)

Base.ndims(equ::SPSDE) = equ.d
Base.axes(equ::SPSDE) = axes(equ.q₀[begin])
GeometricBase.nsamples(equ::SPSDE) = length(equ.q₀)
GeometricBase.periodicity(equ::SPSDE) = equ.periodicity

initial_conditions(equ::SPSDE) = (equ.t₀, equ.q₀, equ.p₀)

hasparameters(::SPSDEPT{<:NullParameters}) = false
hasparameters(::SPSDEPT{<:NamedTuple}) = true

_get_v( equ::SPSDE) = hasparameters(equ) ? (t,q,p,v) -> equ.v( t, q, p, v, equ.parameters) : equ.v
_get_f1(equ::SPSDE) = hasparameters(equ) ? (t,q,p,f) -> equ.f1(t, q, p, f, equ.parameters) : equ.f1
_get_f2(equ::SPSDE) = hasparameters(equ) ? (t,q,p,f) -> equ.f2(t, q, p, f, equ.parameters) : equ.f2
_get_B( equ::SPSDE) = hasparameters(equ) ? (t,q,p,B) -> equ.B( t, q, p, B, equ.parameters) : equ.B
_get_G1(equ::SPSDE) = hasparameters(equ) ? (t,q,p,G) -> equ.G1(t, q, p, G, equ.parameters) : equ.G1
_get_G2(equ::SPSDE) = hasparameters(equ) ? (t,q,p,G) -> equ.G2(t, q, p, G, equ.parameters) : equ.G2

function functions(equ::SPSDE)
    names = (:v,:f1,:f2,:B,:G1,:G2)
    equs  = (_get_v(equ), _get_f1(equ), _get_f2(equ), _get_B(equ), _get_G1(equ), _get_G2(equ))

    NamedTuple{names}(equs)
end


text:
import Oracle
import Oracle.Timestamps

if VERSION < v"0.7-"
    using Base.Test
    using Missings
    using Base.Dates
else
    using Test
    using Dates
end

@testset "Timestamps" begin

    @testset "Timestamp" begin
        ts = Timestamps.Timestamp(2018, 12, 31, 23, 58, 59, 999_200_300)

        @test year(ts) == 2018
        @test month(ts) == 12
        @test day(ts) == 31
        @test hour(ts) == 23
        @test minute(ts) == 58
        @test second(ts) == 59
        @test millisecond(ts) == 999
        @test microsecond(ts) == 200
        @test nanosecond(ts) == 300
    end

    @testset "TimestampTZ" begin
        @test_throws AssertionError Timestamps.TimestampTZ(false, 2018, 12, 31, 23, 58, 59, 999_200_300, -3, 30)

        ts = Timestamps.TimestampTZ(false, 2018, 12, 31, 23, 58, 59, 999_200_300, -3, -30)
        @test Timestamps.is_ltz(ts) == false
        @test year(ts) == 2018
        @test month(ts) == 12
        @test day(ts) == 31
        @test hour(ts) == 23
        @test minute(ts) == 58
        @test second(ts) == 59
        @test millisecond(ts) == 999
        @test microsecond(ts) == 200
        @test nanosecond(ts) == 300
    end

    @testset "TimestampLTZ" begin
        ts = Timestamps.TimestampTZ(true, 2018, 12, 31, 23, 58, 59, 999_200_300, -3, -30)
        @test Timestamps.is_ltz(ts) == true
        @test year(ts) == 2018
        @test month(ts) == 12
        @test day(ts) == 31
        @test hour(ts) == 23
        @test minute(ts) == 58
        @test second(ts) == 59
        @test millisecond(ts) == 999
        @test microsecond(ts) == 200
        @test nanosecond(ts) == 300
    end

    @testset "OraTimestamp" begin
        @testset "Constructors" begin
            ora_ts = Oracle.OraTimestamp(Timestamps.Timestamp(2018, 12, 31, 23, 58, 59, 999_200_300))
            @test ora_ts.year == 2018
            @test ora_ts.month == 12
            @test ora_ts.day == 31
            @test ora_ts.hour == 23
            @test ora_ts.minute == 58
            @test ora_ts.second == 59
            @test ora_ts.fsecond == 999_200_300

            ora_ts_tz = Oracle.OraTimestamp(Timestamps.TimestampTZ(false, 2018, 12, 31, 23, 58, 59, 999_200_400, -3, -30))
            @test ora_ts_tz.year == 2018
            @test ora_ts_tz.month == 12
            @test ora_ts_tz.day == 31
            @test ora_ts_tz.hour == 23
            @test ora_ts_tz.minute == 58
            @test ora_ts_tz.second == 59
            @test ora_ts_tz.fsecond == 999_200_400
            @test ora_ts_tz.tzHourOffset == -3
            @test ora_ts_tz.tzMinuteOffset == -30
        end

        @testset "ts Conversion" begin
            ts = Timestamps.Timestamp(2018, 12, 31, 23, 58, 59, 999_200_300)
            @test Oracle.parse_timestamp(Oracle.OraTimestamp(ts)) == ts
        end

        @testset "DateTime comparison" begin
            dt = now()
            ts = Timestamps.Timestamp(dt)
            @test ts == dt
            @test dt == ts
        end
    end
end