text:
abstract type AbstractTranslocation end
abstract type AbstractDissipativeTranslocation <: AbstractTranslocation end
abstract type AbstractLosslessTranslocation <: AbstractTranslocation end

@mix @columns @flattenable struct Trans{F}
    κtra::F        | 0.6     | _   | Beta(2.0, 2.0)  | [0.0,1.0]      | _ | "Reserve flux allocated to translocation"
end

@mix @columns @flattenable struct MultiTrans{D,P}
    destnames::D   | false | (:leaf,) | _ | _              | _         | _ | "The organ/s translocated to"
    proportions::P | true  | (1.0,)   | _ | Beta(2.0, 2.0) | [0.0,1.0] | _ | "The proportion of translocation sent in the first translocation. Only for inetermediaries. nothing = 100%"
end

@mix @columns @flattenable struct DissTrans{MoMo}
    y_E_ET::MoMo   | true | 0.8       | mol*mol^-1      | Beta(2.0, 2.0)  | [0.0,1.0]    | _ | "Translocated reserve:"
end

# Dissipation types
@Trans @MultiTrans @DissTrans struct DissipativeMultipleTranslocation{} <: AbstractDissipativeTranslocation end
@Trans @DissTrans struct DissipativeTranslocation{} <: AbstractDissipativeTranslocation  end
@Trans @MultiTrans struct LosslessMultipleTranslocation{} <: AbstractLosslessTranslocation end
@Trans struct LosslessTranslocation{} <: AbstractLosslessTranslocation end

abstract type AbstractRejection end

@columns struct DissipativeRejection{MoMo} <: AbstractRejection
    y_EC_ECT::MoMo       | 1.0             | mol*mol^-1      | Beta(2.0, 2.0)  | [0.0,1.0] | _ | "Translocated C-reserve"
    y_EN_ENT::MoMo       | 1.0             | mol*mol^-1      | Beta(2.0, 2.0)  | [0.0,1.0] | _ | "Translocated N-reserve"
end
struct LosslessRejection <: AbstractRejection end

# Recurse through all organs. A loop would not be type-stable.
# translocation!(organs::Tuple, destorgans::Tuple) = begin
#     props = buildprops(organs[1])
#     translocation!(organs[1], destorgans, organs[1].params.translocation.destnames, props)
#     translocation!(tail(organs), destorgans)
# end
# translocation!(organs::Tuple{}, destorgans::Tuple) = nothing
# translocation!(organ::Organ, destorgans::Tuple, destnames::Symbol, props) =
#     translocation!(organ, destorgans, (destnames,), props)
# # Translocate to organs with names in the destnames list
# translocation!(organ::Organ, destorgans::Tuple, destnames, props) = begin
#     for i = 1:length(destnames)
#         if destorgans[1].params.name == destnames[i]
#             reuse_rejected!(organ, destorgans[1], props[i])
#             translocate!(organ, destorgans[1], props[i])
#             break
#         end
#     end
#     translocation!(organ, tail(destorgans), destnames, props)
# end
# translocation!(organ::Organ, destorgans::Tuple{}, destnames, props) = nothing

# Add the last remainder proportion (so that its not a model parameter)
# buildprops(o::Organ) = buildprops(o.params.translocation.proportions)
# buildprops(x::Nothing) = (1.0)
# buildprops(x::Number) = (x, 1 - x)
# buildprops(xs::Tuple) = (xs..., 1 - sum(xs))

"""
Versions for E, CN and CNE reserves.

Translocation is occurs between adjacent organs.
This function is identical both directiono, and ox represents
whichever is not the current organs. Will not run with less than 2 organs.
"""
translocate!(o1, o2, prop) = translocate!(trans_pars(o1), o1, o2, prop)
translocate!(p::Nothing, o1, o2, prop) = nothing
translocate!(p::AbstractDissipativeTranslocation, o1, o2, prop) = begin
    # outgoing translocation
    trans = κtra(o1) * o1.J1[:E,:ctb]
    reserve_drain!(o1, Val(:tra), trans)

    # incoming translocation
    transx = κtra(o2) * o2.J1[:E,:ctb]
    o1.J[:E,:tra] += transx * y_E_ET(o2)

    # loss = transx * (1 - y_E_ET(o2))
    # reserve_loss!(o2, loss)
    # conversion_loss!(o2, transx * y_E_ET(o2), n_N_E(o2))
    nothing
end

translocate!(p::AbstractLosslessTranslocation, o1, o2, prop) = begin
    # outgoing translocation
    trans = κtra(o1) * o1.J1[:E,:ctb]
    reserve_drain!(o1, Val(:tra), trans)

    # incoming translocation
    transx = κtra(o2) * o2.J1[:E,:ctb]
    o1.J[:E,:tra] += transx

    # conversion_loss!(o2, transx, n_N_E(o2))
    nothing
end

"""
Translocation occurs between adjacent organs.
This function is identical both directiono, so on represents
whichever is not the current organs.

Will not run with less than 2 organs.
"""
translocation!(organs::Tuple{T1,T2}) where {T1, T2} = begin
    translocate_rejected!(organs[1], organs[2], 1.0)
    translocate_rejected!(organs[2], organs[1], 1.0)
    translocate!(organs[1], organs[2], 1.0)
    translocate!(organs[2], organs[1], 1.0)
end
# translocation!(organs::Tuple) = translocation!(organs...)
translocation!(organs::Tuple{T}) where T = nothing
translocation!(organs::Tuple{}) = nothing

"""
Reallocate state rejected from synthesizing units.
TODO: add a 1-organs method. How does this interact with assimilation?
"""
translocate_rejected!(source, dest, prop) = translocate_rejected!(rejection_pars(source), source, dest, prop)
translocate_rejected!(f::Nothing, source, dest, prop) = nothing
translocate_rejected!(f::DissipativeRejection, source, dest, prop) = begin
    Js, J1s, Jd = flux(source), flux1(source), flux(dest)
    transC = J1s[:C,:rej] # * (1 - κEC(o))
    transN = J1s[:N,:rej] # * (1 - κEN(o))
    Js[:C,:rej] = -transC
    Js[:N,:rej] = -transN
    Jd[:C,:tra] = f.y_EC_ECT * transC
    Jd[:N,:tra] = f.y_EN_ENT * transN
    # J1[:C,:los] += transC * (1 - y_EC_ECT(o)) + transN * (1 - y_EN_ENT(o))
    # J1[:N,:los] += (transC * (1 - y_EC_ECT(o)), transN * (1 - y_EN_ENT(o))) ⋅ (n_N_EC(o), n_N_EN(o))
    nothing
end
translocate_rejected!(rejected::LosslessRejection, source, dest, prop) = begin
    Js, J1s, Jd = flux(source), flux1(source), flux(dest)
    transC = J1s[:C,:rej] # * (1 - κEC(o))
    transN = J1s[:N,:rej] # * (1 - κEN(o))
    Js[:C,:rej] = -transC
    Js[:N,:rej] = -transN
    Jd[:C,:tra] = transC
    Jd[:N,:tra] = transN
    nothing
end

κtra(trans_pars::AbstractTranslocation) = trans_pars.κtra


text:
module BitHelper
include("bit.jl")
export bit
end # module

ext: jl


text:
#!/usr/bin/env julia

using SphericalFourierBesselDecompositions
SFB = SphericalFourierBesselDecompositions

using Test

@testset "NDIterators" begin
    # test iteration
    N = 4
    Lmin = fill(0, N)
    Lmax = rand(1:9, N)
    ndi = SFB.NDIterator(Lmin, Lmax)
    for L1=0:Lmax[1], L2=0:Lmax[2], L3=0:Lmax[3], L4=0:Lmax[4]
        @test SFB.advance(ndi) == true
        @test ndi[1] == L1
        @test ndi[2] == L2
        @test ndi[3] == L3
        @test ndi[4] == L4
    end
    @test SFB.advance(ndi) == false

    # max < min should give no iteration
    N = 4
    Lmin = [2, 3, 4, 5]
    Lmax = [2, 3, 3, 5]
    L = SFB.NDIterator(Lmin, Lmax)
    @test SFB.advance(L) == false

    # test construction
    N = 4
    Lmin = fill(0, N)
    Lmax = rand(1:9, N)
    L1 = SFB.NDIterator(Lmin, Lmax)
    M1 = SFB.NDIterator(-L1, L1)
    M2 = SFB.NDIterator(0, L1)
    M3 = SFB.NDIterator(L1, 10)
    M4 = SFB.NDIterator(0, 10)
    M5 = SFB.NDIterator(0, 10; N=5)
    @test length(L1) == N
    @test length(M1) == N
    @test length(M2) == N
    @test length(M3) == N
    @test length(M4) == 1
    @test length(M5) == 5

    # test broadcast
    L = SFB.NDIterator(0, 10, N=6)
    SFB.advance(L)  # make L valid
    arr = @. L + 1
    @test arr == fill(1, 6)

    # test array conversion
    @test convert(Array, L) == fill(0,6)
    #Array(L)  # doesn't work

    # firstindex, lastindex
    @test L[begin] == 0
    @test L[end] == 0
end

# vim: set sw=4 et sts=4 :


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

module MolecularGraphUtilTest
    using Test
    using MolecularGraph.MolecularGraphUtil
    include("./util/iterator.jl")
end

module MolecularGraphModelTest
    using Test
    using MolecularGraph.MolecularGraphModel

    include("./graph/dfs.jl")
    include("./graph/shortestpath.jl")
    include("./graph/generator.jl")

    include("./graph/plaingraph.jl")
    include("./graph/plaindigraph.jl")
    include("./graph/plainhypergraph.jl")

    include("./graph/multigraph.jl")
    include("./graph/dag.jl")
    include("./graph/cycle.jl")
    include("./graph/connectivity.jl")
    include("./graph/inducedsubgraph.jl")
    include("./graph/disjointunion.jl")
    include("./graph/linegraph.jl")

    include("./graph/triangle.jl")
    include("./graph/clique.jl")
    include("./graph/bipartite.jl")
    include("./graph/planarity.jl")
    include("./graph/product.jl")
    include("./graph/isomorphism/vf2.jl")
    include("./graph/isomorphism/cliquebased.jl")
end

module MolecularGraphGeometryTest
    using Test
    using LinearAlgebra
    using MolecularGraph.MolecularGraphGeometry

    include("./geometry/coords2d.jl")
    include("./geometry/coords3d.jl")
    include("./geometry/coordsinternal.jl")
end

module MolecularGraphTest
    using Test
    using MolecularGraph
    using MolecularGraph.MolecularGraphGeometry
    using MolecularGraph.MolecularGraphModel

    include("./model/atom.jl")
    include("./model/molgraph.jl")

    include("sdfilereader.jl")
    include("sdfilewriter.jl")

    include("./smarts/base.jl")
    include("./smarts/logicaloperator.jl")
    include("./smarts/atom.jl")
    include("./smarts/bond.jl")
    include("./smarts/molecule.jl")
    include("./smarts/smiles.jl")
    include("./smarts/smarts.jl")

    include("properties.jl")
    include("preprocess.jl")
    include("wclogp.jl")
    include("substructure.jl")
    include("funcgroup.jl")
    include("mcs.jl")

    include("./draw/base.jl")
    include("./draw/svg.jl")
end


text:
export Platte

struct Platte{T<:AbstractFloat, Ti<:Integer}
    xL::T
    yL::T
    nx::Ti
    ny::Ti
    rho::Matrix{T}
    lambda::Matrix{T}
    cp::Matrix{T}
    temp::Matrix{T}
end

function Platte(xL::T, yL::T, nx::Ti, ny::Ti, dens, heat, cap, tv::Union{Number,Function}) where {T,Ti}
    if nx <= 0
        throw(ArgumentError("nx muss > 0 sein, ist aber $nx"))
    end

        function setgrid(f::Number)
            return [f for ix in 1:nx, iy in 1:ny]
        end

        function setgrid(f::Function)
            return [f((ix - 0.5) * xL / nx, (iy - 0.5) * xL / ny) for ix in 1:nx, iy in 1:ny]
        end

    rho = setgrid(dens)
    temp = setgrid(tv)
    lambda= setgrid(heat)
    cp = setgrid(cap)

    return Platte{T, Ti}(xL, yL, nx, ny, rho, lambda, cp, temp)
end

function tempiter(p::Platte, deltat::T) where {T}

    dtemp = [heattransfer(ix, iy, p) for ix in 1:p.nx, iy in 1:p.ny]
    tempn = dtemp * deltat + p.temp

    return Platte(p.xL, p.yL, p.nx, p.ny, p.rho, p.lambda, p.cp, tempn)
end

function heattransfer(ix, iy, p::Platte)

    produkt = [p.lambda[ix, iy] * nabla(p.temp, p)[ix, iy] for ix in 1:p.nx, iy in 1:p.ny]

    tempdiff = nabla(produkt, p)[ix, iy] / (p.cp[ix, iy] * p.rho[ix, iy])

end

function nabla(vm::Matrix{<:Vector}, p::Platte)
    vmx = [vm[ix, iy][1] for ix in 1:p.nx, iy in 1:p.ny]
    vmy = [vm[ix, iy][2] for ix in 1:p.nx, iy in 1:p.ny]
    [nablax(vmx, ix, iy, p) + nablay(vmy, ix, iy, p) for ix in 1:p.nx, iy in 1:p.ny]
end

function nabla(m::Matrix{<:Number}, p::Platte)
    [[nablax(m, ix, iy, p), nablay(m, ix, iy, p)] for ix in 1:p.nx, iy in 1:p.ny]
end

function nablax(m::Matrix, ix, iy, p::Platte)
    dx = p.xL / p.nx
    ix == 1 ? (m[ix + 1, iy] - m[ix, iy]) / dx :
    ix == p.nx ? (m[ix, iy] - m[ix - 1, iy]) / dx :
    (m[ix + 1, iy] - m[ix - 1, iy]) / (2 * dx)
end

function nablay(m::Matrix, ix, iy, p::Platte)
    dy = p.yL / p.ny
    iy == 1 ? (m[ix, iy + 1] - m[ix, iy]) / dy :
    iy == p.ny ? (m[ix, iy] - m[ix, iy - 1]) / dy :
    (m[ix, iy + 1] - m[ix, iy - 1]) / (2 * dy)
end


text:
using Distributed
addprocs(16)
@everywhere include("get_objective.jl")
get_Javg_vs_s(2)

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
module Layers
    using ...HyperSphere
    import ..Initializer
    import ...AbstractObject

    #Initialize all pre allocated arrays
    @Fun(LayerInitializer, Function, InputArray::AbstractArray, OutputArray::AbstractArray, Constants::AbstractArray)

    export LayerDesign

    struct LayerDesign
        sig::Signature
        input_shape
        output_shape
        constant_shape
        const_initializer
        layer_init::LayerInitializer

        LayerDesign(sig, input_shape, output_shape, constant_shape,
            const_initializer, layer_init) = new(sig, input_shape, output_shape, constant_shape, const_initializer, layer_init)
    end

    #Designed for passing around a layer function generator from a given input size
    #Should be able to make multiple copies (ie the copies are independent of one another)
    @Fun(LayerGenerator, LayerDesign, ProgramBuilder::ProgramBuilder, Settings::MemorySettings, Input_Shape::Any)

    include("Dense.jl")
    include("Convolutional.jl")
    include("Softmax.jl")
    include("Reshape.jl")
    include("ElementWise.jl")
    include("MatrixWise.jl")
end


text:
import JuliaDB: tuplesetindex

@testset "Utilities" begin
    @testset "tuplesetindex" begin
        @test tuplesetindex((1,2,3), 4, 2) == (1,4,3)
        @test tuplesetindex((x=1,y=2,z=3), 4, 2) == (x=1,y=4,z=3)
        @test tuplesetindex((x=1,y=2,z=3), 4, :y) == (x=1,y=4,z=3)
    end
end

import JuliaDB: Interval, hasoverlap

@testset "Interval" begin

    @testset "hasoverlap" begin
        @test hasoverlap(Interval(0,2), Interval(1,2))
        @test !(hasoverlap(Interval(2,1), Interval(1,2)))
        @test hasoverlap(Interval(0,2), Interval(0,2))
        @test hasoverlap(Interval(0,2), Interval(-1,2))
        @test hasoverlap(Interval(0,2), Interval(2,3))
        @test !hasoverlap(Interval(0,2), Interval(4,5))
    end

    @testset "in" begin
        @test 1 in Interval(0, 2)
        @test !(1 in Interval(2, 1))
        @test !(3 in Interval(0, 2))
        @test !(-1 in Interval(0, 2))
    end

end

@testset "vcat PooledArray Vector{Union{T,Missing}}" begin
    a = PooledArray(["x"])
    b = Vector{Union{String,Missing}}(["y"])
    c = vcat(a, b)
    @test eltype(c) == Union{String,Missing}
    @test c == ["x", "y"]
end


text:
mutable struct Linear
    W
    b
end

Linear() = Linear(rand(1,1), rand(1))

function linear_cost(m::Linear, x, y)
    sum((m.W * x .+ m.b .- y) .^ 2)
end

function linear_cost_flat(W, b, x, y)
    sum((W * x .+ b .- y) .^ 2)
end

@testset "struct" begin
    W = randn(500, 784)
    b = randn(500)
    m = Linear(W, b)
    x = randn(784, 100)
    y = randn(500)
    inputs = [:m => m, :x => x, :y => y]
    flat_inputs = [:W => W, :b => b, :x => x, :y => y]

    # first check that df with flat model is correct
    test_compare(linear_cost_flat; flat_inputs...)

    # now check that results from flat and structured models are the same
    df_flat = xdiff(linear_cost_flat; flat_inputs...)
    r, dW_, db_, dx_, dy_ = df_flat(W, b, x, y)

    df = xdiff(linear_cost; inputs...)
    r, dm, dx, dy = df(m, x, y)

    @test isapprox(dW_, dm.W)
    @test isapprox(db_, dm.b)
    @test isapprox(dx_, dx)
    @test isapprox(dy_, dy)
end


text:
#xortraintest-3.jl

#trains an xornet using noisy and unreliable data.  Ten data inputs, top two xor'ed to get the
#correct values, 5% of the time.
function xortrain_3()

  println("working on unreliable xor data set with swarm optimization")

  input_matrix = rand(Bool, 10, 500)
  training_results = Array{Bool,2}(1,500)
  training_results[:] = [input_matrix[1, col] $ input_matrix[2, col] $ (rand() < 0.05) for col in 1:size(input_matrix,2)]

  xornet = GenML.MLP.MultilayerPerceptron{Float64,(10,2,1)}()

  GenML.MLP.randomizelayers(xornet, () -> randn())

  GenML.Optimizers.swarmoptimize(xornet, input_matrix, training_results, cost)

  #verify that the optimization has resulted in a good data set.
  wrongcount = 0
  for x = 1:50
    input_vector = rand(Bool, 10)
    wrongcount += (xornet(input_vector)[1] > 0.5) != (input_vector[1] $ input_vector[2])
  end
  println("incorrect responses for unreliable xor: $wrongcount")

  return wrongcount
end

xortrain_3() == 0 || xortrain_3() == 0 || throw(ErrorException("two executions failed."))

max_stars_repo_name: interplanetary-robot/GenML

max_issues_repo_name: interplanetary-robot/GenML

max_forks_repo_name: interplanetary-robot/GenML


text:
"""
    GSDist <: ContinuousUnivariateDistribution

Generalized S-Distribution for approximating univariate distributions.
"""
struct GSDist <: ContinuousUnivariateDistribution
    F₀::Real
    x₀::Real
    α::Real
    g::Real
    k::Real
    γ::Real
    dist::Union{UnivariateDistribution, Nothing} # underlying distribution
    function GSDist(F₀, x₀, α, g, k, γ, dist)
        F₀ < 0 || F₀ > 1 && throw(DomainError(F₀, "F₀ must be between 0 and 1"))
        α ≤ 0 && throw(DomainError(α, "α must be a non-negative real number"))
        g ≤ 0 && throw(DomainError(g, "g must be a non-negative real number"))
        k ≤ 0 && throw(DomainError(k, "k must be a non-negative real number"))
        γ ≤ 0 && throw(DomainError(γ, "γ must be a non-negative real number"))

        new(F₀, x₀, α, g, k, γ, dist)
    end
end

function GSDist(dist::UnivariateDistribution, F₀::Real=0.5;
    n::Int=21, diff::Function=three_point_midpoint, h::Real=1.0)
    F₀ < 0 || F₀ > 1 && throw(DomainError(F₀, "F₀ must be between 0 and 1"))
    F₀ = Float64(F₀)

    l,u = extrema(dist)
    if isinf(l) l = quantile(dist, 1e-6) end
    if isinf(u) u = quantile(dist, 1 - 1e-6) end
    ql, qu = cdf.(dist, (l,u))

    q = range(ql, qu, length=n)
    X = typeof(dist) <: DiscreteUnivariateDistribution ? Set(quantile.(dist, q)) : quantile.(dist, q)

    FX = cdf.(dist, X)
    dF = y->diff(x->cdf(dist,x), y, h)
    fX = typeof(dist) <: DiscreteUnivariateDistribution ? dF.(X) : pdf.(dist, X)

    f = (t, p) -> p[1] * t.^p[2] .* (1.0 .- t.^p[3]).^p[4]

    p0 = [1, 0.5, 1, 0.5]
    fit = curve_fit(f, FX, fX, p0)
    p = coef(fit)

    GSDist(F₀, quantile(dist, F₀), p..., dist)
end

function Base.show(io::IO, ::MIME"text/plain", G::GSDist)
    α, g, k, γ = params(G)
    println(io, "Generalized S-Distribution ($(G.dist))")
    println(io, " α: $α")
    println(io, " g: $g")
    println(io, " k: $k")
    print(io, " γ: $γ")
end

params(G::GSDist) = (G.α, G.g, G.k, G.γ)

function quantile(D::GSDist, p::Real)
    q = _gsd_quantile(p, D.F₀, D.x₀, D.α, D.g, D.k, D.γ)
    if isnan(q)
        @debug "Unable to calculate the quantile of the GSDist. Falling back to the quantile of the underlying distribution."
        quantile(D.dist, p)
    else
        q
    end
end

function mean(G::GSDist)
    m = _gsd_mean(G.F₀, G.x₀, G.α, G.g, G.k, G.γ)
    if isnan(m)
        @debug "Unable to calculate the mean of the GSDist. Falling back to the mean of the underlying distribution."
        mean(G.dist)
    else
        m
    end
end

function var(G::GSDist)
    m1 = mean(G)
    m2 = try
        _gsd_moment(G.F₀, G.x₀, G.α, G.g, G.k, G.γ, moment=2)
    catch y
        if isa(y, DomainError)
            @debug "Unable to calculate the variance of the GSDist. Falling back to the variance of the underlying distribution."
            return var(G.dist)
        end
    end
    m2 - m1^2
end

std(G::GSDist) = sqrt(var(G))


text:
using LinearAlgebra

function sparseApprox(A, b, k, opts...)
    if length(opts) == 0
        ε = .1
        σ = -2.5
        doprint = false
        print_debug = false
        info_frac = zeros(k)
        for i = 1 : k
            info_frac[i] = (1 / (k - (i - 1)))
        end
    elseif length(opts) == 1
        ε = opts[1]
        σ = -2.5
        doprint = false
        print_debug = false
        info_frac = zeros(k)
        for i = 1 : k
            info_frac[i] = (1 / (k - (i - 1)))
        end
    elseif length(opts) == 2
        ε = opts[1]
        σ = opts[2]
        doprint = false
        print_debug = false
        info_frac = zeros(k)
        for i = 1 : k
            info_frac[i] = (1 / (k - (i - 1)))
        end
    elseif length(opts) == 3
        ε = opts[1]
        σ = opts[2]
        doprint = opts[3]
        print_debug = false
        info_frac = zeros(k)
        for i = 1 : k
            info_frac[i] = (1 / (k - (i - 1)))
        end
    elseif length(opts) == 4
        ε = opts[1]
        σ = opts[2]
        doprint = opts[3]
        print_debug = opts[4]
        info_frac = zeros(k)
        for i = 1 : k
            info_frac[i] = (1 / (k - (i - 1)))
        end
    elseif length(opts) > 4
        ε = opts[1]
        σ = opts[2]
        doprint = opts[3]
        print_debug = opts[4]
        info_frac = opts[5]
    end

    n, numtar = size(A)
    info_bound = info_frac[1] * (1 - ε)
    can_use_prelim = falses(numtar)
    max_expl = zeros(Float64, numtar)
    prohib = Array{Any, 1}(nothing, k)
    prohib_counter = zeros(Int64, k)
    improvement = .03
    counter_prelim = 0;
    bestres = 1
    bestsol = zeros(numtar)
    iter = 1
    δ = [0.1, 0.1]
    tol = 0.001
    lambda = 0.01
    maxiter_gd = 100
    pre_al = 100000
    prohib[1] = zeros(Int64, numtar)
    #combinations of size 1 and 2 are treated as base cases
    tar_norm = norm(b, 1)
    if doprint
        println("Checking combinations of size: 1")
        println("Maximum residual allowed: ", 1 - info_bound)
    end
    for i = 1 : numtar
        tmp = ladConstrained_1d(b, A[:, i], σ)
        tmpsol = zeros(numtar)
        tmpsol[i] = tmp[1]
        tmpres = norm(b - A * tmpsol, 1) / norm(b, 1)
        tmprawres = norm(A * tmpsol, 1)
        max_expl[i] = 1 - tmpres
        if tmpres < bestres
            bestres = tmpres
            bestsol .= tmpsol
        end
        if max_expl[i] >= info_bound
            if doprint
                println(i)
                println("Residual = ", tmpres)
            end
            can_use_prelim[i] = true
            counter_prelim = counter_prelim + 1
            prohib_counter[1] += 1
            prohib[1][prohib_counter[1]] = i
        end
    end

    cmblist = zeros(Int64, counter_prelim)
    iter_cmb = 1
    for i = 1 : numtar
        if can_use_prelim[i]
            cmblist[iter_cmb] = i
            iter_cmb += 1
        end
    end
    outer_loop_cap = length(cmblist)
    prohib[2] = zeros(Int64, pre_al, 2)
    can_use_phase2 = Array{Any, 1}(undef, outer_loop_cap)
    for i = 1 : outer_loop_cap
        can_use_phase2[i] = falses(numtar)
    end
    max_expl_phase2 = zeros(Float64, pre_al)

    S = zeros(n, 2)
    info_bound = info_frac[2] * (1 - ε)
    if doprint
        println("# of valid 1-combinations : ", outer_loop_cap)
        println("Checking combinations of size: 2")
        println("Maximum residual allowed: ", 1 - info_bound)
    end
    for i = 1 : outer_loop_cap
        if doprint
            println("Combination# ", i)
        end

        for j = 1 : numtar
            if cmblist[i] != j && !can_use_phase2[i][j] && !isprohib_v4([cmblist[i] j], prohib, 2, prohib_counter[2]) &&
                (max_expl[j] + max_expl[cmblist[i]]) >= info_bound

                S[:, 1] .= A[:, cmblist[i]]
                S[:, 2] .= A[:, j]

                tmp = ladFitGd_constrained_v2(S, b, σ, tol, δ, maxiter_gd)
                tmpsol = zeros(numtar)
                tmpsol[cmblist[i]] = tmp[1]
                tmpsol[j] = tmp[2]
                tmpres = norm(b - A * tmpsol, 1) / norm(b, 1)

                if tmpres < bestres
                    if doprint
                        println([cmblist[i] j])
                        println("Residual = ", tmpres)
                    end
                    bestres = tmpres
                    bestsol .= tmpsol
                end
                if 1 - (tmpres) >= info_bound
                    can_use_phase2[i][j] = true
                    prohib_counter[2] += 1
                    max_expl_phase2[prohib_counter[2]] = 1 - tmpres
                    prohib[2][prohib_counter[2], 1] = cmblist[i]
                    prohib[2][prohib_counter[2], 2] = j
                end
            end
        end
    end

    k_iter = 2

    if k > k_iter
        cmblist_prev = similar(cmblist)
        cmblist_prev[:] .= cmblist[:]
        can_use_prev = can_use_phase2
        max_expl_prev = max_expl_phase2
    end

    prohib_old = prohib[2]
    num_prohibited_old = 0
    while k_iter < k
        k_iter += 1
        counter_loop = Int(0)
        cmblist_rec = zeros(Int64, pre_al, k_iter - 1)
        iter_cmb = Int(1)
        if k_iter == 3
            for i = 1 : length(can_use_prev)
                for j = 1 : numtar
                    if can_use_prev[i][j]
                        cmblist_rec[iter_cmb, :] .= [cmblist_prev[i, :] ; j]
                        iter_cmb += 1
                    end
                end
            end
        else
            c1, c2 = size(can_use_prev)
            for i = 1 : c1
                for j = 1 : c2
                    if can_use_prev[i, j]
                        cmblist_rec[iter_cmb, :] .= [cmblist_prev[i, :] ; j]
                        iter_cmb += 1
                    end
                end
            end
        end
        num_cmb = iter_cmb - 1

        prohib_mat = zeros(Int64, pre_al, k_iter)
        can_use_tmp = falses(num_cmb, numtar)
        max_expl_rec = zeros(Float64, pre_al)
        S = zeros(n, k_iter)
        info_bound = max(info_frac[k_iter] * (1 - ε), (1 - bestres) + (1 - bestres) * improvement)
        if doprint
            println("# of valid ", k_iter - 1, "-combinations : ", num_cmb)
            println("Checking combinations of size: ", k_iter)
            println("Maximum residual allowed: ", 1 - info_bound)
        end

        for j = 1 : num_cmb
            if doprint
                println("Combination# ", j)
            end

            for l = 1 : numtar
                idxisallowed = true
                for q = 1 : k_iter - 1
                    if l == cmblist_rec[j, q]
                        idxisallowed = false
                        break
                    end
                end

                flag = true
                idxcount = 0

                for ii = 1 : counter_loop
                    for jj = 1 : k_iter
                        for kk = 1 : k_iter - 1
                            if cmblist_rec[ii, kk] == prohib_mat[ii, jj]
                                idxcount += 1
                            end
                        end
                        if l == prohib_mat[ii, jj] && idxcount == k_iter - 1
                            flag = false
                            break
                        end
                    end
                end

                if idxisallowed && !can_use_tmp[j, l] && flag && (max_expl[l] + max_expl_prev[j]) >= info_bound
                    S[:, 1] .= A[:, cmblist_rec[j, 1]]
                    for q = 2 : k_iter - 1
                        S[:, q] .= A[:, cmblist_rec[j, q]]
                    end
                    S[:, k_iter] .= A[:, l]
                    tmp = ladFitGd_constrained_v2(S, b, σ, tol, δ, maxiter_gd)
                    tmpsol = zeros(numtar)
                    for q = 1 : k_iter - 1
                        tmpsol[cmblist_rec[j, q]] = tmp[q]
                    end
                    tmpsol[l] =  tmp[end]

                    tmpres = norm(b - A * tmpsol, 1) / norm(b, 1)

                    if tmpres < bestres
                        if doprint
                            for q = 1 : k_iter - 1
                                print(cmblist_rec[j, q], ", ")
                            end
                            println(l)
                            println("Residual = ", tmpres)
                        end
                        bestres = tmpres
                        bestsol .= tmpsol
                    end

                    if 1 - tmpres >= info_bound
                        can_use_tmp[j, l] = true
                        counter_loop += 1
                        if counter_loop > pre_al
                            pre_al *= 2
                            tmpprohib = prohib_mat
                            tmpmaxepl = max_expl_rec
                            prohib_mat = zeros(Int64, pre_al, k_iter)
                            max_expl_rec = zeros(Float64, pre_al)
                            prohib_mat[1:counter_loop - 1, :] .= tmpprohib
                            max_expl_rec[1:counter_loop - 1, :] .= tmpmaxepl
                        end
                        max_expl_rec[counter_loop] = 1 - tmpres
                        prohib_mat[counter_loop, 1:k_iter-1] .= cmblist_rec[j, :]
                        prohib_mat[counter_loop, end] = l
                    end
                else
                    if print_debug && cmblist_rec[j, end] == 219 && l == 307
                        for q = 1 : k_iter - 1
                            print(cmblist_rec[j, q], ", ")
                        end
                        println(l, ", excluded from search")
                        println("Is redundant? ", !idxisallowed)
                        println("Is prohibited? ", isprohib_v4([cmblist_rec[j, :] ; l], prohib_mat, k_iter, counter_loop))
                        println("Info Bound = ", info_bound)
                        println("Upper Bound = ", (1 - max_expl[l]))
                    end
                end

            end
        end
        cmblist_prev = cmblist_rec
        can_use_prev = can_use_tmp
        max_expl_prev = max_expl_rec
        prohib_old = prohib_mat
        num_prohibited_old = counter_loop
    end

    return bestsol, bestres
end;

function ladConstrained_1d(x, y, ε)
    coeff = y \ x
    cy = y * coeff
    r = x - cy

    if length(size(r)) > 1
        r = dropdims(r, dims=2)
    end

    sp = sortperm(r)

    return (x[sp[1]] - ε) / y[sp[1]]
end;

function ladFitGd_constrained_v2(A, b, ε, opts...)
    dim_list = size(A)
    if length(dim_list) > 1
        m = dim_list[1]
        n = dim_list[2]
    else
        m = dim_list[1]
        n = 1
    end
    if length(opts) == 4
        tol = opts[1]
        lambda = opts[2]
        maxiter = opts[3]
        doprint = opts[4]
    elseif length(opts) == 3
        tol = opts[1]
        lambda = opts[2]
        maxiter = opts[3]
        doprint = false
    elseif length(opts) == 2
        tol = opts[1]
        lambda = opts[2]
        maxiter = 1000
        doprint = false
    elseif length(opts) == 1
        tol = opts[1]
        lambda = 0.01
        maxiter = 1000
        doprint = false
    else
        tol = 0.001
        lambda = 0.01
        maxiter = 1000
        doprint = false
    end

    prev_step::Float64 = Float64(1)
    if length(dim_list) > 1
        cur_x = ones(size(A)[2], 1)
    else
        error("use ladConstrained_1d for 1 dimensional x")
    end
    δ = [.1, .1]

    if all(b .- A * cur_x .> ε) && all(cur_x .> 0)
        initres = norm(A*cur_x .- b, 1) + δ[1] * sum(-log.(cur_x)) + δ[2] * sum(-log.(b .- A * cur_x .- ε))
    else
        initres = Inf
    end
    if doprint
        println("initial guess = ", cur_x)
    end
    while isinf(initres)
        cur_x = cur_x .* .5

        if all(b .- A * cur_x .> ε) && all(cur_x .> 0)
            initres = norm(A*cur_x .- b, 1) + δ[1] * sum(-log.(cur_x)) + δ[2] * sum(-log.(b .- A * cur_x .- ε))
        else
            initres = Inf
        end
    end
    if doprint
        println("ε = ", ε)
        println("feasible starting point = ", cur_x)
        println("l1 norm = ", norm(A*cur_x .- b, 1))
        println("2nd term = ", δ[1] * sum(-log.(cur_x)))
        println("3rd term = ", δ[2] * sum(-log.(b .- A * cur_x .- ε)))
        println("initres = ", initres)
    end
    iter = 0
    lambda = .5
    pg = zeros(size(cur_x))
    t1 = zeros(size(cur_x))
    t2 = zeros(size(cur_x))
    t3 = zeros(size(cur_x))
    while iter < maxiter && prev_step > tol
        if all(b .- A * cur_x .> ε) && all(cur_x .> 0)
            res = norm(A*cur_x .- b, 1) + δ[1] * sum(-log.(cur_x)) + δ[2] * sum(-log.(b .- A * cur_x .- ε))
        else
            res = Inf
        end

        t1 .= transpose(A) * sign.(A * cur_x .- b)
        t2[:] .= 0
        if all(cur_x .> 0)
            t2 .= δ[1] * (-1 ./ cur_x)
        elseif !any(cur_x .< 0)
            for i = 1 : size(t2)[1]
                if cur_x[i] > 0
                    t2[i] = δ[1] * (-1 ./ cur_x[i])
                end
            end
        end
        t3[:] .= 0
        for i = 1 : size(A)[1]
            numer = A[i, :]
            denom = b[i] - dot(A[i, :], cur_x) - ε
            if denom > 0
                t3 .+= (numer ./ denom)
            end
        end
        g = t1 + t2 + (δ[2] * t3)
        if doprint
            println("res = ", res)
            println("t1 = ", t1)
            println("t2 = ", t2)
            println("t3 = ", t3)
            println("grad = ", g)
        end

        cur_x = cur_x - (g .* lambda)
        prev_step = norm(g, 2)

        iter = iter + 1
        normdiff = abs(norm(pg + g, 2))
        if doprint
            println("prev_step = ", prev_step)
            println("pg = ", pg)
            println("g = ", g)
            println("diff = ", pg + g)
            println("norm diff = ", normdiff)
            println("norm diff scaled = ", normdiff / lambda)
            println("iter = ", iter)
        end

        if normdiff < .125
             lambda *= .8
        end
        pg .= g
        if doprint
            println("lambda = ", lambda)
            println("cur_x = ", cur_x)
        end
    end

    return dropdims(cur_x, dims = 2)
end;

function isprohib_v4(idx, prohib, level, len)
    l = length(idx)
    if l == 1
        for i = 1 : len
            for j = 1 : level
                if idx == prohib[i, j]
                    return true
                end
            end
        end
    else
        for i = 1 : len
            counter = 0
            for j = 1 : level
                for k = 1 : l
                    if idx[k] == prohib[i, j]
                        counter += 1
                    end
                end
            end
            if counter == level
                return true
            end
        end
    end
    return false
end;

max_stars_repo_name: mushinmohan/NNSA

max_issues_repo_name: mushinmohan/NNSA

max_forks_repo_name: mushinmohan/NNSA


text:
boltzmann(β, E::Number) = exp(-β*E)

function boltzmann(β, E)
    factors = boltzmann.(β, E)
    partition = sum(factors)
    factors ./ partition
end

function boltzmann(β, E, deg)
    factors = map(*, boltzmann.(β, E), deg)
    partition = sum(factors)
    convert(Array{Float64}, factors ./ partition)
end

boltzmann(β, space::StateSpace) = boltzmann(β, space.energy, space.deg)


text:
using  Test
import RecursiveFactorization
import LinearAlgebra
using  LinearAlgebra: norm
using  Random

Random.seed!(12)

const baselu = LinearAlgebra.lu
const mylu = RecursiveFactorization.lu

function testlu(A, MF, BF)
    @test MF.info == BF.info
    @test norm(MF.L*MF.U - A[MF.p, :], Inf) < sqrt(eps(real(first(A))))
    nothing
end

@testset "Test LU factorization" begin
    for p in (Val(true), Val(false)), T in (Float64, Float32, ComplexF64, ComplexF32, Real)
        siz = (50, 100)
        if isconcretetype(T)
            A = rand(T, siz...)
        else
            _A = rand(50, 100)
            A  = Matrix{T}(undef, siz...)
            copyto!(A, _A)
        end
        MF = mylu(A, p)
        BF = baselu(A, p)
        testlu(A, MF, BF)
        for i in 50:7:100 # test `MF.info`
            A[:, i] .= 0
            MF = mylu(A, p, check=false)
            BF = baselu(A, p, check=false)
            testlu(A, MF, BF)
        end
    end
end


text:
module PosteriorSummary

export effective_sample_size,
       hpd_interval,
       hpd_excludes_value

import MCMCDiagnostics.effective_sample_size

function hpd_interval(x::AbstractVector{Float64}; alpha::Float64 = 0.05, sorted::Bool = false)
    n = length(x)
    span = Int(round((1.0 - alpha) * n))
    smallest_rng = Inf
    rng_upper = Inf
    rng_lower = -Inf
    sx::AbstractVector{Float64} = sorted ? x : sort(x)
    lower = 0
    for upper = span:n
        lower = upper - span + 1
        rng = sx[upper] - sx[lower]
        if rng < smallest_rng
            smallest_rng = rng
            rng_lower = sx[lower]
            rng_upper = sx[upper]
        end
    end

    return (lower = rng_lower, upper = rng_upper)
end

function hpd_interval(X::AbstractMatrix{Float64}; args...)
    p = size(X, 2)
    hpds = [hpd_interval(@view X[:, i]; args...) for i = 1:p]
    return hpds
end

function hpd_excludes_value(x::AbstractVector{Float64}, value::Float64; alpha::Float64 = 0.05, sorted::Bool = false)
    hpd = hpd_interval(x, alpha = alpha, sorted = sorted)
    return hpd[1] > value || hpd[2] < value
end

end


text:
# map_many.jl

using MIRT: map_many

using Test: @test, @inferred

N = (4,3)
fun = x -> cat(dims=3, x, 4x)
x1 = ones(N)
x = cat(dims=3, x1, 2x1)

#@inferred map_many(fun, x, N) # fails because fun could be anything!
@test map_many(fun, x, N) == cat(dims=4, fun(x1), fun(2x1))


text:
using BSONify, Test, Dates, NamedTupleTools
import Mongoc

function testReconstitutePrimitive(T, value)
    original = Dict("v" => value)
    # println("original: ", original)
    # println()
    bson = Mongoc.BSON(original)
    # println("bson: ", bson)
    # println()

    reconstituted = as_type(Dict{String, T}, bson)
    # println("reconstituted: ", reconstituted)
    # println()
    @test original["v"] == reconstituted["v"]
end

function testReconstituteStruct(T, value)
    # println(value)
    bson = Mongoc.BSON(value)
    # println(bson)
    reconstituted = as_type(T, bson)
    # println(reconstituted)
    @test value == reconstituted
end

numberTypes = [ Int8, Int16, Int32, Int64, UInt8, UInt16, UInt32, Float16, Float32, Float64 ]

####### Primitives (and dict, given implementation of testReconstitutePrimitive)

testReconstitutePrimitive(Bool, true)
testReconstitutePrimitive(Bool, false)

for T in numberTypes
    testReconstitutePrimitive(T, convert(T, 42))
end
testReconstitutePrimitive(String, "42")
testReconstitutePrimitive(DateTime, now())
testReconstitutePrimitive(DataType, Int16)

# Auto convert string to number (really only used when converting numeric map keys to strings to be bson compatible)
as_type(Int32, "123") == 123
as_type(Float64, "12.3") == 12.3

######## Numbers as dict keys

testReconstitutePrimitive(Dict{Int32, Int32}, Dict(12 => 34, 56 => 78))
testReconstitutePrimitive(Dict{Float32, Int32}, Dict(1.2f0 => 34, 5.6f0 => 78))

######## Arrays of primitives

for T in numberTypes
    testReconstitutePrimitive(Array{T}, collect(0:42))
end
testReconstitutePrimitive(Array{String}, split("It's not a question, but an answer learned in time"))
testReconstitutePrimitive(Array{DateTime}, [ now() ])

####### Tuples of primitives

for T in numberTypes
    testReconstitutePrimitive(Tuple{T}, tuple(collect(0:42)...))
end
testReconstitutePrimitive(Tuple{String}, tuple(split("It's not a question, but an answer learned in time")...))
testReconstitutePrimitive(Tuple{DateTime}, ( now(), ))

######## Named tuple

symbols = Meta.parse.(string.(numberTypes))
vals=[collect(1:length(numberTypes))...]
nt = namedtuple(symbols, vals)
testReconstitutePrimitive(typeof(nt), nt)

####### UNION TYPES

TestUnionType = Union{Int32, String, Nothing}

testReconstitutePrimitive(TestUnionType, 123)
testReconstitutePrimitive(TestUnionType, "blah")
testReconstitutePrimitive(TestUnionType, nothing)

########## STRUCTS #######

# Allow comparing structs with (mutable) arrays: https://stackoverflow.com/questions/62336686/struct-equality-with-arrays
abstract type Comparable end
import Base.==
function ==(a::T, b::T) where T <: Comparable
    f = fieldnames(T)
    getfield.(Ref(a),f) == getfield.(Ref(b),f)
end

####### Simple struct

struct SimpleStruct
    primitive1::Int32
    primitive2::String
end
testReconstituteStruct(SimpleStruct, SimpleStruct(123, "lalala"))

####### Parametric type struct

struct ParametricStruct{T} <: Comparable
    primitive1::Int32
    parametric::T
    array::Array{T}
end
testReconstituteStruct(ParametricStruct{Float32}, ParametricStruct{Float32}(
    123,
    4.56,
    [7, 8, 9]
))

####### Nested structs

struct GrandChildStruct{T, S} <: Comparable
    val1::T
    val2::S
    val3::Dict{T, S}
end

struct ChildStruct{T, S} <: Comparable
    grandChild::GrandChildStruct{T, S}
end

struct ParentStruct{T, S} <: Comparable
    child::ChildStruct{T, S}
end

testReconstituteStruct(
    ParentStruct{String, Int32},
    ParentStruct{String, Int32}(
        ChildStruct{String, Int32}(
            GrandChildStruct{String, Int32}(
                "foo",
                123,
                Dict("bar" => 456, "baz" => 789)
            )
        )
    )
)

mutable struct RecursiveStruct <: Comparable
    val::Int64
    child::Union{RecursiveStruct, Nothing}
end

testReconstituteStruct(
    RecursiveStruct,
    RecursiveStruct(
        1, RecursiveStruct(
            2, RecursiveStruct(
                3, nothing
            )
        )
    )
)

# EXCEPTIONS

@test_throws TypeValueMismatch as_type(Int32, 1.2)
@test_throws TypeValueMismatch as_type(String, 12)
@test_throws TypeValueMismatch as_type(SimpleStruct, 12)
@test_throws TypeValueMismatch testReconstitutePrimitive(TestUnionType, 12.3)

@test_throws FieldsMismatch as_type(SimpleStruct, Dict("primitive1" => 12))
@test_throws FieldsMismatch as_type(SimpleStruct, Dict("primitive1" => 12, "primitive3" => "12"))
@test_throws FieldsMismatch as_type(SimpleStruct, Dict("primitive1" => 12, "primitive2" => "12", "primitive3" => "12"))

@test_throws TypeIsNotAMapping Mongoc.BSON(123)

max_stars_repo_name: OliverColeman/BSONify

max_forks_repo_name: OliverColeman/BSONify


text:
using Test, LinearAlgebra
using QBase

@testset "./src/types/unitaries.jl" begin

@testset "is_unitary()" begin
    @testset "valid unitaries" begin
        valid_unitaries = [
            [0 1;1 0],
            [0 -im 0;im 0 0;0 0 1],
            [1 0 0 0;0 1 0 0;0 0 0 1;0 0 1 0],
        ]

        for U in valid_unitaries
            @test is_unitary(U)
        end
    end

    @testset "invalid unitaries" begin
        invalid_unitaries = [
            [0.5 0;0 0.5],
            [0 im 0;-im 0 0;0 0 0],
            [1 0 0;0 1 0],
        ]

        for U in invalid_unitaries
            @test !is_unitary(U)
        end
    end
end

@testset "Unitary()" begin
    @testset "valid cases" begin
        U = Unitary([0 1;1 0])
        @test U == [0 1;1 0]
        @test U isa Unitary{Int64}
        @test U isa Unitary

        U = Unitary([1/sqrt(2) 1/sqrt(2) 0;-1im/sqrt(2) 1im/sqrt(2) 0;0 0 im])
        @test U == [1/sqrt(2) 1/sqrt(2) 0;-1im/sqrt(2) 1im/sqrt(2) 0;0 0 im]
        @test U isa Unitary{Complex{Float64}}
        @test U isa Unitary
    end

    @testset "invalid unitaries" begin
        @test_throws DomainError Unitary([0 1;1 0;0 0])
        @test_throws DomainError Unitary([0.5 0.5;0.5 0.5])
    end

    @testset "atol" begin
        @test_throws DomainError Unitary([1 0.1;0 1])
        @test Unitary([1 0.1;0 1], atol=0.2) isa Unitary{Float64}
    end
end

@testset "unitary multiplication" begin
    Id = σx*σx
    @test Id == [1 0;0 1]
    @test Id isa Unitary{Int64}

    U_prod = *(random_unitary.([3,3,3,3,3])...)
    @test U_prod isa Unitary{Complex{Float64}}
end

@testset "unitary kron" begin
    σ_xz = kron(σx, σz)
    @test σ_xz == [0 0 1 0;0 0 0 -1;1 0 0 0;0 -1 0 0]
    @test σ_xz isa Unitary{Int64}

    σ_yII = kron(σy,σI,σI)
    @test σ_yII == [zeros(4,4) -im*Matrix(I,4,4);im*Matrix(I,4,4) zeros(4,4)]
    @test σ_yII isa Unitary{Complex{Int64}}
end

@testset "unitary adjoint" begin
    rx = qubit_rotation(π/2)
    @test rx ≈ [1 -im;-im 1]/sqrt(2)

    adj_rx = rx'
    @test adj_rx ≈ [1 im;im 1]/sqrt(2)
    @test adj_rx isa Unitary{Complex{Float64}}

    U = Unitary([1 0;0 2],  atol=3.01)
    @test U' isa Unitary{Int64}
    @test U' == [1 0;0 2]
end

end


text:
using TOML
using Base: open

include("event.jl")

struct File
    adress::String
end

struct Url
    adress::String
end

struct Card
    name::String
    events::Vector{Event}
    description::String
    image::Union{File,Url}
end

"""
    Card_FromFile(file::AbstractString)

# Arguments
- `file::AbstractString`: The file storing the card(s)
                          see `example.toml` for structure
"""
function Card_FromFile(file::AbstractString)
    cards = Card[]
    file = TOML.parsefile(file)
    for name in keys(file)
        # manipulate image to later determine if file or url
        image = lowercase(file[name]["image"])
        if startswith(image, "url:")
            image = Url(image)
        elseif startswith(image, "file:")
            image = File(replace(image))
        else
            @warn "Cannot determine image type of image for card $name"
        end

        # create events for card
        events = Event_CreateEvent.(file[name]["events"])

        description = file[name]["description"]
        push!(cards, Card(name, events, description, image))
    end

    return cards
end

max_stars_repo_name: Flobii/discord-bot

max_issues_repo_name: Flobii/discord-bot

max_forks_repo_name: Flobii/discord-bot


text:
# This file was generated, do not modify it. # hide
#hideall
team = [
  (name="Jane Doe", title="CEO & Founder", vitae="Phasellus eget enim eu lectus faucibus vestibulum", email="example@example.com"),
  (name="Mike Ross", title="Art Director", vitae="Phasellus eget enim eu lectus faucibus vestibulum", email="example@example.com"),
  (name="John Doe", title="Designer", vitae="Phasellus eget enim eu lectus faucibus vestibulum", email="example@example.com")
  ]

"@@cards @@row" |> println
for person in team
  """
  @@column
    \\card{$(person.name)}{$(person.title)}{$(person.vitae)}{$(person.email)}
  @@
  """ |> println
end
println("@@ @@") # end of cards + row

raw"""
~~~
<style>
.column {
  float:left;
  width:30%;
  margin-bottom:16px;
  padding:0 8px;
}
@media (max-width:62rem) {
  .column {
    width:45%;
    display:block;
  }
}
@media (max-width:30rem){
  .column {
    width:95%;
    display:block;
  }
}
.card{
  box-shadow: 0 4px 8px 0 rgba(0,0,0,0.2);
}
.card img {
  padding-left:0;
  width: 100%;
}
.container {
  padding: 0 16px;
}
.container::after, .row::after{
  content:"";
  clear:both;
  display:table;
}
.title {
  color:grey;
}
.vitae {
  margin-top:0.5em;
}
.email {
  font-family:courier;
  margin-top:0.5em;
  margin-bottom:0.5em;
}
.button{
  border:none;
  outline:0;
  display:inline-block;
  padding:8px;
  color:white;
  background-color:#000;
  text-align:center;
  cursor:pointer;
  width:100%;
}
.button:hover{
  background-color:#555;
}
</style>
~~~
""" |> println

max_stars_repo_name: mcognetta/franklindocs

max_issues_repo_name: mcognetta/franklindocs

max_forks_repo_name: mcognetta/franklindocs


text:
using Documenter, SearchLight, Validation, PostgreSQLDatabaseAdapter, Database

push!(LOAD_PATH,  "../../src")
push!(LOAD_PATH,  "../../src/database_adapters")

makedocs()


text:
module InferenceWrapper

using JavaCall

const AlternativeHypothesis = @jimport org.hipparchus.stat.inference.AlternativeHypothesis
const BinomialTest = @jimport org.hipparchus.stat.inference.BinomialTest
const ChiSquareTest = @jimport org.hipparchus.stat.inference.ChiSquareTest
const Class = @jimport java.lang.Class
const Collection = @jimport java.util.Collection
const Enum = @jimport java.lang.Enum
const GTest = @jimport org.hipparchus.stat.inference.GTest
const InferenceTestUtils = @jimport org.hipparchus.stat.inference.InferenceTestUtils
const JString = @jimport java.lang.String
const KolmogorovSmirnovTest = @jimport org.hipparchus.stat.inference.KolmogorovSmirnovTest
const MannWhitneyUTest = @jimport org.hipparchus.stat.inference.MannWhitneyUTest
const NaNStrategy = @jimport org.hipparchus.stat.ranking.NaNStrategy
const Object = @jimport java.lang.Object
const OneWayAnova = @jimport org.hipparchus.stat.inference.OneWayAnova
const Optional = @jimport java.util.Optional
const RealDistribution = @jimport org.hipparchus.distribution.RealDistribution
const StatisticalSummary = @jimport org.hipparchus.stat.descriptive.StatisticalSummary
const TTest = @jimport org.hipparchus.stat.inference.TTest
const TiesStrategy = @jimport org.hipparchus.stat.ranking.TiesStrategy
const WilcoxonSignedRankTest = @jimport org.hipparchus.stat.inference.WilcoxonSignedRankTest

export AlternativeHypothesis
export BinomialTest
export ChiSquareTest
export GTest
export InferenceTestUtils
export KolmogorovSmirnovTest
export MannWhitneyUTest
export OneWayAnova
export TTest
export WilcoxonSignedRankTest
export anova_f_value
export anova_p_value
export anova_test
export approximate_p
export binomial_test
export bootstrap
export cdf
export cdf_exact
export chi_square
export chi_square_data_sets_comparison
export chi_square_test
export chi_square_test_data_sets_comparison
export compare_to
export describe_constable
export equals
export exact_p
export g
export g_data_sets_comparison
export g_test
export g_test_data_sets_comparison
export g_test_intrinsic
export get_class
export get_declaring_class
export hash_code
export homoscedastic_t
export homoscedastic_t_test
export kolmogorov_smirnov_statistic
export kolmogorov_smirnov_test
export ks_sum
export mann_whitney_u
export mann_whitney_u_test
export name
export notify
export notify_all
export one_way_anova_f_value
export one_way_anova_p_value
export one_way_anova_test
export ordinal
export paired_t
export paired_t_test
export pelz_good
export root_log_likelihood_ratio
export t
export t_test
export to_string
export value_of
export values
export wait
export wilcoxon_signed_rank
export wilcoxon_signed_rank_test

include("alternative_hypothesis.jl")
include("binomial_test.jl")
include("chi_square_test.jl")
include("g_test.jl")
include("inference_test_utils.jl")
include("kolmogorov_smirnov_test.jl")
include("mann_whitney_u_test.jl")
include("one_way_anova.jl")
include("t_test.jl")
include("wilcoxon_signed_rank_test.jl")

end


text:
"""`ciid(f)`

`RandVar` that is __c__onditionally __i__ndependent given its parents, but __i__dentically  __d__istributed to `f`

`ciid` is the primary mechanism to construct a `RandVar` from a function.

```julia
function x_(ω)
  x = normal(ω, 0, 1)
end
x = ciid(x_)
rand(x)
```

__Important:__  any parents of `f` are shared.
That is, if `X` is a `RandVar`, then `X(ω)` will return the same value from
whatever context it is called.

Example:

```julia
numflips = poisson(2)

flips_(ω) = [bernoulli(ω, 0.5, Bool) for i = 1:numflips(ω)]
flips = ciid(flips_)

"At least one of numflips is true"
anyheads_(ω) = any(flips(ω))
anyheads = ciid(anyheads_)

"All flips are true"
allheads_(ω) = all(flips(ω))
allheads = ciid(allheads_)  # `allheads` and `anyheads` share `flips`

rand((numflips, flips, anyheads, allheads))
```
"""
ciid(f; id = uid()) = URandVar(f, (), id)

"ciid(x::RandVar) \n\n  `RandVar` identically distributed to `x` but conditionally independent given parents`"
ciid(x::T; id = uid()) where T <: RandVar =  T(params(x)..., id)
@spec equaldist(x, _res)

"""$(SIGNATURES) ciid with arguments

If arguments are random variables they are resolved to values.

Equivalent to:

`ciid(ω -> f(ω, (arg isa RandVar ? arg(ω) : arg for arg in args)...))`

Example:

```julia
function f_(ω, n)
  x = 0.0
  for i = 1:n
    x += uniform(ω, 0, 1)
  end
  x
end

f = ciid(f_, poisson(3))
```
"""
ciid(f, args...; id = uid()) = URandVar(reifyapply, (f, args...), id)

@inline reifyapply(ωπ, f, args...) = f(ωπ, reify(ωπ, args)...)
@inline reifyapply(ωπ, f) = f(ωπ)

Base.:~(x) = ciid(x)
Base.:~(id::ID, x) = ciid(x; id = id)
Base.:~(id, x) = ciid(x; id = toid(id))

## Consisten
# Distributions.jl / OmegaCore interopt
# using Distributions: Normal, Uniform, Bernoulli
# Base.:~(id, x::Normal) = id ~ normal(x.μ, x.σ)


text:
function channelprune(x::AbstractArray{<:Any, 4}, level;
                      magnitude = x -> norm(x, 1))
    sⱼs = [magnitude(x[:, :, :, j]) for j in axes(x, 4)]
    p = sortperm(sⱼs; rev = true)
    n = ceil(Int, (1 - level) * length(p))

    return mask(x, :, :, :, p[1:n])
end
channelprune(m::Conv, level; kwargs...) =
    mappruneable1(w -> channelprune(w, level; kwargs...), m)

struct ChannelPrune{T, S}
    level::T
    magnitude::S
end
ChannelPrune(level; magnitude = x -> norm(x, 1)) = ChannelPrune(level, magnitude)
(p::ChannelPrune)(m) = channelprune(m, p.level; magnitude = p.magnitude)

prune(strategy::ChannelPrune, m) =
    mappruneable(strategy, m; exclude = x -> x isa Conv)


text:
__precompile__()
module Mea

include("Green.jl")
include("Periodize.jl")
include("Sigmadc.jl")
include("Transport.jl")

#using Mea.Green, Mea.Sigmadc, Mea.Periodize
#export

# green.jl

# periodize.jl
#ModelVector, buildmodelvec

#sigmadc.jl
#calc_sigmadc

#transport
#calc_sigmadc

end # module


text:
#
# Correctness Tests
#

module TestCategoricalArrays
    fatalerrors = length(ARGS) > 0 && ARGS[1] == "-f"
    quiet = length(ARGS) > 0 && ARGS[1] == "-q"
    anyerrors = false

    using Test
    using CategoricalArrays

    tests = [
        "01_typedef.jl",
        "02_buildorder.jl",
        "03_buildfields.jl",
        "04_constructors.jl",
        "05_convert.jl",
        "05_copy.jl",
        "06_show.jl",
        "06_length.jl",
        "07_levels.jl",
        "08_equality.jl",
        "08_string.jl",
        "09_hash.jl",
        "10_isless.jl",
        "11_array.jl",
        "12_missingarray.jl",
        "13_arraycommon.jl",
        "14_view.jl",
        "15_extras.jl",
        "16_recode.jl"
    ]

    @testset "$test" for test in tests
        include(test)
    end
end


text:
# versions of the CUDA toolkit
const toolkits = [v"1.0", v"1.1",
                  v"2.0", v"2.1", v"2.2",
                  v"3.0", v"3.1", v"3.2",
                  v"4.0", v"4.1", v"4.2",
                  v"5.0", v"5.5",
                  v"6.0", v"6.5",
                  v"7.0", v"7.5",
                  v"8.0",
                  v"9.0", v"9.1"]

const cudnns = [v"1.0", v"2.0", v"3.0", v"4.0", v"5.0", v"5.1", v"6.0", v"7.0"]

struct VersionRange
    lower::VersionNumber
    upper::VersionNumber
end

Base.in(v::VersionNumber, r::VersionRange) = (v >= r.lower && v < r.upper)

Base.colon(a::VersionNumber, b::VersionNumber) = VersionRange(a, b)

Base.intersect(v::VersionNumber, r::VersionRange) =
    v < r.lower ? (r.lower:v) :
    v > r.upper ? (v:r.upper) : (v:v)

const lowest = v"0"
const highest = v"999"

# GCC compilers supported by the CUDA toolkit

# Source: CUDA/include/host_config.h
const cuda_gcc_db = Dict(
    v"5.5" => lowest:v"4.9-",   # (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ > 8)) && #error
    v"6.0" => lowest:v"4.9-",   # (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ > 8)) && #error
    v"6.5" => lowest:v"4.9-",   # (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ > 8)) && #error
    v"7.0" => lowest:v"4.10-",  # (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ > 9)) && #error
    v"7.5" => lowest:v"4.10-",  # (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ > 9)) && #error
    v"8.0" => lowest:v"6-",     # (__GNUC__ > 5)                                          && #error
    v"9.0" => lowest:v"7-",     # (__GNUC__ > 6)                                          && #error
    v"9.1" => lowest:v"7-",     # (__GNUC__ > 6)                                          && #error
)

function gcc_for_cuda(ver::VersionNumber)
    match_ver = VersionNumber(ver.major, ver.minor)
    return get(cuda_gcc_db, match_ver) do
        error("no support for CUDA $ver")
    end
end

# devices supported by the CUDA toolkit

# Source:
# - https://en.wikipedia.org/wiki/CUDA#GPUs_supported
# - ptxas |& grep -A 10 '\--gpu-name'
const dev_cuda_db = Dict(
    v"1.0" => lowest:v"7.0",
    v"1.1" => lowest:v"7.0",
    v"1.2" => lowest:v"7.0",
    v"1.3" => lowest:v"7.0",
    v"2.0" => lowest:v"9.0",
    v"2.1" => lowest:v"9.0",
    v"3.0" => v"4.2":highest,
    v"3.2" => v"6.0":highest,
    v"3.5" => v"5.0":highest,
    v"3.7" => v"6.5":highest,
    v"5.0" => v"6.0":highest,
    v"5.2" => v"7.0":highest,
    v"5.3" => v"7.5":highest,
    v"6.0" => v"8.0":highest,
    v"6.1" => v"8.0":highest,
    v"6.2" => v"8.0":highest,
    v"7.0" => v"9.0":highest
)

function devices_for_cuda(ver::VersionNumber)
    match_ver = VersionNumber(ver.major, ver.minor)

    caps = Set{VersionNumber}()
    for (cap,r) in dev_cuda_db
        if match_ver in r
            push!(caps, cap)
        end
    end
    return caps
end

# PTX ISAs supported by the CUDA toolkit

# Source:
# - PTX ISA document, Release History table
# NOTE: this table lists e.g. sm_20 being supported on CUDA 9.0, which is wrong?
const isa_cuda_db = Dict(
    v"1.0" => v"1.0":highest,
    v"1.1" => v"1.1":highest,
    v"1.2" => v"2.0":highest,
    v"1.3" => v"2.1":highest,
    v"1.4" => v"2.2":highest,
    v"1.5" => v"2.2":highest,   # driver 190
    v"2.0" => v"3.0":highest,   # driver 195
    v"2.1" => v"3.1":highest,   # driver 256
    v"2.2" => v"3.2":highest,   # driver 260
    v"2.3" => v"4.2":highest,   # driver 295, or driver 285 with 4.1
    v"3.0" => v"4.1":highest,   # driver 285
    v"3.1" => v"5.0":highest,   # driver 302
    v"3.2" => v"5.5":highest,   # driver 319
    v"4.0" => v"6.0":highest,   # driver 331
    v"4.1" => v"6.5":highest,   # driver 340
    v"4.2" => v"7.0":highest,   # driver 346
    v"4.3" => v"7.5":highest,   # driver 351
    v"5.0" => v"8.0":highest,   # driver 361
    v"6.0" => v"9.0":highest    # driver 384
)

function isas_for_cuda(ver::VersionNumber)
    match_ver = VersionNumber(ver.major, ver.minor)

    caps = Set{VersionNumber}()
    for (cap,r) in isa_cuda_db
        if match_ver in r
            push!(caps, cap)
        end
    end
    return caps
end

# devices supported by the LLVM NVPTX back-end

# Source: LLVM/lib/Target/NVPTX/NVPTX.td
const dev_llvm_db = Dict(
    v"2.0" => v"3.2":highest,
    v"2.1" => v"3.2":highest,
    v"3.0" => v"3.2":highest,
    v"3.2" => v"3.7":highest,
    v"3.5" => v"3.2":highest,
    v"3.7" => v"3.7":highest,
    v"5.0" => v"3.5":highest,
    v"5.2" => v"3.7":highest,
    v"5.3" => v"3.7":highest,
    v"6.0" => v"3.9":highest,
    v"6.1" => v"3.9":highest,
    v"6.2" => v"3.9":highest
)

function devices_for_llvm(ver::VersionNumber)
    match_ver = VersionNumber(ver.major, ver.minor)

    caps = Set{VersionNumber}()
    for (cap,r) in dev_llvm_db
        if match_ver in r
            push!(caps, cap)
        end
    end
    return caps
end

# PTX ISAs supported by the LVM NVPTX back-end

# Source: LLVM/lib/Target/NVPTX/NVPTX.td
const isa_llvm_db = Dict(
    v"3.0" => v"3.2":v"3.6",
    v"3.1" => v"3.2":v"3.6",
    v"3.2" => v"3.5":highest,
    v"4.0" => v"3.5":highest,
    v"4.1" => v"3.7":highest,
    v"4.2" => v"3.7":highest,
    v"4.3" => v"3.9":highest,
    v"5.0" => v"3.9":highest,
    v"6.0" => v"6.0":highest
)

function isas_for_llvm(ver::VersionNumber)
    match_ver = VersionNumber(ver.major, ver.minor)

    caps = Set{VersionNumber}()
    for (cap,r) in isa_llvm_db
        if match_ver in r
            push!(caps, cap)
        end
    end
    return caps
end

# other

shader(cap::VersionNumber) = "sm_$(cap.major)$(cap.minor)"


text:
# LibTOAST.jl: interface to the TOAST++ library
# Copyright (C) 2019 Samuel Powell

@testset "Regularisation" begin

  fdiskip = 300     # Finite difference applied to every xth element

  nx, ny = 300, 300
  pixelmap = PixelMap(mesh2D, (nx,ny))
  vtx,  = data(mesh2D)

  # Create some reasonably smooth function
  nprm = sin.(5.0*vtx[:,1]./15) .+ cos.(4.0*vtx[:,2]./15)

  nc = NodalCoeff(mesh2D, nprm)
  rc = SolutionCoeff(pixelmap, nc)
  rcpdx = similar(rc)
  rcmdx = similar(rc)
  rcp1 = copy(rc)
  rcp1 .+= 1
  r0 = SolutionCoeff(pixelmap, zeros(slen(pixelmap)))

  @testset "TK0" begin

    reg = RegulTK0(rc)

    @test val(reg, rc) == 0
    @test val(reg, r0) ≈ norm(rc.data)^2

    @test norm(LibTOAST.grad(reg,rc)) == 0
    @test norm(LibTOAST.grad(reg, r0)) ≈ 2*norm(rc)

  end

  @testset "TK1" begin

    reg = RegulTK1(rc)

    @test val(reg, rc) == 0.0
    @test ≈(val(reg, rcp1), 0.0, atol=1e-20)

    # Finite difference the grad
    g = grad(reg, rc)

    for i = 1:fdiskip:length(g)
      dx = 0.0001*rc[i]
      rcpdx .= rc; rcpdx[i] += dx
      rcmdx .= rc; rcmdx[i] -= dx
      @test abs((val(reg, rcpdx) - val(reg, rcmdx))/(2*dx) - g[i]) < 1e-5
    end

  end

  @testset "TV" begin

    reg = RegulTV(rc)

    @test val(reg, rc) == 0.0
    @test ≈(val(reg, rcp1), 0.0, atol=1e-20)

    # Finite difference the grad
    g = grad(reg, rc)

    for i = 1:fdiskip:length(g)
      dx = 0.0001*rc[i]
      rcpdx .= rc; rcpdx[i] += dx
      rcmdx .= rc; rcmdx[i] -= dx
      @test abs((val(reg, rcpdx) - val(reg, rcmdx))/(2*dx) - g[i]) < 1e-5
    end

  end

  @testset "PM" begin

    reg = RegulPM(rc)

    @test val(reg, rc) == 0.0
    @test ≈(val(reg, rcp1), 0.0, atol=1e-20)

    # Finite difference the grad
    g = grad(reg, rc)

    for i = 1:fdiskip:length(g)
      dx = 0.0001*rc[i]
      rcpdx .= rc; rcpdx[i] += dx
      rcmdx .= rc; rcmdx[i] -= dx
      @test abs((val(reg, rcpdx) - val(reg, rcmdx))/(2*dx) - g[i]) < 1e-5
    end

  end

end


text:
#############################################################################
#Copyright (c) 2016-2019 Charles Le Losq
#
#The MIT License (MIT)
#
#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the #Software without restriction, including without limitation the rights to use, copy, #modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, #and to permit persons to whom the Software is furnished to do so, subject to the #following conditions:
#
#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
#
#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, #INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
#
# baseline.jl contains several "baseline" functions. It is directly dependent on functions.jl.
#
#############################################################################

"""
	baseline(x::Array{Float64},y::Array{Float64},roi::Array{Float64},basetype::AbstractString;polynomial_order=1, s = 1.0, lam = 10^5, p = 0.01, ratio = 0.01, niter = 10, p0_exp = [1.,1.,1.],p0_log =[1.,1.,1.])

Allows subtracting a baseline under a x y spectrum; uses the baseline function from the rampy Python package.

Parameters
----------
	x_input : ndarray
		x values.
	y_input : ndarray
		y values.
	bir : ndarray
		Contain the regions of interest, organised per line.
		For instance, roi = np.array([[100., 200.],[500.,600.]]) will
		define roi between 100 and 200 as well as between 500 and 600.
		Note: This is NOT used by the "als" and "arPLS" algorithms, but still is a requirement when calling the function.
		bir and method probably will become args in a futur iteration of rampy to solve this.
	methods : str
		"poly": polynomial fitting, with splinesmooth the degree of the polynomial.
		"unispline": spline with the UnivariateSpline function of Scipy, splinesmooth is
					 the spline smoothing factor (assume equal weight in the present case);
		"gcvspline": spline with the gcvspl.f algorythm, really robust.
					 Spectra must have x, y, ese in it, and splinesmooth is the smoothing factor;
					 For gcvspline, if ese are not provided we assume ese = sqrt(y).
					 WARNING: Requires the installation of the gcvspline Python package prior to use in the Python ENV used by Julia.
					 See website for install instructions
		"exp": exponential background;
		"log": logarythmic background;
		"rubberband": rubberband baseline fitting;
		"als": automatic least square fitting following Eilers and Boelens 2005;
		"arPLS": automatic baseline fit using the algorithm from Baek et al. 2015
				 Baseline correction using asymmetrically reweighted penalized least squares smoothing, Analyst 140: 250-257.

Options
------
	polynomial_order : Int
		The degree of the polynomial (0 for a constant), default = 1.
	s : Float
		spline smoothing coefficient for the unispline and gcvspline algorithms.
	lam : Float
		float, the lambda smoothness parameter for the ALS and ArPLS algorithms. Typical values are between 10**2 to 10**9, default = 10**5.
	p : Float
		float, for the ALS algorithm, advised value between 0.001 to 0.1, default = 0.01.
	ratio : float
		ratio parameter of the arPLS algorithm. default = 0.01.
	niter : Int
		number of iteration of the ALS algorithm, default = 10.
	p0_exp : List
		containg the starting parameter for the exp baseline fit with curve_fit. Default = [1.,1.,1.].
	p0_log : List
		containg the starting parameter for the log baseline fit with curve_fit. Default = [1.,1.,1.,1.].

Returns
-------
	out1 : ndarray
		Contain the corrected signal.
	out2 : ndarray
		Contain the baseline.

Example
-------

Consider a Y signal sampled at X. We want to fit a polynomial baseline (2nd order)
in the regions of interest comprised between x1 and x2, as well as x3 and x4.

	```julia-repl
	julia> Y_corr, Y_baseline = baseline(X,Y,[x1 x2; x3 x4],"poly",polynomial_order=2)
	```
"""
function baseline(x::Array{Float64},y::Array{Float64},roi::Array{Float64},basetype::AbstractString;polynomial_order=1, s = 1.0, lam = 10^5, p = 0.01, ratio = 0.01, niter = 10, p0_exp = [1.,1.,1.],p0_log =[1.,1.,1.])
	yout, bas = rampy.baseline(x,y,roi,basetype,polynomial_order=polynomial_order, s=s, lam=lam, niter=niter, p0_exp=p0_exp, p0_log=p0_log)
end


text:
"""
    Equation(R, K)

Construct an Equation that represents a PDE of the form

    ∂ϕ/∂t = ∂/∂z (K ∂ϕ/∂z) + R ,

so that K is a diffusivity and R captures 'remaining' terms.

In OceanTurb, K and R are arrays of functions.
The functions K[j](m, i) and R[j](m, i) calculate the
diffusivity K and remaining terms R for the model `m` at
grid point `i`.
"""
struct Equation{TM, TL, TK, TR, U}
    M       :: TM
    L       :: TL
    K       :: TK
    R       :: TR
    update! :: U
end

minuszero(args...) = -0
funkynothing(args...) = nothing

function Equation(;
         K,
         M = Tuple(minuszero for i=1:length(K)),
         L = Tuple(minuszero for i=1:length(K)),
         R = Tuple(minuszero for i=1:length(K)),
    update = funkynothing
    )
    return Equation(M, L, K, R, update)
end


text:
#################### Binary Metropolised Gibbs Sampler ####################

#################### Types and Constructors ####################

const BMGForm = Union{Int, Vector{Vector{Int}}}

mutable struct BMGTune{F<:BMGForm} <: SamplerTune
  logf::Union{Function, Missing}
  k::F

  BMGTune{F}() where {F<:BMGForm} = new{F}()

  BMGTune{F}(x::Vector, k::F) where {F<:BMGForm} = new{F}(missing, k)

  BMGTune{F}(x::Vector, k::F, logf::Function) where {F<:BMGForm} =
    new{F}(logf, k)
end

const BMGIntVariate = SamplerVariate{BMGTune{Int}}
const BMGVecVariate = SamplerVariate{BMGTune{Vector{Vector{Int}}}}

BMGVariate(x::Vector, logf::Function; k::F=1) where {F<:BMGForm} =
  SamplerVariate{BMGTune{F}}(x, k, logf)

function validate(v::BMGIntVariate)
  n = length(v)
  v.tune.k <= n || throw(ArgumentError("k exceeds variate length $n"))
  validatebinary(v)
end

function validate(v::BMGVecVariate)
  n = length(v)
  mapreduce(maximum, max, v.tune.k) <= n ||
    throw(ArgumentError("indices in k exceed variate length $n"))
  validatebinary(v)
end

#################### Sampler Constructor ####################
"""
    BMG(params::ElementOrVector{Symbol}; k::F=1) where {F<:BMGForm}

Construct a `Sampler` object for BMG sampling. Parameters are assumed to have
 binary numerical values (0 or 1).

Returns a `Sampler{BMGTune{typeof(k)}}` type object.

* `params`: stochastic node(s) to be updated with the sampler.

* `k`:  number of parameters or vector of parameter indices to select at random for simultaneous updating in each call of the sampler.
"""
function BMG(params::ElementOrVector{Symbol}; k::F=1) where {F<:BMGForm}
  samplerfx = function(model::Model, block::Integer)
    block = SamplingBlock(model, block)
    v = SamplerVariate(block, k)
    sample!(v, x -> logpdf!(block, x))
    relist(block, v)
  end
  Sampler(params, samplerfx, BMGTune{F}())
end

#################### Sampling Functions ####################

sample!(v::SamplerVariate{BMGTune{F}}) where {F<:BMGForm} = sample!(v, v.tune.logf)
"""
    sample!(v::SamplerVariate{BMGTune{F}}, logf::Function) where {F<:BMGForm}

Draw one sample from a target distribution using the BMG sampler. Parameters
are assumed to have binary numerical values (0 or 1).

Returns `v` updated with simulated values and associated tuning parameters.
"""
function sample!(v::SamplerVariate{BMGTune{F}}, logf::Function) where {F<:BMGForm}
  n = length(v)
  probs = Vector{Float64}(undef, n)
  idx = randind(v)

  pbernoulli! = function(x, probs)
    for i in idx
      x_i = x[i]

      x[i] = 0.0
      logf0 = logf(x)
      x[i] = 1.0
      logf1 = logf(x)

      p = invlogit(logf1 - logf0)
      probs[i] = 0.0 < p < 1.0 ? p : 0.5

      x[i] = x_i
    end
    probs
  end

  x = v[:]
  pbernoulli!(x, probs)
  theta = map(p -> rand() < p, probs[idx])

  if n == 1
    v[idx] = theta
  else
    y = v[:]
    y[idx] = theta
    qy = mapreduce(i -> y[i] == 1.0 ? log(probs[i]) : log(1.0 - probs[i]), +,
                   idx)

    pbernoulli!(y, probs)
    qx = mapreduce(i -> x[i] == 1.0 ? log(probs[i]) : log(1.0 - probs[i]), +,
                   idx)

    if rand() < exp((logf(y) - qy) - (logf(x) - qx))
      v[idx] = theta
    end
  end

  v
end

randind(v::BMGIntVariate) = sample(1:length(v), v.tune.k, replace=false)
randind(v::BMGVecVariate) = sample(v.tune.k)


text:
#####
##### Stepsize and adaptation tests
#####

isinteractive() && include("common.jl")

@testset "stepsize general rootfinding" begin
    Δ = 3.0                   # shift exponential so that ϵ = 1 is not in interval
    A = ϵ -> exp(-ϵ*Δ)
    invA = A -> -log(A) / Δ
    @test 0.1 ≈ invA(A(0.1)) atol = 1e-10 # inverse is OK
    params = InitialStepsizeSearch()
    @test A(params.ϵ₀) ≤ params.a_min || A(params.ϵ₀) ≥ params.a_max # outside interval
    constantA(ϵ) = params.a_max + 0.1
    # parameters and defaults
    @test params.a_min == 0.25
    @test params.a_max ≤ 1.0
    @test params.ϵ₀ == 1.0
    @test params.C == 2.0
    @test params.maxiter_crossing == 400
    @test params.maxiter_bisect == 400
    @test_throws ArgumentError InitialStepsizeSearch(; a_min = 0.9, a_max = 0.1) # not <
    @test_throws ArgumentError InitialStepsizeSearch(; C = 0.5) # not > 1
    @test_throws ArgumentError InitialStepsizeSearch(; maxiter_crossing = 2)
    @test_throws ArgumentError InitialStepsizeSearch(; maxiter_bisect = 2)
    # crossing from below
    ϵ₀, Aϵ₀, ϵ₁, Aϵ₁ = find_crossing_stepsize(params, A, invA(params.a_min-0.1))
    @test ϵ₀ ≥ invA(params.a_min) ≥ ϵ₁
    @test Aϵ₀ ≤ params.a_min ≤ Aϵ₁
    @test_throws ErrorException find_crossing_stepsize(params, constantA, 100.0)
    # crossing from above
    ϵ₀, Aϵ₀, ϵ₁, Aϵ₁ = find_crossing_stepsize(params, A, invA(params.a_max+0.1))
    @test ϵ₀ ≤ invA(params.a_max) ≤ ϵ₁
    @test Aϵ₀ ≥ params.a_max ≥ Aϵ₁
    @test_throws ErrorException find_crossing_stepsize(params, constantA, 0.1)
    # bisection
    ϵ = bisect_stepsize(params, A, invA(params.a_max + 0.1), invA(params.a_min - 0.1))
    @test params.a_min ≤ A(ϵ) ≤ params.a_max
    @test_throws ArgumentError bisect_stepsize(params, A, 0.4, 0.3) # order
    @test_throws ArgumentError bisect_stepsize(params, A, # already in interval
                                               invA(params.a_max - 0.1),
                                               invA(params.a_min + 0.1))
    # combined algorithm - single test
    ϵ = find_initial_stepsize(params, A)
    @test params.a_min ≤ A(ϵ) ≤ params.a_max
    # combined algorithm - random tests
    for _ in 1:1000
        α = abs(randn())
        β = -abs(randn())
        A = ϵ -> exp(-(α*ϵ+β)*ϵ)  # A(0)=1, has a hump, goes to A(∞)→0
        Aϵ₀ = A(params.ϵ₀)
        # uncomment line belows for coverage debugging or inverse
        # println("α = $α, β = $β")
        # println(Aϵ₀ ≤ params.a_min ? "below" : (Aϵ₀ ≥ params.a_max ? "above" : "in"))
        # invA(A) = -(log(A)+β)/A
        ϵ = find_initial_stepsize(params, A)
        @test params.a_min ≤ A(ϵ) ≤ params.a_max
    end
end

"""
$(SIGNATURES)

A parametric random acceptance rate that depends on the stepsize. For unit
testing acceptance rate tuning.
"""
dummy_acceptance_rate(ϵ, σ = 0.05) = min(1/ϵ * exp(randn()*σ - σ^2/2), 1)

mean_dummy_acceptance_rate(ϵ, σ = 0.05) = mean(dummy_acceptance_rate(ϵ, σ) for _ in 1:10000)

@testset "dual averaging far" begin
    ϵ₀ = 100.0                # way off
    δ = 0.65
    dual_averaging = DualAveraging(; δ = δ)
    A = initial_adaptation_state(dual_averaging, ϵ₀)
    @test A.logϵ̄ == 0           # ϵ̄₀ = 1 in Gelman and Hoffman (2014)
    @test A.m == 0
    @test A.H̄ == 0
    for _ in 1:500
        A = adapt_stepsize(dual_averaging, A, dummy_acceptance_rate(current_ϵ(A)))
    end
    @test mean_dummy_acceptance_rate(final_ϵ(A)) ≈ δ atol = 0.02
end

@testset "dual averaging close" begin
    ϵ₀ = 2.0
    δ = 0.65
    dual_averaging = DualAveraging(; δ = δ)
    A = initial_adaptation_state(dual_averaging, ϵ₀)
    for _ in 1:2000
        A = adapt_stepsize(dual_averaging, A, dummy_acceptance_rate(current_ϵ(A)))
    end
    @test mean_dummy_acceptance_rate(final_ϵ(A)) ≈ δ atol = 0.01
end

@testset "dual averaging far and noisy" begin
    ϵ₀ = 20.0
    δ = 0.65
    dual_averaging = DualAveraging(; δ = δ)
    A = initial_adaptation_state(dual_averaging, ϵ₀)
    for _ in 1:2000
        A = adapt_stepsize(dual_averaging, A, dummy_acceptance_rate(current_ϵ(A), 2.0))
    end
    @test mean_dummy_acceptance_rate(final_ϵ(A), 2.0) ≈ δ atol = 0.03
end

@testset "fixed stepsize sanity checks" begin
    fs = FixedStepsize()
    ϵ = 1.0
    A = initial_adaptation_state(fs, ϵ)
    @test A == adapt_stepsize(fs, A, ϵ)
    @test current_ϵ(A) == ϵ
    @test final_ϵ(A) == ϵ
end

@testset "find reasonable stepsize - random H, z" begin
    p = InitialStepsizeSearch()
    for _ in 1:100
        @unpack H, z = rand_Hz(rand(3:5))
        ϵ = find_initial_stepsize(p, H, z)
        logA = logdensity(H, leapfrog(H, z, ϵ)) - logdensity(H, z)
        @test p.a_min ≤ exp(logA) ≤ p.a_max
    end
end

@testset "error for non-finite initial density" begin
    p = InitialStepsizeSearch()
    @unpack H, z = rand_Hz(2)
    z = DynamicHMC.PhasePoint(z.Q, [NaN, NaN])
    @test_throws DomainError find_initial_stepsize(p, H, z)
end


text:
abstract type NewFitness end

abstract type SingleObjectiveFitness <: NewFitness end

abstract type MultiObjectiveFitness <: NewFitness end

# Some evaluators go through several stages of fitness evaluation:
#   1. a candidate is simulated -> *simulation(s)*
#   2. simulations have *characteristics* that are currently in *focus*
#   3. characteristics have *fitnessvalues* in relation to current *goals*
# A certain fitness value summarizes all of this information. By default
# this complex model is not really needed and we set the characteristics == fitnessvalues
# and the evaluatorversion (which determines both the focus and the goals that we active
# for this version of the evaluator) is a constant, i.e. the focus and goals do not change.
# To judge if the versions are current some fitness values saves a reference to the evaluator,
# but the default is not to have such a reference.
# For more advanced schemes you should override these methods but here are the defaults:

focusversion(f::NewFitness) = 0
goalversion(f::NewFitness) = 0
simulations(f::NewFitness) = characteristics(f)
characteristics(f::NewFitness) = fitnessvalues(f)
evaluator(f::NewFitness) = nothing

struct SingleNumberFitness{T <: Real} <: SingleObjectiveFitness
    fvalue::T
end

struct VectorFitness{T <: Real} <: MultiObjectiveFitness
    fvalues::Vector{T}
end

makefitness(fvalue::T) where {T <: Real} = SingleNumberFitness(fvalue)
makefitness(fvalues::Vector{T}) where {T <: Real} = VectorFitness(fvalues)

fitnessvalues(f::VectorFitness{T}) where {T <: Real} = f.fvalues
fitnessvalues(f::SingleNumberFitness{T}) where {T <: Real} = [f.fvalue]


text:
"""
treeknit

# Arguments

- `nwk1`: Newick file for first tree
- `nwk2`: Newick file for second tree

# Options

- `-o, --outdir <arg>`: directory for results; Example `-o=treeknit_results`
- `-g, --gamma <arg>`: value of γ; Example `-g=2`
- `--seq-lengths <arg>`: length of the sequences. Example: `--seq-length "1500 2000"`
- `--n-sa-it <arg>`: number of SA iterations per temperature and per leaf.

# Flags

- `--naive`: Naive inference (overrides `-g`).
- `--no-likelihood`: Do not use branch length likelihood test to sort between different MCCs
- `--no-resolve`: Do not attempt to resolve trees before inferring MCCs.
"""
@main function treeknit(
	nwk1::AbstractString, nwk2::AbstractString;
	# options
	outdir::AbstractString = "treeknit_results",
	gamma::Float64 = 2.,
	seq_lengths::AbstractString = "1 1",
	n_sa_it::Float64 = 1.,
	# flags
	naive::Bool = false,
	no_likelihood::Bool = false,
	no_resolve::Bool = false,
)
	println("Treeknit: ")
	println("Input trees: $nwk1 \t $nwk2")
	println("Results directory: $outdir")
	println("γ: $gamma")

	# Reading trees
	t1 = read_tree(nwk1)
	t2 = read_tree(nwk2)
	if !TreeTools.share_labels(t1, t2)
		error("Trees must share leaves")
	end

	# Setting up OptArgs
	sl = map(i->parse(Int, i), split(seq_lengths, " "))
	oa = OptArgs(;
		γ = gamma,
		likelihood_sort = !no_likelihood,
		resolve = !no_resolve,
		Md = 1 / n_sa_it,
	)

	# Setting up directories
	mkpath(outdir)

	#
	MCCs = computeMCCs(t1, t2, oa; naive, seqlengths = sl)
	resolve!(t1, t2, MCCs)
	arg, rlm, lm1, lm2 = SRG.arg_from_trees(t1, t2, MCCs)

	# Write output
	write_mccs(outdir * "/" * "MCCs.dat", MCCs)
	for (nwk, t) in zip((nwk1,nwk2),(t1,t2))
		fn = basename(nwk)
		name, ext = splitext(fn)
		write_newick(outdir * "/" * name * ".resolved" * ext, t)
	end
	write(outdir * "/" * "arg.nwk", arg)
	write_rlm(outdir * "/" * "nodes.dat", rlm)
end

function write_rlm(filename, rlm)
	open(filename, "w") do io
		for (n, (k,v)) in enumerate(rlm)
			write(io, k)
			write(io, ",")
			for (i,vv) in enumerate(v)
				if isnothing(vv)
					write(io, " ")
				else
					write(io, vv)
				end
				if i != length(v)
					write(io, ",")
				end
			end
			if n != length(rlm)
				write(io, "\n")
			end
		end
	end
end

max_stars_repo_name: PierreBarrat/OldTreeKnit

max_issues_repo_name: PierreBarrat/OldTreeKnit

max_forks_repo_name: PierreBarrat/OldTreeKnit


text:
# Autogenerated wrapper script for ROCmDeviceLibs_jll for x86_64-linux-musl
export bitcode_path

using Zlib_jll
JLLWrappers.@generate_wrapper_header("ROCmDeviceLibs")
JLLWrappers.@declare_file_product(bitcode_path)
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll)
    JLLWrappers.@init_file_product(
        bitcode_path,
        "amdgcn/bitcode",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using MittagLefflerFunctions
using PyPlot
using Polynomials

x = range(0, 10, length=201)
α_vals = [0.2, 0.4, 0.6, 0.8, 1.0]

figure(1)
r = 5
m = r+1
n = r
style = ("-", "--", ":", "-.", "-")
for row = 1:length(α_vals)
    α = α_vals[row]
    Eα_approx = MLPade(α, 1.0, m, n, :svd)
    Eα_ref = MLQuad(α, 1.0, 14, :hyperbola)
    labelstring = latexstring("\$\\alpha=$α\$")
    semilogy(x, abs.(Eα_ref.(-x) - Eα_approx.(x)),
             style[row], label=labelstring)
end
xlabel(L"$x$")
legend()
grid(true)

savefig("figure10.pdf")

α = 1.85
β = 1.0
println("Setting α = $α, β = $β.")
for (m, n) in [(8,7), (10,9)]
    println("\nRoots of q when m = $m, n = $n:")
    E = MLPade(α, β, m, n, :svd)
    q = Polynomial(E.q)
    display(roots(q))
end


text:
include("globals.jl")

module sources
using globals;

function hard_gaussian_source!(ez, pos, time, delay=30., width=100.)
    ez[pos] = exp(-(time - delay) * (time - delay) / width);
end

function additive_gaussian_source!(ez, pos, time, delay=30, width=100.)
    ez[pos] += exp(-(time - delay) * (time - delay) / width);
end

function gaussian_source(pos, time, delay=30, width=100.)
    ez_inc = exp(-(time + 0.5 - (-0.5) - delay) * (time + 0.5 - (-0.5) - delay) / width);
    hy_inc = exp(-(time - delay) * (time - delay) / width) / globals.imp0;

    # ez_inc = exp(-((time - pos - delay) / globals.cdtds / width)^2 );
    return ez_inc, hy_inc
end

# function harmonic_source(pos, time, phase=0., amp=1.)
#    ez

end


text:
abstract type Toroid end

struct DCore <: Toroid

    FlatHeight::Float64
    MaxHeight::Float64 #Height from CENTER LINE, so the actual core is 2 times as tall overall. This is the half-height effectively
    RadiusAtPeak::Float64 #This is how far from the center axis the peak of the geometry is, not the radius of curvature at the apex
    ID::Float64 #Thner diameter of the core
    OD::Float64 #Thner diameter of the core
    Turns::Int64
    Layers::Int64
    WireLength::Float64
    Resistance::Float64
end

struct Circle <: Toroid

    CoreRadius::Float64
    CenterRadius::Float64
    ID::Float64 #Thner diameter of the core
    OD::Float64 #Thner diameter of the core
    Turns::Int64
    Layers::Int64
    WireLength::Float64
    Resistance::Float64
end
struct GeneralParams <: Toroid
    SingleLayerInductance::Float64
    WireDiameter::Float64
    TargetInductance::Float64

end

mutable struct Geom
    DCore::DCore
    Circ::Circle
    General::GeneralParams
end
"""
This function makes an optimal D-Core toroid given the following inputs:
    Dia - wire diamter in meters
    LTarget-  target inductance in Henries

    kwarg:
    NumLayers - Number of windings layers
    CoreMu - permeability scaling factor
    Alpha - ratio of OD/ID use an integer between 2-5
    CuFillFactor - scales resistance value. If using litz it will be <1
    ExportFileName - if you want to export the SVG of the core cross section you can use this. If you enter a file name E.g. "Test.svg", it will save it in the present working directory
    NPts - number of points in the exported svg

"""
function ToroidOptimizer(
    Dia,
    LTarget;
    NumLayers = 2,
    CoreMu = 1,
    Alpha = 2,
    CuFillFactor = 1,
    ExportFileName = nothing,
    NPts = 100
)
    ## Summary: To use this function enter the wire diameter and the target
    #inductance (in meters and Henries). The calculations assume "NumLayers" fully
    #wound layers as it is a good balance between heat dissipation and
    #geometric efficiency (2 is a good option, unless heat is an issue). The output is a struct containing the geometric
    #information for a toroid with a D shaped core, and a circular core.

    ## Constants
    μ_o = 4 * pi * 10^(-7)
    μ = μ_o * CoreMu
    L0 = μ * Dia / (2 * pi)
    L = LTarget / NumLayers^2 #Use the assumption that there will be 'NumLayers' fully wound layers...
    #If there are two fully wound layers, it can be calculated at just one
    #layer with 1/4 the inductance. Then double it when winding
    LperL0 = L / L0 #An intermediate dimensionless parameter that comes from the papers which these formulas are derived.
    SingleLayerInductance = L
    WireDiameter = Dia
    TargetInductance = LTarget
    General =
        GeneralParams(SingleLayerInductance, WireDiameter, TargetInductance)

    ## D shaped core
    #Alpha = 2; #User can change this to an integer between 2 and 5 (inclusive) varies the aspect ratio of the D
    #A high alpha value such as 5 will be a taller/bigger core with fewer
    #turns. It is theoretically more efficient but will probably have more
    #leakage flux.

    # AlphaMat Key. The columns represent [alpha e h p s t]

    AlphaMat = [
        2 0.26 0.645 3.6 0.72 1.77
        3 0.85 1.5 8.0 2.74 4.42
        4 1.6 2.4 12.8 5.76 8.09
        5 2.45 3.4 17.9 9.61 12.6
    ] #Table 1 From  "D-Shaped toroidal cage inductors" P.N. Murgatroyd & D. Belahrache,1989

    T = AlphaMat[Alpha-1, 6] #Picking a value in the table for clarity
    S = AlphaMat[Alpha-1, 5] #Also Selecting value for clarity
    P = AlphaMat[Alpha-1, 4] # As above

    KFunk_D(k) =
        abs.((2 * pi)^0.5 * (S / (P^(3 / 2))) * k .^ (3 / 2) + 0.25 * k - LperL0) #From "The Optimal Form for Coreless Inductor", P. Murgatroyd IEEE TMI, 25 No. 3 1989
    Ktmp = optimize(KFunk_D, 0, 1e5) #solving for ideal K with bisection method
    K = Ktmp.minimizer
    WireLength_D = K * Dia #Equation 6 in "Economic designs for..."
    Turns_D = (2 * pi * K / P)^(1 / 2) # Also from P. Murgatroyd, section 4, Eq. 20 in "Economic designs for single-layer toroidal inductors"
    # TurnsD = sqrt(2*pi*L/(T*μ*B)); #Turns needed in a D shaped core toroid
    # ID_DToroid = 2*B;
    # B = WireLength_D/(Turns_D*P); #B=inner radius, so B = total length/(number of turns*Perimeter of each turn)
    B = Dia / 2 + Turns_D * Dia / (2 * pi) #Making it so that the wires all tough on the inner edge of the core
    ## For the following reference Fig 4 in "D-Shaped toroidal cage inductors" P.N. Murgatroyd & D. Belahrache,1989
    FlatHeight = B * AlphaMat[Alpha-1, 2]
    MaxHeight = B * AlphaMat[Alpha-1, 3] #Height from CENTER LINE, so the actual core is 2 times as tall overall. This is the half-height effectively
    RadiusAtPeak = Alpha^0.5 * B #This is how far from the center axis the peak of the geometry is, not the radius of curvature at the apex
    ID = B * 2#Thner diameter of the core
    OD = 2 * Alpha * B
    Turns = round(Turns_D * NumLayers)
    Layers = NumLayers
    WireLength = WireLength_D * NumLayers
    D_WireR = Length2Resist(WireLength, Dia; FillFac = CuFillFactor)
    DCoreGeom = DCore(
        FlatHeight,
        MaxHeight,
        RadiusAtPeak,
        ID,
        OD,
        round(Turns),
        Layers,
        WireLength,
        D_WireR,
    )

    ## Two-layer circular core

    KFunk(k) = abs.(0.2722 * k .^ (3 / 2) + 0.25 * k - LperL0) #From "The Optimal Form for Coreless Inductor", P. Murgatroyd IEEE TMI, 25 No. 3 1989
    Kopt = optimize(KFunk, 0, 1e5) #solving for ideal K
    K = Kopt.minimizer
    WireLength = K * Dia
    Turns = 0.8165 * K^(1 / 2) # Also from P. Murgatroyd, section 2 in "Economic designs for single-layer toroidal inductors"
    R = WireLength / (2 * pi * Turns)# Equation 2 in "Economic designs for single-layer toroidal inductors",P. Murgatroyd
    T = Dia / (sin(pi / Turns) * 2) + R# Equation 3 in "Economic designs for single-layer toroidal inductors"

    Turns = Turns * NumLayers
    WireLength = WireLength * NumLayers
    Layers = NumLayers
    ID = 2 * (T - R)
    OD = 2 * (T + R)
    CoreRadius = R
    CenterRadius = T
    Circ_WireR = Length2Resist(WireLength, Dia; FillFac = CuFillFactor)
    CircleGeom = Circle(
        CoreRadius,
        CenterRadius,
        ID,
        OD,
        round(Turns),
        Layers,
        WireLength,
        Circ_WireR,
    )

    ## Back-checking the work with two approximations.

    Approx_EnclosedArea_DCore =
        DCoreGeom.MaxHeight * 2 * (DCoreGeom.OD - DCoreGeom.ID) / 2 * 0.75#This ROUGH approximation is assuming the core fills 75# of the circumscribed rectangle around the core (roughly a semi-cirle)
    SanityCheckInductance =
        μ * DCoreGeom.Turns^2 * Approx_EnclosedArea_DCore /
        (pi * (DCoreGeom.OD + DCoreGeom.ID) / 2)

    MeanPathLength = (DCoreGeom.OD + DCoreGeom.ID) * pi

    Laa = μ * CircleGeom.Turns^2 * R^2 / (2 * T)#Another formula for inductance of circular toroid
    #
    # println("Target inductance is: $(round(LTarget*1e6)) μH")
    # println(
    #     "Sanity check inductance is: $(round(SanityCheckInductance*1e6)) μH",
    # )
    # println("
    # FlatHeight     :$(round(DCoreGeom.FlatHeight*1000;sigdigits=2)) mm
    # MaxHeight      :$(round(DCoreGeom.MaxHeight*1000;sigdigits=2)) mm
    # RadiusAtPeak   :$(round(DCoreGeom.RadiusAtPeak*1000;sigdigits=2)) mm
    # ID             :$(round(DCoreGeom.ID*1000;sigdigits=2)) mm
    # OD             :$(round(DCoreGeom.OD*1000;sigdigits=2)) mm
    # Turns          :$(DCoreGeom.Turns) Turns
    # Layers         :$(DCoreGeom.Layers) Turns
    # WireLength     :$(round(DCoreGeom.WireLength;sigdigits=2)) meters
    # WireResist     :$(round(DCoreGeom.Resistance;sigdigits=2)) Ω")
    CoilGeom = Geom(DCoreGeom, CircleGeom, General)

    if ~(ExportFileName===nothing)
        DCoreParams2SVG(DCoreGeom.ID*1000 , DCoreGeom.OD*1000,"DCoreGeom"*ExportFileName;NPts = NPts)
    end
    return CoilGeom
end

function Length2Resist(L, Diam; ρ = 1.68e-8, FillFac = 1)
    if Diam > 0.1
        println(
            "Diameter must be in METERS, it was automatically multiplied by 1e-3 in case you entered it in mm",
        )
        Diam = Diam * 1e-3
    end
    A = pi / 4 .* Diam .^ 2
    return L * ρ / FillFac / A
end

"""
This function takes in an inner and outer radius and makes the ideal "D-shaped" toroid core

It returns a list of coordinates (r,z) where r is the radial location
on the toroid, and z is the height from the center axis. In a cross-section
view this is essentially x and y.

r₁ - inner radius
r₂ - outer radius
dr (kwarg) is the dr between each sampled location

"""
function DCoreGeom(r₁, r₂; dr = 0.0001,NPts = nothing,PlotOn=false,UpsamplePoints = 1e4)
    dZdr(r) = log.(sqrt(r₁ * r₂) ./ r) ./ sqrt.(log.(r ./ r₁) .* log.(r₂ ./ r))
    if ~(NPts===nothing)
        dr = (r₂ - r₁) / UpsamplePoints
    end
    rᵢ = r₁ + dr / 10#starting at a small fraction above zero, to avoid an inf.
    r = [rᵢ]
    zᵢ = 0.0
    z = [zᵢ] # initializing vector
    II = 1 #iterator
    NumIters = floor((r₂ - r₁) / dr)

    while (NumIters > II)
        dz = dZdr(rᵢ)
        zᵢ +=dz*dr
        rᵢ += dr

        push!(z, zᵢ)
        push!(r, rᵢ)
        II += 1
    end
    z .-= z[end] #the boundary condition at the end must be fixed to equal 0

    z[1] = 0 #because of the inf initial dz/dr, the initial boundary must be manually set to zero
    ZFlatHeight = z[2]
    StartPts = Int(ceil(ZFlatHeight/dr))
    for i in 1:StartPts
        insert!(z,1+i,z[1+i]*i/StartPts)
        insert!(r,1+i,r[1])
    end

    z = vcat(z[:],-1 .*reverse(z[2:end-2]))
    r = vcat(r[:],reverse(r[2:end-2]))
    DownsampleIndex = 1:Int(round(UpsamplePoints/NPts)):length(z)
    z = z[DownsampleIndex]
    r = r[DownsampleIndex]
    if PlotOn
        plot(r, z)
    end
    return r, z
end

"""
This function takes in the parameters of a circular Rogowski coil
    N: Num turns,
    ID: ID of circle in meters
    OD: OD of circle in meters
    and kwargs:
    Current: Current in the wire being sensed
    μ: permeability of the material
    ω: frequency (rad/sec) of current
"""
function Rogowski_Calc(N,ID,OD;Current=1,μ=4*π*1e-7,ω = 25e3*2*π)

    M₂₁ =μ*N/2*(ID+OD-2*√(ID*OD)) #From multiple sources.

    V = ω*M₂₁*Current
    return (M₂₁,V)
end

"""
This function takes in the parameters of a rectagular Rogowski coil
    N: Num turns,
    h: height of rectangle in meters
    ID: ID of rectangle in meters
    OD: OD of rectangle in meters
    and kwargs:
    Current: Current in the wire being sensed
    μ: permeability of the material
    ω: frequency (rad/sec) of current
"""
function Rogowski_Calc(N,h,ID,OD;Current=1,μ=4*π*1e-7,ω = 25e3*2*π)

    M₂₁ =μ*N/(2*π)*h*(log(OD/ID)) #From multiple sources.

    V = ω*M₂₁*Current
    return (M₂₁,V)
end

function DCore_DetermineIdealInduct(ID,N,α)
    #Equation 7
     # AlphaMat Key. The columns represent [alpha e h p s t]
     μ₀ = 4*pi*1e-7
     AlphaMat = [
        2 0.26  0.645 3.6   0.72  1.77
        3 0.85  1.5   8.0   2.74  4.42
        4 1.6   2.4   12.8  5.76  8.09
        5 2.45  3.4   17.9  9.61  12.6
        6 3.4   4.5   23.3  14.2  17.8
        7 4.4   5.6   28.9  19.4  23.7
    ] #Table 1 From  "D-Shaped toroidal cage inductors" P.N. Murgatroyd & D. Belahrache,1989

    T = AlphaMat[α-1, 6] #Picking a value in the table for clarity
    S = AlphaMat[α-1, 5] #Also Selecting value for clarity
    #S₁ = AlphaMat[α, 5] #Also Selecting value for clarity
    P = AlphaMat[α-1, 4] # As above
    E = AlphaMat[α-1, 2] # As above
    B = ID/2 #For consistency with the paper
    #L = μ₀*N^2*B / (2*π)*((2*E+1)/4 + 2/3*S + S₁)
    println(T)
    L = T * μ₀*N^2*B / (2*π)
    return L
end

function CircCore_DetermineIdealInduct(IRad,N,α)
    IRad = IRad*100
    ORad = IRad*α
    R = (ORad+IRad)/2
    a = (ORad-IRad)/2
    return 1e-6*(0.01257*N^2*(R-√*(R^2-a^2)))
end


text:
"""
    TimeBoundary(dataSource::DataSource; <keyword arguments>)

Returns a Dict{String, String} with "timestamp", "minTime" and/or "maxTime"
(depending on `bound`) as keys.

A method with all arguments as keyword arguments is also provided.

# Arguments
- bound::String = nothing
- filter::Filter = nothing
- virtualColumns::Vector{<:VirtualColumn} = nothing
- context::Dict = nothing
"""
mutable struct TimeBoundary <: Query
    queryType::String
    dataSource::DataSource
    bound
    filter
    virtualColumns
    context
    function TimeBoundary(dataSource; bound=nothing, filter=nothing, virtualColumns=nothing, context=nothing)
        bound === nothing || bound ∈ ["minTime", "maxTime"] || error("Invalid bound")
        nothing_or_type(filter, Filter)
        nothing_or_type(virtualColumns, Vector{<:VirtualColumn})
        nothing_or_type(context, Dict)
        new("timeBoundary", dataSource, bound, filter, virtualColumns, context)
    end
end
TimeBoundary(; dataSource, bound=nothing, filter=nothing, virtualColumns=nothing, context=nothing) =
    TimeBoundary(dataSource; bound, filter, virtualColumns, context)

function execute(client::Client, query::TimeBoundary; pretty=false)
    res = JSON.parse(execute_native_query(client, query; pretty))
    d = Dict{String, String}()
    d["timestamp"] = res[1]["timestamp"]
    if haskey(res[1]["result"], "minTime")
        d["minTime"] = res[1]["result"]["minTime"]
    end
    if haskey(res[1]["result"], "maxTime")
        d["maxTime"] = res[1]["result"]["maxTime"]
    end
    d
end