text:
module ParameterDistributionStorage

## Imports
using Distributions
using Statistics
using StatsBase
using Random

## Exports

#types

export ParameterDistributionType
#objects
export Parameterized, Samples
export ParameterDistribution
export Constraint

#functions
export get_name, get_distribution, get_total_dimension, get_dimensions, get_all_constraints, get_n_samples
export sample_distribution
export no_constraint, bounded_below, bounded_above, bounded
export transform_constrained_to_unconstrained, transform_unconstrained_to_constrained
export get_logpdf, get_cov, get_var, get_mean, batch

## Objects
# for the Distribution
abstract type ParameterDistributionType end

"""
    struct Parameterized <: ParameterDistributionType

A distribution constructed from a parametrized formula (e.g Julia Distributions.jl)
"""
struct Parameterized <: ParameterDistributionType
    distribution::Distribution
end

"""
    struct Samples{FT<:Real} <: ParameterDistributionType

A distribution comprised of only samples, stored as columns of parameters
"""
struct Samples{FT<:Real} <: ParameterDistributionType
    distribution_samples::Array{FT,2} #parameters are columns
    Samples(distribution_samples::Array{FT,2}; params_are_columns=true) where {FT <: Real} = params_are_columns ? new{FT}(distribution_samples) : new{FT}(permutedims(distribution_samples,(2,1)))
    #Distinguish 1 sample of an ND parameter or N samples of 1D parameter, and store as 2D array
    Samples(distribution_samples::Array{FT,1}; params_are_columns=true) where {FT <: Real} = params_are_columns ? new{FT}(reshape(distribution_samples,1,:)) : new{FT}(reshape(distribution_samples,:,1))

end

# For the transforms
abstract type ConstraintType end

"""
    struct Constraint <: ConstraintType

Contains two functions to map between constrained and unconstrained spaces.
"""
struct Constraint <: ConstraintType
    constrained_to_unconstrained::Function
    unconstrained_to_constrained::Function
end

"""
    function no_constraint()

Constructs a Constraint with no constraints, enforced by maps x -> x and x -> x.
"""
function no_constraint()
    c_to_u = (x -> x)
    u_to_c = (x -> x)
    return Constraint(c_to_u,u_to_c)
end

"""
    function bounded_below(lower_bound::FT) where {FT <: Real}

Constructs a Constraint with provided lower bound, enforced by maps x -> log(x - lower_bound) and x -> exp(x) + lower_bound.
"""
function bounded_below(lower_bound::FT) where {FT <: Real}
    c_to_u = ( x -> log(x - lower_bound) )
    u_to_c = ( x -> exp(x) + lower_bound )
    return Constraint(c_to_u, u_to_c)
end

"""
    function bounded_above(upper_bound::FT) where {FT <: Real}

Constructs a Constraint with provided upper bound, enforced by maps x -> log(upper_bound - x) and x -> upper_bound - exp(x).
"""
function bounded_above(upper_bound::FT) where {FT<:Real}
    c_to_u = ( x -> log(upper_bound - x) )
    u_to_c = ( x -> upper_bound - exp(x) )
    return Constraint(c_to_u, u_to_c)
end

"""
    function bounded(lower_bound::FT, upper_bound::FT) where {FT <: Real}

Constructs a Constraint with provided upper and lower bounds, enforced by maps
x -> log((x - lower_bound) / (upper_bound - x))
and
x -> (upper_bound * exp(x) + lower_bound) / (exp(x) + 1)

"""
function bounded(lower_bound::FT, upper_bound::FT) where {FT <: Real}
    if (upper_bound <= lower_bound)
        throw(DomainError("upper bound must be greater than lower bound"))
    end
    c_to_u = ( x -> log( (x - lower_bound) / (upper_bound - x)) )
    u_to_c = ( x -> (upper_bound * exp(x) + lower_bound) / (exp(x) + 1) )
    return Constraint(c_to_u, u_to_c)
end

"""
    function len(c::Array{CType})

The number of constraints, each constraint has length 1.
"""
function len(c::CType) where {CType <: ConstraintType}
    return 1
end

function len(carray::Array{CType}) where {CType <: ConstraintType}
    return size(carray)[1]
end

"""
    function get_dimensions(d<:ParametrizedDistributionType)

The number of dimensions of the parameter space
"""
function dimension(d::Parameterized)
    return length(d.distribution)
end

function dimension(d::Samples)
    return size(d.distribution_samples)[1]
end

"""
    function n_samples(d::Samples)

The number of samples in the array
"""
function n_samples(d::Samples)
    return size(d.distribution_samples)[2]
end
function n_samples(d::Parameterized)
    return "Distribution stored in Parameterized form, draw samples using `sample_distribution` function"
end

"""
    struct ParameterDistribution

Structure to hold a parameter distribution, always stored as an array of distributions
"""
struct ParameterDistribution{PDType <: ParameterDistributionType, CType <: ConstraintType, ST <: AbstractString}
    distributions::Array{PDType}
    constraints::Array{CType}
    names::Array{ST}

    function ParameterDistribution(parameter_distributions::Union{PDType,Array{PDType}},
                                   constraints::Union{CType, Array{CType}, Array},
                                   names::Union{ST, Array{ST}}) where {PDType <: ParameterDistributionType,
                                                                       CType <: ConstraintType,
                                                                       ST <: AbstractString}

        parameter_distributions = isa(parameter_distributions, PDType) ? [parameter_distributions] : parameter_distributions
        n_parameters_per_dist = [dimension(pd) for pd in parameter_distributions]
        constraints = isa(constraints, Union{<:ConstraintType,Array{<:ConstraintType}}) ? [constraints] : constraints #to calc n_constraints_per_dist
        names = isa(names, ST) ? [names] : names

        n_constraints_per_dist = [len(c) for c in constraints]
        n_dists = length(parameter_distributions)
        n_names = length(names)
        if !(n_parameters_per_dist == n_constraints_per_dist)
            throw(DimensionMismatch("There must be one constraint per parameter in a distribution, use NoConstraint() type if no constraint is required"))
        elseif !(n_dists == n_names)
            throw(DimensionMismatch("There must be one name per parameter distribution"))
        else
            constraints = cat(constraints...,dims=1)

            new{PDType,ConstraintType,ST}(parameter_distributions, constraints, names)
        end
    end

end

## Functions

"""
    function get_name(pd::ParameterDistribution)

Returns a list of ParameterDistribution names
"""
function get_name(pd::ParameterDistribution)
    return pd.names
end

"""
    function get_dimensions(pd::ParameterDistribution)

The number of dimensions of the parameter space
"""
function get_dimensions(pd::ParameterDistribution)
    return [dimension(d) for d in pd.distributions]
end
function get_total_dimension(pd::ParameterDistribution)
    return sum(dimension(d) for d in pd.distributions)
end

"""
    function get_n_samples(pd::ParameterDistribution)

The number of samples in a Samples distribution
"""
function get_n_samples(pd::ParameterDistribution)
    return Dict{String,Any}(pd.names[i] => n_samples(d) for (i,d) in enumerate(pd.distributions))
end

"""
    function get_all_constraints(pd::ParameterDistribution)

returns the (flattened) array of constraints of the parameter distribution
"""
function get_all_constraints(pd::ParameterDistribution)
    return pd.constraints
end

"""
    function batch(pd:ParameterDistribution)

Returns a list of contiguous [collect(1:i), collect(i+1:j),... ] used to split parameter arrays by distribution dimensions
"""
function batch(pd::ParameterDistribution)
    #chunk xarray to give to the different distributions.
    d_dim = get_dimensions(pd) #e.g [4,1,2]
    d_dim_tmp = Array{Int64}(undef,size(d_dim)[1]+1)
    d_dim_tmp[1] = 0
    for i = 2:size(d_dim)[1]+1
        d_dim_tmp[i] = sum(d_dim[1:i-1]) # e.g [0,4,5,7]
    end

    return [collect(d_dim_tmp[i]+1:d_dim_tmp[i+1]) for i = 1:size(d_dim)[1]] # e.g [1:4, 5:5, 6:7]
end

"""
    function get_distribution(pd::ParameterDistribution)

Returns a `Dict` of `ParameterDistribution` distributions, with the parameter names as dictionary keys. For parameters represented by `Samples`, the samples are returned as a 2D (parameter_dimension x n_samples) array
"""
function get_distribution(pd::ParameterDistribution)
    return Dict{String,Any}(pd.names[i] => get_distribution(d) for (i,d) in enumerate(pd.distributions))
end

function get_distribution(d::Samples)
    return d.distribution_samples
end
function get_distribution(d::Parameterized)
    return d.distribution
end

"""
    function sample_distribution(pd::ParameterDistribution)

Draws samples from the parameter distributions returns an array, with parameters as columns
"""
function sample_distribution(pd::ParameterDistribution)
    return sample_distribution(pd,1)
end

function sample_distribution(pd::ParameterDistribution, n_draws::IT) where {IT <: Integer}
    return cat([sample_distribution(d,n_draws) for d in pd.distributions]...,dims=1)
end

function sample_distribution(d::Samples, n_draws::IT)  where {IT <: Integer}
    n_stored_samples = n_samples(d)
    samples_idx = StatsBase.sample(collect(1:n_stored_samples), n_draws)
    if dimension(d) == 1
        return reshape(d.distribution_samples[:,samples_idx],:,n_draws) #columns are parameters
    else
        return d.distribution_samples[:,samples_idx]
    end
end

function sample_distribution(d::Parameterized, n_draws::IT) where {IT <: Integer}
    if dimension(d) == 1
        return reshape(rand(d.distribution,n_draws), :, n_draws) #columns are parameters
    else
        return rand(d.distribution, n_draws)
    end
end

"""
    function logpdf(pd::ParameterDistribution, xarray::Array{<:Real,1})

Obtains the independent logpdfs of the parameter distributions at xarray (non-Samples Distributions only), and returns their sum.
"""
function get_logpdf(d::Parameterized, xarray::Array{FT,1}) where {FT <: Real}
    return logpdf.(d.distribution, xarray)
end

function get_logpdf(pd::ParameterDistribution, xarray::Array{FT,1}) where {FT <: Real}
    #first check we don't have sampled distribution
    for d in pd.distributions
        if typeof(d) <: Samples
            throw(ErrorException("Cannot compute get_logpdf of Samples distribution. Consider using a Parameterized type for your prior."))
        end
    end
    #assert xarray correct dim/length
    if size(xarray)[1] != get_total_dimension(pd)
        throw(DimensionMismatch("xarray must have dimension equal to the parameter space"))
    end

    # get the index of xarray chunks to give to the different distributions.
    batches = batch(pd)

    # perform the logpdf of each of the distributions, and returns their sum
    return sum(cat([get_logpdf(d, xarray[batches[i]]) for (i,d) in enumerate(pd.distributions)]...,dims=1))
end

"""
    function get_cov(pd::ParameterDistribution)

returns a blocked covariance of the distributions
"""
function get_cov(d::Parameterized)
    return cov(d.distribution)
end

function get_cov(d::Samples)
    return cov(d.distribution_samples, dims=2) #parameters are columns
end

function get_var(d::Parameterized)
    return var(d.distribution)
end

function get_var(d::Samples)
    return var(d.distribution_samples)
end

function get_cov(pd::ParameterDistribution)
    #first check we don't have sampled distribution

    d_dims = get_dimensions(pd)

    # create each block (co)variance
    block_cov = Array{Any}(undef,size(d_dims)[1])
    for (i,dimension) in enumerate(d_dims)
        if dimension == 1
            block_cov[i] = get_var(pd.distributions[i])
        else
            block_cov[i] = get_cov(pd.distributions[i])
        end
    end

    return cat(block_cov...,dims=(1,2)) #build the block diagonal (dense) matrix

end

"""
    get_mean(pd::Parameterized)

returns a mean of parameterized distribution
"""
function get_mean(d::Parameterized)
    return mean(d.distribution)
end

"""
    get_mean(pd::Samples)

returns a mean of the samples
"""
function get_mean(d::Samples)
    return mean(d.distribution_samples, dims=2) #parameters are columns
end

"""
    function get_mean(pd::ParameterDistribution)

returns a mean of the distributions
"""
function get_mean(pd::ParameterDistribution)
    return cat([get_mean(d) for d in pd.distributions]...,dims=1)
end

#apply transforms

"""
    function transform_constrained_to_unconstrained(pd::ParameterDistribution, x::Array{<:Real})

Apply the transformation to map (possibly constrained) parameters `xarray` into the unconstrained space
"""
function transform_constrained_to_unconstrained(pd::ParameterDistribution, xarray::Array{FT,1}) where {FT <: Real}
    return cat([c.constrained_to_unconstrained(xarray[i]) for (i,c) in enumerate(pd.constraints)]...,dims=1)
end

"""
    function transform_constrained_to_unconstrained(pd::ParameterDistribution, xarray::Array{FT,2})

Apply the transformation to map (possibly constrained) parameter samples `xarray` into the unconstrained space.
Here, `xarray` contains parameters as columns and samples as rows.
"""
function transform_constrained_to_unconstrained(pd::ParameterDistribution, xarray::Array{FT,2}) where {FT <: Real}
    return Array(hcat([c.constrained_to_unconstrained.(xarray[i,:]) for (i,c) in enumerate(pd.constraints)]...)')
end

"""
    function transform_unconstrained_to_constrained(pd::ParameterDistribution, xarray::Array{<:Real,1})

Apply the transformation to map parameters `xarray` from the unconstrained space into (possibly constrained) space
"""
function transform_unconstrained_to_constrained(pd::ParameterDistribution, xarray::Array{FT,1}) where {FT <: Real}
    return cat([c.unconstrained_to_constrained(xarray[i]) for (i,c) in enumerate(pd.constraints)]...,dims=1)
end

"""
    function transform_unconstrained_to_constrained(pd::ParameterDistribution, xarray::Array{FT,2})

Apply the transformation to map parameter samples `xarray` from the unconstrained space into (possibly constrained) space.
Here, `xarray` contains parameters as columns and samples as rows.
"""
function transform_unconstrained_to_constrained(pd::ParameterDistribution, xarray::Array{FT,2}) where {FT <: Real}
    return Array(hcat([c.unconstrained_to_constrained.(xarray[i,:]) for (i,c) in enumerate(pd.constraints)]...)')
end

"""
    function transform_unconstrained_to_constrained(pd::ParameterDistribution, xarray::Array{Array{FT,2},1})

Apply the transformation to map parameter sample ensembles `xarray` from the unconstrained space into (possibly constrained) space.
Here, `xarray` contains parameters sample ensembles for different EKP iterations.
"""
function transform_unconstrained_to_constrained(pd::ParameterDistribution, xarray::Array{Array{FT,2},1}) where {FT <: Real}
    transf_xarray = []
    for elem in xarray
           push!(transf_xarray, transform_unconstrained_to_constrained(pd, elem))
       end
    return transf_xarray
end

end # of module


text:
# ------------------------------------------------------------------------------------------
# # Playing with matrices
# (and if you get the subliminal message about abstractions, we'll be thrilled!)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Julia is a dynamic language.  You don't need type declarations, and can change variable
# types dynamically and interactively.
#
# For working with simple numbers, arrays, and strings, its syntax is *superficially*
# similar to Matlab, Python, and other popular languages.
#
# In order to execute the `In` cells, select the cell and press `Shift-Enter`, or press the
# `Play` button above. To run the entire notebook, navigate to the `Cell` menu and then `Run
# All`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # We tell beginners, you don't need types!
# ------------------------------------------------------------------------------------------

typeof(1.0)

typeof(1)

S = "Hello Julia Class"
typeof(S)

#  Exercise fix my spelling in the cell above
S

# ------------------------------------------------------------------------------------------
# # Now forget all that (for now): Julia is not so different from your favorite dynamic
# language
# ------------------------------------------------------------------------------------------

1 + 1  # shift + enter to run

A = rand(5, 5)

using LinearAlgebra # A LinearAlgebra standard package contains structured matrices
A = SymTridiagonal(rand(6), rand(5))

b = rand(6)
A \ b

A = fill(3.15, 5, 5) # Fill a 5x5 array with 3.15's

# ------------------------------------------------------------------------------------------
# # Let's create some addition tables
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# First let's use a nested for loop to fill in a matrix that's initially zero:
# ------------------------------------------------------------------------------------------

A = zeros(5, 5)

for i in 1:5
    for j in 1:5
        A[i, j] = i+j  # Square brackets for indices.  Also: indices start at 1, not 0.
    end
end

A

# ------------------------------------------------------------------------------------------
# We can abbreviate this using a double `for` loop:
# ------------------------------------------------------------------------------------------

for i in 1:5, j in 1:5
   A[i, j] = i+j  # Square brackets for indices.  Also: indices start at 1, not 0.
end

A

# ------------------------------------------------------------------------------------------
# The Julia way would be to use a so-called **array comprehension**:
# ------------------------------------------------------------------------------------------

[i+j for i in 1:5, j in 1:5]

# Equivalently,
[i+j for i = 1:5, j = 1:5]

# ------------------------------------------------------------------------------------------
# **Explore**: What does the following do?
# ------------------------------------------------------------------------------------------

[i for i in (1:7).^2]

# What happens when  we remove the dot syntax?
[i for i in (1:7)^2]

[i^2 for i in 1:7]

# ------------------------------------------------------------------------------------------
# **Explore**: What does the following do?
# ------------------------------------------------------------------------------------------

sort(unique(x^2 + y^2 for x in 1:5, y in 1:5))  # The inner parentheses define a **generator**

# ------------------------------------------------------------------------------------------
# Suppose we want to see $n \times n$ multiplication tables for $n=1,2,3,4,5$:
# ------------------------------------------------------------------------------------------

for n in 1:5
    display([i*j for i=1:n, j=1:n])
end

# ------------------------------------------------------------------------------------------
# # `Interact.jl` is a Julia *package* for interacting with data
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# It's way more fun to **interact** with our data.
# We install the `Interact.jl` package as follows; this needs to be executed only once for
# any given Julia installation:
# ------------------------------------------------------------------------------------------

# using Pkg; Pkg.add("Interact")

# ------------------------------------------------------------------------------------------
# Now we load the package with the following `using` command, in each Julia session:
# ------------------------------------------------------------------------------------------

using Interact

# ------------------------------------------------------------------------------------------
# The package contains a `@manipulate` macro, that is wrapped around a `for` loop:
# ------------------------------------------------------------------------------------------

@manipulate for n in 1:1000
    n
end

@manipulate for n in 1:20
    [i*j for i in 1:n, j in 1:n]
end

# ------------------------------------------------------------------------------------------
# We use a double `for` loop to get a double slider!
# ------------------------------------------------------------------------------------------

@manipulate for n in 3:10, i in 1:9
   A = fill(0, n, n)
   A[1:3, 1:3] .= i    # fill a sub-block
A
end

# ------------------------------------------------------------------------------------------
# # Functions
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Julia is built around functions: all "commands" or "operations" in Julia are functions:
# ------------------------------------------------------------------------------------------

# verbose form:
function f(x)
    x^2
end

# one-line form:
f2(x) = x^2

# anonymous form:
f3 = x -> x^2;

f(10)

# ------------------------------------------------------------------------------------------
# Functions just work, as long as they make sense:
# ------------------------------------------------------------------------------------------

# The square of a matrix is unambiguously defined
f(rand(3, 3))

# What the 'square of a vector' means is ambiguous
f(rand(3))

# In the definition below, `a` and `power` are optional arguments to `f`, supplied with default values.
function f(x, a=1, power=2)
    a*x^power
end

# `a` defaults to 1 and `power` to 2
f(7)

# The first optional argument passed is assigned to the local variable `a`
# `power` defaults to 2
f(10, 3)

f(10, 3, 3)

# ------------------------------------------------------------------------------------------
# Let's define a function to insert a block in a matrix:
# ------------------------------------------------------------------------------------------

function insert_block(A, i, j, what=7)
    B = A[:,:]        # B is a copy of A
    B[i:i+2, j:j+2] = fill(what, 3, 3)

    return B          # the `return` keyword is optional
end

A = fill(0, 9, 9)
insert_block(A, 3, 5)  # this returns the new matrix

A = fill(0, 9, 9)
insert_block(A, 3, 5, 2)  # Use 2 instead of 7

# ------------------------------------------------------------------------------------------
# We can move the block around:
# ------------------------------------------------------------------------------------------

A = fill(0, 10, 10)
n = size(A, 1)

@manipulate for i in 1:n-2, j in 1:n-2
    insert_block(A, i, j)
end

# ------------------------------------------------------------------------------------------
# # Strings
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Julia can manipulate strings easily:
# ------------------------------------------------------------------------------------------

S = "Hello"

replace(S, "H" => "J")

a = 3

string(S, " ", S, " ", "Julia; a = ", a)  # build a string by concatenating things

# ------------------------------------------------------------------------------------------
# More about strings:  Julia
# Doc on Strings
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Functions in Julia try to be **generic**, i.e. to work with as many kinds of object as
# possible:
# ------------------------------------------------------------------------------------------

A = fill("Julia", 5, 5)

# ------------------------------------------------------------------------------------------
# Julia allows us to display objects in different ways. For example, the following code
# displays a matrix of strings
# in the notebook using an HTML representation:
# ------------------------------------------------------------------------------------------

function Base.show(io::IO, ::MIME"text/html", M::Matrix{T}) where {T<:String}
    max_length = maximum(length.(M))
    dv="
"
    print(io, dv*join([join("
".*M[i,:].*"
", " ") for i in 1:size(M, 1)]
            , "
$dv")*"
")
end

A

# Remember this ????
A = fill(0, 10, 10)
n = size(A, 1)

@manipulate for i in 1:n-2, j in 1:n-2
    insert_block(A, i,j)
end

# ------------------------------------------------------------------------------------------
# Let's use the **same code**, but now with strings:
# ------------------------------------------------------------------------------------------

A = fill("Julia", 10, 10)
n = size(A, 1)

@manipulate for i in 1:n-2, j in 1:n-2
    insert_block(A, i,j, "[FUN]")
end

@manipulate for i in 1:n-2, j in 1:n-2
    insert_block(A, i, j, "π")
end

@manipulate for i in 1:n-2, j in 1:n-2
    insert_block(A, i, j, "♡")
end

airplane = "✈"
alien = "👽"
rand([airplane, alien], 5, 5)

A = fill(airplane, 9, 9)
n = size(A, 1)

@manipulate for i in 1:n-2, j in 1:n-2
    insert_block(A, i, j, alien)
end

# ------------------------------------------------------------------------------------------
# # Colors
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# The `Colors` package provides objects representing colours:
# ------------------------------------------------------------------------------------------

# using Pkg; Pkg.add("Colors")
using Colors

distinguishable_colors(12)

@manipulate for n in 1:80
    distinguishable_colors(n)
end

colors = distinguishable_colors(100)

# Remember this ????
A = fill(0, 10, 10)
n = size(A, 1)

@manipulate for i in 1:n-2, j in 1:n-2
    insert_block(A, i, j)
end

# ------------------------------------------------------------------------------------------
# What happens if we use colors instead?
# ------------------------------------------------------------------------------------------

A = fill(colors[1], 10, 10)
n = size(A, 1)

@manipulate for i in 1:n-2, j in 1:n-2
    insert_block(A, i, j, colors[4])
end

# Exercise: Create Tetris Pieces, have them fall from the top

ext: jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using Behavior
using Niancat.Persistence
using Niancat.Instances
using Niancat.Games.NiancatGames
using Niancat.Languages
using Niancat.Scores

@then("the solution board lists {String} as solving {String}") do context, username, word
    game = context[:game]
    db = context[:db]
    dictionary = context[:dictionary]

    board = getsolutionboard(db, gameinstanceid(game), gameround(game))

    users = board.solutions[normalizedword(dictionary, word)]
    expecteduser = getuser(db, username, "defaultteam")

    @expect expecteduser in users
end

@then("the solution board does not list {String} as solving {String}") do context, username, word
    game = context[:game]
    db = context[:db]
    dictionary = context[:dictionary]

    board = getsolutionboard(db, gameinstanceid(game), gameround(game))

    users = board.solutions[normalizedword(dictionary, word)]
    expecteduser = getuser(db, username, "defaultteam")

    @expect !(expecteduser in users)
end


text:
using Test
using DelimitedFiles: readdlm

import Seis
using Beamforming

@testset "Stack" begin
    # These test data are generated with the script test/SAC/sac_stack.jl
    let t = Seis.sample_data(:array)
        # Time limits of stack, relative to A marker
        t1, t2 = -5, 10

        # Stacking weights
        weights = readdlm(joinpath(@__DIR__, "SAC", "stack_weights.txt"))

        for stack_type in ("linear", "nthroot", "phaseweight"), n in 1:2:5
            stack_file = joinpath(@__DIR__, "SAC", "$(stack_type)")
            if stack_type == "linear"
                n > 1 && continue
                nstring = ""
            else
                stack_file *= "-$n"
                nstring = "$n"
            end
            stack_file *= ".SHZ"

            sacstack = Seis.read_sac(stack_file)
            s = Beamforming.stack(t, (t1, t2), :A,
                method=Symbol(stack_type), n=n)

            @test Seis.trace(Seis.normalise(s)) ≈
                Seis.trace(Seis.normalise(sacstack)) rtol=0.2
        end
    end
end

@testset "Utilities" begin
    @testset "s/° ↔ s/km" begin
        @test Beamforming.s_per_degree(Beamforming.s_per_km(7, 1), 1) ≈ 7
        @test Beamforming.s_per_degree(1, 1) ≈ 2π/360
        @test Beamforming.s_per_degree(1) ≈ 111.19492664455873
    end
end

max_stars_repo_name: anowacki/Beamforming

max_issues_repo_name: anowacki/Beamforming

max_forks_repo_name: anowacki/Beamforming


text:
module profile

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
export flow, flowdθ, flowdσ, flowdψ, STARTING_σ_ψ, STARTING_log_ogip, STARTING_t

const STARTING_σ_ψ      = 0x1.baddbb87af68ap-2 # = 0.432
const STARTING_log_ogip = 0x1.670bf3d5b282dp-1 # = 0.701
const STARTING_t = 2*0.042/(2016-2003)

# functions in case we have months to expiration
@inline flow(  FF::Type, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T,             d::Integer, d1::Integer, Dgt0::Bool, sgn_ext::Bool, τ::Real, geoid::Real, roy::T) where {N,T} = flow(  FF, θ, σ, z, ψ,    d, d1, Dgt0, sgn_ext, geoid, roy)
@inline flowdθ(FF::Type, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T, k::Integer, d::Integer, d1::Integer, Dgt0::Bool, sgn_ext::Bool, τ::Real, geoid::Real, roy::T) where {N,T} = flowdθ(FF, θ, σ, z, ψ, k, d, d1, Dgt0, sgn_ext, geoid, roy)
@inline flowdσ(FF::Type, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T,             d::Integer, d1::Integer, Dgt0::Bool, sgn_ext::Bool, τ::Real, geoid::Real, roy::T) where {N,T} = flowdσ(FF, θ, σ, z, ψ,    d,                    geoid, roy)
@inline flowdψ(FF::Type, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T,             d::Integer,                          sgn_ext::Bool, τ::Real, geoid::Real, roy::T) where {N,T} = flowdψ(FF, θ, σ, z, ψ,    d,           sgn_ext, geoid, roy)
@inline flowdψ(FF::Type, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T,             d::Integer, d1::Integer, Dgt0::Bool, sgn_ext::Bool, τ::Real, geoid::Real, roy::T) where {N,T} = flowdψ(FF, θ, σ, z, ψ,    d,           sgn_ext, geoid, roy)

@inline flowdσ(FF::Type, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T,             d::Integer,                                                  geoid::Real, roy::T) where {N,T} = flowdσ(FF, θ, σ, z, ψ,    d,                    geoid, roy)

# functions in case we have months to expiration
@inline flow(  FF::Type, wp::AbstractUnitProblem, i::Integer, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T,             d::Integer, geoid::Real, roy::T) where {N,T} = flow(  FF, θ, σ, z, ψ,    d, _d1(wp,i), _Dgt0(wp,i), _sgnext(wp,i), geoid, roy)
@inline flowdθ(FF::Type, wp::AbstractUnitProblem, i::Integer, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T, k::Integer, d::Integer, geoid::Real, roy::T) where {N,T} = flowdθ(FF, θ, σ, z, ψ, k, d, _d1(wp,i), _Dgt0(wp,i), _sgnext(wp,i), geoid, roy)
@inline flowdσ(FF::Type, wp::AbstractUnitProblem, i::Integer, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T,             d::Integer, geoid::Real, roy::T) where {N,T} = flowdσ(FF, θ, σ, z, ψ,    d,                                        geoid, roy)
@inline flowdψ(FF::Type, wp::AbstractUnitProblem, i::Integer, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T,             d::Integer, geoid::Real, roy::T) where {N,T} = flowdψ(FF, θ, σ, z, ψ,    d,                         _sgnext(wp,i), geoid, roy)

# @inline function flowdθ(::Type, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T, k::Integer, d::Integer, d1::Integer, Dgt0::Bool, sgn_ext::Bool, geoid::Real, roy::T)::T where {N,T}

# --------------------------- common revenue functions & derivatives  --------------------------------------

@inline    rev_exp_restricted(θ1::T, σ::T, logp::Real,          ψ::Real, Dgt0::Bool, geoid::Real, roy::Real) where {T} = rev_exp(   1, θ1, 1, STARTING_log_ogip, STARTING_σ_ψ,             σ, logp,    ψ, Dgt0, geoid, roy)
@inline drevdσ_exp_restricted(θ1::T, σ::T, logp::Real,          ψ::Real,             geoid::Real, roy::Real) where {T} = drevdσ_exp(1, θ1, 1, STARTING_log_ogip, STARTING_σ_ψ,             σ, logp,    ψ,       geoid, roy)
@inline drevdψ_exp_restricted(θ1::T, σ::T, logp::Real,          ψ::Real,             geoid::Real, roy::Real) where {T} = drevdψ_exp(1, θ1, 1, STARTING_log_ogip, STARTING_σ_ψ,             σ, logp,    ψ,       geoid, roy)

@inline    rev_exp_restricted(θ1::T, σ::T, logp::Real, t::Real, ψ::Real, Dgt0::Bool, geoid::Real, roy::Real) where {T} = rev_exp(   1, θ1, 1, STARTING_log_ogip, STARTING_σ_ψ, STARTING_t, σ, logp, t, ψ, Dgt0, geoid, roy)
@inline drevdσ_exp_restricted(θ1::T, σ::T, logp::Real, t::Real, ψ::Real,             geoid::Real, roy::Real) where {T} = drevdσ_exp(1, θ1, 1, STARTING_log_ogip, STARTING_σ_ψ, STARTING_t, σ, logp, t, ψ,       geoid, roy)
@inline drevdψ_exp_restricted(θ1::T, σ::T, logp::Real, t::Real, ψ::Real,             geoid::Real, roy::Real) where {T} = drevdψ_exp(1, θ1, 1, STARTING_log_ogip, STARTING_σ_ψ, STARTING_t, σ, logp, t, ψ,       geoid, roy)

# chebshev polynomials
# See http://www.aip.de/groups/soe/local/numres/bookcpdf/c5-8.pdf
@inline checkinterval(x::Real,min::Real,max::Real) =  min <= x <= max || throw(DomainError("x = $x must be in [$min,$max]"))
@inline checkinterval(x::Real) = checkinterval(x,-1,1)
@inline cheb0(x::Real) = (checkinterval(x); return one(Float64))
@inline cheb1(x::Real) = (checkinterval(x); return x)
@inline cheb2(x::Real) = (checkinterval(x); return 2*x^2 - 1)
@inline cheb3(x::Real) = (checkinterval(x); return 4*x^3 - 3*x)
@inline cheb4(x::Real) = (checkinterval(x); return 8*(x^4 - x^2) + 1)

# -----------------------------------------
# Flows
# -----------------------------------------

include("flow-payoffs-more.jl")

@inline function flow(::Type{Val{:cheb3_cost_tech}}, wp::AbstractUnitProblem, i::Integer, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T, d::Integer, geoid::Real, roy::T) where {N,T}
# @inline function flow(::Type{Val{:cheb3_cost_tech}}, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T, d::Integer, d1::Integer, Dgt0::Bool, sgn_ext::Bool, geoid::Real, roy::T) where {N,T}
    if d == 0
        _sgnext(wp,i) && return θ[11]
        return zero(T)
    end
    logp, logc, t = z
    u = rev_exp(1,θ[1],1,θ[2],θ[3],θ[4],σ, logp, t, ψ, _Dgt0(wp,i), geoid,roy) + (d==1 ? θ[5] : θ[6] )*cheb0(t) + θ[7]*cheb1(t) + θ[8]*cheb2(t) + θ[9]*cheb3(t) + θ[10]*exp(logc)
    d>1 && (u *= d)
    return u::T
end

@inline function flow(::Type{Val{:cheb3_cost_tech_restr}}, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T, d::Integer, d1::Integer, Dgt0::Bool, sgn_ext::Bool, geoid::Real, roy::T) where {N,T}
    if d == 0
        sgn_ext && return θ[8]
        return zero(T)
    end
    logp, logc, t = z
    u = rev_exp_restricted(θ[1], σ, logp, t, ψ, Dgt0, geoid, roy) + (d==1 ? θ[2] : θ[3] )*cheb0(t) + θ[4]*cheb1(t) + θ[5]*cheb2(t) + θ[6]*cheb3(t) + θ[7]*exp(logc)
    d>1 && (u *= d)
    return u::T
end

# -----------------------------------------
# number of parms
# -----------------------------------------

function number_of_model_parms(FF::Symbol)::Int
    FF ∈ (:one_restr,)                                                           && return  3
    FF ∈ (:dgt1_restr,:Dgt0_restr,)                                              && return  4
    FF ∈ (:one,:dgt1_ext_restr, :dgt1_d1_restr,:dgt1_cost_restr,:cheb2_restr)   && return  5
    FF ∈ (:dgt1,:Dgt0,:dgt1_cost_Dgt0_restr,:dgt1_pricecost_restr,:cheb3_restr) && return  6
    FF ∈ (:dgt1_ext,:dgt1_d1,:dgt1_cost,:dgt1_pricebreak_restr,:cheb2, :cheb3_dgt1_restr,)  && return  7
    FF ∈ (:dgt1_cost_Dgt0,:dgt1_pricecost,:cheb3,:cheb3_cost_restr,:ttbuild_cost_restr,:cheb3_cost_tech_restr,)  && return  8
    FF ∈ (:dgt1_pricebreak,:cheb3_dgt1,)                                         && return  9
    FF ∈ (:cheb3_cost,:ttbuild_cost,)                                            && return 10
    FF ∈ (:cheb3_cost_tech,)                                                     && return 11
    throw(error("FF = $(FF) not recognized"))
end

# -----------------------------------------
# dθ
# -----------------------------------------

@inline function flowdθ(::Type{Val{:cheb3_cost_tech}}, wp::AbstractUnitProblem, i::Integer, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T, k::Integer, d::Integer, geoid::Real, roy::T) where {N,T}
# @inline function flowdθ(::Type{Val{:cheb3_cost_tech}}, θ::AbstractVector{T}, σ::T,     z::NTuple{N,T}, ψ::T, k::Integer,d::Integer,           d1::Integer, Dgt0::Bool, sgn_ext::Bool, geoid::Real, roy::T)::T where {N,T}
    d == 0 && !_sgnext(wp,i) && return zero(T)
    Dgt0 = _Dgt0(wp,i)
    logp, logc, t = z

    # revenue
    k == 1  && return   d * rev_exp(1,θ[1],1,θ[2],θ[3],θ[4],σ,logp,t,ψ,Dgt0,geoid,roy)
    k == 2  && return   d * rev_exp(1,θ[1],1,θ[2],θ[3],θ[4],σ,logp,t,ψ,Dgt0,geoid,roy) * geoid
    k == 3  && return   d * rev_exp(1,θ[1],1,θ[2],θ[3],θ[4],σ,logp,t,ψ,Dgt0,geoid,roy) * ( Dgt0 ? ψ : ψ*_ρ(σ) + θ[k]*(1-_ρ2(σ)))
    k == 4  && return   d * rev_exp(1,θ[1],1,θ[2],θ[3],θ[4],σ,logp,t,ψ,Dgt0,geoid,roy) * t

    # drilling cost
    k == 5   && return d != 1 ? zero(T) : T(   cheb0(t) )
    k == 6   && return d <= 1 ? zero(T) : T( d*cheb0(t) )
    k == 7   && return d == 0 ? zero(T) : T( d*cheb1(t) )
    k == 8   && return d == 0 ? zero(T) : T( d*cheb2(t) )
    k == 9   && return d == 0 ? zero(T) : T( d*cheb3(t) )
    k == 10  && return d == 0 ? zero(T) : T( d*exp(logc) )

    # extension cost
    k == 11  && return d == 0 && _sgnext(wp,i) ? one(T) : zero(T)

    throw(error("$k out of bounds"))
end

@inline function flowdθ(::Type{Val{:cheb3_cost_tech_restr}}, θ::AbstractVector{T}, σ::T,     z::NTuple{N,T}, ψ::T, k::Integer,d::Integer,           d1::Integer, Dgt0::Bool, sgn_ext::Bool, geoid::Real, roy::T)::T where {N,T}
    d == 0 && !sgn_ext && return zero(T)

    logp, logc, t = z

    # revenue
    k == 1  && return   d * rev_exp_restricted(θ[1], σ, logp, t, ψ, Dgt0, geoid, roy)

    # drilling cost
    k == 2  && return d != 1 ? zero(T) : T(   cheb0(t) )
    k == 3  && return d <= 1 ? zero(T) : T( d*cheb0(t) )
    k == 4  && return d == 0 ? zero(T) : T( d*cheb1(t) )
    k == 5  && return d == 0 ? zero(T) : T( d*cheb2(t) )
    k == 6  && return d == 0 ? zero(T) : T( d*cheb3(t) )
    k == 7  && return d == 0 ? zero(T) : T( d*exp(logc) )
    # extension cost
    k == 8  && return d == 0 && sgn_ext ? one(T) : zero(T)

    throw(error("$k out of bounds"))
end

# -----------------------------------------
# dσ
# -----------------------------------------

@inline function flowdσ(::FF, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T, d::Integer, geoid::Real, roy::T)::T where {FF <: Union{ Type{Val{:one}}, Type{Val{:dgt1}}, Type{Val{:Dgt0}}, Type{Val{:dgt1_ext}}, Type{Val{:dgt1_d1}}, Type{Val{:dgt1_cost}}, Type{Val{:dgt1_cost_Dgt0}}, Type{Val{:cheb2}}, Type{Val{:cheb3}}, Type{Val{:cheb3_dgt1}}, Type{Val{:cheb3_cost}}, Type{Val{:ttbuild_cost}} }, N, T}
    d == 0 && return zero(T)
    return d * drevdσ_exp(1,θ[1],1,θ[2],θ[3],σ,z[1],ψ,geoid,roy)
end

@inline function flowdσ(::FF, θ::AbstractVector{T}, σ::T, z::NTuple{N,T},  ψ::T, d::Integer, geoid::Real, roy::T)::T where {FF <: Union{ Type{Val{:one_restr}}, Type{Val{:dgt1_restr}}, Type{Val{:Dgt0_restr}}, Type{Val{:dgt1_ext_restr}}, Type{Val{:dgt1_d1_restr}}, Type{Val{:dgt1_cost_restr}}, Type{Val{:dgt1_cost_Dgt0_restr}}, Type{Val{:cheb2_restr}}, Type{Val{:cheb3_restr}}, Type{Val{:cheb3_dgt1_restr}}, Type{Val{:cheb3_cost_restr}}, Type{Val{:ttbuild_cost_restr}} },N, T}
    d == 0 && return zero(T)
    return d * drevdσ_exp_restricted(θ[1],σ,z[1],ψ,geoid,roy)
end

@inline function flowdσ(::FF, wp::AbstractUnitProblem, i::Integer, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T, d::Integer, geoid::Real, roy::T)::T where {FF <: Union{ Type{Val{:cheb3_cost_tech}} }, N, T}
    d == 0 && return zero(T)
    return d * drevdσ_exp(1,θ[1],1,θ[2],θ[3],θ[4],σ,first(z),last(z),ψ,geoid,roy)
end

@inline function flowdσ(::FF, θ::AbstractVector{T}, σ::T, z::NTuple{N,T},  ψ::T, d::Integer, geoid::Real, roy::T)::T where {FF <: Union{ Type{Val{:cheb3_cost_tech_restr}}  },N, T}
    d == 0 && return zero(T)
    return d * drevdσ_exp_restricted(θ[1],σ,first(z),last(z),ψ,geoid,roy)
end

# -----------------------------------------
# dψ
# -----------------------------------------

@inline function flowdψ(::FF, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T, d::Integer, sgn_ext::Bool, geoid::Real, roy::T) where {FF <: Union{ Type{Val{:one}}, Type{Val{:dgt1}}, Type{Val{:Dgt0}}, Type{Val{:dgt1_d1}}, Type{Val{:dgt1_cost}}, Type{Val{:dgt1_cost_Dgt0}}, Type{Val{:cheb2}}, Type{Val{:cheb3}}, Type{Val{:cheb3_dgt1}}, Type{Val{:cheb3_cost}}, Type{Val{:ttbuild_cost}} },N, T}
    d == 0  && return zero(T) # sgn_ext ? θ[10] : zero(T)
    return (d * drevdψ_exp(1,θ[1],1,θ[2],θ[3],σ,z[1],ψ,geoid,roy))::T
end

@inline function flowdψ(::FF, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T, d::Integer, sgn_ext::Bool, geoid::Real, roy::T) where {FF <: Union{ Type{Val{:one_restr}}, Type{Val{:dgt1_restr}}, Type{Val{:Dgt0_restr}}, Type{Val{:dgt1_d1_restr}}, Type{Val{:dgt1_cost_restr}}, Type{Val{:dgt1_cost_Dgt0_restr}}, Type{Val{:cheb2_restr}}, Type{Val{:cheb3_restr}}, Type{Val{:cheb3_dgt1_restr}}, Type{Val{:cheb3_cost_restr}}, Type{Val{:ttbuild_cost_restr}} },N, T}
    d == 0  && return zero(T) # sgn_ext ? θ[10] : zero(T)
    return (d * drevdψ_exp_restricted(θ[1],σ,z[1],ψ,geoid,roy))::T
end

# @inline function flow(  FF::Type, wp::AbstractUnitProblem, i::Integer, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T,             d::Integer, geoid::Real, roy::T)
# @inline function flowdθ(FF::Type, wp::AbstractUnitProblem, i::Integer, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T, k::Integer, d::Integer, geoid::Real, roy::T)
# @inline function flowdσ(FF::Type, wp::AbstractUnitProblem, i::Integer, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T,             d::Integer, geoid::Real, roy::T)
# @inline function flowdψ(FF::Type, wp::AbstractUnitProblem, i::Integer, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T,             d::Integer, geoid::Real, roy::T)

@inline function flowdψ(::FF, wp::AbstractUnitProblem, i::Integer, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T, d::Integer, geoid::Real, roy::T) where {FF <: Union{ Type{Val{:cheb3_cost_tech}} }, N, T}
    d == 0  && return zero(T) # _sgnext(wp,i) ? θ[10] : zero(T)
    return (d * drevdψ_exp(1,θ[1],1,θ[2],θ[3],θ[4],σ,first(z),last(z),ψ,geoid,roy))::T
end

@inline function flowdψ(::FF, θ::AbstractVector{T}, σ::T, z::NTuple{N,T}, ψ::T, d::Integer, sgn_ext::Bool, geoid::Real, roy::T) where {FF <: Union{ Type{Val{:cheb3_cost_tech_restr}} }, N, T}
    d == 0  && return zero(T) # sgn_ext ? θ[10] : zero(T)
    return (d * drevdψ_exp_restricted(θ[1],σ,first(z),last(z),ψ,geoid,roy))::T
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 = "marble"
version = v"20.12.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/KDE/marble/archive/v20.12.0.tar.gz", "20ea52f071bd255109d723565fc46eb46ee4415acb402f5d45f6c6dbd6b312b7"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd marble-20.12.0/
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/win32.patch"
mkdir build
cd build/
apk del ninja
apk add g++
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release ../.
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "windows"),
    Platform("x86_64", "linux"; libc="glibc"),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libastro", :astro),
    LibraryProduct("libCycleStreetsPlugin", :CycleStreetsPlugin, "plugins"),
    LibraryProduct("libNominatimSearchPlugin", :NominatimSearchPlugin, "plugins"),
    LibraryProduct("libPn2Plugin", :Pn2Plugin, "plugins"),
    LibraryProduct("libHostipPlugin", :HostipPlugin, "plugins"),
    LibraryProduct("libStarsPlugin", :StarsPlugin, "plugins"),
    LibraryProduct("libOpenLocationCodeSearchPlugin", :OpenLocationCodeSearchPlugin, "plugins"),
    LibraryProduct("libLicense", :License, "plugins"),
    LibraryProduct("libCachePlugin", :CachePlugin, "plugins"),
    LibraryProduct("libYoursPlugin", :YoursPlugin, "plugins"),
    LibraryProduct("libGosmoreRoutingPlugin", :GosmoreRoutingPlugin, "plugins"),
    LibraryProduct("libElevationProfileFloatItem", :ElevationProfileFloatItem, "plugins"),
    LibraryProduct("libLocalOsmSearchPlugin", :LocalOsmSearchPlugin, "plugins"),
    LibraryProduct("libMapQuestPlugin", :MapQuestPlugin, "plugins"),
    LibraryProduct("libGeoUriPlugin", :GeoUriPlugin, "plugins"),
    LibraryProduct("libRoutinoPlugin", :RoutinoPlugin, "plugins"),
    LibraryProduct("libJsonPlugin", :JsonPlugin, "plugins"),
    LibraryProduct("libMeasureTool", :MeasureTool, "plugins"),
    LibraryProduct("libFoursquarePlugin", :FoursquarePlugin, "plugins"),
    LibraryProduct("libPostalCode", :PostalCode, "plugins"),
    LibraryProduct("libOsmPlugin", :OsmPlugin, "plugins"),
    LibraryProduct("libmarblewidget-qt5", :marblewidget),
    LibraryProduct("libNavigationFloatItem", :NavigationFloatItem, "plugins"),
    LibraryProduct("libOverviewMap", :OverviewMap, "plugins"),
    LibraryProduct("libEclipsesPlugin", :EclipsesPlugin, "plugins"),
    LibraryProduct("libLocalDatabasePlugin", :LocalDatabasePlugin, "plugins"),
    ExecutableProduct("marble-qt", :marble_qt),
    LibraryProduct("libPntPlugin", :PntPlugin, "plugins"),
    LibraryProduct("libSpeedometer", :Speedometer, "plugins"),
    LibraryProduct("libLogPlugin", :LogPlugin, "plugins"),
    LibraryProduct("libProgressFloatItem", :ProgressFloatItem, "plugins"),
    LibraryProduct("libNotesPlugin", :NotesPlugin, "plugins"),
    LibraryProduct("libAtmospherePlugin", :AtmospherePlugin, "plugins"),
    LibraryProduct("libCompassFloatItem", :CompassFloatItem, "plugins"),
    LibraryProduct("libOpenRouteServicePlugin", :OpenRouteServicePlugin, "plugins"),
    LibraryProduct("libMapScaleFloatItem", :MapScaleFloatItem, "plugins"),
    LibraryProduct("libmarbledeclarative", :marbledeclarative),
    LibraryProduct("libNominatimReverseGeocodingPlugin", :NominatimReverseGeocodingPlugin, "plugins"),
    LibraryProduct("libMonavPlugin", :MonavPlugin, "plugins"),
    LibraryProduct("libPositionMarker", :PositionMarker, "plugins"),
    LibraryProduct("libEarthquakePlugin", :EarthquakePlugin, "plugins"),
    LibraryProduct("libFlightGearPositionProviderPlugin", :FlightGearPositionProviderPlugin, "plugins"),
    LibraryProduct("libGpsInfo", :GpsInfo, "plugins"),
    LibraryProduct("libCrosshairsPlugin", :CrosshairsPlugin, "plugins"),
    LibraryProduct("libGosmoreReverseGeocodingPlugin", :GosmoreReverseGeocodingPlugin, "plugins"),
    LibraryProduct("libRoutingPlugin", :RoutingPlugin, "plugins"),
    LibraryProduct("libGpsbabelPlugin", :GpsbabelPlugin, "plugins"),
    LibraryProduct("libAprsPlugin", :AprsPlugin, "plugins"),
    LibraryProduct("libGraticulePlugin", :GraticulePlugin, "plugins"),
    LibraryProduct("libSatellitesPlugin", :SatellitesPlugin, "plugins"),
    LibraryProduct("libLatLonPlugin", :LatLonPlugin, "plugins"),
    LibraryProduct("libAnnotatePlugin", :AnnotatePlugin, "plugins"),
    LibraryProduct("libGpxPlugin", :GpxPlugin, "plugins"),
    LibraryProduct("libKmlPlugin", :KmlPlugin, "plugins"),
    LibraryProduct("libSunPlugin", :SunPlugin, "plugins"),
    LibraryProduct("libOSRMPlugin", :OSRMPlugin, "plugins"),
    LibraryProduct("libElevationProfileMarker", :ElevationProfileMarker, "plugins")
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Qt_jll", uuid="ede63266-ebff-546c-83e0-1c6fb6d0efc8"))
    Dependency(PackageSpec(name="Libglvnd_jll", uuid="7e76a0d4-f3c7-5321-8279-8d96eeed0f29"))
]

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


text:
module Bonsai

using HTTP, JSON3, StructTypes, Dates, LRUCache, FilePaths
import HTTP: Request, Response, Stream
using Base.Meta: show_sexpr

using ExprManipulation

export ResponseCodes

# we need keep a reference to the original atdoc function
# because atdoc! will overwrite it but i need i later
const expander = Core.atdoc

__init__() = (Core.atdoc!(doc_hook); nothing)

doc_hook(args...) = expander(args...)
# I don't know if this second one is needed but it's here because I adapted the code
# from  https://github.com/JuliaDocs/DocStringExtensions.jl/blob/master/src/templates.jl
# which also had it :/
doc_hook() = Core.atdoc!(expander)

_expr = Ref{Any}(nothing)
_doc_expr = Ref{Any}(nothing)

# this hack let's use document expression that should normally be allowed
# such as app.get("/index") do stream end
function doc_hook(a::LineNumberNode, b::Module, c::String, expr::Expr)

	# @info "args" a=a b=b c=c d=expr
	# @info "typeof" a=typeof(a) b=typeof(b) c=typeof(c) d=typeof(expr)

	# _expr[] = expr

	# TODO: Still need to handle this case so you can add a summary doc string

	# "Some doc string for this"
	# app.summary("/pets/")

	# this won't hanlde the case of assigment although it's unlikely you need that
	# app.summary("/pets/")
	# mexpr = MExpr(:call, MExpr(:. , Capture(:app),  MExpr(:_)), Capture(:route))

	# this won't hanlde the case of assigment although it's unlikely you need that
	# f = app.get("/") do stream end
	mexpr = MExpr(:do, MExpr(:call, Capture(:_), Capture(:route)), Capture(:_))

	m = match(mexpr, expr)

	if !isnothing(m)
		app = expr.args[1].args[1].args[1]
		app = eval(getfield(b, app))
		if app isa App

			method = expr.args[1].args[1].args[2].value
			route = m.route

			# handles cases where route is passed in as a variable
			if !(route isa AbstractString)
				route = eval(getfield(b, route))
			end

			# @info  "route" route=route app= app
			# show_sexpr(expr)
			_expr[] = expr

			v = gensym()

			# we need to modify the expression to assign to a variable
			# so we are able to associate documentation with the handler
			expr = quote
				$v = $expr;
			end

			app.paths_docs[method][route] = c

		end

	end
 	return expander(a,b,c,expr)
end

include("ResponseCodes.jl")
include("mime_types.jl")
include("http.jl")
include("exceptions.jl")
include("json_schema.jl")
include("handlers.jl")
include("io.jl")
include("cancel_token.jl")
include("router.jl")
include("middleware.jl")
include("openapi.jl")
include("app.jl")
include("server.jl")
include("docs.jl")

end


text:
"""
    oocmcca(;input::AbstractArray=[], outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractArray=[], pseudocount::AbstractArray=[], colmeanlist::AbstractArray="", colvarlist::Union{Nothing,AbstractString,AbstractArray}="", dim::Number=3, lower::Number=0, upper::Number=1.0f+38, offsetStoch::Number=1f-20, initW::Union{Nothing,AbstractArray}=nothing, perm::Bool=false)
Online CCA solved by Out-of-core CCA.
Input Arguments
---------
- `input` : Julia Binary file generated by `OnlineCCA.csv2bin` function.
- `outdir` : The directory specified the directory you want to save the result.
- `scale` : {log,ftt,raw}-scaling of the value.
- `pseudocount` : The number specified to avoid NaN by log10(0) and used when `Feature_LogMeans.csv` <log10(mean+pseudocount) value of each feature> is generated.
- `colmeanlist` : The mean of each column of matrix. The CSV file is generated by `OnlineCCA.sumr` functions.
- `colvarlist` : The variance of each column of matrix. The CSV file is generated by `OnlineCCA.sumr` functions.
- `dim` : The number of dimension of CCA.
- `lower` : Stopping Criteria (When the relative change of error is below this value, the calculation is terminated)
- `upper` : Stopping Criteria (When the relative change of error is above this value, the calculation is terminated)
- `offsetStoch` : Off set value for avoding overflow when calculating stochastic gradient
- `initW` : The CSV file saving the initial values of eigenvectors.
- `perm` : Whether the data matrix is shuffled at random.
Output Arguments
---------
- `W` : Eigen vectors of covariance matrix (No. columns of the data matrix × dim)
- `λ` : Eigen values (dim × dim)
- `V` : Latent variables of covariance matrix (No. rows of the data matrix × dim)
- `ExpCorVar` : Explained correlation/variance by the eigenvectors
- `TotalCorVar` : Total correlation/variance of the data matrix
- `Objective` : The summation of object function (e.g. w_i' X_i' X_j w_j)
- `AvgLatentCor` : The mean of correlation efficient of latent variables (e.g. cor(X_i w_i, X_j w_j))
"""
function oocmcca(;input::AbstractArray=[], outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractArray=[], pseudocount::AbstractArray=[], colmeanlist::AbstractArray="", colvarlist::Union{Nothing,AbstractString,AbstractArray}="", dim::Number=3, lower::Number=0, upper::Number=1.0f+38, offsetStoch::Number=1f-20, initW::Union{Nothing,AbstractArray}=nothing, perm::Bool=false)
    # Initial Setting
    cca = OOCMCCA()
    pseudocount, W, v, D, colmeanlist, colmeanvec, colvarlist, colvarvec, N, M, TotalCorVar, offsetStoch = init(input, pseudocount, dim, colmeanlist, colvarlist, initW, cca, offsetStoch, scale)
    # Perform CCA
    out = oocmcca(input, outdir, scale, pseudocount, colmeanlist, colvarlist, dim, cca, W, v, D, colmeanvec, colvarvec, N, M, TotalCorVar, offsetStoch, perm)
    if outdir isa String
        output(outdir, out, cca)
    end
    return out
end

function oocmcca(input, outdir, scale, pseudocount, colmeanlist, colvarlist, dim, cca, W, v, D, colmeanvec, colvarvec, N, M, TotalCorVar, offsetStoch, perm)
    N, M = nm(input)
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    l = length(input)
    p = dim + 5
    its = 3
    for i=1:l
        @assert 0 < dim ≤ l ≤ min(N[i], M[i])
    end
    x = Any[]
    normx = Any[]
    Ω = Any[]
    for i=1:l
        push!(x, zeros(UInt32, M[i]))
        push!(normx, zeros(Float32, M[i]))
        push!(Ω, rand(Float32, M[i], p) * offsetStoch)
    end
    Y = Array{Any,2}(undef, l, l)
    for i = 1:l
        for j = 1:l
            Y[i,j] = zeros(Float32, M[i], p)
        end
    end
    stream = Any[]
    #
    # Random Projection
    #
    println("Random Projection : Y_ij = X_i' * X_j * Ω_j")
    progress = Progress(N[1])
    for n = 1:N[1]
        next!(progress)
        if n == 1
            stream = Any[]
            # Stream
            for i=1:l
                push!(stream, ZstdDecompressorStream(open(input[i])))
            end
            # tmpN, tmpN
            for i=1:l
                read!(stream[i], tmpN)
                read!(stream[i], tmpM)
            end
        else
            # x, norm
            for i=1:l
                read!(stream[i], x[i])
                normx[i] = normalizex(x[i], n, stream[i], scale[i], pseudocount[i], colmeanlist[i], colmeanvec[i], colvarlist[i], colvarvec[i])
                if perm
                    normx[i] .= normx[randperm(length(normx[i]))]
                end
            end
            # Random Projection
            for i=1:l
                for j=1:l
                    Y[i,j] .+= normx[i] * (normx[j]'Ω[j])
                end
            end
            # close
            if n == N[1]
                for i=1:l
                    close(stream[i])
                end
            end
        end
    end

    #
    # Renormalize with LU factorization
    #
    F = Array{Any,2}(undef, l, l)
    for i = 1:l
        for j = 1:l
            tmpF = lu!(Y[i,j])
            F[i,j] = Matrix(tmpF.L)
        end
    end

    for s in 1:its
        println("Renormalization : X_i' * X_j * L_ij")
        progress = Progress(N[1])
        for n = 1:N[1]
            next!(progress)
            if n == 1
                stream = Any[]
                # Stream
                for i=1:l
                    push!(stream, ZstdDecompressorStream(open(input[i])))
                end
                # tmpN, tmpN
                for i=1:l
                    read!(stream[i], tmpN)
                    read!(stream[i], tmpM)
                end
            else
                # x, norm
                for i=1:l
                    read!(stream[i], x[i])
                    normx[i] = normalizex(x[i], n, stream[i], scale[i], pseudocount[i], colmeanlist[i], colmeanvec[i], colvarlist[i], colvarvec[i])
                    if perm
                        normx[i] .= normx[randperm(length(normx[i]))]
                    end
                end
                # Renormalize with LU factorization
                for i=1:l
                    for j=1:l
                        Y[i,j] .+= normx[i]*normx[j]'F[j,i]
                    end
                end
                # close
                if n == N[1]
                    for i=1:l
                        close(stream[i])
                    end
                end
            end
        end
        if s < its
            # Renormalize with LU factorization
            for i = 1:l
                for j = 1:l
                    tmpF = lu!(Y[i,j])
                    F[i,j] = Matrix(tmpF.L[:,1:p])
                end
            end
        else
            # Renormalize with QR factorization
            for i = 1:l
                for j = 1:l
                    tmpF = qr!(Y[i,j])
                    F[i,j] = Matrix(tmpF.Q[:,1:p])
                end
            end
        end
    end
    #
    # Small Matrix
    #
    println("Calculation of small matrix : R")
    R = Array{Any,2}(undef, l, l)
    for i = 1:l
        for j = 1:l
            R[i,j] = ones(Float32, p, p)
        end
    end
    progress = Progress(N[1])
    for n = 1:N[1]
        next!(progress)
        if n == 1
            stream = Any[]
            # Stream
            for i=1:l
                push!(stream, ZstdDecompressorStream(open(input[i])))
            end
            # tmpN, tmpN
            for i=1:l
                read!(stream[i], tmpN)
                read!(stream[i], tmpM)
            end
        else
            # x, norm
            for i=1:l
                read!(stream[i], x[i])
                normx[i] = normalizex(x[i], n, stream[i], scale[i], pseudocount[i], colmeanlist[i], colmeanvec[i], colvarlist[i], colvarvec[i])
                if perm
                    normx[i] .= normx[randperm(length(normx[i]))]
                end
            end
            # Small Matrix
            for i=1:l
                for j=1:l
                    if i != j
                        R[i,j] .+= (F[i,j]'normx[i]) * (normx[j]'F[j,i])
                    end
                end
            end
            # close
            if n == N[1]
                for i=1:l
                    close(stream[i])
                end
            end
        end
    end

    # EVD with small matrix
    println("EVD with small matrix : R W = W Λ")
    newR = zeros(Float32, l*p, l*p)
    xstart = 1
    xend = p
    ystart = 1
    yend = p
    for i=1:l
        for j=1:l
            newR[xstart:xend, ystart:yend] = R[i,j]
            ystart = ystart + p
            yend = yend + p
        end
        ystart = 1
        yend = p
        xstart = xstart + p
        xend = xend + p
    end
    smallW = eigvecs(newR)

    # W[i] = M[i] * dim
    wstart = 1
    wend = p
    for i=1:l
        W[i] = F[i,i] * smallW[wstart:wend, 1:dim]
        wstart = wstart + p
        wend = wend + p
    end

    # Return, W, λ, V
    out = WλV(W, input, dim, scale, pseudocount, colmeanlist, colmeanvec, colvarlist, colvarvec, TotalCorVar)
    return (out[1], out[2], out[3], out[4], out[5], out[6], out[7])
end


text:
@with_kw struct Neuro1DRealSignal
    signal::Matrix{Float64}
    sample_rate::Float64
end

@with_kw struct Neuro1DAnalyticSignal
    signal::Matrix{ComplexF64}
    sample_rate::Float64
end

function to_1ds(x::Matrix{T} where T <: Real; sample_rate::Float64)
    check(x)
    Neuro1DRealSignal(
        signal = convert(Matrix{Float64}, x),
        sample_rate = sample_rate
    )
end


text:
using PLINK
using IHT

# problem dimensions
n = 5000    # number of cases
p = 23999   # number of predictors
k = 10      # number of true predictors
s = 0.1     # standard deviation of noise
T = Float64 # bits type used for SharedArrays

# make simulated float data
srand(2016)         # fix random seed, reproducible x_temp
x_temp = randn(n,p) # simulated data

# make simulated model
b = zeros(p)      # vector to house statistical model
b[1:k] = randn(k) # k random nonzero coefficients drawn from standard normal distribution
shuffle!(b)       # random indices for nonzero coefficients
bidx = find(b)    # save locations of nonzero coefficients

# simulate a noisy response variable
# this makes response vector with noise level drawn from N(0,0.01)
y = x_temp*b + s*randn(n)

# add grand mean to data
X = zeros(n,p+1)             # array to house x_temp + grand mean
setindex!(X, x_temp, :, 1:p) # (efficiently) assign x_temp to correct location in x
X[:,end] = 1.0               # rightmost column (24,000th one) now contains grand mean

# make noisier response variable
s  = 5.0                   # very noisy!
y2 = x_temp*b + s*randn(n) # y2 has noise distribution N(0,25), substantially noiser than y

# simulate (or recover from PLINK.jl) some GWAS data
fpath = expanduser("~/.julia/v0.6/PLINK/data/")         # path to simulated data from PLINK module
genopath  = fpath * "x_test.bed"                        # path to original BED file
tgenopath = fpath * "xt_test.bed"                       # path to transposed BED file
covpath   = fpath * "covfile.txt"                       # path to covariate file
xbed  = BEDFile(genopath, tgenopath, covpath, pids=[1]) # load the data to master process

n,p    = size(xbed)        # dimensions of data
fill!(xbed.covar.x, 1)     # this fills the grand mean with ones
fill!(xbed.covar.xt, 1)    # this fills the *transposed* grand mean with ones; remember to always change both!
mean!(xbed)                # compute means in-place
prec!(xbed)                # compute precisions in-place
xbed.means[end] = 0        # index "end" substitutes for position of grand mean in x.means!
xbed.precs[end] = 1        # same as above

# now simulate model, response with GWAS data
bbed      = SharedArray{T}((p,), pids=[1])   # a model b to use with the BEDFile
bbed[1:k] = randn(k)                         # random coefficients
shuffle!(bbed)                               # random model
bidxbed   = find(bbed)                       # store locations of nonzero coefficients
idx       = bbed .!= 0                       # need BitArray indices of nonzeroes in b for A_mul_B
xb        = A_mul_B(xbed, bbed, idx, k)      # compute x*b
ybed2     = xb + 0.1*randn(n)                # yields a Vector, so we must convert it to SharedVector
ybed      = SharedVector{T}((n,), pids=[1])  # our response variable with the BEDFile
copy!(ybed, ybed2)
ypath     = ENV["TMPDIR"] * "y.bin"          # ybed uses path to Julia TMP directory
write(open(ypath, "w"), ybed)                # "w"rite ybed to file

# crossvalidation parameters
nfolds = 5     # number of crossvalidation folds
srand(2016)    # reset seed before crossvalidation

println("Simulation complete.")


text:
function golden_method(f::Function, a::AbstractVector, b::AbstractVector;
                       tol=eps()*10, maxit=1000)
    α1 = (3 - sqrt(5)) / 2
    α2 = 1 - α1
    d = b - a
    x1 = a + α1*d
    x2 = a + α2*d
    s = ones(a)
    f1 = f(x1)
    f2 = f(x2)

    d = α1*α2*d
    it = 0

    while any(d .> tol) && it < maxit
        it += 1
        i = f2 .> f1
        x1[i] = x2[i]
        f1[i] = f2[i]
        d *= α2
        x2 = x1 + s.*(i- map(!, i)).*d
        s = sign.(x2 .- x1)
        f2 = f(x2)
    end

    it >= maxit && warn("`golden_method`: maximum iterations exceeded")

    i = f2 .> f1
    x1[i] = x2[i]
    f1[i] = f2[i]

    x1, f1
end

# function golden_method(f::Function, a::AbstractVector, b::AbstractVector;
#                        tol=eps()*10, maxit=1000)
#     out_x = similar(a)
#     out_f = similar(a)
#     N = length(a)
#     @inbounds for i=1:N
#         out_x[i], out_f[i] = golden_method(f, a[i], b[i]; tol=tol, maxit=maxit)
#     end
#     out_x, out_f
# end

"""
Applies Golden-section search to search for the _maximum_ of a function in
the interval (a, b)

https://en.wikipedia.org/wiki/Golden-section_search
"""
function golden_method(f::Function, a::Real, b::Real; tol::Float64=10*eps(),
                       maxit::Int=1000)
    α1 = (3 - sqrt(5)) / 2
    α2 = 1 - α1
    d = b - a
    x1 = a + α1*d
    x2 = a + α2*d
    s = 1.0
    f1 = f(x1)::Float64
    f2 = f(x2)::Float64

    d = α1*α2*d

    it = 0

    while d > tol && it < maxit
        it += 1

        d *= α2

        if f2 > f1
            x1 = x2
            f1 = f2
            x2 = x1 + s*d
        else
            x2 = x1 - s*d
        end

        s = sign(x2 - x1)

        f2 = f(x2)::Float64
    end

    it >= maxit && warn("`golden_method`: maximum iterations exceeded")

    if f2 > f1
        x1 = x2
        f1 = f2
    end

    x1, f1::Float64
end


text:
using Random, Distributions

@testset "test MCMC sampler on CIR" begin
    function simulate_CIR_data(;Nsteps_CIR = 200)
        Random.seed!(1)

        δ = 3.
        γ = 2.5
        σ = 2.
        Nobs = 5
        dt_CIR = 0.011
        # max_time = .1
        # max_time = 0.001
        λ = 1.

        time_grid_CIR = [k*dt_CIR for k in 0:(Nsteps_CIR-1)]
        X_CIR = DualOptimalFiltering.generate_CIR_trajectory(time_grid_CIR, 3, δ, γ, σ)
        Y_CIR = map(λ -> rand(Poisson(λ), Nobs), X_CIR);
        data_CIR = Dict(zip(time_grid_CIR, Y_CIR))
        return data_CIR, X_CIR,  δ, γ, σ, λ
    end
    data_CIR, X_CIR,  δ, γ, σ, λ = simulate_CIR_data(;Nsteps_CIR = 40)
    prior_δ = truncated(Normal(5, 4), 0, Inf)
    prior_γ = truncated(Normal(5, 4), 0, Inf)
    prior_σ = truncated(Normal(5, 4), 0, Inf)
    prior_logpdf(δi, γi, σi) = logpdf(prior_δ, δi) + logpdf(prior_γ, γi) + logpdf(prior_σ, σi)

    function unnormalised_logposterior(δ, γ, σ)::Float64
    prior_contribution = prior_logpdf(δ, γ, σ)
        if isinf(prior_contribution)
            return prior_contribution
        else
            return prior_contribution + full_loglikelihood(δ, γ, σ, λ, data_CIR)::Float64
        end
    end

    full_loglikelihood(δ, γ, σ, λ, data, fun::Function)::Float64 = rand(Normal())
    full_loglikelihood(δ, γ, σ, λ, data)::Float64 = full_loglikelihood(δ, γ, σ, λ, data, x -> x)

    unnormalised_logposterior_vec(v) = unnormalised_logposterior(v...)

    @test_nowarn DualOptimalFiltering.estimate_jumping_rule(50, DualOptimalFiltering.Jtnorm_create, [1.,1.,1.], unnormalised_logposterior_vec)
    @test_nowarn DualOptimalFiltering.get_mcmc_samples(15, [1.,1.,1.], DualOptimalFiltering.Jtnorm_create, unnormalised_logposterior_vec; print_acceptance_rate = true)

end;

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
# This file contains code that was formerly part of Julia. License is MIT: http://julialang.org/license

using QuadGK, Test

@testset "quadgk" begin
    @test quadgk(cos, 0,0.7,1)[1] ≈ sin(1)
    @test quadgk(x -> exp(im*x), 0,0.7,1)[1] ≈ (exp(1im)-1)/im
    @test quadgk(x -> exp(im*x), 0,1im)[1] ≈ -1im*expm1(-1)
    @test isapprox(quadgk(cos, 0,BigFloat(1),order=40)[1], sin(BigFloat(1)),
                   atol=1000*eps(BigFloat))
    @test quadgk(x -> exp(-x), 0,0.7,Inf)[1] ≈ 1.0
    @test quadgk(x -> exp(x), -Inf,0)[1] ≈ 1.0
    @test quadgk(x -> exp(-x^2), -Inf,Inf)[1] ≈ sqrt(pi)
    @test quadgk(x -> [exp(-x), exp(-2x)], 0, Inf)[1] ≈ [1,0.5]
    @test quadgk(cos, 0,0.7,1, norm=abs)[1] ≈ sin(1)
end

module Test19626
    using QuadGK
    using Test

    # Define a mock physical quantity type
    struct MockQuantity <: Number
        val::Float64
    end

    # Following definitions needed for quadgk to work with MockQuantity
    import Base: +, -, *, abs, isnan, isinf, isless
    +(a::MockQuantity, b::MockQuantity) = MockQuantity(a.val+b.val)
    -(a::MockQuantity, b::MockQuantity) = MockQuantity(a.val-b.val)
    *(a::MockQuantity, b::Number) = MockQuantity(a.val*b)
    abs(a::MockQuantity) = MockQuantity(abs(a.val))
    isnan(a::MockQuantity) = isnan(a.val)
    isinf(a::MockQuantity) = isinf(a.val)
    isless(a::MockQuantity, b::MockQuantity) = isless(a.val, b.val)

    # isapprox only needed for test purposes
    Base.isapprox(a::MockQuantity, b::MockQuantity) = isapprox(a.val, b.val)

    # Test physical quantity-valued functions
    @test QuadGK.quadgk(x->MockQuantity(x), 0.0, 1.0, atol=MockQuantity(0.0))[1] ≈
        MockQuantity(0.5)
end


text:
include("../required.jl")

##Change PredTSN to newID starting from 1
# Sort by ID
sort!(edgelist_T, :PredTSN)
# Get unique IDs
unique_ids = unique(edgelist_T.PredTSN)
#create new IDs
newid_pred = collect(1:length(unique_ids))
# Create Dict where id => newid
id_dict = Dict(unique_ids[i] => newid_pred[i] for i in 1:length(unique_ids))
# Add new ids in dataframe
edgelist_T.newPredID = [id_dict[id] for id in edgelist_T.PredTSN]

### Do same for prey
## Change PreyTSN to new ID starting from end of newPredID
# Sort by ID
sort!(edgelist_T, :PreyTSN)
# Get unique IDs
unique_ids_prey = unique(edgelist_T.PreyTSN)
# Create new IDs
newid_prey = collect(length(id_dict)+1:length(id_dict)+length(unique_ids_prey))
# Create Dict
id_dict2 = Dict(unique_ids_prey[i] => newid_prey[i] for i in 1:length(unique_ids_prey))
# Add new ids in dataframe
edgelist_T.newPreyID = [id_dict2[id] for id in edgelist_T.PreyTSN]

transect = edgelist_T[:, [:newPredID, :newPreyID]]

#sort
sort!(transect, :newPredID)

#get unique values
unique!(transect)

max_stars_repo_name: kiristern/reseaux

max_issues_repo_name: kiristern/reseaux

max_forks_repo_name: kiristern/reseaux


text:
export Configuration,ReadConfig

function ReadConfig(filename::String)
    input = open(filename, "r")
    inputs=JSON.parse(input)
    close(input)
    return inputs
end

mutable struct Molecule
    type::Array{Int64,1}
    mass::Array{Float64,1}
    coords::Array{Float64,2}
    bonds::Array{Int64,2}
    angles::Array{Int64,2}
    vels::Array{Float64,2}
end

#function Molecule(type::Array{Int64,1},
#    mass::Array{Float64,1},
#    coords::Array{Float64,2},
#    bonds::Array{Array{Int64,2},1},
#    angles::Array{Array{Int64,2},1},
#    vels::Array{Float64,2})
#new(
#type::Array{Int64,1},
#mass::Array{Float64,1},
#coords::Array{Float64,2},
#bonds::Array{Array{Int64,2},1},
#angles::Array{Array{Int64,2},1},
#vels::Array{Float64,2}
#)
#end

struct Configuration
    start::Int64
    box::Array{Float64,1}
    moles::Array{Molecule,1}
    moles_length::Array{Int64,1}
    function Configuration(data::Dict)
        start=1
        box=zeros(3)
        if "box" in keys(data)
            box=convert(Array{Float64,1}, data["box"])
            @info "Box Lx Ly Lz: $(box[1]) $(box[2]) $(box[3])"
        end

        moles=Vector{Molecule}(undef,0)
        moles_length=zeros(Int64,0)

        if "molecules" in keys(data)
            molecules=data["molecules"]
            num_moles=length(keys(molecules))

            moles=Vector{Molecule}(undef,num_moles)

            moles_length=zeros(Int64,num_moles)
            for molei in 1:num_moles
                moles_length[molei]=length(keys(molecules[string(molei)]["atoms"]))
                #molecule=Molecule(moles_length[molei])
                type=zeros(Int64,moles_length[molei])
                mass=zeros(Int64,moles_length[molei])
                bonds=zeros(Int64,0,3)
                angles=zeros(Int64,moles_length[molei],4)
                coords=zeros(Float64,moles_length[molei],3)
                vels=zeros(Float64,moles_length[molei],3)
                for atomi in 1:moles_length[molei]
                    type[atomi]=molecules[string(molei)]["atoms"][string(atomi)]["type"]
                    mass[atomi]=molecules[string(molei)]["atoms"][string(atomi)]["mass"]
                    coords[atomi,1:3]=convert(Array{Float64,1},molecules[string(molei)]["atoms"][string(atomi)]["coords"])
                    if "vels" in keys(molecules[string(molei)]["atoms"][string(atomi)])
                        vels[atomi,1:3]=convert(Array{Float64,1},molecules[string(molei)]["atoms"][string(atomi)]["vels"])
                    end
                end
                if "bonds" in keys(molecules[string(molei)])
                    num_bonds=length(keys(molecules[string(molei)]["bonds"]))
                    bonds=zeros(Int64,num_bonds,3)
                    for bondi in 1:num_bonds
                        bonds[bondi,1:3]=convert(Array{Int64,1},molecules[string(molei)]["bonds"][string(bondi)])
                    end
                end
                moles[molei]=Molecule(type,mass,coords,bonds,angles,vels)
            end
        end

        new(start,box,
            moles,
            moles_length)
    end
end


text:
# we collect the list of unit test files here

include("TestHybridQuasiGaussian.jl")
include("TestHybridPerformance.jl")

max_stars_repo_name: sschlenkrich/HybridMonteCarlo

max_issues_repo_name: sschlenkrich/HybridMonteCarlo

max_forks_repo_name: sschlenkrich/HybridMonteCarlo


text:
using CSV, TimeSeries, Dates, Statistics, Measures, Plots, StatsPlots; pyplot()

df = CSV.read("../data/oneOnEpsilonBlogs.csv",copycols = true)

tsA = TimeArray(Date.(df.Day,Dates.DateFormat("m/d/y")),df.Users)
tsB = moving(mean,tsA,7,padding = true)
tsC = TimeArray(timestamp(tsA), values(tsA) - values(tsB))

dow = dayofweek.(timestamp(tsA));
dayDiv = [filter((x)->!isnan(x) && x >= -50 && x <= 50, values(tsC)[dow .== d])
            for d in 1:7]

start2020 = findfirst((d)->Year(d)==Year(20),timestamp(tsA))
indexLast = length(tsA)
indexes2020 = start2020:indexLast

dayGroup1 = [7] #Sun
dayGroup2 = [1,6] #Mon,Sat
dayGroup3 = [2,3,4,5] #Tue, Wed, Thu, Fri
dayGroups = [dayGroup1, dayGroup2, dayGroup3]
groupDivs = [vcat(dayDiv[g]...) for g in dayGroups]

default(legend = :topleft)
labels = ["Daily" "7 Day Average"]
p1 = plot(1:indexLast, [values(tsA) values(tsB)], c=[:blue :red],label = labels,
        xlabel = "Day", ylabel = "Daily Users", ylim = (0,200))

p2 = plot(indexes2020, [values(tsA)[indexes2020] values(tsB)[indexes2020]],
        c=[:blue :red], label=labels, xlabel="Day in 2020",
        ylabel="Daily Users",ylim=(0,200))

p3 = plot(1:indexLast, values(tsC),label = labels, c=:black,
        xlabel="Day",ylabel="Variation",ylim=(-50,50),legend=false)

p4 = plot(indexes2020, values(tsC)[indexes2020],label = labels, c=:black,
         xlabel="Day in 2020",ylabel="Variation",ylim=(-50,50),legend=false)

dayNames = dayname.(timestamp(tsA)[4:10])
p5 = density(dayDiv, label = hcat(dayNames...),legend = :topright,
                xlabel = "Variation", ylabel = "Frequency",)

dayGroupNames = ["Sun", "Mon+Sat", "Tue+Wed+Thu+Fri"]
p6 = density(groupDivs,label = hcat(dayGroupNames...), legend = :topright,
                xlabel = "Variation", ylabel = "Frequency",c=[:blue :red :green])

plot(p1, p2, p3, p4, p5, p6, layout=(3,2), size = (900,600), margin = 5mm)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
# https://fenicsproject.org/docs/dolfin/latest/python/demos/mixed-poisson/demo_mixed-poisson.py.html

module mixed_poisson

using FenicsPy

# Create mesh
mesh = UnitSquareMesh(32, 32)

# Define finite elements spaces and build mixed space
BDM = FiniteElement("BDM", mesh.ufl_cell(), 1)
DG  = FiniteElement("DG", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, BDM * DG)

# Define trial and test functions
(sigma, u) = TrialFunctions(W)
(tau, v) = TestFunctions(W)

# Define source function
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
g = Expression("sin(5*x[0])", degree=2)

# Define variational form
a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx
L = - f*v*dx

# Define function G such that G \cdot n = g
py"""
from math import sin
from dolfin import Cell, UserExpression
class BoundarySource(UserExpression):
    def __init__(self, mesh, **kwargs):
        self.mesh = mesh
        super().__init__(**kwargs)
    def eval_cell(self, values, x, ufc_cell):
        cell = Cell(self.mesh, ufc_cell.index)
        n = cell.normal(ufc_cell.local_facet)
        g = sin(5*x[0])
        values[0] = g*n[0]
        values[1] = g*n[1]
    def value_shape(self):
        return (2,)
"""

G = Expression(py"BoundarySource"(mesh.pyobject, degree=2))

# Define essential boundary
#def boundary(x):
#    return x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS
boundary = "x[1] < $DOLFIN_EPS || x[1] > 1.0 - $DOLFIN_EPS"

bc = DirichletBC(W.sub(0), G, boundary)

# Compute solution
w = FeFunction(W)
solve(a == L, w, bc)
(sigma, u) = w.split()

# Save solution to file in VTK format
vtkfile = File("mixed_poisson/sigma.pvd")
vtkfile << sigma

vtkfile = File("mixed_poisson/u.pvd")
vtkfile << u

# Plot sigma and u
#plot(sigma)
#plot(u)

end # module mixed_poisson


text:
## BFS.m is a best first search of the ColorGraph problem.
#Begin the choices list
choices=A[:,:].^2+B[:,:].^2+C[:,:].^2
options=choices[7,:]
if n==1
    ## if this is the first node added
    #Generate scores for each first choice
    best=find(options == minimum(options))
    # Put it all together and pass back as ADDED
    ADDED=[n,0,best]
else()
    ## add later nodes
    # What are the possible node placements
    OriginOptions=unique(rainbow[:,3])
    OriginOptions=OriginOptions[find(OriginOptions~=0)]
    for L=1:length(OriginOptions)
        options=[options;choices[OriginOptions[L],:]]
    end
    # which find row and column of best
    [row,col]=find(choices == minimum(options))
    if row==7
        ADDED=[n,0,col]
    else()
        ADDED=[n,row,col]
    end
end
#add to rainbow
rainbow[n,:]=ADDED
clear seeded L OriginOptions ADDED best seeded row col

max_stars_repo_name: DesignEngrLab/ColorGraph

max_issues_repo_name: DesignEngrLab/ColorGraph

max_forks_repo_name: DesignEngrLab/ColorGraph


text:
abstract type AbstractStep <: AbstractPolicy end

# function stepsize(step::AbstractStep,k::Integer,fval::Real,x::AbstractVector,g::AbstractVector)
#     error("No defined step function for step policy ", typeof(step))
# end

function stepsize() end

function step_wrapper(klocal::Cint, kglobal::Cint, fval::Cdouble, xbegin::Ptr{Cdouble},
                      xend::Ptr{Cdouble}, gbegin::Ptr{Cdouble},
                      step_data::Ptr{Cvoid})
    step_policy = unsafe_pointer_to_objref(step_data)::AbstractStep
    ptrdiff = Int(xend - xbegin)
    N = divrem(ptrdiff, sizeof(Cdouble))[1]
    x = unsafe_wrap(Array, xbegin, N)
    g = unsafe_wrap(Array, gbegin, N)
    return stepsize(step_policy, klocal, kglobal, fval, x, g)
end

module Step

using LinearAlgebra
using Parameters
using POLO: AbstractStep, AbstractPolicyParameters
import POLO: initialize!, stepsize

function initialize!(policy::AbstractStep,x₀::Vector{Float64})
    nothing
end

include("constant.jl")
include("decreasing.jl")
include("bb.jl")

end


text:
#=
    Variable bounds
        SingleVariable -in- LessThan
        SingleVariable -in- GreaterThan
        SingleVariable -in- EqualTo
        SingleVariable -in- Interval

        SingleVariable -in- ZeroOne
        SingleVariable -in- Integer
        SingleVariable -in- Semiinteger
        SingleVariable -in- Semicontinuous
=#
constrdict(model::LinQuadOptimizer, ::SVCI{LE}) = cmap(model).upper_bound
constrdict(model::LinQuadOptimizer, ::SVCI{GE}) = cmap(model).lower_bound
constrdict(model::LinQuadOptimizer, ::SVCI{EQ}) = cmap(model).fixed_bound
constrdict(model::LinQuadOptimizer, ::SVCI{IV}) = cmap(model).interval_bound

constrdict(model::LinQuadOptimizer, ::SVCI{MOI.ZeroOne}) = cmap(model).binary
constrdict(model::LinQuadOptimizer, ::SVCI{MOI.Integer}) = cmap(model).integer

constrdict(model::LinQuadOptimizer, ::SVCI{MOI.Semicontinuous{Float64}}) = cmap(model).semicontinuous
constrdict(model::LinQuadOptimizer, ::SVCI{MOI.Semiinteger{Float64}}) = cmap(model).semiinteger

"""
    change_both_variable_bounds!(model::LinQuadOptimizer, columns::Vector{Int},
        lower_bounds::Vector{Float64}, upper_bounds::Vector{Float64})

Set the lower bound of column `columns[i]` to `lower_bounds[i]` and the upper
bound to `upper_bounds[i]`. Alternatively, the lower or upper bound can be left
blank by passing an array of length 0 instead.

Examples:
    change_both_variable_bounds!(model, [1, 2], [-0.5, 0.0], [1.0, 2.0])
    change_both_variable_bounds!(model, [1, 2], [-0.5, 0.0], Float64[])
    change_both_variable_bounds!(model, [1, 2], Float64[], [1.0, 2.0])
"""
function change_both_variable_bounds!(
        model::LinQuadOptimizer,
        columns::Vector{Int},
        lower_bounds::Vector{Float64},
        upper_bounds::Vector{Float64})
    if length(lower_bounds) > 0 && length(upper_bounds) > 0
        columns = vcat(columns, columns)
    end
    change_variable_bounds!(model,
        columns,
        vcat(lower_bounds, upper_bounds),
        vcat(
            fill(backend_type(model, Val{:Lowerbound}()), length(lower_bounds)),
            fill(backend_type(model, Val{:Upperbound}()), length(upper_bounds))
        )
    )
end

function set_variable_bound(model::LinQuadOptimizer, v::SinVar, set::LE)
    change_both_variable_bounds!(
        model, [get_column(model, v)], Float64[], [set.upper])
end

function set_variable_bound(model::LinQuadOptimizer, v::SinVar, set::GE)
    change_both_variable_bounds!(
        model, [get_column(model, v)], [set.lower], Float64[])
end

function set_variable_bound(model::LinQuadOptimizer, v::SinVar, set::EQ)
    change_both_variable_bounds!(
        model, [get_column(model, v)], [set.value], [set.value])
end

function set_variable_bound(model::LinQuadOptimizer, v::SinVar, set::IV)
    change_both_variable_bounds!(
        model, [get_column(model, v)], [set.lower], [set.upper])
end

function set_variable_bounds(
        model::LinQuadOptimizer, vars::Vector{SinVar}, sets::Vector{LE})
    change_both_variable_bounds!(model, get_column.(model, vars), Float64[],
        [set.upper for set in sets])
end

function set_variable_bounds(
        model::LinQuadOptimizer, vars::Vector{SinVar}, sets::Vector{GE})
    change_both_variable_bounds!(model, get_column.(model, vars),
        [set.lower for set in sets], Float64[])
end

function set_variable_bounds(
        model::LinQuadOptimizer, vars::Vector{SinVar}, sets::Vector{EQ})
    values = [set.value for set in sets]
    change_both_variable_bounds!(model, get_column.(model, vars), values, values)
end

function set_variable_bounds(
        model::LinQuadOptimizer, vars::Vector{SinVar}, sets::Vector{IV})
    change_both_variable_bounds!(model, get_column.(model, vars),
        [set.lower for set in sets], [set.upper for set in sets])
end

"""
    has_value(dict::Dict{K, V}, value::V) where {K, V}

Return true if `dict` has a `key=>value` pair with `value`.
"""
function has_value(dict::Dict{K, V}, value::V) where {K, V}
    return value in values(dict)
end

function MOI.add_constraint(model::LinQuadOptimizer, variable::SinVar, set::S) where S <: LinSets
    __assert_supported_constraint__(model, SinVar, S)
    variable_type = model.variable_type[variable.variable]
    if variable_type == SEMI_CONTINUOUS || variable_type == SEMI_INTEGER
        error("Cannot set bounds because variable is of type: $(variable_type).")
    end
    set_variable_bound(model, variable, set)
    model.last_constraint_reference += 1
    index = SVCI{S}(model.last_constraint_reference)
    dict = constrdict(model, index)
    dict[index] = variable.variable
    return index
end

function MOI.add_constraints(model::LinQuadOptimizer, variables::Vector{SinVar},
                             sets::Vector{S}) where S <: LinSets
    __assert_supported_constraint__(model, SinVar, S)
    for variable in variables
        variable_type = model.variable_type[variable.variable]
        if variable_type == SEMI_CONTINUOUS || variable_type == SEMI_INTEGER
            error("Cannot set bounds because variable is of type: $(variable_type).")
        end
    end
    set_variable_bounds(model, variables, sets)
    indices = SVCI{S}[]
    for variable in variables
        model.last_constraint_reference += 1
        index = SVCI{S}(model.last_constraint_reference)
        dict = constrdict(model, index)
        dict[index] = variable.variable
        push!(indices, index)
    end
    return indices
end

function MOI.delete(model::LinQuadOptimizer, index::SVCI{S}) where S <: LinSets
    __assert_valid__(model, index)
    delete_constraint_name(model, index)
    dict = constrdict(model, index)
    variable = dict[index]
    model.variable_type[variable] = CONTINUOUS
    set_variable_bound(model, SinVar(variable), IV(-Inf, Inf))
    delete!(dict, index)
    return
end

# constraint set
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintSet, index::SVCI{LE})
    MOI.LessThan{Float64}(
        get_variable_upperbound(model, get_column(model, model[index]))
    )
end
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintSet, index::SVCI{GE})
    MOI.GreaterThan{Float64}(
        get_variable_lowerbound(model, get_column(model, model[index]))
    )
end
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintSet, index::SVCI{EQ})
    MOI.EqualTo{Float64}(
        get_variable_lowerbound(model, get_column(model, model[index]))
    )
end
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintSet, index::SVCI{IV})
    column = get_column(model, model[index])
    return MOI.Interval{Float64}(get_variable_lowerbound(model, column),
                                 get_variable_upperbound(model, column))
end

# constraint function
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintFunction, index::SVCI{<: LinSets})
    return SinVar(model[index])
end

# modify
MOI.supports(::LinQuadOptimizer, ::MOI.ConstraintSet, ::Type{SVCI{S}}) where S<:LinSets = true
function MOI.set(model::LinQuadOptimizer, ::MOI.ConstraintSet, index::SVCI{S}, newset::S) where S <: LinSets
    set_variable_bound(model, SinVar(model[index]), newset)
end

#=
    Binary constraints

For some reason CPLEX doesn't respect bounds on a binary variable, so we should
store the previous bounds so that if we delete the binary constraint we can
revert to the old bounds.

Xpress is worse, once binary, the bounds are changed independently of what the
user does.
=#
function MOI.add_constraint(model::LinQuadOptimizer, variable::SinVar, set::MOI.ZeroOne)
    __assert_supported_constraint__(model, SinVar, MOI.ZeroOne)
    variable_type = model.variable_type[variable.variable]
    if variable_type != CONTINUOUS
        error("Cannot make variable binary because it is $(variable_type).")
    end
    model.variable_type[variable.variable] = BINARY
    model.last_constraint_reference += 1
    index = SVCI{MOI.ZeroOne}(model.last_constraint_reference)
    column = get_column(model, variable)
    dict = constrdict(model, index)
    dict[index] = (variable.variable, get_variable_lowerbound(model, column),
                   get_variable_upperbound(model, column))
    change_variable_types!(model, [column], [backend_type(model, set)])
    change_variable_bounds!(model,
        [column, column],
        [0.0, 1.0],
        [backend_type(model, Val{:Lowerbound}()),
         backend_type(model, Val{:Upperbound}())]
    )
    make_problem_type_integer(model)
    return index
end

function MOI.delete(model::LinQuadOptimizer, index::SVCI{MOI.ZeroOne})
    __assert_valid__(model, index)
    delete_constraint_name(model, index)
    dict = constrdict(model, index)
    (variable, lower, upper) = dict[index]
    model.variable_type[variable] = CONTINUOUS
    column = get_column(model, variable)
    change_variable_types!(
        model, [column], [backend_type(model, Val{:Continuous}())])
    change_variable_bounds!(model,
       [column, column],
       [lower, upper],
       [backend_type(model, Val{:Lowerbound}()),
        backend_type(model, Val{:Upperbound}())]
    )
    delete!(dict, index)
    if !has_integer(model)
        make_problem_type_continuous(model)
    end
end

function MOI.get(::LinQuadOptimizer, ::MOI.ConstraintSet, ::SVCI{MOI.ZeroOne})
    return MOI.ZeroOne()
end

function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintFunction, index::SVCI{MOI.ZeroOne})
    return SinVar(model[index][1])
end

#=
    Integer constraints
=#

function MOI.add_constraint(model::LinQuadOptimizer, variable::SinVar, set::MOI.Integer)
    __assert_supported_constraint__(model, SinVar, MOI.Integer)
    variable_type = model.variable_type[variable.variable]
    if variable_type != CONTINUOUS
        error("Cannot make variable integer because it is $(variable_type).")
    end
    model.variable_type[variable.variable] = INTEGER
    change_variable_types!(model, [get_column(model, variable)],
                           [backend_type(model, set)])
    model.last_constraint_reference += 1
    index = SVCI{MOI.Integer}(model.last_constraint_reference)
    dict = constrdict(model, index)
    dict[index] = variable.variable
    make_problem_type_integer(model)
    return index
end

function MOI.delete(model::LinQuadOptimizer, index::SVCI{MOI.Integer})
    __assert_valid__(model, index)
    delete_constraint_name(model, index)
    dict = constrdict(model, index)
    variable = dict[index]
    model.variable_type[variable] = CONTINUOUS
    change_variable_types!(model, [get_column(model, variable)],
                           [backend_type(model, Val{:Continuous}())])
    delete!(dict, index)
    if !has_integer(model)
        make_problem_type_continuous(model)
    end
end

function MOI.get(::LinQuadOptimizer, ::MOI.ConstraintSet, ::SVCI{MOI.Integer})
    return MOI.Integer()
end

function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintFunction,
                 index::SVCI{MOI.Integer})
    return SinVar(model[index])
end

#=
    Semicontinuous / Semiinteger constraints
=#
const SEMI_TYPES = Union{MOI.Semicontinuous{Float64}, MOI.Semiinteger{Float64}}
variable_type_(::Type{<:MOI.Semicontinuous}) = SEMI_CONTINUOUS
variable_type_(::Type{<:MOI.Semiinteger}) = SEMI_INTEGER
function MOI.add_constraint(model::LinQuadOptimizer, variable::SinVar, set::S) where S <: SEMI_TYPES
    __assert_supported_constraint__(model, SinVar, S)
    variable_type = model.variable_type[variable.variable]
    if variable_type != CONTINUOUS
        error("Cannot make variable $(S) because it is $(variable_type).")
    end
    model.variable_type[variable.variable] = variable_type_(S)
    column = get_column(model, variable)
    change_variable_types!(model, [column], [backend_type(model, set)])
    change_variable_bounds!(model,
        [column, column],
        [set.lower, set.upper],
        [backend_type(model, Val{:Lowerbound}()),
         backend_type(model, Val{:Upperbound}())]
    )
    model.last_constraint_reference += 1
    index = SVCI{S}(model.last_constraint_reference)
    dict = constrdict(model, index)
    dict[index] = variable.variable
    make_problem_type_integer(model)
    return index
end

function MOI.delete(model::LinQuadOptimizer, index::SVCI{<:SEMI_TYPES})
    __assert_valid__(model, index)
    delete_constraint_name(model, index)
    dict = constrdict(model, index)
    variable = dict[index]
    column = get_column(model, variable)
    model.variable_type[variable] = CONTINUOUS
    change_variable_types!(model, [column], [backend_type(model, Val{:Continuous}())])
    change_variable_bounds!(model,
        [column, column],
        [-Inf, Inf],
        [backend_type(model, Val{:Lowerbound}()),
         backend_type(model, Val{:Upperbound}())]
    )
    delete!(dict, index)
    if !has_integer(model)
        make_problem_type_continuous(model)
    end
end

function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintSet, index::SVCI{S}) where S <: SEMI_TYPES
    dict = constrdict(model, index)
    column = get_column(model, dict[index])
    return S(get_variable_lowerbound(model, column),
             get_variable_upperbound(model, column))
end

function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintFunction, index::SVCI{<:SEMI_TYPES})
    return SinVar(model[index])
end


text:
@inline function get_interpolation(wp::Interpolations.Extrapolation{T, 3}, pt::CartesianPoint{T}, S::Val{:cartesian})::T where {T <: SSDFloat}
    return wp(pt.x, pt.y, pt.z)::T
end
@inline function get_interpolation(wp::Interpolations.Extrapolation{T, 3}, pt::CartesianPoint{T}, S::Val{:cylindrical})::T where {T <: SSDFloat}
    p::CylindricalPoint{T} = CylindricalPoint(pt)
    return wp(p.r, p.φ, p.z)::T
end

function add_signal!(signal::Vector{T}, path::Vector{CartesianPoint{T}}, charge::T, wp::Interpolations.Extrapolation{T, 3}, S::Union{Val{:cylindrical}, Val{:cartesian}})::Nothing where {T <: SSDFloat}
    @inbounds for i in eachindex(signal)
        signal[i] += get_interpolation(wp, path[i], S)::T * charge
    end
    nothing
end
function add_signal!(signal::Vector{T}, path::DriftPath{T}, charge::T, wp::Interpolations.Extrapolation{T, 3}, S::Union{Val{:cylindrical}, Val{:cartesian}})::Nothing where {T <: SSDFloat}
    add_signal!(signal, path.e_path, -charge, wp, S) # electrons induce negative charge
    add_signal!(signal, path.h_path,  charge, wp, S)
    nothing
end
function add_signal!(signal::Vector{T}, paths::Vector{DriftPath{T}}, charges::Vector{T}, wp::Interpolations.Extrapolation{T, 3}, S::Union{Val{:cylindrical}, Val{:cartesian}})::Nothing where {T <: SSDFloat}
    for ipath in eachindex(paths)
        add_signal!(signal, paths[ipath], charges[ipath], wp, S)
    end
    nothing
end

function get_signal(path::Vector{CartesianPoint{T}}, charge::T, wp::Interpolations.Extrapolation{T, 3}, S::Union{Val{:cylindrical}, Val{:cartesian}})::Vector{T} where {T <: SSDFloat}
    signal::Vector{T} = zeros(T, length(path))
    add_signal!(signal, path, charge, wp, S)
    return signal
end