text:
using Test
using StaticArrays, IntervalSets
import ClimaCore.DataLayouts: IJFH
import ClimaCore:
    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry
using StaticArrays, IntervalSets, LinearAlgebra

using OrdinaryDiffEq

@testset "Scalar Poisson problem - ∇⋅∇ = f in 2D" begin
    # Poisson equation
    # - ∇⋅(∇ u(x,y)) = f(x,y)

    # True solution (eigenfunction): u(x,y) = sin(c₁ + k₁ x) * sin(c₂ + k₂ y)
    # => - ∇⋅(∇ u(x,y)) = f(x,y) = (k₁^2 + k₂^2) * u(x,y)

    FT = Float64

    domain = Domains.RectangleDomain(
        Geometry.XPoint{FT}(-π) .. Geometry.XPoint{FT}(π),
        Geometry.YPoint{FT}(-π) .. Geometry.YPoint{FT}(π),
        x1periodic = true,
        x2periodic = true,
    )

    mesh = Meshes.RectilinearMesh(domain, 10, 10)
    grid_topology = Topologies.Topology2D(mesh)

    Nq = 6
    quad = Spaces.Quadratures.GLL{Nq}()
    points, weights = Spaces.Quadratures.quadrature_points(Float64, quad)
    space = Spaces.SpectralElementSpace2D(grid_topology, quad)

    # Define eigensolution
    u = map(Fields.coordinate_field(space)) do coord
        c₁ = 0.0
        k₁ = 1.0
        c₂ = 2.0
        k₂ = 3.0
        sin(c₁ + k₁ * coord.x) * sin(c₂ + k₂ * coord.y)
    end

    function laplacian_u(space)
        coords = Fields.coordinate_field(space)
        laplacian_u = map(coords) do coord
            c₁ = 0.0
            k₁ = 1.0
            c₂ = 2.0
            k₂ = 3.0
            (k₁^2 + k₂^2) * sin(c₁ + k₁ * coord.x) * sin(c₂ + k₂ * coord.y)
        end

        return laplacian_u
    end

    function diffusion(u)
        grad = Operators.Gradient()
        wdiv = Operators.WeakDivergence()
        diff = @. -wdiv(grad(u))
        Spaces.weighted_dss!(diff)
        return diff
    end

    # Call the diffusion operator
    diff = diffusion(u)

    exact_solution = laplacian_u(space)

    @test norm(diff .- exact_solution) ≤ 5e-3
end


text:
"""
Format to which the shader conforms. It can be a text format (with a language such as FormatGLSL or FormatHLSL),
or a binary format (such as SPIR-V).
"""
abstract type ShaderFormat end

Base.broadcastable(x::ShaderFormat) = Ref(x)

"""
Human-readable shader format, usually written in a shading language
such as FormatGLSL or FormatHLSL.
"""
abstract type TextFormat <: ShaderFormat end

"""
Efficient format that requires disassembly by a specialized program to be human-readable.
"""
abstract type BinaryFormat <: ShaderFormat end

"""
OpenGL Shading Language format.
"""
struct FormatGLSL <: TextFormat end

"""
High Level Shading Language format, used by DirectX.
"""
struct FormatHLSL <: TextFormat end

"""
Standard Portable Intermediate Representation format, from The Khronos Group.
"""
struct FormatSPIRV <: BinaryFormat end

"""
    ShaderFile(file, format)
    ShaderFile(file, stage)

Shader resource as a file.
If one of the format or stage is not provided but the other is, then the former
will be retrieved from the file extension, if possible, or an error will be thrown.

## Examples

```julia
julia> ShaderFile("my_shader.geom", FormatGLSL())
ShaderFile{FormatGLSL}

julia> ShaderFile("my_shader.FormatGLSL", GeometryStage())
ShaderFile{FormatGLSL}
```
"""
struct ShaderFile{F<:ShaderFormat}
    file::String
    format::F
    stage::ShaderStageFlag
end

ShaderFile(file, format::TextFormat) = ShaderFile(file, format, shader_stage(file, format))
ShaderFile(file, stage::ShaderStageFlag) = ShaderFile(file, shader_format(file), stage)

bytes(shader::ShaderFile) = open(io -> readavailable(io), shader.file)

Base.convert(::Type{ShaderFile{F1}}, shader::ShaderFile{F2}) where {F1<:TextFormat, F2<:TextFormat} = ShaderFile(shader.file, F1(), shader.stage)

function file_ext(::Type{<:Union{FormatGLSL,FormatHLSL}}, stage::ShaderStageFlag)
    '.' * @match stage begin
        &SHADER_STAGE_VERTEX_BIT => "vert"
        &SHADER_STAGE_FRAGMENT_BIT => "frag"
        &SHADER_STAGE_TESSELLATION_CONTROL_BIT => "tesc"
        &SHADER_STAGE_TESSELLATION_EVALUATION_BIT => "tese"
        &SHADER_STAGE_GEOMETRY_BIT => "geom"
        &SHADER_STAGE_COMPUTE_BIT => "comp"
        &SHADER_STAGE_RAYGEN_BIT_KHR => "rgen"
        &SHADER_STAGE_INTERSECTION_BIT_KHR => "rint"
        &SHADER_STAGE_ANY_HIT_BIT_KHR => "rahit"
        &SHADER_STAGE_CLOSEST_HIT_BIT_KHR => "rchit"
        &SHADER_STAGE_MISS_BIT_KHR => "rmiss"
        &SHADER_STAGE_CALLABLE_BIT_KHR => "rcall"
        &SHADER_STAGE_MESH_BIT_NV => "mesh"
        &SHADER_STAGE_TASK_BIT_NV => "task"
        _ => error("Unknown stage $stage")
    end
end

"""
    shader_stage(file_ext)

Automatically detect a shader stage from a file extension.
Can only be used with [`FormatGLSL`](@ref) and [`FormatHLSL`](@ref).

## Examples

```julia
julia> shader_stage("my_shader.frag", FormatGLSL()) == SHADER_STAGE_FRAGMENT_BIT
true

julia> shader_stage("my_shader.geom", FormatHLSL()) == SHADER_STAGE_GEOMETRY_BIT
true
```
"""
function shader_stage(file::AbstractString, ::Union{FormatHLSL,FormatGLSL})
    _, file_ext = splitext(file)
    @match file_ext begin
        ".vert" => SHADER_STAGE_VERTEX_BIT
        ".frag" => SHADER_STAGE_FRAGMENT_BIT
        ".tesc" => SHADER_STAGE_TESSELLATION_CONTROL_BIT
        ".tese" => SHADER_STAGE_TESSELLATION_EVALUATION_BIT
        ".geom" => SHADER_STAGE_GEOMETRY_BIT
        ".comp" => SHADER_STAGE_COMPUTE_BIT
        ".rgen" => SHADER_STAGE_RAYGEN_BIT_KHR
        ".rint" => SHADER_STAGE_INTERSECTION_BIT_KHR
        ".rahit" => SHADER_STAGE_ANY_HIT_BIT_KHR
        ".rchit" => SHADER_STAGE_CLOSEST_HIT_BIT_KHR
        ".rmiss" => SHADER_STAGE_MISS_BIT_KHR
        ".rcall" => SHADER_STAGE_CALLABLE_BIT_KHR
        ".mesh" => SHADER_STAGE_MESH_BIT_NV
        ".task" => SHADER_STAGE_TASK_BIT_NV
        _ => error("Unknown file extension $file_ext")
    end
end

"""
    shader_format(file_ext)

Automatically detect a [`ShaderFormat`](@ref) from the file extension.
Currently, only .spv, .FormatGLSL and .FormatHLSL are recognized, mapping to
`FormatSPIRV()`, `FormatGLSL()` and `FormatHLSL()`.

## Examples

```julia
julia> shader_format("my_shader.FormatGLSL")
FormatGLSL()
julia> shader_format("my_shader.spv")
FormatSPIRV()
```
"""
function shader_format(file::AbstractString)
    _, file_ext = splitext(file)
    @match file_ext begin
        ".spv" => FormatSPIRV()
        ".glsl" => FormatGLSL()
        ".hlsl" => FormatHLSL()
        _ => error("Unknown file extension $file_ext")
    end
end

ext: jl


text:
using Random;
using DelimitedFiles;
using Random: seed!

include("DoME.jl");

seed!(1); # Make the experiment repeatable

numData = 1000;

inputs = Array{Float64,2}(undef, numData, 3);
targets = Array{Float64,1}(undef, numData);
G = 6.67392e-11; # N * (m/kg)^2

minimumMass = 1e23; maximumMass = 1e25;
minimumDist = 1e8;  maximumDist = 1e12;

planet1Mass = 10. .^(log10(minimumMass) .+ rand(numData).*(log10(maximumMass)-log10(minimumMass)));
planet2Mass = 10. .^(log10(minimumMass) .+ rand(numData).*(log10(maximumMass)-log10(minimumMass)));
distance    = 10. .^(log10(maximumDist) .+ rand(numData).*(log10(minimumDist)-log10(maximumDist)));

inputs = [planet1Mass planet2Mass distance];
targets     = G.*planet1Mass.*planet2Mass./(distance.^2);

# As explained in the paper, with so high differences in the variables ranges,
#  in irder to avoid rounding errors, the goal is the mean of the targets
(trainingMSE, _, _, bestTree) = dome(inputs, targets;
    goalMSE = mean(targets) ,
    strategy = Strategy1 ,
    showText = true
    );

expression = string(bestTree);
expression = replace(expression, "X1" => "M1");
expression = replace(expression, "X2" => "M2");
expression = replace(expression, "X3" => "d");
println("Best expression found: ", expression);

########################################################################################################
# A different execution: without so big differences in inputs and targets (not documented in the paper)

println("---------------------------------------------------------------------------------------------")
println("Additional exucution with lower values in masses and distances (not documented in the paper):")

maximumMass/minimumMass
planet1Mass ./= minimumMass;
planet2Mass ./= minimumMass;
distance ./= minimumDist;
G = G*minimumMass/minimumDist;

inputs = [planet1Mass planet2Mass distance];
targets     = G.*planet1Mass.*planet2Mass./(distance.^2);

# In this case, since no so big rounding errors are expected, the MSE goal can be 0
# However, some small rounding erros might happen. For this reason, it is better to have a small MSE goal (1e-10)

(trainingMSE, _, _, bestTree) = dome(inputs, targets;
    goalMSE = 1e-10 ,
    strategy = Strategy1 ,
    showText = true
    );

expression = string(bestTree);
expression = replace(expression, "X1" => "M1");
expression = replace(expression, "X2" => "M2");
expression = replace(expression, "X3" => "d");
println("Best expression found: ", expression);


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# https://cometscome.github.io/JuliaFromFortran/build/chapter2/01/

# %% [markdown]
# https://docs.julialang.org/en/v1/manual/variables-and-scoping/#Loops-and-Comprehensions

# %%
function f()
    i = -99
    for i in 1:3 end
    i
end

@show f();

# %%
function g()
    local i
    for outer i in 1:3 end
    i
end

@show g();

# %%
function mandel()
    nx = 61; ny = 31; maxiter = 90
    x0 = -2.0; x1 = 2.0; y0 = -2.0; y1 = 2.0
    mandelbrot = zeros(Int64,nx, ny)
    mandelbrot .= maxiter
    for iy=1:ny
        y = y0 + (y1 - y0) * (iy - 1) / real(ny - 1)
        for ix=1:nx
            x = x0 + (x1 - x0) * (ix - 1) / real(nx - 1)
            c = x + im*y
            z = 0im
            for iter=0:maxiter
                z = z * z + c
                if abs(z) > 2
                    mandelbrot[ix,iy] = iter
                    break
                end
            end

        end
    end
    for iy=1:ny
        for ix=1:nx
            if mandelbrot[ix,iy] == maxiter
                print("*")
            else
                print((mandelbrot[ix,iy]+9) ÷ 10)
            end
        end
        println("\t")
    end
end
mandel()

# %%
function mandel_outer()
    nx = 61; ny = 31; maxiter = 90
    x0 = -2.0; x1 = 2.0; y0 = -2.0; y1 = 2.0
    mandelbrot = zeros(Int64,nx, ny)
    for iy=1:ny
        y = y0 + (y1 - y0) * (iy - 1) / real(ny - 1)
        for ix=1:nx
            x = x0 + (x1 - x0) * (ix - 1) / real(nx - 1)
            c = x + im*y
            z = 0im
            local iter; for outer iter=0:maxiter
                z = z * z + c
                abs(z) > 2 && break
            end
            mandelbrot[ix,iy] = iter
        end
    end
    for iy=1:ny
        for ix=1:nx
            if mandelbrot[ix,iy] == maxiter
                print("*")
            else
                print((mandelbrot[ix,iy]+9) ÷ 10)
            end
        end
        println("\t")
    end
end
mandel_outer()

# %%

max_stars_repo_path: 0019/for outer.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0019/for outer.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0019/for outer.jl

max_forks_repo_name: genkuroki/public


text:
module TestOpenCL
    using Base.Test

    using OpenCL

    function create_test_buffer()
        ctx = cl.create_some_context()
        queue = cl.CmdQueue(ctx)
        testarray = zeros(Float32, 1000)
        buf = cl.Buffer(Float32, ctx, (:rw, :copy), hostbuf=testarray)
        return (queue, buf, testarray, ctx)
    end

    include("test_platform.jl")
    include("test_context.jl")
    include("test_device.jl")
    include("test_cmdqueue.jl")
    include("test_macros.jl")
    include("test_event.jl")
    include("test_program.jl")
    include("test_kernel.jl")
    include("test_behaviour.jl")
    include("test_memory.jl")
    include("test_buffer.jl")
    include("test_array.jl")

end # module


text:
#####
##### primitive distributions
#####

"Standard multivariate normal with the given dimension `K`."
struct StandardMultivariateNormal <: SamplingLogDensity
    K::Int
end

dimension(ℓ::StandardMultivariateNormal) = ℓ.K

logdensity(ℓ::StandardMultivariateNormal, x) = -0.5 * sum(abs2, x) - ℓ.K/2 * log(2*π)

logdensity_and_gradient(ℓ::StandardMultivariateNormal, x) = logdensity(ℓ, x), -x

hypercube_transform(ℓ::StandardMultivariateNormal, x) = norminvcdf.(x)


text:
"""
Module for interactive testing of highlighters in the REPL. Crayons is not a dependency and
must be installed separately for this module to work.
"""
module Console

using Crayons, Highlights

function Format.render(io::IO, ::MIME"text/ansi", tokens::Format.TokenIterator)
    for (str, id, style) in tokens
        fg = style.fg.active ? map(Int, (style.fg.r, style.fg.g, style.fg.b)) : :nothing
        bg = style.bg.active ? map(Int, (style.bg.r, style.bg.g, style.bg.b)) : :nothing
        crayon = Crayon(
            foreground = fg,
            background = bg,
            bold       = style.bold,
            italics    = style.italic,
            underline  = style.underline,
        )
        print(io, crayon, str, inv(crayon))
    end
end

test(text, lexer) = Highlights.highlight(stdout, MIME("text/ansi"), text, lexer)

end


text:
abstract type AbstractIonizationHistory{T, IT<:AbstractInterpolation{T}} end

abstract type IonizationIntegrator end
struct Peebles <: IonizationIntegrator end

struct IonizationHistory{T, IT} <: AbstractIonizationHistory{T, IT}
    Xₑ::IT
    τ::IT
    τ′::IT
    τ′′::IT
    g̃::IT
    g̃′::IT
    g̃′′::IT
end

# this Peebles history comes from Callin+06, peep the plots from examples/
# which match that paper perfectly
function IonizationHistory(integrator::Peebles, par::ACP, bg::AB) where
                           {T, ACP<:AbstractCosmoParams{T}, AB<:AbstractBackground}
    x_grid = bg.x_grid
    Xₑ_function = Bolt.saha_peebles_recombination(par)
    τ, τ′ = τ_functions(x_grid, Xₑ_function, par)
    g̃ = g̃_function(τ, τ′)

    Xₑ_ = spline(x_grid, Xₑ_function.(x_grid))
    τ_ = spline(x_grid, τ.(x_grid))
    g̃_ = spline(x_grid, g̃.(x_grid))
    IT = typeof(Xₑ_)

    # TO FIX, WHY DOES THIS CONSTRUCTOR REQUIRE {I, IT}???
    return IonizationHistory{T, IT}(
        Xₑ_,
        τ_,
        spline_∂ₓ(τ_, x_grid),
        spline_∂ₓ²(τ_, x_grid),
        g̃_,
        spline_∂ₓ(g̃_, x_grid),
        spline_∂ₓ²(g̃_, x_grid),
    )
end

# Saha Equation
# Useful for high ionization fractions.

# auxillary equations for saha_rhs
const T₀ = ustrip(natural(2.725u"K"))  # CMB temperature [K]  # TODO: make this a parameter of the ionization
n_b(a, par) = par.Ω_b * ρ_crit(par) / (m_H * a^3)
n_H(a, par) = n_b(a, par)  # ignoring helium for now
T_b(a, par) = T₀ / a
saha_rhs(a, par) = (m_e * T_b(a, par) / 2π)^(3/2) / n_H(a, par) *
    exp(-ε₀_H / T_b(a, par))  # rhs of Callin06 eq. 12

function saha_Xₑ(x, par::AbstractCosmoParams)
    rhs = saha_rhs(x2a(x), par)
    return  (√(rhs^2 + 4rhs) - rhs) / 2  # solve Xₑ² / (1-Xₑ) = RHS, it's a polynomial
end
saha_Xₑ(par) = (x -> saha_Xₑ(x, par))

# Peebles Equation
# Use this for Xₑ < 0.99, i.e. z < 1587.4

# recombination parameters for Saha/Peebles
const Λ_2s_to_1s = ustrip(natural(8.227u"s^-1"))  # rate of hydrogen double transition from 2s to 1s
const ε₀_H = ustrip(natural(13.605698u"eV"))  # ionization energy of hydrogen
const m_e = ustrip(natural(float(ElectronMass)))
const m_H = ustrip(natural(float(ProtonMass)))
const α = ustrip(natural(float(FineStructureConstant)))
const σ_T = ustrip(natural(float(ThomsonCrossSection)))

# auxillary equations
ϕ₂(T_b) = 0.448 * log(ε₀_H / T_b)
α⁽²⁾(T_b) = (64π / √(27π)) * (α^2 / m_e^2) * √(ε₀_H / T_b) * ϕ₂(T_b)
β(T_b) = α⁽²⁾(T_b) * (m_e * T_b / (2π))^(3/2) * exp(-ε₀_H / T_b)
β⁽²⁾(T_b) = β(T_b) * exp(3ε₀_H / 4T_b)
n₁ₛ(a, Xₑ, par) = (1 - Xₑ) * n_H(a, par)
Λ_α(a, Xₑ, par) = H_a(a, par) * (3ε₀_H)^3 / ((8π)^2 * n₁ₛ(a, Xₑ, par))
Cᵣ(a, Xₑ, T_b, par) = (Λ_2s_to_1s + Λ_α(a, Xₑ, par)) / (
    Λ_2s_to_1s + Λ_α(a, Xₑ, par) + β⁽²⁾(T_b))

# RHS of Callin06 eq. 13
function peebles_Xₑ′(Xₑ, par, x)
    a = exp(x)
    T_b_a = BigFloat(T_b(a, par))  # handle overflows by switching to bigfloat
    return float(Cᵣ(a, Xₑ, T_b_a, par) / H_a(a, par) * (
        β(T_b_a) * (1 - Xₑ) - n_H(a, par) * α⁽²⁾(T_b_a) * Xₑ^2))
end

"""
    peebles_Xₑ(par, Xₑ₀, x_start, x_end)

Solve the Peebles equation over a span of scale factors, and then
construct an interpolator mapping scale factor to the resulting
ionization fraction.

# Arguments:
- `par`: cosmological parameters
- ` Xₑ₀`: initial ionization fraction
- `x_start`: scale factor to begin integration
- `x_end`: scale factor to end integration

# Returns:
- `generic function`: interpolator for Xₑ(x)
"""
function peebles_Xₑ(par, Xₑ₀, x_start, x_end)
    # set up problem and integrate dXₑ/dx = peebles_Xₑ′
    prob = ODEProblem(peebles_Xₑ′, Xₑ₀, (x_start, x_end), par)
    sol = solve(prob, Tsit5(), reltol=1e-11, abstol=1e-11, dense=true)
    return sol  # ode solutions work as interpolator
end

"""
    saha_peebles_recombination(par::AbstractCosmoParams)

Utility function for generating a decent approximation to Xₑ in ΛCDM recombination,
using the Saha equation until z=1587.4 and then the Peebles equation for the rest.
"""
function saha_peebles_recombination(par::AbstractCosmoParams{T}) where {T}
    z_transition = 1587.4
    x_transition = z2x(z_transition)
    saha_z_grid = 1800:-10:z_transition
    peebles_z_grid = z_transition:-10:100
    early_time_Xₑ = Bolt.saha_Xₑ(par)
    late_time_Xₑ = Bolt.peebles_Xₑ(
        par, early_time_Xₑ(x_transition), x_transition, 0.0)
    Xₑ = x -> (x < x_transition) ? early_time_Xₑ(x) : late_time_Xₑ(x)
    return Xₑ
end

function τ_functions(x, Xₑ_function, par::AbstractCosmoParams)
    @assert x[2] > x[1]  # CONVENTION: x increasing always
    # do a reverse cumulative integrate
    rx = reverse(x)
    τ_primes = [τ′(x_, Xₑ_function, par) for x_ in x]
    τ_integrated = reverse(cumul_integrate(rx, reverse(τ_primes)))

    τ̂ = interpolate((x,),τ_integrated,Gridded(Linear()))
    τ̂′ = interpolate((x,),τ_primes,Gridded(Linear()))
    return τ̂, τ̂′
end

function τ̇(x, Xₑ_function, par)
    a = x2a(x)
    return Xₑ_function(x) * n_H(a, par) * a
end

function τ′(x, Xₑ_function, par)
    a = x2a(x)
    return -Xₑ_function(x) * n_H(a, par) * a * σ_T / ℋ_a(a, par)
end

function g̃_function(τ_x_function, τ′_x_function)
    return x -> -τ′_x_function(x) * exp(-τ_x_function(x))
end


text:
# adapted from http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_hw3a.py

import GDALUtils; const GU = GDALUtils

GU.registerdrivers() do
    GU.read("ospy/data3/sites.shp") do sitesDS
        GU.read("ospy/data3/cache_towns.shp") do townsDS
            siteslayer = GU.fetchlayer(sitesDS, 0)
            townslayer = GU.fetchlayer(townsDS, 0)

            # use an attribute filter to restrict cache_towns.shp to "Nibley"
            GU.setattributefilter(townslayer, "NAME = 'Nibley'")

            GU.fetchfeature(townslayer, 0) do nibleyFeature
                # get the Nibley geometry and buffer it by 1500
                nibleyGeom = GU.getgeom(nibleyFeature)
                bufferGeom = GU.buffer(nibleyGeom, 1500)

                # use bufferGeom as a spatial filter on sites.shp to get all sites
                # within 1500 meters of Nibley
                GU.setspatialfilter(siteslayer, bufferGeom)

                # loop through the remaining features in sites.shp and print their
                # id values
                for sitefeature in siteslayer
                    id_index = GU.getfieldindex(sitefeature, "ID")
                    println(GU.fetchfield(sitefeature, id_index))
                end
            end
        end
    end
end


text:
module Shrike

    using DataStructures
    using Distances
    using LightGraphs
    using LinearAlgebra
    using StatsBase: median
    using SparseArrays
    using ThreadsX

    include("shrike_index.jl")
    include("nearest_neighbors.jl")

    export ShrikeIndex, traverse_to_leaves, ann, knngraph, allknn

end  # module Shrike


text:
export LRNModeType, LRNMode

@compat abstract type LRNModeType end
module LRNMode
import ..LRNModeType
type AcrossChannel <: LRNModeType end
type WithinChannel <: LRNModeType end
end # module LRNMode

################################################################################
# Local Response Normalization Layer
################################################################################
@defstruct LRNLayer Layer (
  name :: AbstractString = "lrn",
  (kernel :: Int = 5, kernel > 0),
  (scale :: Number = 1, isreal(scale)),
  (shift :: Number = 1, isreal(shift)),
  (power :: Number = 0.75, isreal(power)),
  (channel_dim :: Int = -2, channel_dim != 0),
  (tops :: Vector{Symbol} = Symbol[], length(tops) == 1),
  (bottoms :: Vector{Symbol} = Symbol[], length(bottoms) == 1),
  mode :: LRNModeType = LRNMode.AcrossChannel(),
)
@characterize_layer(LRNLayer,
  can_do_bp => true
)

type LRNLayerState <: LayerState
  layer      :: LRNLayer
  blobs      :: Vector{Blob}
  blobs_diff :: Vector{Blob}

  do_split   :: LayerState
  do_square  :: LayerState
  do_pool    :: LayerState
  do_power   :: LayerState
  do_div     :: LayerState
end

function setup(backend::Backend, layer::LRNLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  for i = 1:length(inputs)
    # currently we only handle 4D-tensor
    @assert ndims(inputs[i]) == 4
  end

  split_layer = SplitLayer(no_copy=true, tops=Array{Symbol}(2), bottoms=Array{Symbol}(1))
  do_split = setup(backend, split_layer, inputs, diffs)

  square_layer = PowerLayer(power=2, tops=Array{Symbol}(1), bottoms=Array{Symbol}(1))
  do_square = setup(backend, square_layer,
      Blob[do_split.blobs[1]], Blob[do_split.blobs_diff[1]])

  pre_pad = div(layer.kernel-1,2)
  if isa(layer.mode, LRNMode.AcrossChannel)
    post_pad = layer.kernel - pre_pad - 1
    pool_layer = ChannelPoolingLayer(tops=Array{Symbol}(1), bottoms=Array{Symbol}(1),
        kernel=layer.kernel, stride=1, pad=(pre_pad,post_pad),
        pooling=Pooling.Mean(), channel_dim=layer.channel_dim)
  elseif isa(layer.mode, LRNMode.WithinChannel)
    pool_layer = PoolingLayer(tops=Array{Symbol}(1), bottoms=Array{Symbol}(1),
        kernel=(layer.kernel,layer.kernel), stride=(1,1), pad=(pre_pad,pre_pad),
        pooling=Pooling.Mean())
  end
  do_pool = setup(backend, pool_layer, do_square.blobs, do_square.blobs_diff)

  power_layer = PowerLayer(tops=Array{Symbol}(1), bottoms=Array{Symbol}(1),
      power=layer.power, scale=layer.scale, shift=layer.shift)
  do_power = setup(backend, power_layer, do_pool.blobs, do_pool.blobs_diff)

  div_layer = ElementWiseLayer(tops=Array{Symbol}(1), bottoms=Array{Symbol}(2),
      operation = ElementWiseFunctors.Divide())
  do_div = setup(backend, div_layer,
      Blob[do_split.blobs[2],do_power.blobs[1]],
      Blob[do_split.blobs_diff[2],do_power.blobs_diff[1]])

  state = LRNLayerState(layer, do_div.blobs, do_div.blobs_diff,
      do_split, do_square, do_pool, do_power, do_div)
end
function shutdown(backend::Backend, state::LRNLayerState)
  shutdown(backend, state.do_split)
  shutdown(backend, state.do_square)
  shutdown(backend, state.do_pool)
  shutdown(backend, state.do_power)
  shutdown(backend, state.do_div)
end

function forward(backend::Backend, state::LRNLayerState, inputs::Vector{Blob})
  forward(backend, state.do_split, inputs)
  forward(backend, state.do_square, Blob[state.do_split.blobs[1]])
  forward(backend, state.do_pool, state.do_square.blobs)
  forward(backend, state.do_power, state.do_pool.blobs)
  forward(backend, state.do_div, Blob[state.do_split.blobs[2],state.do_power.blobs[1]])
end

function backward(backend::Backend, state::LRNLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  if !isa(diffs[1], NullBlob)
    backward(backend, state.do_div,
        Blob[state.do_split.blobs[2],state.do_power.blobs[1]],
        Blob[state.do_split.blobs_diff[2],state.do_power.blobs_diff[1]])
    backward(backend, state.do_power, state.do_pool.blobs, state.do_pool.blobs_diff)
    backward(backend, state.do_pool, state.do_square.blobs, state.do_square.blobs_diff)
    backward(backend, state.do_square,
        Blob[state.do_split.blobs[1]],
        Blob[state.do_split.blobs_diff[1]])
    backward(backend, state.do_split, inputs, diffs)
  end
end


text:
Base.@kwdef struct Plane{T₁ <: Vec3, T₂ <: Real}
    normal::T₁
    distance::T₂
end

function get_normal(plane::Plane)
    plane.normal
end

function get_distance(plane::Plane)
    plane.distance
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule QLever_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("QLever")
JLLWrappers.@generate_main_file("QLever", UUID("6a1e3932-4b28-5fd5-abac-1dd504afb24e"))
end  # module QLever_jll


text:
import ModelingToolkit
import ModelingToolkit: build_function, ReactionSystem, substitute
import Catalyst
using StaticArrays
using Transducers
import Base.show

struct DirectAggregator{U}
    sumrate::Float64
    rates::Vector{Float64}
    # a mapping from the trajectory's recorded reaction ids to the ReactionSet indices of the jump aggregator
    update_map::U
    tspan::Tuple{Float64, Float64}
    tprev::Float64
    weight::Float64
end

struct ReactionSet
    rates::Vector{Float64}
    rstoich::Vector{Vector{Pair{Int64, Int64}}}
    nstoich::Vector{Vector{Pair{Int64, Int64}}}
end

function ReactionSet(js::ModelingToolkit.JumpSystem, p)
    parammap  = map(Pair, ModelingToolkit.parameters(js), p)
    statetoid = Dict(ModelingToolkit.value(state) => i for (i,state) in enumerate(ModelingToolkit.states(js)))

    rates = Float64[]
    rstoich_vec = Vector{Pair{Int64, Int64}}[]
    nstoich_vec = Vector{Pair{Int64, Int64}}[]

    for eq in ModelingToolkit.equations(js)
        rate = ModelingToolkit.value(ModelingToolkit.substitute(eq.scaled_rates, parammap))
        rstoich = sort!([statetoid[ModelingToolkit.value(spec)] => stoich for (spec, stoich) in eq.reactant_stoch])
        nstoich = sort!([statetoid[ModelingToolkit.value(spec)] => stoich for (spec, stoich) in eq.net_stoch])

        push!(rates, rate)
        push!(rstoich_vec, rstoich)
        push!(nstoich_vec, nstoich)
    end

    ReactionSet(rates, rstoich_vec, nstoich_vec)
end

add_weight(agg::DirectAggregator, Δweight::Float64, t::Float64) = DirectAggregator(agg.sumrate, agg.rates, agg.update_map, agg.tspan, t, agg.weight + Δweight)

function get_update_index(agg::DirectAggregator, i::Int)
    if checkbounds(Bool, agg.update_map, i)
        j = @inbounds agg.update_map[i]
        if j > 0
            j
        else
            nothing
        end
    else
        nothing
    end
end

struct TrajectoryDistribution{U}
    reactions::ReactionSet
    aggregator::DirectAggregator{U}
end

function TrajectoryDistribution(reactions, update_map = 1:num_reactions)
    num_clusters = maximum(update_map)
    agg = DirectAggregator(0.0, zeros(num_clusters), update_map, (0.0, 0.0), 0.0, 0.0)
    TrajectoryDistribution(reactions, agg)
end

myzero_fn(x) = 0.0
distribution(rn::ReactionSystem, p; update_map=1:Catalyst.numreactions(rn)) = TrajectoryDistribution(ReactionSet(convert(ModelingToolkit.JumpSystem, rn), p), update_map)

@fastmath function Distributions.logpdf(dist::TrajectoryDistribution, trajectory)::Float64
    first = iterate(trajectory)
    if first === nothing
        return 0.0
    end
    tprev = 0.0
    result = 0.0

    agg = dist.aggregator
    for (u, t, i) in trajectory
        agg = update_rates(agg, u, dist.reactions)

        dt = t - tprev
        result -= dt * agg.sumrate
        agg_i = get_update_index(agg, i)
        if agg_i !== nothing
            @inbounds result += log(agg.rates[agg_i])
        end

        tprev = t
    end

    result
end

@fastmath @inline function fold_logpdf(dist::TrajectoryDistribution, agg::DirectAggregator, (u, t, i))
    agg = update_rates(agg, u, dist.reactions)
    agg_i = get_update_index(agg, i)
    dt = t - agg.tprev
    if agg_i !== nothing
        @inbounds log_jump_prob = log(agg.rates[agg_i])
    else
        log_jump_prob = 0.0
    end
    log_surv_prob = - dt * agg.sumrate
    add_weight(agg, log_surv_prob + log_jump_prob, t)
end

function trajectory_energy(dist::TrajectoryDistribution, traj; tspan=(0.0, Inf64))
    agg = dist.aggregator
    agg = DirectAggregator(0.0, agg.rates, agg.update_map, tspan, tspan[1], 0.0)

    f = let dist=dist
        function(agg, (u, t, i))
            if t <= agg.tspan[1]
                return agg
            end
            if t > agg.tspan[2]
                agg = fold_logpdf(dist, agg, (u, agg.tspan[2], 0))
                return Transducers.reduced(agg)
            else
                return fold_logpdf(dist, agg, (u, t, i))
            end
        end
    end

    result = foldxl(f, traj; init=agg)
    result.weight
end

function cumulative_logpdf!(result::AbstractVector, dist::TrajectoryDistribution, traj, dtimes::AbstractVector)
    agg = dist.aggregator
    tspan = (first(dtimes), last(dtimes))
    result[1] = zero(eltype(result))
    agg = DirectAggregator(0.0, agg.rates, agg.update_map, tspan, tspan[1], 0.0)
    foldxl(traj; init=(agg, 1)) do (agg, k), (u, t, i)
        if t <= agg.tspan[1]
            return agg, k
        end

        t = min(t, agg.tspan[2])
        agg = update_rates(agg, u, dist.reactions)

        tprev = agg.tprev
        while k <= length(dtimes) && dtimes[k] < t
            result[k] -= (dtimes[k] - tprev) * agg.sumrate
            tprev = dtimes[k]
            k += 1
            result[k] = result[k - 1]
        end
        result[k] -= (t - tprev) * agg.sumrate

        agg_i = get_update_index(agg, i)
        if agg_i !== nothing
            @inbounds log_jump_prob = log(agg.rates[agg_i])
        else
            log_jump_prob = 0.0
        end
        result[k] += log_jump_prob

        agg = add_weight(agg, -(t - agg.tprev) * agg.sumrate + log_jump_prob, t)

        if t >= agg.tspan[2]
            return Transducers.reduced((agg, k))
        else
            return agg, k
        end
    end

    result
end

cumulative_logpdf(dist::TrajectoryDistribution, trajectory, dtimes::AbstractVector) = cumulative_logpdf!(zeros(length(dtimes)), dist, trajectory, dtimes)

@inline @fastmath function evalrxrate(speciesvec::AbstractVector, rxidx::Int64, rs::ReactionSet)
    val = 1.0
    @inbounds for specstoch in rs.rstoich[rxidx]
        specpop = speciesvec[specstoch[1]]
        val    *= specpop
        @inbounds for k = 2:specstoch[2]
            specpop -= one(specpop)
            val     *= specpop
        end
    end

    @inbounds return val * rs.rates[rxidx]
end

@fastmath function update_rates(aggregator::DirectAggregator, speciesvec::AbstractVector, reactions::ReactionSet)
    sum = 0.0
    aggregator.rates .= 0.0
    for i in eachindex(reactions.rates)
        agg_i = get_update_index(aggregator, i)
        if agg_i !== nothing
            rate = evalrxrate(speciesvec, i, reactions)
            @inbounds aggregator.rates[agg_i] += rate
            sum += rate
        end
    end
    DirectAggregator(sum, aggregator.rates, aggregator.update_map, aggregator.tspan, aggregator.tprev, aggregator.weight)
end


text:
# NOTE: Bolt's background functions are in terms of x ≡ ln(a), the log scale factor

# derived quantities (I've chosen natural units, possibly the wrong choice)
const km_s_Mpc_100 = ustrip(natural(100.0u"km/s/Mpc"))  # [eV]
const G_natural = ustrip(natural(float(NewtonianConstantOfGravitation))) #[eV^-2]
const ζ = 1.2020569 #Riemann ζ(3) for phase space integrals
#In the natural units used here, ħ=kb=c=1, G as above

H₀(par::AbstractCosmoParams) = par.h * km_s_Mpc_100
ρ_crit(par::AbstractCosmoParams) = (3 / 8π) * H₀(par)^2 / G_natural  # [eV⁴]
function Ω_Λ(par::AbstractCosmoParams)
    #Below can definitely be more streamlined, I am just making it work for now
    Tγ = (15/ π^2 *ρ_crit(par) *par.Ω_r)^(1/4)
    νfac = (90 * ζ /(11 * π^4)) * (par.Ω_r * par.h^2 / Tγ) *((par.N_ν/3)^(3/4))
    #^the factor that goes into nr approx to neutrino energy density, plus equal sharing ΔN_eff factor for single massive neutrino
    Ω_ν = par.Σm_ν*νfac/par.h^2
    return 1 - (par.Ω_r*(1+(2/3)*(7par.N_ν/8)*(4/11)^(4/3))  # dark energy density
                                         + par.Ω_b + par.Ω_m
                                         + Ω_ν) #assume massive nus are non-rel today
end

# Hubble parameter ȧ/a in Friedmann background
function H_a(a, par::AbstractCosmoParams)
    ρ_ν,_ = ρP_0(a,par) # we don't atually need pressure?
    return H₀(par) * √((par.Ω_m + par.Ω_b ) * a^(-3)
                        + ρ_ν/ρ_crit(par)
                        + par.Ω_r*(1+(2/3)*(7par.N_ν/8)*(4/11)^(4/3)) * a^(-4)
                        + Ω_Λ(par))
end
# conformal time Hubble parameter, aH
ℋ_a(a, par::AbstractCosmoParams) = a * H_a(a, par)

# functions in terms of x
H(x, par::AbstractCosmoParams) = H_a(x2a(x),par)
ℋ(x, par::AbstractCosmoParams) = ℋ_a(x2a(x), par)

# conformal time
function η(x, par::AbstractCosmoParams)
    return quadgk(a -> 1.0 / (a * ℋ_a(a, par)), 0.0, x2a(x),rtol=1e-6)[1]
end

#background FD phase space
function f0(q,par::AbstractCosmoParams)
    Tν =  (par.N_ν/3)^(1/4) *(4/11)^(1/3) * (15/ π^2 *ρ_crit(par) *par.Ω_r)^(1/4)
    gs =  2 #should be 2 for EACH neutrino family (mass eigenstate)
    return gs / (2π)^3 / ( exp(q/Tν) +1)
end

function dlnf0dlnq(q,par::AbstractCosmoParams) #this is actually only used in perts
    Tν =  (par.N_ν/3)^(1/4) * (4/11)^(1/3) * (15/ π^2 *ρ_crit(par) *par.Ω_r)^(1/4)
    return -q / Tν /(1 + exp(-q/Tν))
end

#FIXME use FastGaussQuadrature instead of slower adaptive 9as in perts)
function ρP_0(a,par::AbstractCosmoParams)
    #Background phase space energy density and pressure
    Tν =  (par.N_ν/3)^(1/4) * (4/11)^(1/3) * (15/ π^2 *ρ_crit(par) *par.Ω_r)^(1/4)
    m = par.Σm_ν
    qmin=1e-18
    qmax=1e1
    ρ = 4π * a^(-4) * quadgk(q ->  q^2 * √( q^2 + (a*m)^2 ) * f0(q,par) ,qmin, qmax,rtol=1e-6)[1]
    P = 4π/3 * a^(-4) * quadgk(q -> q^2 * q^2 /√( q^2 + (a*m)^2) * f0(q,par), qmin, qmax,rtol=1e-6)[1]
    return ρ,P#,norm
end
# now build a Background with these functions

# a background is parametrized on the scalar type T, the interpolator type IT,
# and a type for the grid GT
abstract type AbstractBackground{T, IT<:AbstractInterpolation{T,1}, GT} end

struct Background{T, IT, GT} <: AbstractBackground{T, IT, GT}
    H₀::T
    η₀::T
    ρ_crit::T
    Ω_Λ::T

    x_grid::GT
    quad_pts::Array{T,1}
    quad_wts::Array{T,1}

    ℋ::IT
    ℋ′::IT
    ℋ′′::IT
    η::IT
    η′::IT
    η′′::IT
end

function Background(par::AbstractCosmoParams{T}; x_grid=-20.0:0.01:0.0, nq=15) where T
    ℋ_  = spline([ℋ(x, par) for x in x_grid], x_grid)
    η_   = spline([η(x, par) for x in x_grid], x_grid)
    quad_pts, quad_wts =  gausslegendre( nq ) #12 should get 1e-3, 15 conservative
    return Background(
        T(H₀(par)),
        T(η(0.0, par)),
        T(ρ_crit(par)),
        T(Ω_Λ(par)),

        x_grid,
        convert(Array{T,1},quad_pts), #explicit call to convert instead of constructor for arrays
        convert(Array{T,1},quad_wts),

        ℋ_,
        spline_∂ₓ(ℋ_, x_grid),
        spline_∂ₓ²(ℋ_, x_grid),

        η_,
        spline_∂ₓ(η_, x_grid),
        spline_∂ₓ²(η_, x_grid),

    )
end


text:
# produce the alpha matrix for cell i from equation 2
# given:
# 1) a collection of ortholog ids L
# 2) a collection of ortholog ids R
# 3) single cell data in the form of a 3d matrix
# 4)
function getalpha(df,Lig,Rec,ϵ)
    Li = filteravg(df,Lig)
    Ri = filteravg(df,Rec)
    LR = Li .* Ri'
    a = similar(LR)
    for aa in eachindex(a)
        if LR[aa] <= 1/10^ϵ
            a[aa] = 0
        else
            a[aa] = exp(-1/LR[aa])
        end
    end
    round.(a,digits=ϵ)
end

function getbeta(df,uptarg,ϵ)
    ys = filtersum(df,uptarg)
    m = length(uptarg)
    for i in eachindex(ys)
        if ys[i] <= 1/10^ϵ
            ys[i] = 0
        else
            ys[i] = exp(-m/ys[i])
        end
    end
    round.(ys,digits=ϵ)
end

function getkappa(alpha,beta,ϵ)
    k = alpha .+ beta
    for i in eachindex(k)
        if k[i] <= 1/10^ϵ
            k[i] = 0
        else
            k[i] = alpha[i] / k[i]
        end
    end
    round.(k,digits=ϵ)
end

# all up targets
# alpha * beta * kappa
function getup(df,l,r,ut,ϵ)
    # make sure the genes are expressed
    alpha = getalpha(df,l,r,ϵ)
    beta = getbeta(df,ut,ϵ)
    kappa = getkappa(alpha,beta,ϵ)
    p = alpha .* kappa .* beta
    round.(p,digits=ϵ)
end

function getupavg(df,lrt,ϵ;verbose=false)
    P = zeros(length(unique(df.barcode)),length(unique(df.barcode)))
    n = size(lrt)[1]
    for t in 1:n
        l = filterexpressed(df,lrt[t,:].ligand)
        r = filterexpressed(df,lrt[t,:].receptor)
        ut = filterexpressed(df,lrt[t,:].target)
        if all((!isnothing).([l,r,ut]))
            verbose ? println("adding triple $t of $n") : nothing
            p = getup(df,l,r,ut,ϵ)
            P += round.(p,digits=ϵ)
        else
            verbose ? println("triple $t not detected") : nothing
        end
    end

    if any(P .> 0)
        return P ./ n
    else
        verbose ? println("no lrt expression") : nothing
        return P ./ n
    end
end

function filterexpressed(df,gnames)
    nms = names(df)[findall((!isnothing).([findfirst(feat->feat==g,gnames) for g in names(df)]))]
    if length(nms) > 0
        return nms
    else
        return nothing
    end
end

# get the average
function groupavg(p,groups)
    n = size(p,1)
    n_g = length(groups)
    avg1 = Array{Float64,2}(undef,n,n_g)
    for i in 1:length(groups)
        grp = groups[i]
        avg1[:,i] = mean(p[:,grp],dims=2)
    end
    avg2 = Array{Float64,2}(undef,n_g,n_g)
    for i in 1:length(groups)
        grp = groups[i]
        avg2[i,:] = mean(avg1[grp,:],dims=1)
    end
    avg2
end

# get the groups
# find the indices of each unique value in the provided columns
# return the indices and the values
function getgroups(df,cnames...)
    nms = DataFrame([String,Vector{Int64}],[:name,:inds])
    for cn in cnames
        dfcol = df[!,cn]
        unames = unique(skipmissing(dfcol))
        for nm in unames
            nmind = []
            for i in 1:length(dfcol)
                if !ismissing(dfcol[i])
                    if dfcol[i] == nm
                        push!(nmind,i)
                    end
                end
            end
            push!(nms,(name=nm,inds=nmind))
        end
    end
    nms
end


text:
import Base: -, \, +, /, *
import Base: copy
import Dolo: maxabs

function copy(m::MSM)
    dd = deepcopy(m.data)
    return MSM(dd, m.sizes)
end

function -(a::MSM)
    N = length(a.data)
    sizes = [length(e) for e in a.views]
    # c = zeros_like(b)
    data = -a.data
    return MSM(data, sizes)
end

function -(a::MSM, b::MSM)
    N = length(a.data)
    sizes = [length(e) for e in a.views]
    # c = zeros_like(b)
    data = a.data .- b.data
    return MSM(data, sizes)
end

function +(a::MSM, b::MSM)
    N = length(a.data)
    sizes = [length(e) for e in a.views]
    # c = zeros_like(b)
    data = a.data .+ b.data
    return MSM(data, sizes)
end

function \(a::MSM, b::MSM)
    N = length(a.data)
    sizes = [length(e) for e in a.views]
    # c = zeros_like(b)
    data = a.data .\ b.data
    return MSM(data, sizes)
end

function /(a::MSM, b::MSM)
    N = length(a.data)
    sizes = [length(e) for e in a.views]
    # c = zeros_like(b)
    data = a.data ./ b.data
    return MSM(data, sizes)
end

function /(a::MSM, b::Number)
    N = length(a.data)
    sizes = [length(e) for e in a.views]
    # c = zeros_like(b)
    data = a.data ./ b
    return MSM(data, sizes)
end

function *(a::MSM, b::Number)
    N = length(a.data)
    sizes = [length(e) for e in a.views]
    # c = zeros_like(b)
    data = a.data .* b
    return MSM(data, sizes)
end

norm(a::MSM) = maximum( u-> maximum(abs,u), a.data )
maxabs(a::MSM) = maximum( u-> maximum(abs,u), a.data )

function *(x::MSM, y::MSM)
    data = x.data .* y.data
    return MSM(data, x.sizes)
end

function invert!(x::MSM)
    for i=1:length(x.data)
        x.data[i] = inv(x.data[i] )
    end
end

function ldiv(a::MSM{S}, b::MSM{T}) where S where T
    r = a.data .\ b.data
    MSM(r, a.sizes)
end


text:
# For reading information & metadata about the file

function _get_fc(file::String) # convenience function for `get_duration` and `get_start_time`
    v = open(file)
    return unsafe_load(v.apFormatContext[1])
end
get_duration(fc::AVFormatContext) = fc.duration / 1e6

"""
    get_duration(file::String) -> Float64

Return the duration of the video `file` in seconds (float).
"""
get_duration(file::String) = get_duration(_get_fc(file))

get_start_time(fc::AVFormatContext) = Dates.unix2datetime(fc.start_time_realtime)

"""
    get_start_time(file::String) -> DateTime

Return the starting date & time of the video `file`. Note that if the starting date & time are missing, this function will return the Unix epoch (00:00 1st January 1970).
"""
get_start_time(file::String) = get_start_time(_get_fc(file))

get_time_duration(fc::AVFormatContext) = (get_start_time(fc), get_duration(fc))

"""
    get_time_duration(file::String) -> (DateTime, Microsecond)

Return the starting date & time as well as the duration of the video `file`. Note that if the starting date & time are missing, this function will return the Unix epoch (00:00 1st January 1970).
"""
get_time_duration(file::String) = get_time_duration(_get_fc(file))


text:
const TRAILING_BYTE_ON_NEGATIVE_NUMBERS = 0x66

@inline isnull(num::OraNumber) = num.len == 0xff
@inline isnegative(num::OraNumber) = num.ex & 0x80 == 0x00
@inline ispositive(num::OraNumber) = !isnegative(num)
@inline iszero(num::OraNumber) = sizeof_mantissa(num) == 0x00 && ispositive(num)
@inline is_negative_1e126(num::OraNumber) = sizeof_mantissa(num) == 0x00 && isnegative(num)

@inline function sizeof_mantissa(num::OraNumber) :: UInt8
    @assert !isnull(num) "null number"

    len_byte = num.len

    if isnegative(num) && num.mantissa[len_byte-UInt(1)] == TRAILING_BYTE_ON_NEGATIVE_NUMBERS
        # ignore trailing byte that contains the number 102 on negative numbers
        return len_byte - UInt(2)
    else
        return len_byte - UInt(1)
    end
end

@inline function decode_exponent_byte(num::OraNumber) :: Int8
    @assert !isnull(num) "null value"
    iszero(num) && return 0x00
    is_negative_1e126(num) && return UInt8(126)

    byte = num.ex

    if isnegative(num)
        # bits are inverted for negative numbers
        byte = ~byte
    end

    # the byte has an offset of 11000001 (bit string)
    return reinterpret(Int8, byte - 0xc1)
end

# inverse of `decode_exponent_byte`. Should not be used for zero/special number
@inline function encode_exponent_byte(ex::Int8, is_negative::Bool) :: UInt8
    local byte::UInt8 = reinterpret(UInt8, ex) + 0xc1

    if is_negative
        byte = ~byte
    end

    return byte
end

@inline function encode_mantissa_byte(byte::UInt8, is_negative::Bool) :: UInt8
    if is_negative
        return 0x65 - byte
    else
        return byte + 0x01
    end
end

# returns a mantissa byte as a 100-base number
@inline function decode_mantissa_byte(num::OraNumber, index::Integer, is_negative::Bool=isnegative(num)) :: UInt8
    @assert index <= MAX_SIZEOF_MANTISSA "Out of bounds."

    # len is exponent + mantissa bytes, so we subtract one
    if index > num.len - UInt(1)
        # unused byte
        return 0x00
    end

    if is_negative
        if index == num.len - UInt(1) && num.mantissa[index] == TRAILING_BYTE_ON_NEGATIVE_NUMBERS
            # trailing byte on negative numbers
            return 0x00
        else
            # negative numbers are subtracted from the value 101
            return 0x65 - num.mantissa[index]
        end
    else
        # positive numbers have 1 added to them
        return num.mantissa[index] - 0x01
    end
end


text:
using Core: MethodInstance
using Core.Compiler: CodeInfo, IRCode

include("compilerutils.jl")
using .CompilerUtils: ircode_from_codeinfo, methodinstance_from_expr_new_opaque_closure

struct ChildTaskInfo
    mi::MethodInstance
    inst::Expr
    position::Int
end

"""
    child_task_info(ir::IRCode) -> info::Vector{ChildTaskInfo}

Extract `MethodInstance`s from `:new_opaque_closure` instructions.
"""
function child_task_info(ir::IRCode)
    info = ChildTaskInfo[]
    for i in 1:length(ir.stmts.inst)
        inst = ir.stmts.inst[i]
        Meta.isexpr(inst, :new_opaque_closure) || continue
        mi = methodinstance_from_expr_new_opaque_closure(inst)
        push!(info, ChildTaskInfo(mi, inst, i))
    end
    return info
end

function Core.CodeInfo(info::ChildTaskInfo)
    ci = Base.uncompressed_ir(info.mi.def)::CodeInfo
    ci = Core.Compiler.copy(ci)  # not sure about the ownership
    if ci.parent === nothing
        ci.parent = info.mi
    end
    return ci
end

ci_child_tasks(ir::IRCode) = map(CodeInfo, child_task_info(ir))
ircode_child_tasks(ir::IRCode) = map(ircode_from_codeinfo, ci_child_tasks(ir))

InteractiveUtils.code_typed(info::ChildTaskInfo) = CodeInfo(info) => Any

InteractiveUtils.code_llvm(info::ChildTaskInfo; kwargs...) =
    InteractiveUtils.code_llvm(stdout, info; kwargs...)

function InteractiveUtils.code_llvm(
    io::IO,
    info::ChildTaskInfo;
    world::UInt = typemax(UInt),
    raw::Bool = false,
    dump_module::Bool = false,
    optimize::Bool = true,
    debuginfo::Symbol = :default,
)
    wrapper::Bool = false
    strip_ir_metadata::Bool = !raw
    str = InteractiveUtils._dump_function_linfo_llvm(
        info.mi::MethodInstance,
        world,
        wrapper,
        strip_ir_metadata,
        dump_module,
        optimize,
        debuginfo,
        Base.CodegenParams(),
    )
    if InteractiveUtils.highlighting[:llvm] && get(io, :color, false)
        InteractiveUtils.print_llvm(io, str)
    else
        print(io, str)
    end
end


text:
mutable struct SwitchProposeState{T}
    choices::DynamicChoiceMap
    weight::Float64
    retval::T
    SwitchProposeState{T}(choices, weight) where T = new{T}(choices, weight)
end

function process!(gen_fn::Switch{C, N, K, T},
                  index::Int,
                  args::Tuple,
                  state::SwitchProposeState{T}) where {C, N, K, T}

    (submap, weight, retval) = propose(getindex(gen_fn.branches, index), args)
    state.choices = submap
    state.weight += weight
    state.retval = retval
end

@inline process!(gen_fn::Switch{C, N, K, T}, index::C, args::Tuple, state::SwitchProposeState{T}) where {C, N, K, T} = process!(gen_fn, getindex(gen_fn.cases, index), args, state)

function propose(gen_fn::Switch{C, N, K, T},
                 args::Tuple) where {C, N, K, T}

    index = args[1]
    choices = choicemap()
    state = SwitchProposeState{T}(choices, 0.0)
    process!(gen_fn, index, args[2:end], state)
    return state.choices, state.weight, state.retval
end

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
"""
    CorrBijector <: Bijector{2}

A bijector implementation of Stan's parametrization method for Correlation matrix:
https://mc-stan.org/docs/2_23/reference-manual/correlation-matrix-transform-section.html

Basically, a unconstrained strictly upper triangular matrix `y` is transformed to
a correlation matrix by following readable but not that efficient form:

```
K = size(y, 1)
z = tanh.(y)

for j=1:K, i=1:K
    if i>j
        w[i,j] = 0
    elseif 1==i==j
        w[i,j] = 1
    elseif 1<i==j
        w[i,j] = prod(sqrt(1 .- z[1:i-1, j].^2))
    elseif 1==i<j
        w[i,j] = z[i,j]
    elseif 1<i<j
        w[i,j] = z[i,j] * prod(sqrt(1 .- z[1:i-1, j].^2))
    end
end
```

It is easy to see that every column is a unit vector, for example:

```
w3' w3 ==
w[1,3]^2 + w[2,3]^2 + w[3,3]^2 ==
z[1,3]^2 + (z[2,3] * sqrt(1 - z[1,3]^2))^2 + (sqrt(1-z[1,3]^2) * sqrt(1-z[2,3]^2))^2 ==
z[1,3]^2 + z[2,3]^2 * (1-z[1,3]^2) + (1-z[1,3]^2) * (1-z[2,3]^2) ==
z[1,3]^2 + z[2,3]^2 - z[2,3]^2 * z[1,3]^2 + 1 -z[1,3]^2 - z[2,3]^2 + z[1,3]^2 * z[2,3]^2 ==
1
```

And diagonal elements are positive, so `w` is a cholesky factor for a positive matrix.

```
x = w' * w
```

Consider block matrix representation for `x`

```
x = [w1'; w2'; ... wn'] * [w1 w2 ... wn] ==
[w1'w1 w1'w2 ... w1'wn;
 w2'w1 w2'w2 ... w2'wn;
 ...
]
```

The diagonal elements are given by `wk'wk = 1`, thus `x` is a correlation matrix.

Every step is invertible, so this is a bijection(bijector).

Note: The implementation doesn't follow their "manageable expression" directly,
because their equation seems wrong (7/30/2020). Insteadly it follows definition
above the "manageable expression" directly, which is also described in above doc.
"""
struct CorrBijector <: Bijector{2} end

function (b::CorrBijector)(x::AbstractMatrix{<:Real})
    w = cholesky(x).U  # keep LowerTriangular until here can avoid some computation
    r = _link_chol_lkj(w)
    return r + zero(x)
    # This dense format itself is required by a test, though I can't get the point.
    # https://github.com/TuringLang/Bijectors.jl/blob/b0aaa98f90958a167a0b86c8e8eca9b95502c42d/test/transform.jl#L67
end

(b::CorrBijector)(X::AbstractArray{<:AbstractMatrix{<:Real}}) = map(b, X)

function (ib::Inverse{<:CorrBijector})(y::AbstractMatrix{<:Real})
    w = _inv_link_chol_lkj(y)
    return w' * w
end
(ib::Inverse{<:CorrBijector})(Y::AbstractArray{<:AbstractMatrix{<:Real}}) = map(ib, Y)

function logabsdetjac(::Inverse{CorrBijector}, y::AbstractMatrix{<:Real})
    K = LinearAlgebra.checksquare(y)

    result = float(zero(eltype(y)))
    for j in 2:K, i in 1:(j - 1)
        @inbounds abs_y_i_j = abs(y[i, j])
        result += (K - i + 1) * (
            IrrationalConstants.logtwo - (abs_y_i_j + LogExpFunctions.log1pexp(-2 * abs_y_i_j))
        )
    end

    return result
end
function logabsdetjac(b::CorrBijector, X::AbstractMatrix{<:Real})
    #=
    It may be more efficient if we can use un-contraint value to prevent call of b
    It's recommended to directly call
    `logabsdetjac(::Inverse{CorrBijector}, y::AbstractMatrix{<:Real})`
    if possible.
    =#
    return -logabsdetjac(inv(b), (b(X)))
end
function logabsdetjac(b::CorrBijector, X::AbstractArray{<:AbstractMatrix{<:Real}})
    return mapvcat(X) do x
        logabsdetjac(b, x)
    end
end

function _inv_link_chol_lkj(y)
    K = LinearAlgebra.checksquare(y)

    w = similar(y)

    @inbounds for j in 1:K
        w[1, j] = 1
        for i in 2:j
            z = tanh(y[i-1, j])
            tmp = w[i-1, j]
            w[i-1, j] = z * tmp
            w[i, j] = tmp * sqrt(1 - z^2)
        end
        for i in (j+1):K
            w[i, j] = 0
        end
    end

    return w
end

"""
    function _link_chol_lkj(w)

Link function for cholesky factor.

An alternative and maybe more efficient implementation was considered:

```
for i=2:K, j=(i+1):K
    z[i, j] = (w[i, j] / w[i-1, j]) * (z[i-1, j] / sqrt(1 - z[i-1, j]^2))
end
```

But this implementation will not work when w[i-1, j] = 0.
Though it is a zero measure set, unit matrix initialization will not work.

For equivelence, following explanations is given by @torfjelde:

For `(i, j)` in the loop below, we define

    z₍ᵢ₋₁, ⱼ₎ = w₍ᵢ₋₁,ⱼ₎ * ∏ₖ₌₁ⁱ⁻² (1 / √(1 - z₍ₖ,ⱼ₎²))

and so

    z₍ᵢ,ⱼ₎ = w₍ᵢ,ⱼ₎ * ∏ₖ₌₁ⁱ⁻¹ (1 / √(1 - z₍ₖ,ⱼ₎²))
            = (w₍ᵢ,ⱼ₎ * / √(1 - z₍ᵢ₋₁,ⱼ₎²)) * (∏ₖ₌₁ⁱ⁻² 1 / √(1 - z₍ₖ,ⱼ₎²))
            = (w₍ᵢ,ⱼ₎ * / √(1 - z₍ᵢ₋₁,ⱼ₎²)) * (w₍ᵢ₋₁,ⱼ₎ * ∏ₖ₌₁ⁱ⁻² 1 / √(1 - z₍ₖ,ⱼ₎²)) / w₍ᵢ₋₁,ⱼ₎
            = (w₍ᵢ,ⱼ₎ * / √(1 - z₍ᵢ₋₁,ⱼ₎²)) * (z₍ᵢ₋₁,ⱼ₎ / w₍ᵢ₋₁,ⱼ₎)
            = (w₍ᵢ,ⱼ₎ / w₍ᵢ₋₁,ⱼ₎) * (z₍ᵢ₋₁,ⱼ₎ / √(1 - z₍ᵢ₋₁,ⱼ₎²))

which is the above implementation.
"""
function _link_chol_lkj(w)
    K = LinearAlgebra.checksquare(w)

    z = similar(w) # z is also UpperTriangular.
    # Some zero filling can be avoided. Though diagnoal is still needed to be filled with zero.

    # This block can't be integrated with loop below, because w[1,1] != 0.
    @inbounds z[1, 1] = 0

    @inbounds for j=2:K
        z[1, j] = atanh(w[1, j])
        tmp = sqrt(1 - w[1, j]^2)
        for i in 2:(j - 1)
            p = w[i, j] / tmp
            tmp *= sqrt(1 - p^2)
            z[i, j] = atanh(p)
        end
        z[j, j] = 0
    end

    return z
end


text:
struct MLS{AT} <: IRMeasurement
    sig::AT
    reps::Int
    gain::Float64
    prepad::Int
end

MLS(sig::AT, reps, gain, prepad) where AT = MLS{AT}(sig, reps, gain, prepad)

prepadding(m::MLS) = m.prepad

"""
    mls(L, reps=2; options...)

Generate a maximum-length sequence (MLS), also known as Schroeder's method.
Actual length is `L-1`. The sequence will be repeated `reps` times, which should
be at least 2, as only the last `reps-1` repetitions are used.

## Options:
$optiondoc_gain
$optiondoc_prepad
"""
function mls(len, reps=2; gain=mls_gain, prepad=default_prepad, samplerate=1)
    N = Int(log2(nextpow(2, durnorm(len, samplerate))))
    1 <= N <= length(mlspolys) || throw(ArgumentError("N must be positive and <= $(length(mlspolys))"))
    poly = mlspolys[N]
    # where we'll keep the state used to compute the next value. The
    # rightmost bit is the most recent value. Nothing super special about this
    # initialization value, as long as there's at least 1 nonzero bit in the
    # filter span
    register = 0x0000000000000001
    # bit mask for the filter taps
    taps = reduce(|, (0x0000000000000001 << (N-1-shift) for shift in poly))
    out = zeros(2^N-1)
    # initialize our output array to match the register above - with only the
    # most recent sample as 1.
    out[1:N-1] .= -1
    out[N] = 1
    for i in N+1:length(out)
        nextval = UInt64(count_ones(register & taps)) & 0x0000000000000001
        register = (register << 1) | nextval
        out[i] = float(nextval) * 2 - 1
    end

    MLS(out, reps, gain, durnorm(prepad, samplerate))
end

stimulus(m::MLS) = [zeros(m.prepad); repeat(m.sig, m.reps)*m.gain]

"""
Perform circular cross-correlation, assuming `u` and `v` are equal length, but
allowing for the possibility that their length is not amenable to FFT
processing.

Uses algorithm from:

Daigle, J. N., & Xiang, N. (2006). A specialized fast cross-correlation for
acoustical measurements using coded sequences. The Journal of the Acoustical
Society of America, 119(1), 330–335. https://doi.org/10.1121/1.2141236

"""
function circxcorr(u, v)
    L = length(u)
    @assert L == length(v)

    if L == nextfastfft(L)
        _circxcorr(u,v)
    else
        N = nextfastfft(2L)
        u = DSP._zeropad(u, N)
        v = DSP._zeropad(v, N)
        xc = _circxcorr(u,v)
        xc[1:L] + xc[end-L+1:end]
    end
end

"""
Perform circular cross-correlation, assuming `u` and `v` are vectors of equal
length and that length is amenable to FFT processing.
"""
_circxcorr(u, v) = irfft(rfft(u) .* conj.(rfft(v)), length(u))

function _analyze(m::MLS, response::AbstractArray)
    if m.reps < 2
        throw(ArgumentError("Can't perform MLS analysis with only one repetition"))
    end
    L = length(m.sig)
    avgresp = zeros(L, Base.tail(size(response))...)
    for rep in 1:m.reps-1
        avgresp .+= timeslice(response, (1:L) .+ (m.prepad + rep*L))
    end
    avgresp ./= (m.reps-1)
    mapslices(avgresp, dims=1) do v
        # compensate for the amplitude drop
        circxcorr(v, m.sig) ./ L / m.gain
    end
end

# these polynomial coefficients represent the tap indices of the filter, where
# zero is the oldest value (delay of N)
# they are from:
# Stahnke, Wayne 1973. Primitive Binary Polynomials.
# Mathematics of Computation 27(124): 977–980.
# the paper has higher orders, but this seems like more than enough for
# reasonable audio purposes. Feel free to submit a PR with more if they are
# useful to you.

const mlspolys = [
    [0],            # 01
    [1,  0],        # 02
    [1,  0],        # 03
    [1,  0],        # 04
    [2,  0],        # 05
    [1,  0],        # 06
    [1,  0],        # 07
    [6,  5,  1, 0], # 08
    [4,  0],        # 09
    [3,  0],        # 10
    [2,  0],        # 11
    [7,  4,  3, 0], # 12
    [4,  3,  1, 0], # 13
    [12, 11, 1, 0], # 14
    [1,  0],        # 15
    [5,  3,  2, 0], # 16
    [3,  0],        # 17
    [7,  0],        # 18
    [6,  5,  1, 0], # 19
    [3,  0],        # 20
    [2,  0],        # 21
    [1,  0],        # 22
    [5,  0],        # 23
    [4,  3,  1, 0], # 24
    [3,  0],        # 25
    [8,  7,  1, 0], # 26
    [8,  7,  1, 0], # 27
    [3,  0],        # 28
    [2,  0],        # 29
    [16, 15, 1, 0]  # 30
]


text:
using Revise
using ZigZagBoomerang
using DataStructures
using LinearAlgebra
using Random
using SparseArrays
using Test
using FileIO
using CairoMakie, AbstractPlotting

include("./../gridlaplace.jl")

# Define precision operator of a Gaussian random field (sparse matrix operating on `vec`s of `n*n` matrices)
n = 100
Γ = 0.01I + gridlaplacian(Float64, n, n)
mat(x) = reshape(x, (n, n)) # vector to matrix

# Γ is very sparse
@show nnz(Γ)/length(Γ) # 0.000496

# Corresponding Gaussian potential
# ϕ(x', Γ) = 0.5*x'*Γ*x  # not needed

# Define ∇ϕ(x, i, Γ) giving the partial derivative of ϕ(x) with respect to x[i]
∇ϕ(x, i, Γ) = ZigZagBoomerang.idot(Γ, i, x) # more efficient that dot(Γ[:, i], x)

# Random initial values
t0 = 0.0
x0 = randn(n*n)
θ0 = rand([-1.0,1.0], n*n)

# Rejection bounds
c = [norm(Γ[:, i], 2) for i in 1:n*n]

# Define ZigZag
Z = ZigZag(Γ, x0*0)
# or try the FactBoomerang
#Z = FactBoomerang(Γ, x0*0, 0.1)

# Run sparse ZigZag for T time units and collect trajectory
T = 20.0
@time trace, (tT, xT, θT), (acc, num) = spdmp(∇ϕ, t0, x0, θ0, T, c, Z, Γ)
su = false
adapt = true
κ = [1.0 for i in eachindex(x0)]
trace, (t, x, θ), (acc, num), c = @time sspdmp(∇ϕ, t0, x0, θ0, T, c, Z, κ, Γ;
                                                strong_upperbounds = su ,
                                                adapt = adapt)

@time traj = collect(discretize(trace, 0.1))

# Prepare surface plot
M = Node(mat(traj[1].second))
p1 = surface(M; shading=false, show_axis=false, colormap = :deep)
scale!(p1, 1.0, 1.0, 7.5)
zlims!(p1, -2.0, 2.0)

# Movie frames
mkpath(joinpath(@__DIR__, "output"))
for i in Iterators.take(eachindex(traj), 100)
    M[] = mat(traj[i].second)
    FileIO.save(joinpath(@__DIR__, "output", "surf$i.png"), p1)
end

# Make video
dir = joinpath(@__DIR__, "output")
run(`ffmpeg -y -r 40 -f image2 -i $dir/surf%d.png -vcodec libx264 -crf 25  -pix_fmt yuv420p $(typeof(Z).name)field.mp4`)


text:
using Downloads

# WARNING: THIS FILE IS WORK-IN-PROGRESS

#--------------------------------------------------------------------
# Metadata info: dimensions
#--------------------------------------------------------------------
function listdimensions(apiurl::String, dataflow::String)

    io = IOBuffer()
    resp = Downloads.download(apiurl * "data/" * dataflow * "?detail=serieskeysonly&lastNObservations=1&format=jsondata", io) |> take!

    ds = JSON3.read(resp).structure.dimensions.series

    dim = OrderedDict()

    for dsᵢ in ds
        dim[dsᵢ.name] = NamedTuple(Symbol(j.id) => j.name for j in dsᵢ.values)
    end

    dim
end

#--------------------------------------------------------------------
# Basic Data download
#--------------------------------------------------------------------

"""Generates a SDMX API compliant url
#kwargs
    filter::Union{nothing,NamedTuple},
    updatedAfter::DateTime,
    firstNObservations::Int,
    lastNObservations::Int,
    dimensionAtObservation,
    attributes = "dsd",
    measures = "all",
    includeHistory = false
"""
function generateurl(;context = "*", agencyID = "*", resourceID = "*", version = "*", key::String="*", kwargs...) # TODO: substitute query args by kwargs...
    baseurl = join(["https://ws-entry-point/data/",context,agencyID,resourceID,version,key],"/")
    query = "/?"*join([String(k)*"="*kwargs[k] for k in keys(kwargs)],"&")
end

""" Generates key part of url from a dictionary"""
function generatekey(dims)
    key = ""
    for v in dims
        key = key*join[v,"+"]*"." # FIXME = sobraría un punto??
    end
end

""" Fetch data and creates a SDMX.Datatable"""
function getseries(url)
    io = IOBuffer()
    return Downloads.download(url, io) |> take! |> SDMX.read(alldims = false)
end

max_stars_repo_name: lsaenzt/SDMX

max_issues_repo_name: lsaenzt/SDMX

max_forks_repo_name: lsaenzt/SDMX


text:
tagrng(ω, rng::AbstractRNG) = tag(ω, (rng = rng,))
rng(ω::Ω) = Random.GLOBAL_RNG
rng(ω::TaggedΩ) = haskey(ω.tags, :rng) ? ω.tags.rng : Random.GLOBAL_RNG

function combinetags(::Type{Val{:rng}}, a, b)
  # @assert false
  a
end


text:
# Growth screens seem to be highly susceptible to the distribution of
# genes that promote growth versus that hurt it. This experiment tests
# the behavior of growth screens under several different library setups

function growth_sensitivity_library(filepath; debug=false, quiet=false)

    if !debug

        inc_space = linspace(0.02, 0.1, 5)
        cent_space = linspace(0.1, 0.55, 5)

        libs = Library[]
        for (inc, cent, crisprtype) in IterTools.product(inc_space, cent_space, [CRISPRi(), CRISPRn()])

            max_phenotype_dists = Dict{Symbol, Tuple{Float64, Sampleable}}(
                :inactive => (1-inc-0.1-0.05, Delta(0.0)),
                :negcontrol => (0.05, Delta(0.0)),
                :increasing => (inc, TruncatedNormal(cent, 0.1, 0.025, 1)),
                :decreasing => (0.1, TruncatedNormal(-0.55, 0.2, -1, -0.1))
            )
            push!(libs, Library(max_phenotype_dists, crisprtype))
        end

        parameters = Dict{Symbol, Vector}(
                :representation => [100],
                :bottleneck_representation => [100],
                :seq_depth => [100],
                :num_bottlenecks => collect(1:20)
            )
        num_runs = 25
    else
        libs = Library[Library(CRISPRi()), Library(CRISPRn())]

        parameters = Dict{Symbol, Vector}(
                :representation => [100],
                :bottleneck_representation => [100],
                :seq_depth => [100],
                :num_bottlenecks => [7, 20]
            )
        num_runs = 1
    end

    runs = grouped_param_space(GrowthScreen(), parameters, libs, [:num_bottlenecks], num_runs);

    function compute_auprc(barcodes, genes)
        i = auprc(genes[:pvalmeanprod_bin2_div_bin1], genes[:class], Set([:increasing]))[1]
        d = auprc(genes[:pvalmeanprod_bin2_div_bin1], genes[:class], Set([:decreasing]), rev=false)[1]
        (i, d)
    end

    before = time()
    results = pmap(args -> run_exp(args[1], args[2], compute_auprc;
                   run_idx=args[3], flatten_func=flatten_both), runs)
    (!quiet) && println("$(time() - before) seconds")

    data = DataFrame(permutedims(hcat(results...), [2, 1]))
    hierarchy = reshape([:inc; :dec], (2, 1))
    new_names = [[:measure, :score]...; fieldnames(GrowthScreen)...; array_names(Library)...; :run_idx]

    data = construct_hierarchical_label(hierarchy, data, new_names)

    CSV.write(filepath, data)
end


text:
module SerializeTables

export serialize_table_rows, deserialize_table_rows

using ArgCheck: @argcheck
using DocStringExtensions: FUNCTIONNAME, SIGNATURES
using Serialization: serialize, deserialize, Serializer, Serialization
using Tables
using TranscodingStreams: TranscodingStream, Codec

# preamble

const SIGNATURE = "SerializedTables.jl table data"

struct Version{V}
    function Version{V}() where V
        @argcheck V isa Int64
        new{V}()
    end
end

const VERSION = Version{1}()

function write_preamble(io, ::Version{V}) where V
    write(io, SIGNATURE)
    write(io, V)
end

function checked_preamble_version(io)
    for c in Vector{UInt8}(SIGNATURE)
        @argcheck read(io, UInt8) == c "invalid signature"
    end
    Version{read(io, Int64)}()
end

# generic serialize and deserialize

"""
$(FUNCTIONNAME)([codec], filename, table; kwargs...)

Write representation of `table` (which conforms to the `Tables` interface and provides row
access) to `filename`, by rows, using Julia's serializer from the standard library.

When `codec` is specified, it is used with `TranscodingStreams` for compression.
"""
serialize_table_rows(filename::AbstractString, table; kwargs...) =
    open(io -> serialize_table_rows(io, table; kwargs...), filename, "w")

function serialize_table_rows(codec::Codec, filename::AbstractString, table; kwargs...)
    open(filename, "w") do io
        compressed_io = TranscodingStream(codec, io)
        serialize_table_rows(compressed_io, table; kwargs...)
        close(compressed_io)
    end
end

function serialize_table_rows(io::IO, table; kwargs...)
    write_preamble(io, VERSION)
    _serialize_table_rows(VERSION, io, table; kwargs...)
end

"""
$(SIGNATURES)

Return an iterable for the table rows in filename, that supports the rows interface of `Tables`.

When `codec` is specified, it is used with `TranscodingStreams` for decompression.
"""
deserialize_table_rows(filename::AbstractString) = # FIXME finalizer?
    deserialize_table_rows(open(filename, "r"))

deserialize_table_rows(codec::Codec, filename::AbstractString) =
    deserialize_table_rows(TranscodingStream(codec, open(filename, "r")))

function deserialize_table_rows(io::IO)
    V = checked_preamble_version(io)
    # FIXME finalizer?
    _deserialize_table_rows(V, io)
end

# Version 1

"""
$(SIGNATURES)

# Layout

A serialized stream of

1. the schema (of rows)
2. contents of rows, directly as named tuples.
"""
function _serialize_table_rows(V::Version{1}, io::IO, table)
    @argcheck Tables.istable(table)
    s = Serializer(io)
    Serialization.writeheader(s)
    rows = Tables.rows(table)
    serialize(s, Tables.schema(rows))
    for row in rows
        serialize(s, row)
    end
end

mutable struct SerializedRows{S,I,T}
    schema::S
    io::I
    serializer::T
    atstart::Bool
end

_close(sr::SerializedRows) = close(sr.io)

function _deserialize_table_rows(::Version{1}, io::IO)
    serializer = Serializer(io)
    schema = deserialize(serializer)
    SerializedRows(schema, io, serializer, true)
end

Tables.istable(::Type{<:SerializedRows}) = true

Tables.rowaccess(::Type{<:SerializedRows}) = true

Tables.schema(sr::SerializedRows) = sr.schema

Tables.rows(sr::SerializedRows) = sr

Base.IteratorSize(::Type{<:SerializedRows}) = Base.SizeUnknown()

function Base.iterate(sr::SerializedRows, atstart = true)
    if atstart
        @argcheck sr.atstart "Iteration in progress, open file again."
        sr.atstart = false
    end
    eof(sr.io) && (_close(sr); return nothing)
    deserialize(sr.serializer), false
end

end # module


text:
module BuildBlink
# put into module b/c some globals are defined in install.jl
using RelocatableFolders

include(joinpath(@__DIR__, "../src/AtomShell/install.jl"))

function get_installed_version()
    _path = Sys.isapple() ?
        joinpath(deps_folder(), "version") :
        joinpath(deps_folder(), "atom", "version")
    strip(read(_path, String), 'v')
end

if !isinstalled() || version != get_installed_version()
    install()
end

end


text:
# open the following 3 lines in case you need to use stand-alone substitute generation
using Knet, JLD
include("model2.jl")
include("preprocess.jl")

function gensub(parameters, states, sequence; lp=false)
    result = Array(Any, length(sequence))
    hlayers = length(states) / 4
    hlayers = convert(Int, hlayers)

    # forward lstm
    fhiddens = Array(Any, length(sequence))
    t1 = states[1:2*hlayers]
    for i=1:length(sequence)-1
        input = oftype(parameters[1], sequence[i])
        x = input * parameters[end-1]
        fhiddens[i+1] = forward(parameters[1:2*hlayers], t1, x)
    end
    fhiddens[1] = oftype(fhiddens[2], zeros(size(fhiddens[2])))

    # backward lstm
    bhiddens = Array(Any, length(sequence))
    t2 = states[2*hlayers+1:4*hlayers]
    for i=length(sequence):-1:2
        input = oftype(parameters[1], sequence[i])
        x = input * parameters[end]
        bhiddens[i-1] = forward(parameters[2*hlayers+1:4*hlayers], t2, x)
    end
    bhiddens[end] = oftype(bhiddens[end-1], zeros(size(bhiddens[2])))

    # merge layer
    for i=1:length(fhiddens)
        ypred = hcat(fhiddens[i], bhiddens[i]) * parameters[end-2] .+ parameters[end-3]
        ynorm = logp(ypred, 2)
        result[i] = (lp ? ynorm : exp(ynorm))
    end
    return result
end

function printsub(parameters, states, sequence, index_to_word; lp=true)
    pbvec = gensub(parameters, states, sequence; lp=lp)
    for i=1:length(sequence)
        subs = zip(index_to_word, pbvec[i])
        index = find(x->x==true, sequence[i])
        word = index_to_word[index][1]
        print(word, "\t")
        for (sub, lprob) in subs
            print(sub," ", lprob,"\t")
        end
        println()
    end
    # return pbvec
end

function test()
    x = load("initialmodel.jld")
    param = x["model"]
    state = x["initstates"]
    tdata = Data("readme_data.txt";batchsize=1)
    sequences = Any[]
    for item in tdata; push!(sequences, item);end;
    printsub(param, state, sequences[1], tdata.index_to_word)
end


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

name = "OpenCV"
version = v"4.5.2"
julia_version = v"1.6.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/opencv/opencv.git", "39d25787f16c4dd6435b9fe0a8253394ac51e7fb"),
    GitSource("https://github.com/opencv/opencv_contrib.git", "f5d7f6712d4ff229ba4f45cf79dfd11c557d56fd"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
# Apply patch for BB specific CMake changes
cd opencv_contrib
git apply ../patches/opencv-julia.patch
cd ..
mkdir build && cd build
export USE_QT="ON"
if [[ "${target}" == *-apple-* ]]; then
    # We want to use OpenBLAS over Accelerate framework...
    export OpenBLAS_HOME=${prefix}
    export CXXFLAGS=""
    # ...but we also need to rename quite a few symbols
    for symbol in sgemm dgemm cgemm zgemm; do
        # Rename CBLAS symbols for ILP64
        CXXFLAGS="${CXXFLAGS} -Dcblas_${symbol}=cblas_${symbol}64_"
    done
    for symbol in sgesv_ sposv_ spotrf_ sgesdd_ sgeqrf_ sgels_ dgeqrf_ dgesdd_ sgetrf_ dgesv_ dposv_ dgels_ dgetrf_ dpotrf_ dgeev_; do
        # Rename LAPACK symbols for ILP64
        CXXFLAGS="${CXXFLAGS} -D${symbol}=${symbol}64_"
    done
    # Apply patch to help CMake find our 64-bit OpenBLAS
    atomic_patch -p1 -d../opencv ../patches/find-openblas64.patch
    # Disable QT
    export USE_QT="OFF"
elif [[ "${target}" == *-w64-* ]]; then
    # Needed for mingw compilation of big files
    export CXXFLAGS="-Wa,-mbig-obj"
    export USE_QT="OFF"
fi
cmake -DCMAKE_FIND_ROOT_PATH=${prefix} \
      -DJulia_PREFIX=${prefix} \
      -DWITH_JULIA=ON \
      -DCMAKE_INSTALL_PREFIX=${prefix} \
      -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
      -DCMAKE_BUILD_TYPE=Release \
      -DWITH_QT=${USE_QT} \
      -DOPENCV_EXTRA_MODULES_PATH=../opencv_contrib/modules \
      -DBUILD_LIST=core,imgproc,imgcodecs,highgui,videoio,dnn,features2d,objdetect,calib3d,video,gapi,stitching,julia \
      ../opencv/
if [[ "${target}" == *-freebsd* ]]; then
    atomic_patch -p1 ../patches/freebsd-malloc-h.patch
fi
make -j${nproc}
make install

# Install also libopencv_julia
cp lib/libopencv_julia.* ${libdir}/.
# Move julia bindings to the prefix
cp -R OpenCV ${prefix}

install_license ../opencv/{LICENSE,COPYRIGHT}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
# We don't have Qt5 for Musl platforms
filter!(p -> libc(p) != "musl", platforms)
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct(["libopencv_calib3d", "libopencv_calib3d452"], :libopencv_calib3d),
    LibraryProduct(["libopencv_objdetect", "libopencv_objdetect452"], :libopencv_objdetect),
    LibraryProduct(["libopencv_core", "libopencv_core452"], :libopencv_core),
    LibraryProduct(["libopencv_dnn", "libopencv_dnn452"], :libopencv_dnn),
    LibraryProduct(["libopencv_imgcodecs", "libopencv_imgcodecs452"], :libopencv_imgcodecs),
    LibraryProduct(["libopencv_highgui", "libopencv_highgui452"], :libopencv_highgui),
    LibraryProduct(["libopencv_flann", "libopencv_flann452"], :libopencv_flann),
    LibraryProduct(["libopencv_gapi", "libopencv_gapi452"], :libopencv_gapi),
    LibraryProduct(["libopencv_imgproc", "libopencv_imgproc452"], :libopencv_imgproc),
    LibraryProduct(["libopencv_features2d", "libopencv_features2d452"], :libopencv_features2d),
    LibraryProduct(["libopencv_stitching", "libopencv_stitching452"], :libopencv_stitching),
    LibraryProduct(["libopencv_video", "libopencv_video452"], :libopencv_video),
    LibraryProduct(["libopencv_videoio", "libopencv_videoio452"], :libopencv_videoio),
    LibraryProduct("libopencv_julia", :libopencv_julia)#,
    # FileProduct("OpenCV.jl.tar", :OpenCV_jl)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Qt5Base_jll", uuid="ea2cea3b-5b76-57ae-a6ef-0a8af62496e1"))
    Dependency(PackageSpec(name="Libglvnd_jll", uuid="7e76a0d4-f3c7-5321-8279-8d96eeed0f29"))
    BuildDependency(PackageSpec(name="libjulia_jll", version=julia_version))
    Dependency(PackageSpec(name="libcxxwrap_julia_jll"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version = v"7")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
    Structure{GeometryType}

Defines a physical object with a shape and medium.

# Fields
- `geometry`: Defines the geometry of the structure.
- `medium`: Defines the material of the structure.
"""
@with_kw struct Structure{GeometryType}
    geometry::GeometryType
    medium::Medium
end

medium(s::Structure) = s.medium
geometry(s::Structure) = s.geometry

epsilon(s::Structure, freq = Inf) = epsilon(medium(s), freq)
mu(s::Structure, freq = Inf) = mu(medium(s), freq)


text:
@enum OrderingOption lexigraphic edge soporder
"""
     findSplitOrdering(model::SequencingModel,type::OrderingOption=edge)

Find a split ordering for the nodes

# Arguments
- `model::SequencingModel`: the sequencing problem to be evaluated
- `type::SOPOrderingOption`: The desired method to use for finding a split ordering, default is edge
"""
function findSplitOrdering(model::SequencingModel,type::Union{OrderingOption,Nothing}=edge)
    if type == edge || isnothing(type)
        return orderEdgeWeight(model::SequencingModel)
    elseif type == lexigraphic
        return orderLexigraphic(model)
    elseif type == soporder
        return orderSOPPrecedenceToEdge(model)
    end
end

"""
     orderLexigraphic(model::SequencingModel)

Find a variable ordering for the model based on the order of the jobs in the problem definition
Return the variable ordering

# Arguments
- `model::SequencingModel`: the sequencing problem to be evaluated
"""
function orderLexigraphic(model::SequencingModel)
    objType = valtype(model.objective.f.components)
    assignmentOrdering = Vector{objType}()
    for i in 1:length(model)
        push!(assignmentOrdering, i)
    end
    return assignmentOrdering
end

"""
     orderEdgeWeight(model::SequencingModel)

Find a variable ordering for the model based on the distances between nodes
Return the variable ordering

# Arguments
- `model::SequencingModel`: the sequencing problem to be evaluated
"""
function orderEdgeWeight(model::SequencingModel)
    #get objective values of the arcs
    objType = valtype(model.objective.f.components)
    objectiveValues = Dict{objType, Vector{objType}}()
    for (key,value) in model.objective.f.components
        for i in 1:2
            if key[i] in keys(objectiveValues)
                push!(objectiveValues[key[i]], value)
            else
                objectiveValues[key[i]] = [value]
            end
        end
    end

    #get (and sort by) averages of objective values for each node
    averages = Dict{objType, Float64}()
    for (key, value) in objectiveValues
        averages[key] = mean(value)
    end
    return sort(collect(keys(averages)), by=x->averages[x], rev=true)
end


text:
using jubart
using Test

@testset "jubart.jl" begin
    # Testing the friedman_simulation_function
    @test isa(sim_friedman_simple(10,5,0.01),Dict)
    @test isa(sim_friedman_simple(10,5,0.01)["x"],Array{Float64,2})
    @test isa(sim_friedman_simple(10,5,0.01)["y"],Array{Float64,1})
    @test isa(sim_friedman_simple(10,5,0.01)["true_mean"],Array{Float64,1})

    # Testing the Julia model
    data_test=sim_friedman_simple(100,5,0.001)
    min_node_size_test=10
    number_tree=5
    model_test=juBart(data_test["x"],data_test["y"],number_tree,min_node_size_test)

    # Testing just the output
    @test isa(model_test,Dict)

    # Test for check the min node size
    node_size_count=0 # Counter of min_node_size

    for mh_tree in model_test["trees"]
        for m in 1:number_tree
            tree_nodes=mh_tree[m].nodes
            for run_tree_nodes in tree_nodes
                if(length(run_tree_nodes.observations_index)<min_node_size_test)
                    node_size_bol+=1
                end
            end
        end
    end

    if(node_size_count==0) node_size_bol=true else node_size_bol=false end

    # Veryfing the node boolean
    @test node_size_bol

    # Veryfing NA in test model
    @test !any(isnan,model_test["predictions"])

    # Testing the predict function
    prediction_test_model_all=predict_juBart(model_test,rand(100,10),"all") # Creating the all matrix
    prediction_test_model_mean=predict_juBart(model_test,rand(100,10)) # Creating the prediction vector

    # Verfying size
    @test isa(prediction_test_model_all,Array{Union{Nothing,Float64},2})
    @test isa(prediction_test_model_mean,Array{Float64,2})

    # Veryfing NA
    @test !any(isnan,prediction_test_model_all)
    @test !any(isnan,prediction_test_model_mean)

end


text:
# Deprecated
# use various test_*.jl files
using Test

using Dates
using Logging

# see https://discourse.julialang.org/t/writing-tests-in-vs-code-workflow-autocomplete-and-tooltips/57488
# see https://github.com/julia-vscode/julia-vscode/issues/800
if isdefined(@__MODULE__, :LanguageServer)
    # invoked by VS Code
    include("../src/OilData.jl")
    using .OilData
else
    # invoked during test
    using OilData
end

od = OilData

logger = SimpleLogger(stdout, Logging.Info)
global_logger(logger)

@test true == true
@test false == false

@testset "add_day" begin
    d = Date(2021, 1, 1)
    t = Time(15, 32, 33)
    dt = DateTime(d,t)
    dt_1 = DateTime(Date(2021, 1, 2), Time(15,32,33))

    @test dt_1 == add_day(dt, 1.0)
    @test dt_1 == add_day(dt, 1)

    @test_throws MethodError add_day(1, 1.0)
    @test_throws MethodError add_day(1.0, 1.0)

    @info pwd()
end

@testset "find_prt_start_date" begin
   @test basename(data_dir()) == "data"
   path = joinpath(data_dir(), "Project-2.PRT")
   start_datetime = od.find_prt_start_date(path)
   #@info start_datetime
   expected = DateTime(Date(2019, 2, 1), Time(00,00,00))
   @test start_datetime == expected
end


text:
###############################################################################
#
# Subscript
#
###############################################################################
"""
    subscript(ss::Number; option="mathrm", standalone=true)
    subscript(ss::String; option="mathrm", standalone=true)

Return a string of subscript, given
"""
function subscript(ss::Number; option="mathrm", standalone=true)
    return subscript(string(ss); option=option, standalone=standalone)
end

function subscript(ss::String; option="mathrm", standalone=true)
    if ss==""
        return ""
    elseif standalone
        return "\$_\\$(option){$(ss)}\$"
    else
        return "_\\$(option){$(ss)} "
    end
end

###############################################################################
#
# Superscript
#
###############################################################################
"""
    superscript(ss::Number; option="mathrm", standalone=true)
    superscript(ss::String; option="mathrm", standalone=true)

Return a string of superscript, given
"""
function superscript(ss::Number; option="mathrm", standalone=true)
    return superscript(string(ss); option=option, standalone=standalone)
end

function superscript(ss::String; option="mathrm", standalone=true)
    if ss==""
        return ""
    elseif standalone
        return "\$^\\$(option){$(ss)}\$"
    else
        return "^\\$(option){$(ss)} "
    end
end


text:
module EngEconomics

using Dates

export simpleInterest, findFractionalPeriod

# Simple Interest (F = P(1 + i x n))
# F = Future amount
# P = Principal amount
# i = simple interest rate
# n = period of time
function simpleInterest(P::Number, i::Number, n::Number)::Number
	F = P * (1 + i * n)
end

# Find Fraction Period (360 days a year, 12 months, 60 days)
function findFractionalPeriod(date1::Date, date2::Date)::Number
	monthDays = 30 * (Dates.month(date2) - Dates.month(date1))
	d = Dates.day(date2) - Dates.day(date1) + 1
	return (monthDays + d) / 360
end

export learningCurveModel, learningCurveModelCoeff, learningCurvePercentage
function learningCurveModel(T_init, N, b)
	return T_init * N^b
end

function learningCurveModelCoeff(T1, N1, T2, N2)
	b = log(T1 / T2) / log(N1 / N2)
end

function learningCurvePercentage(b)
	return exp(b * log(2))
end

export compoundAmountFactor, presentWorthFactor
# Compound Amount Factor: (F/P, i, N)
# - gives future amount, F, that is equivalent to a present amount, P, when the
#   interest rate is i, and the number of periods is N
# - How much money will be in the bank after 15 years if we invest 100 dollars
#   now at 8% interest compounded semi-annually
function compoundAmountFactor(i, N)
	(1 + i)^N
end

# Present Worth Factor: (P/F, i, N)
# - gives the present amount, P, that is equivalent to the future amount, F,
#   when the interest rate is i, and the number of periods is N
# - How much should be invested now so that in 4 years we will have 1360.5
#   dollars in the bank if the interest rate is 8%
function presentWorthFactor(i, N)
	1 / ((1 + i)^N)
end

export sinkingFundFactor, uniformSeriesCompoundAmoundFactor
# Sinking Fund Factor: (A/F, i, N)
# - gives the size of A of a repeated receipt or disbursement that is equivalent
#   to a future amount F if the interest rate is i and the number of periods is
#   N
# - If you want to accumulate 50,000 in 5 years, how much should we invest
#   yearly if the interest rate is 10%? F * (A/F,10%, 5)
function sinkingFundFactor(i, N)
	i / ((1 + i)^N - 1)
end

# Uniform Series Compound Amount Factor: (F/A, i, N)
# - gives the future value F that is equivalent to a series of equal size
#   receipts or disbursements of A when the interest rate is i and the number of
#   periods is N
# - If you deposit $200 yearly for 45 years if the interest rate is 12% how much
#   will be in your savings account? 200 (A/F, 12%, 45)
function uniformSeriesCompoundAmoundFactor(i, N)
	((1+i)^N - 1) / i
end

export capitalRecoveryFactor, seriesPresentAmountFactor
# Capital Recovery Factor: (A/P, i, N)
# - gives the value of A of the equal periodic payments or receipts that are
#   equivalent to the present amount P when the interest is i and the number of
#   period is N
# - how much should I pay per month to be equal to just paying $100 dollars now
#   for the thing?
function capitalRecoveryFactor(i, N)
	(i * (1 + i)^N) / ((i+1)^N - 1)
end

# Series Present Worth Factor: (P/A,i,N)
# - gives the present amount, P, that is equivalent to an annuity with
#   disbursements or receipts in the amount of A, where the interest rate is i
#   and the number of periods is N
# - how much should I pay now to get out of paying $200 a month for 40 years if
#   the interest is 2%
function seriesPresentAmountFactor(i, N)
	((1 + i)^N - 1) / (i * (1 + i)^N)
end

export singlePaymentPresentWorthArithmetic, sinkingFundPaymentArithmetic
# Single Payment Present Worth for Arithmetic Gradient Series:
# - gives the present worth of paying a gradient cost for N periods at an
#   interest rate of i
# - (P/G, i, N)
function singlePaymentPresentWorthArithmetic(i, N)
	((1+i)^N - i * N - 1) / (i^2 * (1+i)^N)
end

# Sinking Fund Payment for Arithmetic Gradient Series
# - (A/G, i, N)
function sinkingFundPaymentArithmetic(i, N)
	(1 / i - N / ((1+i)^N - 1))
end

export presentWorthConversionFactorGeometric
# Conversion Factor for Geometric Gradient Series (P/A, g, i, N)
# - gives the present worth P that is equivalent to a geometric gradient series
#   where the base receipt or disbursement is A, the growth rate is g, and the
#   interest rate is i, and the number of periods is N
function presentWorthConversionFactorGeometric(g, i, N)
	i0 = (1 + i) / (1 + g) - 1
	(((1 + i0)^N - 1) / (i0 * (1+i0)^N)) * 1 / (1 + g)
end

export paybackPeriod
function paybackPeriod(firstCost, annualSavings)
	return firstCost / annualSavings
end

end


text:
using SHIPs.SphericalHarmonics: SHBasis, index_y
using StaticArrays
using JuLIP: AbstractCalculator, Atoms, JVec
using JuLIP.Potentials: MSitePotential, SZList, ZList
using NeighbourLists: max_neigs, neigs

import JuLIP, JuLIP.MLIPs
import JuLIP: cutoff, alloc_temp, alloc_temp_d
import JuLIP.Potentials: evaluate!, evaluate_d!
import JuLIP.MLIPs: combine
import Base: Dict, convert, ==

export SHIP

struct SHIP{T, NZ, TJ} <: MSitePotential
   J::TJ
   SH::SHBasis{T}
   zlist::SZList{NZ}
   # -------------- A and AA datastructures + coefficients
   alists::NTuple{NZ, AList}
   aalists::NTuple{NZ, AAList}
   coeffs::NTuple{NZ, Vector{Complex{T}}}
end

cutoff(ship::SHIP) = cutoff(ship.J)

==(S1::SHIP, S2::SHIP) =
      all( getfield(S1, i) == getfield(S2, i)
           for i = 1:fieldcount(SHIP) )

Base.length(ship::SHIP) = sum(length.(ship.aalists))

# BO + 1 because BO is the number of neighbours not the actual body-order
bodyorder(ship::SHIP{BO}) where {BO} = BO + 1

# ------------------------------------------------------------
#   Initialisation code
# ------------------------------------------------------------

combine(basis::SHIPBasis, coeffs) = SHIP(basis, coeffs)

function SHIP(basis::SHIPBasis{T, NZ}, coeffs::AbstractVector{T}
              ) where {T, NZ}
   AAcoeffs = ntuple(iz0 -> (coeffs[_get_I_iz0(basis, iz0)]' * basis.A2B[iz0])[:], NZ)
   return SHIP( basis.J, basis.SH, basis.zlist,
                basis.alists, basis.aalists, AAcoeffs )
end

# ------------------------------------------------------------
#   FIO code
# ------------------------------------------------------------

Dict(ship::SHIP{T,NZ}) where {T, NZ} = Dict(
      "__id__" => "SHIPs_SHIP_v2",
      "J" => Dict(ship.J),
      "SH_maxL" => ship.SH.maxL,   # TODO: replace this with Dict(SH)
      "T" => string(eltype(ship.SH)),
      "zlist" => Dict(ship.zlist),
      "alists" => [Dict.(ship.alists)...],
      "aalists" => [Dict.(ship.aalists)...],
      "coeffs_re" => [ real.(ship.coeffs[i]) for i = 1:NZ  ],
      "coeffs_im" => [ imag.(ship.coeffs[i]) for i = 1:NZ  ]
   )

convert(::Val{:SHIPs_SHIP_v2}, D::Dict) = SHIP(D)

# bodyorder - 1 is because BO is the number of neighbours
# not the actual body-order
function SHIP(D::Dict)
   T = Meta.eval(Meta.parse(D["T"]))
   J = TransformedJacobi(D["J"])
   SH = SHBasis(D["SH_maxL"], T)
   zlist = decode_dict(D["zlist"])
   NZ = length(zlist)
   alists = ntuple(i -> AList(D["alists"][i]), NZ)
   aalists = ntuple(i -> AAList(D["aalists"][i], alists[i]), NZ)
   coeffs = ntuple(i -> T.(D["coeffs_re"][i]) + im * T.(D["coeffs_im"][i]), NZ)
   return  SHIP(J, SH, zlist, alists, aalists, coeffs)
end

# ------------------------------------------------------------
#   Evaluation code
# ------------------------------------------------------------

alloc_temp(ship::SHIP{T,NZ}, N::Integer) where {T, NZ} =
   (     J = alloc_B(ship.J),
         Y = alloc_B(ship.SH),
         A = [ zeros(Complex{T}, length(ship.alists[iz])) for iz=1:NZ ],
      tmpJ = alloc_temp(ship.J),
      tmpY = alloc_temp(ship.SH),
         R = zeros(JVec{T}, N),
         Z = zeros(Int16, N)
           )

# compute one site energy
function evaluate!(tmp, ship::SHIP{T},
                   Rs::AbstractVector{JVec{T}},
                   Zs::AbstractVector{<:Integer},
                   z0::Integer) where {T}
   iz0 = z2i(ship, z0)
   precompute_A!(tmp.A[iz0], tmp, ship.alists[iz0], Rs, Zs, ship)
   aalist = ship.aalists[iz0]
   A = tmp.A[iz0]
   c = ship.coeffs[iz0]
   Es = zero(T)
   for iAA = 1:length(aalist)
      Esi = c[iAA] # one(Complex{T})
      for α = 1:aalist.len[iAA]
         Esi *= A[aalist.i2Aidx[iAA, α]]
      end
      Es += real(Esi)
   end
   return Es
end

alloc_temp_d(ship::SHIP{T, NZ}, N::Integer) where {T, NZ} =
      ( J = alloc_B(ship.J),
       dJ = alloc_dB(ship.J),
        Y = alloc_B(ship.SH),
       dY = alloc_dB(ship.SH),
        A = [ zeros(Complex{T}, length(ship.alists[iz])) for iz=1:NZ ],
     dAco = [ zeros(Complex{T}, length(ship.alists[iz])) for iz=1:NZ ],
     tmpJ = alloc_temp_d(ship.J),
     tmpY = alloc_temp_d(ship.SH),
       dV = zeros(JVec{T}, N),
        R = zeros(JVec{T}, N),
        Z = zeros(Int16, N)
      )

# compute one site energy
function evaluate_d!(dEs, tmp, ship::SHIP{T},
                     Rs::AbstractVector{JVec{T}},
                     Zs::AbstractVector{<:Integer},
                     z0::Integer
                     ) where {T}
   iz0 = z2i(ship, z0)
   alist = ship.alists[iz0]
   aalist = ship.aalists[iz0]

   # stage 1: precompute all the A values
   precompute_A!(tmp.A[iz0], tmp, alist, Rs, Zs, ship)

   # stage 2: compute the coefficients for the ∇A_{klm} = ∇ϕ_{klm}
   A = tmp.A[iz0]
   dAco = tmp.dAco[iz0]
   _evaluate_d_stage2!(dAco, A, aalist, ship.coeffs[iz0], ship)

   # stage 3: get the gradients
   fill!(dEs, zero(JVec{T}))
   for (iR, (R, Z)) in enumerate(zip(Rs, Zs))
      R̂ = R / norm(R)
      eval_basis_d!(tmp.J, tmp.dJ, tmp.tmpJ, ship.J, norm(R))
      eval_basis_d!(tmp.Y, tmp.dY, tmp.tmpY, ship.SH, R)
      iz = z2i(ship, Z)
      for iA = alist.firstz[iz]:(alist.firstz[iz+1]-1)
         zklm = alist[iA]
         ik, iy = zklm.k+1, index_y(zklm.l, zklm.m)
         dEs[iR] += real(dAco[iA] * (
                         tmp.J[ik] * tmp.dY[iy] + tmp.dJ[ik] * tmp.Y[iy] * R̂) )
      end
   end
   return dEs
end

function _evaluate_d_stage2!(dAco::AbstractVector{CT}, A, aalist, c, ship
                            ) where {CT}
   fill!(dAco, 0)
   for iAA = 1:length(aalist)
      for α = 1:aalist.len[iAA]
         CxA_α = CT(c[iAA])
         for β = 1:aalist.len[iAA]
            if β != α
               iAβ = aalist.i2Aidx[iAA, β]
               CxA_α *= A[iAβ]
            end
         end
         iAα = aalist.i2Aidx[iAA, α]
         dAco[iAα] += CxA_α
      end
   end
end

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
using SpectralMethodFFT
using Test

@testset "SpectralMethodFFT.jl" begin

    # +++++ 1-dimensional +++++
    c = ConfigFFT((64,), ((0., 2π),))
    x_X = x_Xgen(c)
    x_f, x_g = exp(cos(x_X)), exp(sin(x_X))
    x_fg_exact = exp(cos(x_X) + sin(x_X))

    k_f, k_g = k_x.((x_f, x_g))
    x_fg = x_k(k_f ⊙ k_g)
    @test norm_x(x_fg - x_fg_exact) / (2π) < 2eps()
    x_fg = x_k(k_f ⊗ k_g)
    @test norm_x(x_fg - x_fg_exact) / (2π) < 2eps()

    # +++++ 2-dimensional +++++
    c = ConfigFFT((64, 64), ((0., 2π), (0., 2π)))
    xy_X, xy_Y = xy_Xgen(c), xy_Ygen(c)
    xy_f, xy_g = exp(cos(xy_X)), exp(cos(xy_Y))
    xy_fg_exact = exp(cos(xy_X) + cos(xy_Y))

    kl_f, kl_g = kl_xy.((xy_f, xy_g))
    xy_fg = xy_kl(kl_f ⊙ kl_g)
    @test norm_xy(xy_fg - xy_fg_exact) / ((2π)^2) < 2eps()
    xy_fg = xy_kl(kl_f ⊗ kl_g)
    @test norm_xy(xy_fg - xy_fg_exact) / ((2π)^2) < 2eps()

    # +++++ 3-dimensional +++++
    c = ConfigFFT((64, 64, 64), ((0., 2π), (0., 2π), (0., 2π)))
    xyz_X, xyz_Y, xyz_Z = xyz_Xgen(c), xyz_Ygen(c), xyz_Zgen(c)
    xyz_f = exp(cos(xyz_X) + cos(xyz_Y))
    xyz_g = exp(cos(xyz_Y) + cos(xyz_Z))
    xyz_fg_exact = exp(cos(xyz_X) + 2cos(xyz_Y) + cos(xyz_Z))

    klm_f, klm_g = klm_xyz.((xyz_f, xyz_g))
    xyz_fg = xyz_klm(klm_f ⊙ klm_g)
    @test norm_xyz(xyz_fg - xyz_fg_exact) / ((2π)^3) < 2eps()
    xyz_fg = xyz_klm(klm_f ⊗ klm_g)
    @test norm_xyz(xyz_fg - xyz_fg_exact) / ((2π)^3) < 2eps()

end

max_stars_repo_name: mat-der-D/SpectralMethodFFTW

max_issues_repo_name: mat-der-D/SpectralMethodFFTW

max_forks_repo_name: mat-der-D/SpectralMethodFFTW


text:
# zlib inflate implementation similar to: https://zlib.net/zpipe.c
mutable struct z_stream
    next_in::Ptr{UInt8}
    avail_in::Cuint
    total_in::Culong
    next_out::Ptr{UInt8}
    avail_out::Cuint
    total_out::Culong
    msg::Ptr{UInt8}
    state::Ptr{Cvoid}
    zalloc::Ptr{Cvoid}
    zfree::Ptr{Cvoid}
    opaque::Ptr{Cvoid}
    data_type::Cint
    adler::Culong
    reserved::Culong

    function z_stream()
        new(
            C_NULL,  # next_in
            0,       # avail_in
            0,       # total_in
            C_NULL,  # next_out
            0,       # avail_out
            0,       # total_out
            C_NULL,  # msg
            C_NULL,  # state
            C_NULL,  # zalloc
            C_NULL,  # zfree
            C_NULL,  # opaque
            0,       # data_type
            0,       # adler
            0,       # reserved
        )
    end
end

const Z_OK            = 0
const Z_STREAM_END    = 1
const Z_NEED_DICT     = 2
const Z_ERRNO         = -1
const Z_STREAM_ERROR  = -2
const Z_DATA_ERROR    = -3
const Z_MEM_ERROR     = -4
const Z_BUF_ERROR     = -5
const Z_VERSION_ERROR = -6

@static if Base.VERSION > v"1.3-" && !(Sys.iswindows() && Sys.WORD_SIZE == 32)
    using Pkg
    Pkg.add("Zlib_jll")
    using Zlib_jll: libz
else
    isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl")) ||
        error("PDFIO not properly installed. Please run Pkg.build(\"PDFIO\")")
    include("../deps/deps.jl")
end

_zlibVersion() = ccall((:zlibVersion, libz), Ptr{Cstring}, ())

_inflateInit(stm::z_stream) =
    ccall((:inflateInit2_, libz), Cint,
          (Ptr{Cvoid}, Cint, Ptr{Cstring}, Cint),
          Ref(stm), 47, _zlibVersion(), sizeof(z_stream))

_inflateEnd(stm::z_stream) =
    ccall((:inflateEnd, libz), Cint, (Ref{z_stream},), stm)

_inflate(stm::z_stream) =
    ccall((:inflate, libz), Cint, (Ref{z_stream}, Cint), stm, 0)

function inflate(io::IO)
    CHUNK = 16384

    iob = IOBuffer()

    strm = z_stream()
    ret = _inflateInit(strm)
    ret != Z_OK && _zerror(ret)

    inb = zeros(UInt8, CHUNK)
    oub = zeros(UInt8, CHUNK)

    try
        while ret != Z_STREAM_END
            strm.avail_in = readbytes!(io, inb, CHUNK)
            strm.avail_in == 0 && break
            strm.next_in = pointer(inb)

            strm.avail_out = 0
            while strm.avail_out == 0
                strm.avail_out = CHUNK
                strm.next_out = pointer(oub)
                ret = _inflate(strm)
                ret == Z_STREAM_ERROR && error("zlib stream state clobbered")
                ret != Z_OK && ret != Z_STREAM_END && _zerror(ret)
                have = CHUNK - strm.avail_out
                resize!(oub, have)
                write(iob, oub)
                resize!(oub, CHUNK)
            end
        end
    finally
        _inflateEnd(strm)
    end
    return seekstart(iob)
end

function _zerror(ret::Cint)
    msg =
        ret == Z_STREAM_ERROR  ? "invalid compression level" :
        ret == Z_DATA_ERROR    ? "invalid or incomplete deflate data" :
        ret == Z_MEM_ERROR     ? "out of memory" :
        ret == Z_VERSION_ERROR ? "zlib version mismatch!" :
                                 "zlib internal error"

    error(msg)
end

init_table(size::Int) =
    ([(i <= 256) ? [ UInt8(i-1)] :
      i == 257  ? [ UInt8(0)]   :
      i == 258  ? [ UInt8(0)] : UInt8[] for i = 1:size], 9, 258)
@inline function next_number(bin, cl, sy, si)
    sy -= 1
    si -= 1
    n = 0
    while cl > 0
        @inbounds b::Int = bin[sy + 1]
        m = UInt8(0xff) >> si
        b = b & m
        cl -= (8 - si)
        b <<= cl
        n += b
        si, sy = 0, sy + 1
        1 <= cl <= 7 || continue
        si = cl
        @inbounds b = bin[sy + 1]
        m = UInt8(0xff) << (8 - si)
        b = (b & m) >> (8 - si)
        n += b
        cl = 0
    end
    return n, sy + 1, si + 1
end

function decode_lzw(io::IO, earlyChange::Int = 1)
    bin = read(io)
    util_close(io)
    len = length(bin)
    t, cl, it = init_table(4096)
    hy, hi = 1, 1
    iob = IOBuffer()
    old_set = false
    old, c = 1, 0x0
    s = UInt8[]
    while hy < len && hi <= 8
        n, hy, hi = next_number(bin, cl, hy, hi)
        hy > len && break
        n == 257 && break
        if n == 256
            t, cl, it = init_table(4096)
            old_set = false
        else
            if n < it
                if !old_set
                    old = n
                    write(iob, t[n+1])
                    c = t[n + 1][1]
                    old_set = true
                    continue
                end
                s = t[n + 1]
            elseif n == it
                s = copy(t[old + 1])
                push!(s, c)
            else
                error(E_FAILED_COMPRESSION*" $n : $it")
            end
            write(iob, s)
            c = s[1]
            append!(t[it+1], t[old+1])
            push!(t[it+1], c)
            old = n
            it == 4095 && continue
            if it == ((1 << cl) - 1 - earlyChange)
                cl += 1
            end
            it += 1
        end
    end
    return seekstart(iob)
end