text:
function updateliver!(system::CVSystem,n::Int64,hemoflag="no")
    inflow = Float64(0);

    # sum liver inflows
    for i = 1:length(system.branches.ID)
        if isempty(system.branches.children[i])
            if i == 12 || i == 13 || i == 15 # portal vein inputs
                inflow += system.branches.term[i].Q[n+1,5];
            end
        end
    end

    # Euler forward update to liver
    system.liver.V[n+2,1] = (system.liver.V[n+1,1] + system.solverparams.h*(
        inflow - system.liver.Q[n+1,1]));
    if hemoflag == "no" || (hemoflag == "yes" && system.hemo.injured[14] == false) # hemorrhage not at HA
        system.liver.Q[n+2,1] = ((1 - system.liver.R[1]/system.liver.L*
            system.solverparams.h)*system.liver.Q[n+1,1]+
            system.solverparams.h/system.liver.L*(system.liver.P[n+1,1]-
            system.liver.P[n+1,2]));
        system.liver.V[n+2,2] = (system.liver.V[n+1,2]+
            system.solverparams.h*(system.liver.Q[n+1,1] - system.liver.Q[n+1,2]));
    elseif hemoflag == "yes" && system.hemo.injured[14] == true # sever PV connection to sinusoids if HA severed
        system.liver.Q[n+2,1] = ((1 - system.liver.R[1]/system.liver.L*
            system.solverparams.h)*system.liver.Q[n+1,1]+
            system.solverparams.h/system.liver.L*(system.liver.P[n+1,1]-
            system.hemo.Ph));
        system.liver.V[n+2,2] = (system.liver.V[n+1,2]-
            system.solverparams.h*system.liver.Q[n+1,2]);
        system.hemo.Vloss += system.solverparams.h*system.liver.Q[n+2,1];
    end
    system.liver.V[n+2,3] = (system.liver.V[n+1,3]+
        system.solverparams.h*(system.liver.Q[n+1,2] - system.liver.Q[n+1,3]));
    system.liver.V[n+2,4] = (system.liver.V[n+1,4]+
        system.solverparams.h*(system.branches.term[14].Q[n+1,5] + system.liver.Q[n+1,3] - system.liver.Q[n+1,4]));

    # note pressure, flowrate updated algebraically except portal vein compartment
    system.liver.P[n+2,1] = (1/system.liver.C[1]*
        (system.liver.V[n+2,1] - system.liver.V0[1]));
    system.liver.P[n+2,2] = (1/system.liver.C[2]*
        (system.liver.V[n+2,2] - system.liver.V0[2]));
    system.liver.P[n+2,3] = (1/system.liver.C[3]*
        (system.liver.V[n+2,3] - system.liver.V0[3]));
    system.liver.P[n+2,4] = (1/system.liver.C[4]*
        (system.liver.V[n+2,4] - system.liver.V0[4]));
    system.liver.Q[n+2,2] = ((system.liver.P[n+2,2] - system.liver.P[n+2,3])/
        system.liver.R[2]);
    system.liver.Q[n+2,3] = ((system.liver.P[n+2,3] - system.liver.P[n+2,4])/
        system.liver.R[3]);
    system.liver.Q[n+2,4] = ((system.liver.P[n+2,4] - system.ivc.P[n+2])/
        system.liver.R[4]);
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
using ClimateMachine
const clima_dir = dirname(dirname(pathof(ClimateMachine)));

if parse(Bool, get(ENV, "CLIMATEMACHINE_PLOT_EDMF_COMPARISON", "false"))
    plot_dir = joinpath(clima_dir, "output", "bomex_edmf", "pycles_comparison")
else
    plot_dir = nothing
end

include(joinpath(@__DIR__, "compute_mse.jl"))

data_file = Dataset(joinpath(PyCLES_output_dataset_path, "Bomex.nc"), "r")

#! format: off
best_mse = OrderedDict()
best_mse["prog_ρ"] = 3.4917662870525452e-02
best_mse["prog_ρu_1"] = 3.0715053983126868e+03
best_mse["prog_ρu_2"] = 1.2895738234435595e-03
best_mse["prog_moisture_ρq_tot"] = 4.1331426262587470e-02
best_mse["prog_turbconv_environment_ρatke"] = 6.7642719033845833e+02
best_mse["prog_turbconv_environment_ρaθ_liq_cv"] = 8.5667226792770009e+01
best_mse["prog_turbconv_environment_ρaq_tot_cv"] = 1.6435640302009438e+02
best_mse["prog_turbconv_updraft_1_ρa"] = 8.0846149716812874e+01
best_mse["prog_turbconv_updraft_1_ρaw"] = 8.5071202823418249e-02
best_mse["prog_turbconv_updraft_1_ρaθ_liq"] = 9.0386633006301569e+00
best_mse["prog_turbconv_updraft_1_ρaq_tot"] = 1.0803376910477068e+01
#! format: on

computed_mse = compute_mse(
    solver_config.dg.grid,
    solver_config.dg.balance_law,
    time_data,
    dons_arr,
    data_file,
    "Bomex",
    best_mse,
    400,
    plot_dir,
)

@testset "BOMEX EDMF Solution Quality Assurance (QA) tests" begin
    #! format: off
    test_mse(computed_mse, best_mse, "prog_ρ")
    test_mse(computed_mse, best_mse, "prog_ρu_1")
    test_mse(computed_mse, best_mse, "prog_ρu_2")
    test_mse(computed_mse, best_mse, "prog_moisture_ρq_tot")
    test_mse(computed_mse, best_mse, "prog_turbconv_environment_ρatke")
    test_mse(computed_mse, best_mse, "prog_turbconv_environment_ρaθ_liq_cv")
    test_mse(computed_mse, best_mse, "prog_turbconv_environment_ρaq_tot_cv")
    test_mse(computed_mse, best_mse, "prog_turbconv_updraft_1_ρa")
    test_mse(computed_mse, best_mse, "prog_turbconv_updraft_1_ρaw")
    test_mse(computed_mse, best_mse, "prog_turbconv_updraft_1_ρaθ_liq")
    test_mse(computed_mse, best_mse, "prog_turbconv_updraft_1_ρaq_tot")
    #! format: on
end

ext: jl


text:
# import Base.copy
import Base.-
import Base./

export
    Vector2D,
    set!

mutable struct Vector2D{T <: AbstractFloat}
    x::T
    y::T

    function Vector2D{T}() where {T <: AbstractFloat}
        new(0.0, 0.0)
    end

    function Vector2D{T}(x::T, y::T) where {T <: AbstractFloat}
        new(x, y)
    end
end

# copy specializations
copy(v::Vector2D) = Vector2D{Float64}(v.x, v.y)

# setters/getters set_point
function set!(v::Vector2D{T}, x::T, y::T) where {T <: AbstractFloat}
    v.x = x
    v.y = y;
end

function set!(v::Vector2D{T}, from::Vector2D{T}) where {T <: AbstractFloat}
    v.x = from.x
    v.y = from.y;
end

function length(x::T, y::T) where {T <: AbstractFloat}
    sqrt(x * x + y * y)
end

function length(v::Vector2D{T}) where {T <: AbstractFloat}
    sqrt(v.x * v.x + v.y * v.y)
end

function length_sq(x::T, y::T) where {T <: AbstractFloat}
    x * x + y * y
end

# Scratch vectors
scv = Vector2D{Float64}()
scv2 = Vector2D{Float64}()

function add!(v1::Vector2D{T}, v2::Vector2D{T}, out::Vector2D{T}) where {T <: AbstractFloat}
    set!(out, v1.x + v2.x, v1.y + v2.y)
end

function sub!(v1::Vector2D{T}, v2::Vector2D{T}, out::Vector2D{T}) where {T <: AbstractFloat}
    set!(out, v1.x - v2.x, v1.y - v2.y)
end

function scale!(v::Vector2D{T}, s::T) where {T <: AbstractFloat}
    set!(v, v.x * s, v.y * s)
end

function scale!(v::Vector2D{T}, s::T, out::Vector2D{T}) where {T <: AbstractFloat}
    set!(out, v.x * s, v.y * s)
end

function div!(v::Vector2D{T}, value::T) where {T <: AbstractFloat}
    set!(v, v.x / value, v.y / value)
end

# distance between two vectors
function distance(v1::Vector2D{T}, v2::Vector2D{T}) where {T <: AbstractFloat}
    sub!(v1, v2, scv)
    length(scv)
end

# returns the angle in radians between this vector and the x axis
function angle(v::Vector2D{T}) where {T <: AbstractFloat}
    atan(v.y, v.x)
end

# Returns  multiplied to a length of 1.
# If the point is 0, it returns (1, 0)
function normalize!(v::Vector2D{T}) where {T <: AbstractFloat}
    len = length(v)
    if length(v) == 0.0
        set!(scv, 1.0, 0.0)
        return scv
    end

    div!(v, len);
end

# Calculates dot product of two points.
function dot(v1::Vector2D{T}, v2::Vector2D{T}) where {T <: AbstractFloat}
    v1.x * v2.x + v1.y * v2.y
end

# Calculates cross product of two points.
function cross(v1::Vector2D{T}, v2::Vector2D{T}) where {T <: AbstractFloat}
    v1.x * v2.y - v1.y * v2.x
end

# returns the angle in radians between two vector directions
function angle_between(v1::Vector2D{T}, v2::Vector2D{T}) where {T <: AbstractFloat}
    set!(scv, v1)
    set!(scv2, v2)

    normalize!(scv) #a2
    normalize!(scv2) #b2

    angle = atan(cross(scv, scv2), dot(scv, scv2))

    if abs(angle) < EPSILON
        return 0.0
    else
        return angle
    end
end

# Calculates perpendicular of v, rotated 90 degrees counter-clockwise -- cross(v, perp(v)) >= 0
function perpindicular(x::T, y::T) where {T <: AbstractFloat}
    Vector2D(-y, x)
end

function ccw_perpindicular!(v::Vector2D{T}) where {T <: AbstractFloat}
    set!(v, -v.y, v.x);
end

# Calculates perpendicular of v, rotated 90 degrees clockwise -- cross(v, rperp(v)) <= 0
function cw_perpindicular!(v::Vector2D{T}) where {T <: AbstractFloat}
    set!(v, v.y, -v.x);
end

function set_direction!(v::Vector2D{T}, degrees::T) where {T <: AbstractFloat}
    set!(v, cos(deg2rad(degrees)), sin(deg2rad(degrees)));
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
# ===========================
# Homogeneous case
# ===========================

function reach_homog_VREP!(F::Vector{ReachSet{N, VP}},
                           Ω0::VP,
                           Φ::AbstractMatrix,
                           NSTEPS::Integer,
                           δ::Float64,
                           X::Universe,
                           Δt0::TimeInterval) where {N, VN, VP<:VPOLY{N, VN}}

    # initial reach set
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = ReachSet(Ω0, Δt)

    k = 2
    @inbounds while k <= NSTEPS
        Rₖ = linear_map(Φ, set(F[k-1]))
        Δt += δ
        F[k] = ReachSet(Rₖ, Δt)
        k += 1
    end
    return F
end


text:
# Vector

function _pyconvert_rule_iterable(ans::Vector{T0}, it::Py, ::Type{T1}) where {T0,T1}
    @label again
    x_ = pynext(it)
    if ispynull(x_)
        pydel!(it)
        return pyconvert_return(ans)
    end
    x = @pyconvert_and_del(T1, x_)
    if x isa T0
        push!(ans, x)
        @goto again
    end
    T2 = Utils._promote_type_bounded(T0, typeof(x), T1)
    ans2 = Vector{T2}(ans)
    push!(ans2, x)
    return _pyconvert_rule_iterable(ans2, it, T1)
end

function pyconvert_rule_iterable(::Type{R}, x::Py, ::Type{Vector{T0}}=Utils._type_lb(R), ::Type{Vector{T1}}=Utils._type_ub(R)) where {R<:Vector,T0,T1}
    it = pyiter(x)
    ans = Vector{T0}()
    return _pyconvert_rule_iterable(ans, it, T1)
end

# Set

function _pyconvert_rule_iterable(ans::Set{T0}, it::Py, ::Type{T1}) where {T0,T1}
    @label again
    x_ = pynext(it)
    if ispynull(x_)
        pydel!(it)
        return pyconvert_return(ans)
    end
    x = @pyconvert_and_del(T1, x_)
    if x isa T0
        push!(ans, x)
        @goto again
    end
    T2 = Utils._promote_type_bounded(T0, typeof(x), T1)
    ans2 = Set{T2}(ans)
    push!(ans2, x)
    return _pyconvert_rule_iterable(ans2, it, T1)
end

function pyconvert_rule_iterable(::Type{R}, x::Py, ::Type{Set{T0}}=Utils._type_lb(R), ::Type{Set{T1}}=Utils._type_ub(R)) where {R<:Set,T0,T1}
    it = pyiter(x)
    ans = Set{T0}()
    return _pyconvert_rule_iterable(ans, it, T1)
end

# Dict

function _pyconvert_rule_mapping(ans::Dict{K0,V0}, x::Py, it::Py, ::Type{K1}, ::Type{V1}) where {K0,V0,K1,V1}
    @label again
    k_ = pynext(it)
    if ispynull(k_)
        pydel!(it)
        return pyconvert_return(ans)
    end
    v_ = pygetitem(x, k_)
    k = @pyconvert_and_del(K1, k_)
    v = @pyconvert_and_del(V1, v_)
    if k isa K0 && v isa V0
        push!(ans, k => v)
        @goto again
    end
    K2 = Utils._promote_type_bounded(K0, typeof(k), K1)
    V2 = Utils._promote_type_bounded(V0, typeof(v), V1)
    ans2 = Dict{K2,V2}(ans)
    push!(ans2, k => v)
    return _pyconvert_rule_mapping(ans2, x, it, K1, V1)
end

function pyconvert_rule_mapping(::Type{R}, x::Py, ::Type{Dict{K0,V0}}=Utils._type_lb(R), ::Type{Dict{K1,V1}}=Utils._type_ub(R)) where {R<:Dict,K0,V0,K1,V1}
    it = pyiter(x)
    ans = Dict{K0,V0}()
    return _pyconvert_rule_mapping(ans, x, it, K1, V1)
end

# Tuple

function pyconvert_rule_iterable(::Type{T}, xs::Py) where {T<:Tuple}
    T isa DataType || return pyconvert_unconverted()
    ts = collect(T.parameters)
    if !isempty(ts) && Base.isvarargtype(ts[end])
        isvararg = true
        vartype = ts[end].body.parameters[1]::Type
        pop!(ts)
    else
        isvararg = false
        vartype = Union{}
    end
    zs = Any[]
    for x in xs
        if length(zs) < length(ts)
            t = ts[length(zs) + 1]
        elseif isvararg
            t = vartype
        else
            return pyconvert_unconverted()
        end
        z = @pyconvert_and_del(t, x)
        push!(zs, z)
    end
    return length(zs) < length(ts) ? pyconvert_unconverted() : pyconvert_return(T(zs))
end

# N-Tuple for N up to 16
# TODO: Vararg

for N in 0:16
    Ts = [Symbol("T", n) for n in 1:N]
    zs = [Symbol("z", n) for n in 1:N]
    @eval function pyconvert_rule_iterable(::Type{Tuple{$(Ts...)}}, xs::Py) where {$(Ts...)}
        pylen(xs) == $N || return pyconvert_unconverted()
        $((
            :($z = @pyconvert_and_del($T, pytuple_getitem(xs, $(i-1))))
            for (i, T, z) in zip(1:N, Ts, zs)
        )...)
        return pyconvert_return(($(zs...),))
    end
end

# Pair

function pyconvert_rule_iterable(::Type{R}, x::Py, ::Type{Pair{K0,V0}}=Utils._type_lb(R), ::Type{Pair{K1,V1}}=Utils._type_ub(R)) where {R<:Pair,K0,V0,K1,V1}
    it = pyiter(x)
    k_ = pynext(it)
    if ispynull(k_)
        pydel!(it)
        pydel!(k_)
        return pyconvert_unconverted()
    end
    k = @pyconvert_and_del(K1, k_)
    v_ = pynext(it)
    if ispynull(v_)
        pydel!(it)
        pydel!(v_)
        return pyconvert_unconverted()
    end
    v = @pyconvert_and_del(V1, v_)
    z_ = pynext(it)
    pydel!(it)
    if ispynull(z_)
        pydel!(z_)
    else
        pydel!(z_)
        return pyconvert_unconverted()
    end
    K2 = Utils._promote_type_bounded(K0, typeof(k), K1)
    V2 = Utils._promote_type_bounded(V0, typeof(v), V1)
    return pyconvert_return(Pair{K2,V2}(k, v))
end


text:
module ReferringExpressions

using PyCall, JSON, Iterators
@pyimport cPickle as pickle
const DATADIR = Pkg.dir("ReferringExpressions","data")

struct RefExpData
    # inputs
    datadir
    dataset
    datasplit

    # storage
    imgdir
    refs
    annotations
    images
    categories
    sentences
    img2ref
    img2ann
    ref2ann
    ann2ref
    cat2ref
    sent2ref
    sent2tok

    function RefExpData(datadir, dataset, datasplit)
        dicts = (:data, :refs, :annotations, :images, :categories, :sentences,
                 :img2ref, :img2ann, :ref2ann, :ann2ref, :cat2ref, :sent2ref,
                 :sent2tok)
        for dict in dicts
            @eval $dict = Dict()
        end

        imgdir = _get_image_dir(datadir, dataset)
        data["dataset"] = dataset
        data["refs"] = _load_refs(datadir, dataset, datasplit)
        instances = _get_instances(datadir, dataset)
        data["images"] = instances["images"]
        data["annotations"] = instances["annotations"]
        data["categories"] = instances["categories"]

        for ann in data["annotations"]
            id, iid = ann["id"], ann["image_id"]
            annotations[id] = ann
            img2ann[iid] = vcat(get(img2ann, iid, []), ann)
        end

        for img in data["images"]
            images[img["id"]] = img
        end

        for cat in data["categories"]
            categories[cat["id"]] = cat["name"]
        end

        for ref in data["refs"]
            ref_id = ref["ref_id"]
            ann_id = ref["ann_id"]
            category_id = ref["category_id"]
            image_id = ref["image_id"]

            refs[ref_id] = ref
            img2ref[image_id] = vcat(get(img2ref, image_id, []), ref)
            cat2ref[category_id] = vcat(get(cat2ref, category_id, []), ref)
            ref2ann[ref_id] = annotations[ann_id]
            ann2ref[ann_id] = ref

            for sent in ref["sentences"]
                sent_id = sent["sent_id"]
                sentences[sent_id] = sent
                sent2ref[sent_id] = ref
                sent2tok[sent_id] = sent["tokens"]
            end
        end

        return new(datadir,dataset,datasplit,imgdir,refs,annotations,
                   images,categories,sentences,img2ref,img2ann,ref2ann,
                   ann2ref,cat2ref,sent2ref,sent2tok)
    end

    RefExpData(dataset, datasplit) = RefExpData(DATADIR, dataset, datasplit)
end
export RefExpData

# load methods
for obj in (:refs, :annotations, :images, :categories)
    fun = Symbol(:load_, obj)
    @eval begin
        $fun{T<:Int}(d::RefExpData, id::T) = [d.$(obj)[id]]
        $fun{T<:Array}(d::RefExpData, ids::T) = map(id->d.$(obj)[id], ids)
        export $fun
    end
end

# filter methods
function filter_refs_by_split(d::RefExpData, split="")
    filter_refs_by_split(d.refs, split)
end
function filter_refs_by_split(refs, split="")
    isempty(split) && return refs
    if in(split, ("testA","testB","testC"))
        refs = filter((k,v) -> in(string(split[end]), v["split"]), refs)
    elseif in(split, ("testAB","testBC","testAC"))
        refs = filter((k,v) -> v["split"] == split, refs)
    elseif in(split, ("train","valid","test"))
        refs = filter((k,v) -> startswith(v["split"], split), refs)
    end
    return refs
end
export filter_refs_by_split

_filters = (
    (:categories, "category_id"), (:refs, "ref_id"), (:images, "image_id"))
_arrays = (:ref, :image, :annotation)
for ((_filter,_id), _array) in product(_filters, _arrays)
    F = Symbol(:filter_, _array, :s_by_, _filter)
    @eval $F(d::RefExpData, ids) = filter(
        (k,v) -> in(v[$_id], ids) , d.$_filter)
    @eval $F(data, ids) = filter((k,v) -> in(v[$_id], ids) , data)
    @eval export $F
end

# get methods
get_category_ids(d::RefExpData) = keys(d.categories)
export get_category_ids

function get_ref_ids(
    d::RefExpData; ref_ids=[], cat_ids=[], image_ids=[], split="")
    refs = d.refs
    refs = filter_refs_by_refs(refs, ref_ids)
    refs = filter_refs_by_categories(refs, cat_ids)
    refs = filter_refs_by_images(refs, image_ids)
    refs = filter_refs_by_split(refs, split)
end
export get_ref_ids

function get_annotations_ids(
    d::RefExpData; ref_ids=[], cat_ids=[], image_ids=[])
    data = d.annotations
    data = filter_annotations_by_refs(data, ref_ids)
    data = filter_annotations_by_categories(data, cat_ids)
    data = filter_annotations_by_images(data, image_ids)
end
export get_annotations_ids

# utils
function _get_image_dir(datadir, dataset)::String
    subdir = ""
    if in(dataset, ("refcoco", "refcoco+", "refcocog"))
        subdir = "images/mscoco/images/train2014"
    elseif dataset == "refclef"
        subdir = "images/saiapr_tc-12"
    else
        error("There's no such dataset called as $dataset.")
    end
    return abspath(joinpath(datadir, subdir))
end

function _load_refs(datadir, dataset, datasplit)
    refsfile = abspath(joinpath(datadir, dataset, "refs("*datasplit*").p"))
    refs = open(refsfile, "r") do f
        pickle.loads(pybytes(read(f)))
    end
    return refs
end

function _get_instances(datadir, dataset)
    jsonfile = abspath(joinpath(datadir, dataset, "instances.json"))
    return JSON.parsefile(jsonfile)
end

end # module


text:
mutable struct MCIntegrationVI{T<:Real} <: NumericalVI{T}
    nMC::Int64 #Number of samples for MC Integrations
    clipping::T
    ϵ::T #Convergence criteria
    n_iter::Integer #Number of steps performed
    stoch::Bool # Flag for stochastic optimization
    batchsize::Int #Size of mini-batches
    ρ::T # Scaling coeff. for stoch. opt.
    NaturalGradient::Bool
    HyperParametersUpdated::Bool # Flag for updating kernel matrix
    vi_opt::NVIOptimizer
    function MCIntegrationVI{T}(
        ϵ::T,
        nMC::Int,
        optimiser,
        Stochastic::Bool,
        clipping::Real,
        batchsize::Int,
        natural::Bool,
    ) where {T<:Real}
        return new{T}(
            nMC,
            clipping,
            ϵ,
            0,
            Stochastic,
            batchsize,
            one(T),
            natural,
            true,
            NVIOptimizer(optimiser),
        )
    end
end

"""
    MCIntegrationVI(;ϵ::T=1e-5, nMC::Integer=1000, clipping::Real=Inf, natural::Bool=true, optimiser=Momentum(0.001))

Variational Inference solver by approximating gradients via MC Integration.
It means the expectation `E[log p(y|f)]` as well as its gradients is computed
by sampling from q(f).

## Keyword arguments
- `ϵ::Real` : convergence criteria, which can be user defined
- `nMC::Int` : Number of samples per data point for the integral evaluation
- `clipping::Real` : Limit the gradients values to avoid overshooting
- `natural::Bool` : Use natural gradients
- `optimiser` : Optimiser used for the variational updates. Should be an Optimiser object from the [Flux.jl](https://github.com/FluxML/Flux.jl) library, see list here [Optimisers](https://fluxml.ai/Flux.jl/stable/training/optimisers/) and on [this list](https://github.com/theogf/AugmentedGaussianProcesses.jl/tree/master/src/inference/optimisers.jl). Default is `Momentum(0.01)`
"""
MCIntegrationVI

function MCIntegrationVI(;
    ϵ::T=1e-5,
    nMC::Integer=1000,
    optimiser=Momentum(0.01),
    clipping::Real=Inf,
    natural::Bool=true,
) where {T<:Real}
    return MCIntegrationVI{T}(ϵ, nMC, optimiser, false, clipping, 1, natural)
end

"""
    MCIntegrationSVI(batchsize::Int; ϵ::Real=1e-5, nMC::Integer=1000, clipping=Inf, natural=true, optimiser=Momentum(0.0001))

Stochastic Variational Inference solver by approximating gradients via Monte Carlo integration when using minibatches
See [`MCIntegrationVI`](@ref) for more explanations.

## Argument

-`batchsize::Integer` : Number of samples per mini-batches

## Keyword arguments

- `ϵ::T` : convergence criteria, which can be user defined
- `nMC::Int` : Number of samples per data point for the integral evaluation
- `clipping::Real` : Limit the gradients values to avoid overshooting
- `natural::Bool` : Use natural gradients
- `optimiser` : Optimiser used for the variational updates. Should be an Optimiser object from the [Flux.jl](https://github.com/FluxML/Flux.jl) library, see list here [Optimisers](https://fluxml.ai/Flux.jl/stable/training/optimisers/) and on [this list](https://github.com/theogf/AugmentedGaussianProcesses.jl/tree/master/src/inference/optimisers.jl). Default is `Momentum()` (ρ=(τ+iter)^-κ)
"""
MCIntegrationSVI

function MCIntegrationSVI(
    batchsize::Integer;
    ϵ::T=1e-5,
    nMC::Integer=200,
    optimiser=Momentum(0.001),
    clipping::Real=0.0,
    natural::Bool=true,
) where {T<:Real}
    return MCIntegrationVI{T}(ϵ, nMC, optimiser, true, clipping, batchsize, natural)
end

function grad_expectations!(
    m::AbstractGPModel{T,L,<:MCIntegrationVI{T}}, state, y
) where {T,L}
    raw_samples = randn(T, inference(m).nMC, n_latent(m))
    samples = similar(raw_samples)
    μ = mean_f(m, state.kernel_matrices)
    σ² = var_f(m, state.kernel_matrices)
    num_sample = batchsize(m)
    opt_state = state.opt_state
    for j in 1:num_sample # Loop over every data point
        samples .=
            raw_samples .* sqrt.([σ²[k][j] for k in 1:n_latent(m)])' .+ [μ[k][j] for k in 1:n_latent(m)]'
        grad_samples!(m, samples, opt_state, y[j, :], j) # Compute the gradient for data point j
    end
end

function expec_loglikelihood(
    l::AbstractLikelihood, i::MCIntegrationVI{T}, y, μ, σ²
) where {T} # μ and σ² are tuples of vectors
    num_latent = n_latent(l)
    raw_samples = randn(T, i.nMC, num_latent) # dimension nMC x nLatent
    # samples = similar(raw_samples)
    num_sample = batchsize(i)
    tot = 0.0
    for j in 1:num_sample # Loop over every data point
        samples =
            raw_samples .* sqrt.([σ²[k][j] for k in 1:num_latent])' .+ [μ[k][j] for k in 1:num_latent]'
        # samples is of dimension nMC x nLatent again
        y_j = y[j, :] # Obtain the label for data point j
        for f in eachrow(samples)
            # We now compute the loglikelihood over every sample
            tot += loglikelihood(l, y_j, f)
        end
    end
    return tot / i.nMC
end


text:
## Parser.jl --- functions for parsing PDB files

function Base.read(filename::AbstractString, ::Type{StructureFrame}; kwargs...)
    @assert endswith(filename, ".pdb") "only `.pdb` files are supported for read"
    open(filename) do io
        read(io, StructureFrame; kwargs...)
    end
end

function Base.read(io::IO, ::Type{StructureFrame}; ensemble=false, sort=false)
    ## for each in IO
    # check if it has MODEL record and create frame components
    # populate the frame components accordingly
    # read only the first model if ensemble is false
    ensemble && throw(ArgumentError("ensemble = true is not yet supported"))

    ## create frame components: step = 0, at_pos, at_list, res_list
    at_pos   = SVector{3,Float64}[]
    at_list  = JAtom{Float64}[]
    res_list = JResidue[]

    # read the records
    supposed_to_be_ensemble = false
    finished_with_atoms = false
    ## residue markers
    curr_residue = :o; curr_chain = :k; curr_seqnum = 0; curr_pdb = false
    res_atom_dict = Dict{Symbol,UInt}()
    for (i, line) in enumerate(eachline(io))
        supposed_to_be_ensemble |= startswith(line, "MODEL")

        ## if starts with any of HETATM or ATOM create a JAtom
        ## at the same time take note of residues; when traversed, create a residue
        if startswith(line, "ATOM") || startswith(line, "HETATM")
            jatom, pos, res_, chain_, seqnum_, pdb_ = parseatomline(line)

            ## case for "disordered atom"; use only default one
            isnothing(jatom) && continue
            push!(at_pos, pos); push!(at_list, jatom)

            ## check both to make sure
            if res_ != curr_residue || chain_ != curr_chain || curr_seqnum != seqnum_
                # do the necessary updates
                if !isempty(res_atom_dict)
                    ch_str = curr_chain |> String
                    jresidue = JResidue(curr_residue, ch_str, res_atom_dict, curr_pdb)
                    push!(res_list, jresidue)
                end
                curr_residue = res_; curr_chain = chain_
                curr_seqnum = seqnum_; curr_pdb = pdb_
                res_atom_dict = Dict{Symbol,UInt}() ## allocate new dict
            end

            ## add atom to its residue
            res_atom_dict[jatom.name] = length(at_pos)
        end

        ## small hack: in reaching end, push last recorded residue
        if !finished_with_atoms && startswith(line, "MASTER")
            ch_str = curr_chain |> String
            jresidue = JResidue(curr_residue, ch_str, res_atom_dict, curr_pdb)
            push!(res_list, jresidue)
            finished_with_atoms = true
        end

        if supposed_to_be_ensemble && !ensemble
            ## for now terminate early if found endmdl
            startswith(line, "ENDMDL") && break
        end
    end

    if sort
        ## apparently false < true; if protein residues are first, then rev=true
        sort!(res_list, lt=(x,y)->isless(x.standard_pdb,y.standard_pdb), rev=true)
    end

    return StructureFrame(0, at_pos, at_list, res_list)
end

function parseatomline(line)
    ## create an atom
    ## also return information about the chain and the residue
    (line[17] != 'A' && line[17] != ' ') && return ntuple(i->nothing, 5)

    pdb_ = startswith(line, "ATOM")
    name_ = line[13:16] |> strip |> q->replace(q, '\''=>'p') |> Symbol

    res_name = line[18:20] |> strip |> Symbol
    chainid = line[22] |> Symbol
    seqnum = line[23:26] |> strip |> q->parse(Int, q)

    pos_str = line[31:38], line[39:46], line[47:54]
    pos_vec = map(pos_str) do xistr
        ## try to parse
        xistr |> strip |> q->parse(Float64, q)
    end |> SVector{3,Float64}

    ## for now initialize mass and charge to 0; they're not needed anyway
    ## TODO: as for type; it's not needed but should be inferred directly
    jatom = JAtom(name_, "", zero(Float64), zero(Float64))
    jatom, pos_vec, res_name, chainid, seqnum, pdb_
end


text:
# # Extend Base operations for my types or existing concrete types only
# # else, fallback to methods AbstractMatrix{Float64}

array_types = (:(SparseArrays.SparseMatrixCSC{Float64,Int}),:(Matrix{Float64}))
my_types = (:LazyGenerator, :FullGenerator, :SFMDistribution)
for mt in my_types
    for at in (array_types..., :Float64, :Int)
        @eval *(B::$mt,A::$at) = fast_mul(B,A)
        @eval *(A::$at,B::$mt) = fast_mul(A,B)
    end
end

for mt in (:LazyGenerator, :FullGenerator)
    @eval *(A::$mt,B::$mt) = fast_mul(A,B)
end

+(f::SFMDistribution,g::SFMDistribution) =
    throw(DomainError("cannot add SFMDistributions with different DiscretisedFluidQueues"))
function +(f::SFMDistribution{T},g::SFMDistribution{T}) where T<:Mesh
    !(f.dq==g.dq)&&throw(DomainError("SFMDistributions need the same model & mesh"))
    return SFMDistribution{T}(f.coeffs+g.coeffs,f.dq)
end


text:
export MemoryDataSet

"In Memory Data Set"
struct MemoryDataSet{InputType, OutputType} <: AbstractDataSet{InputType, OutputType}
    data::Vector{DataEntry{InputType, OutputType}}

    function MemoryDataSet(set::AbstractDataSet{I, O}) where {I, O}
        (set isa MemoryDataSet) && return set
        newSet = MemoryDataSet{I, O}(length(set))
        for i in eachindex(set)
            newSet[i] = set[i]
        end
        newSet
    end

    MemoryDataSet{I, O}(length = 0) where {I, O} = new{I, O}(Vector{DataEntry{I, O}}(undef, length))
    Base.append!(set::MemoryDataSet{I, O}, set2::Vector{DataEntry{I, O}}) where {I, O} = append!(set.data, set2)
    Base.push!(set::MemoryDataSet{I, O}, entry::DataEntry{I, O}) where {I, O} = push!(set, entry)
    Base.getindex(set::MemoryDataSet, row) = set.data[row]
    Base.setindex!(set::MemoryDataSet{I, O}, value::DataEntry{I, O}, row) where {I, O} = set.data[row] = value
    Base.deleteat!(set::MemoryDataSet, idxes) = deleteat!(set, idxes)
    Base.length(set::MemoryDataSet)::Int64 = length(set.data)
end


text:
import Base: start, done, next
using LibArchive

struct TickRawRecordType
  Date::UInt32
  Ask::UInt32
  Bid::UInt32
  AskVolume::Float32
  BidVolume::Float32
end

struct TickRecordType
  Date::DateTime
  Ask::AbstractFloat
  Bid::AbstractFloat
  AskVolume::AbstractFloat
  BidVolume::AbstractFloat
end

function convert(raw_rec::TickRawRecordType, dt_chunk, p_digits)
    TickRecordType(
        Base.Dates.Millisecond(raw_rec.Date) + dt_chunk,
        Float64(raw_rec.Ask) / 10^p_digits,
        Float64(raw_rec.Bid) / 10^p_digits,
        Float64(raw_rec.AskVolume),
        Float64(raw_rec.BidVolume),
    )
end

struct TickIter{S<:IO}
    stream::S
    ondone::Function
    dt_chunk::DateTime
    p_digits::UInt
end
function TickIter(stream::S, year, month, day, hour, symb; ondone=()->nothing) where S<:IO
    TickIter{S}(stream, ondone, DateTime(year, month, day, hour), price_digits(symb))
end

function start(itr::TickIter)
    seek(itr.stream, 0)
    nothing
end

function next(itr::TickIter, nada)
    fh = itr.stream
    raw_rec = TickRawRecordType(
        ntoh(read(fh, UInt32)),
        ntoh(read(fh, UInt32)),
        ntoh(read(fh, UInt32)),
        ntoh(read(fh, Float32)),
        ntoh(read(fh, Float32))
    )
    rec = convert(raw_rec, itr.dt_chunk, itr.p_digits)
    rec, nothing
end

function done(itr::TickIter, nada)
    if !eof(itr.stream)
        return false
    end
    itr.ondone()
    true
end

struct TickReader
    filename::AbstractString
    reader::LibArchive.Reader
    itr::TickIter

    function TickReader(dt::DateTime, ticker, filename)
        reader = LibArchive.Reader(filename)

        LibArchive.support_format_raw(reader)
        LibArchive.support_filter_all(reader)
        entry = LibArchive.next_header(reader)
        arr = read(reader)
        close(reader)
        stream = IOBuffer(arr)
        itr = TickIter(stream, Dates.year(dt), Dates.month(dt), Dates.day(dt), Dates.hour(dt), ticker)
        new(filename, reader, itr)
    end
end

function to_arrays(reader::TickReader)
    itr = reader.itr
    seek(itr.stream, 0)
    #arr = collect(itr)  # MethodError: no method matching length(...) with Julia 0.6 (and probably 0.5)
    #a_date = Array{DateTime}(map(rec->rec.Date, arr))
    #a_ask = Array{Float64}(map(rec->rec.Ask, arr))
    #a_bid = Array{Float64}(map(rec->rec.Bid, arr))
    #a_ask_vol = Array{Float64}(map(rec->rec.AskVolume, arr))
    #a_bid_vol = Array{Float64}(map(rec->rec.BidVolume, arr))
    a_date = Vector{DateTime}()
    a_ask = Vector{Float64}()
    a_bid = Vector{Float64}()
    a_ask_vol = Vector{Float64}()
    a_bid_vol = Vector{Float64}()
    state = start(itr);
    while !done(itr, state)
        rec, state = next(itr, state)
        push!(a_date, rec.Date)
        push!(a_ask, rec.Ask)
        push!(a_bid, rec.Bid)
        push!(a_ask_vol, rec.AskVolume)
        push!(a_bid_vol, rec.BidVolume)
    end
    a_date, a_ask, a_bid, a_ask_vol, a_bid_vol
end

using DataFrames
function to_dataframe(reader::TickReader)
    a_date, a_ask, a_bid, a_ask_vol, a_bid_vol = to_arrays(reader)
    columns = [:Date, :Ask, :Bid, :AskVolume, :BidVolume]
    df = DataFrame([a_date a_ask a_bid a_ask_vol a_bid_vol])
    names!(df, columns)
    df
end

using TimeSeries: TimeArray
function to_timearray(reader::TickReader)
    a_date, a_ask, a_bid, a_ask_vol, a_bid_vol = to_arrays(reader)
    columns = ["Ask", "Bid", "AskVolume", "BidVolume"]
    dat = [a_ask a_bid a_ask_vol a_bid_vol]
    TimeArray(a_date, dat, columns)
end


text:
#=
Type and methods to simplify data filtering of tradable assets
=#

using CSV
using DataFrames
using Statistics
using Dates: Date, DateFormat, today, format

curDate = format(today(), DateFormat("yyyymmdd"))

fileName = "$curDate.filtered"
#using HDF5, JLD

#function writer(file::String, sep::Char)
#    d = load("data/computed/data.jld")["data"]
#    open(file, "a") do io
#        cols = join(["coin", "long", "short", "pnl"], sep)
#        println(io, cols)
#        for (i, (k, v)) in enumerate(d)
#            t = join([k, v["long"], v["short"], v["pnl"]], sep)
#            println(io, t)
#        end
#    end
#end

#fileo = "data/computed/pnl.csv"
#writer(fileo, ',')

#df = CSV.File(fileo) |> DataFrame
#max_trades = maximum(df.trades_days_idle)
#mask = df.trades_days_idle ./ max_trades * 100 .> 10
#df[mask, :]
#coin = "AION-ETH.csv"
#df = CSV.File(joinpath(pathCoins, coin)) |> DataFrame
#using Dates
#insertcols!(df, 1, :open_datetime => map(x->(unix2datetime(x/1000)), df.open_time))
#
#
#using CSV
#using DataFrames
#using Statistics
##pathCoins = "/home/bane/ext_devs/hdd/candlestick_retriever/data/"
#pathCoins = "../../../python/candlestick_retriever/data/"
#function filterTradesNum(file::String)
#    open(file, "a") do io
#        cols = join(["coin", "trades_num_total","trades_days_idle","trades_num_mean"], ',')
#        println(io, cols)
#        files = readdir(pathCoins)
#        for f in files
#            df = CSV.File(joinpath(pathCoins, f)) |> DataFrame
#            coin = split(f, '.')[1]
#            trades_num_total = sum(df.number_of_trades)
#            trades_days_idle = sum(df.number_of_trades .== 0)
#            trades_num_mean  = Statistics.mean(df.number_of_trades)
#            row = join([coin, trades_num_total, trades_days_idle, trades_num_mean], ',')
#            println(io, row)
#        end
#    end
#end
#
#fileName = "2011111.filtered"
#
#filterTradesNum(fileName)
#df = CSV.File(fileName)
#max_trades = maximum(df.trades_days_idle)
#println(max_trades)
#mask = df.trades_days_idle ./ max_trades * 100 .> 10
#df[mask, :]
#using HDF5, JLD

max_stars_repo_name: lithiferous/Strathreads

max_issues_repo_name: lithiferous/Strathreads

max_forks_repo_name: lithiferous/Strathreads


text:
# This file is part of the Julia package ModularForms.jl
#
# Copyright (c) 2018-2019: Lynn Engelberts and Alexandru Ghitza.

# (Incomplete) list of eta quotients that are Hecke eigenforms,
# from the table on pages 4852--4854 of Yves Martin, Multiplicative
# eta-quotients, Transactions of the AMS, vol 348, no 12, 1996.

using Nemo

R, q = PowerSeriesRing(ZZ, 20, "q")

ETA_QUOTIENT_EIGENFORM = Dict(
  [[1, 24]] => q - 24*q^2 + 252*q^3 - 1472*q^4 + 4830*q^5 - 6048*q^6 - 16744*q^7 + 84480*q^8 - 113643*q^9 - 115920*q^10 + 534612*q^11 - 370944*q^12 - 577738*q^13 + 401856*q^14 + 1217160*q^15 + 987136*q^16 - 6905934*q^17 + 2727432*q^18 + 10661420*q^19 + O(q^20),
  [[1, 8], [2, 8]] => q - 8*q^2 + 12*q^3 + 64*q^4 - 210*q^5 - 96*q^6 + 1016*q^7 - 512*q^8 - 2043*q^9 + 1680*q^10 + 1092*q^11 + 768*q^12 + 1382*q^13 - 8128*q^14 - 2520*q^15 + 4096*q^16 + 14706*q^17 + 16344*q^18 - 39940*q^19 + O(q^20),
  [[1, 6], [3, 6]] => q - 6*q^2 + 9*q^3 + 4*q^4 + 6*q^5 - 54*q^6 - 40*q^7 + 168*q^8 + 81*q^9 - 36*q^10 - 564*q^11 + 36*q^12 + 638*q^13 + 240*q^14 + 54*q^15 - 1136*q^16 + 882*q^17 - 486*q^18 - 556*q^19 + O(q^20),
  [[2, 12]] => q - 12*q^3 + 54*q^5 - 88*q^7 - 99*q^9 + 540*q^11 - 418*q^13 - 648*q^15 + 594*q^17 + 836*q^19 + O(q^20),
  [[1, 4], [2, 2], [4, 4]] => q - 4*q^2 + 16*q^4 - 14*q^5 - 64*q^8 + 81*q^9 + 56*q^10 - 238*q^13 + 256*q^16 + 322*q^17 - 324*q^18 + O(q^20),
  [[1, 4], [5, 4]] => q - 4*q^2 + 2*q^3 + 8*q^4 - 5*q^5 - 8*q^6 + 6*q^7 - 23*q^9 + 20*q^10 + 32*q^11 + 16*q^12 - 38*q^13 - 24*q^14 - 10*q^15 - 64*q^16 + 26*q^17 + 92*q^18 + 100*q^19 + O(q^20),
  [[1, 2], [2, 2], [3, 2], [6, 2]] => q - 2*q^2 - 3*q^3 + 4*q^4 + 6*q^5 + 6*q^6 - 16*q^7 - 8*q^8 + 9*q^9 - 12*q^10 + 12*q^11 - 12*q^12 + 38*q^13 + 32*q^14 - 18*q^15 + 16*q^16 - 126*q^17 - 18*q^18 + 20*q^19 + O(q^20),
  [[1, 3], [7, 3]] => q - 3*q^2 + 5*q^4 - 7*q^7 - 3*q^8 + 9*q^9 - 6*q^11 + 21*q^14 - 11*q^16 - 27*q^18 + O(q^20),
  [[2, 4], [4, 4]] => q - 4*q^3 - 2*q^5 + 24*q^7 - 11*q^9 - 44*q^11 + 22*q^13 + 8*q^15 + 50*q^17 + 44*q^19 + O(q^20),
  [[1, 2], [2, 1], [4, 1], [8, 2]] => q - 2*q^2 - 2*q^3 + 4*q^4 + 4*q^6 - 8*q^8 - 5*q^9 + 14*q^11 - 8*q^12 + 16*q^16 + 2*q^17 + 10*q^18 - 34*q^19 + O(q^20),
  [[3, 8]] => q - 8*q^4 + 20*q^7 - 70*q^13 + 64*q^16 + 56*q^19 + O(q^20),
  [[1, 2], [11, 2]] => q - 2*q^2 - q^3 + 2*q^4 + q^5 + 2*q^6 - 2*q^7 - 2*q^9 - 2*q^10 + q^11 - 2*q^12 + 4*q^13 + 4*q^14 - q^15 - 4*q^16 - 2*q^17 + 4*q^18 + O(q^20),
  [[2, 3], [6, 3]] => q - 3*q^3 + 2*q^7 + 9*q^9 - 22*q^13 + 26*q^19 + O(q^20),
  [[1, 1], [2, 1], [7, 1], [14, 1]] => q - q^2 - 2*q^3 + q^4 + 2*q^6 + q^7 - q^8 + q^9 - 2*q^12 - 4*q^13 - q^14 + q^16 + 6*q^17 - q^18 + 2*q^19 + O(q^20),
  [[1, 1], [3, 1], [5, 1], [15, 1]] => q - q^2 - q^3 - q^4 + q^5 + q^6 + 3*q^8 + q^9 - q^10 - 4*q^11 + q^12 - 2*q^13 - q^15 - q^16 + 2*q^17 - q^18 + 4*q^19 + O(q^20),
  [[2, -4], [4, 16], [8, -4]] => q + 4*q^3 - 2*q^5 - 24*q^7 - 11*q^9 + 44*q^11 + 22*q^13 - 8*q^15 + 50*q^17 - 44*q^19 + O(q^20),
  [[2, -12], [4, 36], [8, -12]] => q + 12*q^3 + 54*q^5 + 88*q^7 - 99*q^9 - 540*q^11 - 418*q^13 + 648*q^15 + 594*q^17 - 836*q^19 + O(q^20),
  [[4, 6]] => q - 6*q^5 + 9*q^9 + 10*q^13 - 30*q^17 + O(q^20),
  [[2, 2], [10, 2]] => q - 2*q^3 - q^5 + 2*q^7 + q^9 + 2*q^13 + 2*q^15 - 6*q^17 - 4*q^19 + O(q^20),
  [[1, 1], [23, 1]] => q - q^2 - q^3 + q^6 + q^8 - q^13 - q^16 + O(q^20),
  [[2, 1], [4, 1], [6, 1], [12, 1]] => q - q^3 - 2*q^5 + q^9 + 4*q^11 - 2*q^13 + 2*q^15 + 2*q^17 - 4*q^19 + O(q^20),
  [[3, 2], [9, 2]] => q - 2*q^4 - q^7 + 5*q^13 + 4*q^16 - 7*q^19 + O(q^20),
  [[2, -2], [4, 5], [8, 5], [16, -2]] => q + 2*q^3 - 5*q^9 - 14*q^11 + 2*q^17 + 34*q^19 + O(q^20),
  [[4, 2], [8, 2]] => q - 2*q^5 - 3*q^9 + 6*q^13 + 2*q^17 + O(q^20),
  [[6, 4]] => q - 4*q^7 + 2*q^13 + 8*q^19 + O(q^20),
  [[2, 1], [22, 1]] => q - q^3 - q^5 + q^11 + q^15 + O(q^20),
  [[2, -1], [4, 4], [6, -1], [8, -1], [12, 4], [24, -1]] => q + q^3 - 2*q^5 + q^9 - 4*q^11 - 2*q^13 - 2*q^15 + 2*q^17 + 4*q^19 + O(q^20),
  [[2, -3], [4, 9], [6, -3], [8, -3], [12, 9], [24, -3]] => q + 3*q^3 - 2*q^7 + 9*q^9 - 22*q^13 - 26*q^19 + O(q^20),
  [[4, -2], [8, 8], [16, -2]] => q + 2*q^5 - 3*q^9 - 6*q^13 + 2*q^17 + O(q^20),
  [[4, -6], [8, 18], [16, -6]] => q + 6*q^5 + 9*q^9 - 10*q^13 - 30*q^17 + O(q^20),
  [[4, -14], [8, 38], [16, -14]] => q + 14*q^5 + 81*q^9 + 238*q^13 + 322*q^17 + O(q^20))


text:
using DifferentialEquations, ModelingToolkit

###############################################################################
#                               Lead Lag example
# from: http://web.engr.oregonstate.edu/~webbky/ESE499_files/Section%209%20Root-Locus%20Design.pdf
# page 93
###############################################################################
function plant1(;name)
    @variables t y(t) dy(t) ddy(t) dddy(t) x(t)
    D = Differential(t)
    eqs = [
        D(y) ~ dy
        D(dy) ~ ddy
        D(ddy) ~ dddy
        15x ~ dddy+10*ddy+27*dy+18*y
    ]
    ODESystem(eqs, t;name, defaults=[x=>0.0, y=>0.0,
                                     dy=>0.0,ddy=>0.0, dddy=>0.0])
end

function ex_lead_lag()
    @variables t
    @named p = plant1()
    @named lead = lead_lag_factory(k=130.0, zero=4.85, pole=100.0)
    @named lag = lead_lag_factory(k=1.0, zero=0.8, pole=0.08)
    connections = [lead.input ~ 1/(1+exp(-100*(t-1)))-p.y
                   lag.input ~ lead.output
                   p.x ~ lag.output]
    @named connected = ODESystem(connections ; systems=[p, lead, lag])
    sys = structural_simplify(connected)
    prob = ODEProblem(sys,[],(0.0,10.0))
    #plot(solve(prob),vars=[p.y])
end
# end lead lag example

###############################################################################
#                                 PID example
# from: https://ctms.engin.umich.edu/CTMS/index.php?example=Introduction&section=ControlPID
###############################################################################
function plant2(;name)
    @variables t F(t) x(t) dx(t)
    D = Differential(t)
    eqs = [
        D(x) ~ dx
        D(dx) ~ F - 10*dx - 20*x
    ]
    ODESystem(eqs, t;name, defaults=[x=>0.0, dx=>0.0])
end

function ex_pid()
    @variables t
    @named pid1 = PID_factory(kp=350.0,ki=300.0,kd=50.0)
    @named p = plant2()
    connections = [pid1.e ~ 1/(1+exp(-100*(t-1)))-p.x
                   p.F ~ pid1.u]
    @named connected = ODESystem(connections ; systems=[pid1,p])
    sys = structural_simplify(connected)
    prob = ODEProblem(sys,[],(0.0,10.0))
    return prob
end
# end PID example


text:
using PowerModels, Ipopt, Gurobi
using CSV;
using DataFrames;
using JuMP

list_file = readdir("./data/hard_case/Julia_Solvable")
solver = optimizer_with_attributes(Ipopt.Optimizer, "tol" => 1e-6)
case_name = AbstractString[]
case_obj = Float64[]
case_solve_time = Float64[]
case_ter_status = TerminationStatusCode[]
case_dual_status = ResultStatusCode[]
case_primal_status = ResultStatusCode[]
case_obj_lb = Float64[]
case_solution = []

for itema in list_file
    if occursin("case", itema)
        # result = run_opf(data_path * itema, ACPPowerModel, solver)
        result = run_ac_opf("./data/hard_case/Julia_Solvable/" * itema, Ipopt.Optimizer)
        push!(case_name, itema)
        push!(case_solve_time, result["solve_time"])
        push!(case_ter_status, result["termination_status"])
        push!(case_primal_status, result["primal_status"])
        push!(case_dual_status, result["dual_status"])
        push!(case_obj, result["objective"])
        push!(case_obj_lb, result["objective_lb"])
        push!(case_solution, result["solution"])
    end
end

#for Gurobi solver
# # obtain the numerical issue cases
# for iKey in keys(result_ter_cat)
#     if result_ter_cat[iKey] == NUMERICAL_ERROR
#         push!(NumList, iKey)
#     end
# end
# # numerical issues are mostly caused by infeasibility as well

# for iKey in keys(result_ter_cat)
#     if result_ter_cat[iKey] == OPTIMAL
#         push!(OptList, iKey)
#     end
# end

# for iKey in keys(result_ter_cat)
#     if result_ter_cat[iKey] == INFEASIBLE
#         push!(InfeasList, iKey)
#     end
# end

dict_cat = Dict("case_name" => case_name, "solve_time" => case_solve_time, "termination_status" => case_ter_status, "dual_status" => case_dual_status,
    "primal_status" => case_primal_status, "objective" => case_obj, "objective_lb" => case_obj_lb, "solution" => case_solution)

df_result = DataFrame(dict_cat)
CSV.write("result.csv", df_result, bufsize = 2^26)


text:
using Aqua
using Base64
using CompatHelper
using CompatHelper: DepInfo, EntryType
using Dates
using GitForge
using GitForge: GitForge, GitHub, GitLab
using Mocking
using Pkg
using Random
using SHA
using Test
using TimeZones
using TOML
using UUIDs

Mocking.activate()

Aqua.test_all(CompatHelper; ambiguities=false)

@testset "`version =` line in the workflow file" begin
    root_directory = dirname(dirname(@__FILE__))
    project_file = joinpath(root_directory, "Project.toml")
    version = Base.VersionNumber(TOML.parsefile(project_file)["version"])
    major_version = version.major
    @test major_version >= 1
    workflow_dir = joinpath(root_directory, ".github", "workflows")
    workflow_filename = joinpath(workflow_dir, "CompatHelper.yml")
    workflow_filecontents = read(workflow_filename, String)
    @test occursin(Regex("\\sversion = \"$(major_version)\"\n"), workflow_filecontents)
    @test length(findall(r"version[\s]*?=", workflow_filecontents)) == 1
end
include("patches.jl")

@testset "CompatHelper.jl" begin
    include(joinpath("utilities", "ci.jl"))
    include(joinpath("utilities", "git.jl"))
    include(joinpath("utilities", "new_versions.jl"))
    include(joinpath("utilities", "ssh.jl"))
    include(joinpath("utilities", "types.jl"))
    include(joinpath("utilities", "utilities.jl"))
    include(joinpath("utilities", "rate_limiting.jl"))
    include(joinpath("utilities", "timestamps.jl"))

    include("dependencies.jl")
    include("exceptions.jl")
    include("pull_requests.jl")
    include("main.jl")

    COMPATHELPER_RUN_INTEGRATION_TESTS = get(ENV, "COMPATHELPER_RUN_INTEGRATION_TESTS", "")
    if COMPATHELPER_RUN_INTEGRATION_TESTS == "true"
        @testset "Integration Tests" begin
            @info "Running Integration Tests"
            include(joinpath("utilities", "integration_tests.jl"))
            include("integration_tests.jl")
        end
    end
end


text:
using Base.Test

include("grains.jl")

@testset "On squares" begin
    @testset "On square $s" for s = UInt64(1):64
        @test on_square(s) == 2^(s-1)
        @test total_after(s) == 2^s - 1
    end
end

@testset "Invalid values" begin
    @testset "Zero" begin
        @test_throws DomainError on_square(0)
        @test_throws DomainError total_after(0)
    end

    @testset "Negative" begin
        @test_throws DomainError on_square(-1)
        @test_throws DomainError total_after(-1)
    end

    @testset "Greater than 64" begin
        @test_throws DomainError on_square(65)
        @test_throws DomainError total_after(65)
    end
end


text:
@enum RobotCmd::Cint begin
    NO_COMMAND = -1

    # Commands to increase/decrease the value of scalar parameters:
    # velocity, yaw, height, pitch, and roll
    INCREASE_VELOCITY   # Increase forward trot
    DECREASE_VELOCITY   # Decrease forward trot (reverse)
    INCREASE_YAW        # Veer right
    DECREASE_YAW        # Veer left
    INCREASE_HEIGHT     # Increase ride height
    DECREASE_HEIGHT     # Decrease ride height
    PITCH_NOSE_UP       # Nose up, tail down
    PITCH_NOSE_DOWN     # Nose down, tail up
    ROLL_LEFT           # Roll body to the left
    ROLL_RIGHT          # Roll body to the right

    # Commands to change the behavior state
    TOGGLE_ACTIVATION   # DEACTIVATED->REST, REST->DEACTIVATED
    TOGGLE_TROT         # REST->TROT, TROT->REST, HOP->TROT, FINISHHOP->TROT
    CYCLE_HOP           # REST->HOP, HOP->FINISHHOP, FINISHHOP->REST, TROT->HOP

    # Synthesized commands - not currently directly accessible via joystick
    TURN_LEFT           # Turn left while trotting in place
    TURN_RIGHT          # Turn right while trotting in place
end

function execute_robotcmd(s::mjSim, robotcmd::RobotCmd)
    c = s.robot.command # shorthand
    # set_height!, adjust_height!

    # Velocity PgUp / PgDn
    if robotcmd == INCREASE_VELOCITY        c.horizontal_velocity[1] += 0.01
    elseif robotcmd == DECREASE_VELOCITY    c.horizontal_velocity[1] -= 0.01

    # Height Home / End
    elseif robotcmd == INCREASE_HEIGHT      adjust_height!(s.robot, +1)
    elseif robotcmd == DECREASE_HEIGHT      adjust_height!(s.robot, -1)

    # Yaw left / right arrow
    elseif robotcmd == INCREASE_YAW         c.yaw_rate -= 0.02
    elseif robotcmd == DECREASE_YAW         c.yaw_rate += 0.02

    # Pitch up / down arrow
    elseif robotcmd == PITCH_NOSE_UP        set_pitch!(s.robot, -1)
    elseif robotcmd == PITCH_NOSE_DOWN      set_pitch!(s.robot, +1)

    # Roll left (/) / right (+)
    elseif robotcmd == ROLL_LEFT            adjust_roll!(s.robot, -1)
    elseif robotcmd == ROLL_RIGHT           adjust_roll!(s.robot, +1)

    # Toggle activate (-) / trot (+) / hop (Enter)
    elseif robotcmd == TOGGLE_ACTIVATION    toggle_activate(s.robot)
    elseif robotcmd == TOGGLE_TROT          toggle_trot(s.robot)
    elseif robotcmd == CYCLE_HOP            toggle_hop(s.robot)

    # Turn left/right (0/.)
    elseif robotcmd == TURN_LEFT            turn_left(s.robot)
    elseif robotcmd == TURN_RIGHT           turn_right(s.robot)
    end
end


text:
using Formatting
using NCDatasets

function loadData(
    datafiles :: AbstractDict,
    varnames :: AbstractArray{String},
    y :: Int64,
    m :: Int64,
    d :: Int64;
    layers :: Union{UnitRange, Colon} = Colon(),
)

    key = format("{:04d}-{:02d}", y, m)

    full_path = datafiles[key][:fp]

    ds = Dataset(full_path, "r")

    data = Dict()

    for varname in varnames

        v = ds[varname]
        dim = length(size(v))

        if dim == 4
            rng = (:, :, layers, d)
        elseif dim == 3
            rng = (:, :, d)
        else
            throw(ErrorException("Unknown dimension: $(dim) of variable $(varname) in file $(full_path)"))
        end

        data[varname] = nomissing(ds[varname][rng...], NaN)

        if dim == 4
            data[varname] = permutedims(data[varname], [3, 1, 2])
        end
    end
    close(ds)

    return data

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
# ------------ Choice sites ------------ #

@inline function (ctx::ProposeContext)(call::typeof(trace),
                                       addr::T,
                                       d::Distribution{K}) where {T <: Address, K}
    visit!(ctx, addr)
    s = rand(d)
    score = logpdf(d, s)
    add_value!(ctx, addr, score, s)
    return s
end

# ------------ Learnable ------------ #

@inline function (ctx::ProposeContext)(fn::typeof(learnable), addr::Address)
    visit!(ctx, addr)
    has_sub(ctx.params, addr) && return get_sub(ctx.params, addr)
    error("Parameter not provided at address $addr.")
end

# ------------ Black box call sites ------------ #

@inline function (ctx::ProposeContext)(c::typeof(trace),
                                       addr::T,
                                       call::Function,
                                       args...) where T <: Address
    visit!(ctx, addr)
    ps = get_sub(ctx.params, addr)
    ret, submap, w = propose(ps, call, args...)
    set_sub!(ctx.map, addr, submap)
    return ret
end

@inline function (ctx::ProposeContext)(c::typeof(trace),
                                       addr::T,
                                       call::G,
                                       args...) where {G <: GenerativeFunction, T <: Address}
    visit!(ctx, addr)
    ps = get_sub(ctx.params, addr)
    ret, submap, w = propose(ps, call.fn, args...)
    set_sub!(ctx.map, addr, submap)
    return ret
end

# ------------ Convenience ------------ #

function propose(opt::J, params, fn::Function, args...) where J <: CompilationOptions
    ctx = Propose(opt, DynamicMap{Value}(), params)
    ret = ctx(fn, args...)
    return ret, ctx.map, ctx.score
end

function propose(params, fn::Function, args...) where J <: CompilationOptions
    ctx = Propose(DefaultPipeline(), DynamicMap{Value}(), params)
    ret = ctx(fn, args...)
    return ret, ctx.map, ctx.score
end

function propose(fn::Function, args...) where J <: CompilationOptions
    ctx = Propose(DefaultPipeline(), DynamicMap{Value}(), Empty())
    ret = ctx(fn, args...)
    return ret, ctx.map, ctx.score
end

# ------------ Documentation ------------ #

@doc(
"""
```julia
mutable struct ProposeContext{T <: AddressMap, P <: AddressMap} <: ExecutionContext
    tr::T
    score::Float64
    visited::Visitor
    params::P
end
```

`ProposeContext` is used to propose traces for inference algorithms which use custom proposals. `ProposeContext` instances can be passed sets of `AddressMap` to configure the propose with parameters which have been learned by differentiable programming.

Inner constructors:

```julia
ProposeContext(tr::T) where T <: AddressMap = new{T}(tr, 0.0, AddressMap())
```

Outer constructors:

```julia
Propose() = ProposeContext(AddressMap())
```
""", ProposeContext)

@doc(
"""
```julia
ret, g_cl, w = propose(fn::Function, args...)
ret, cs, w = propose(fn::typeof(rand), d::Distribution{K}) where K
```

`propose` provides an API to the `ProposeContext` execution context. You can use this function on any of the matching signatures above - it will return the return value `ret`, a `RecordSite` instance specialized to the call, and the score `w`.
""", propose)


text:
using CI_Test
using Test

@testset "CI_Test.jl" begin
    # Write your own tests here.
    alphabet = "abcdefghijklmnopqrstuvwxyz"
    d = CI_Test.TerribleGreek()
    @test typeof(d) == Dict{Char,Char}
    for c in keys(d)
        @test c in alphabet
    end
end


text:
#
# This file is part of the GenServers.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

id(m::Module) = m

function check(m::Module)
end

"""
```
genserver(m::Module, args...; name=nothing,
    pid=myid(), thrd=false,
    sticky=false, taskref=nothing)
```

Create an actor in `:genserver` mode. It uses the
callbacks in a user provided module `m` when processing
messages.

# Arguments

- `m::Module`: a user provided module in current scope,
- `args...`: arguments to the user provided `init` callback.

# Keyword Arguments
- `name=nothing`: if a `name::Symbol` is provided the server
    is registered and the name is returned,
- `pid=myid()`: worker pid to create the actor on,
- `thrd=false`: thread to create the actor on,
- `sticky=false`: if `true`, the actor is created in
    the same thread,
- `taskref=nothing`: if a `Ref{Task}` variable is
    provided, it gets the created `Task`.
"""
function genserver(m::Module, args...; name=nothing,
                pid=myid(), thrd=false,
                sticky=false, taskref=nothing)
    s = spawn(Bhv(id, m); mode=:genserver, pid=pid,
            thrd=thrd, sticky=sticky, taskref=taskref)
    isnothing(name) || register(name, s)
    init!(s, m.init, args...)
    return isnothing(name) ? s : name
end


text:
sess = Session(); init(sess)
@testset "pcl_impose_Dirichlet_boundary_conditions" begin
using Revise; using AdFem; using PyPlot; using SparseArrays
sess = Session(); init(sess)

    A = sprand(10,10,0.6)
    II, JJ, VV = findnz(A)
    nz = length(VV)
    indices = SparseTensor(A).o.indices

    IND = run(sess, indices).+1

    x = placeholder(rand(nz))
    B = RawSparseTensor(indices, x, 10, 10)
    rhs = rand(10)
    bdnode = [1;4;6]
    B1, rhs1 = impose_Dirichlet_boundary_conditions(B, rhs, bdnode, zeros(3))
    n = length(run(sess, B1.o.values))

    loss = sum(B1.o.values^2)
    g = gradients(loss, x)
    J = pcl_impose_Dirichlet_boundary_conditions(IND, bdnode, n)

    function test_f(x0)
        G = run(sess, g, x=>x0)
        H = pcl_linear_op(J, pcl_square_sum(n))
        return G, H
    end

    err1, err2 = test_hessian(test_f, rand(nz); showfig = true)
    @test all(err1.<err2)
end

@testset "pcl_compute_fem_laplace_matrix1" begin
    mmesh = Mesh(10,10,0.1)
    κ = placeholder(rand(get_ngauss(mmesh)))
    A = compute_fem_laplace_matrix1(κ, mmesh)
    loss = sum((A.o.values - rand(mmesh.elem_ndof^2 * get_ngauss(mmesh)))^2)
    g = gradients(loss, κ)

    function test_f(x0)
        G = run(sess, g, κ=>x0)
        J = pcl_compute_fem_laplace_matrix1(mmesh)
        H = pcl_linear_op(J, pcl_square_sum(mmesh.elem_ndof^2 * get_ngauss(mmesh)))
        return G, H
    end

    err1, err2 = test_hessian(test_f, rand(get_ngauss(mmesh)); showfig = true)
    @test all(err1.<err2)
end


text:
# parse
a0 = parse(Int, "32")
println(typeof(a0))

# trunc
a1 = trunc(Int64, 3.444)
a2 = trunc(Int64, 3.999)
println(typeof(a1))

b1 = floor(Int64, 3.444)
b2 = floor(Int64, 3.999)

c1 = round(Int64, 3.444)
c2 = round(Int64, 3.999)

d1 = ceil(Int64, 3.444)
d2 = ceil(Int64, 3.444)

# Math Functions
e = log10(100)
println(log10(100))
e2 = log(2.718281828459045)

f1 = sind(90)
f2 = sin(pi/2)


text:
# To be executed from CAGD.jl

todisplay = VERSION <= VersionNumber("1.2") ? true : false

using CAGD
using SparseArrays
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
if todisplay  using ViewerGL; GL = ViewerGL  end

include("./CAGD.jl/paper/examples/cellular_complex.jl");
include("./CAGD.jl/paper/examples/CSG_example.jl");
include("./CAGD.jl/paper/examples/cycles_boundaries.jl");
include("./CAGD.jl/paper/examples/planar_algebra_2D.jl");
include("./CAGD.jl/paper/examples/solid_algebra_3D.jl");
include("./CAGD.jl/paper/examples/space_arrangement_2d.jl");
include("./CAGD.jl/paper/examples/space_arrangement.jl");
include("./CAGD.jl/paper/examples/subdivision_of_2_cells.jl");


text:
#---------------------------------------------
#  GP Ensemble - multiple GPs blended together
#---------------------------------------------
export GPEnsemble

type GPEnsemble{N}
    models::Vector{GaussianProcess{N}}
    cntrs::Vector{Vector{Float64}}

    w_gen::Float64
    w_qry::Float64
    sig_n2::Float64
    function GPEnsemble{N}(wg,wq,sn2) where N
        new(GaussianProcess{N}[], Vector{Float64}[], wg, wq, sn2)
    end
end
function Base.show{N}(io::IO, gpm::GPEnsemble{N})
    print(io,string(
        "GPEnsemble{$N}(Nmdls = $(length(gpm.models)), ",
        "w_gen = $(gpm.w_gen)):"
        ))
    for i=1:length(gpm.models)
        @printf(io, "\n\t%3u : ", i)
        _printvec(io, gpm.cntrs[i])
        print(io, " -> ", gpm.models[i])
    end
end

function update!{N}(gpm::GPEnsemble{N}, xn, fn)
    mdl_dists = length(gpm.models) == 0 ? [0.0] :
        [gpm.models[i].K(xn, gpm.cntrs[i]) for i=1:length(gpm.models)]

    max_simil, max_idx = findmax(mdl_dists)
    if max_simil < gpm.w_gen
        #Generate new model
        push!(gpm.cntrs, xn)
        push!(gpm.models, SEGP(SECov{N}(1.0,0.5*ones(N)), gpm.sig_n2, N))
        update!(gpm.models[end],xn,fn)
    else
        #Put into existing model
        lgp = gpm.models[max_idx]
        Xnew = hcat(lgp.Xtrn, xn)
        fnew = [lgp.ftrn; fn]
        popt = parameteropt(lgp, Xnew, fnew)

        gpm.models[max_idx] = SEGP(SECov{N}(popt[1],popt[2:end]),lgp.sig_n2,N)
        train!(gpm.models[max_idx], Xnew, fnew)
        gpm.cntrs[max_idx] = vec(mean(Xnew, 2))
    end
    gpm
end

function train!{N}(gpm::GPEnsemble{N}, X::Matrix, f::Vector)
    for i=1:size(X,2)
        update!(gpm, X[:,i], f[i])
    end
    gpm
end

function query{N}(xq, gpm::GPEnsemble{N})
    @assert length(gpm.models) > 0
    ws = [gpm.models[i].K(xq, gpm.cntrs[i]) for i=1:length(gpm.models)]
    closest_idx = find(d -> d > gpm.w_qry, ws)
    if length(closest_idx) == 0
        closest_idx = [findmax(ws)[2]]
    end

    fmdls = Float64[query(xq, gpm.models[i])[1] for i in closest_idx]

    #Weighted average
    ws_close = ws[closest_idx]
    dot(ws_close, fmdls) / sum(ws_close)
end

function search{N}(ftrgt, tol, gpm::GPEnsemble{N}, startx)
    function optf(x)
        fq = query(x, gpm)
        abs(fq - ftrgt)
    end

    opt_pts = zeros(N, size(startx,2))
    opt_vals = zeros(size(startx,2))
    for i=1:size(startx,2)
        opt_result = Optim.optimize(optf, startx[:,i])
        opt_pts[:,i] = opt_result.minimum
        opt_vals[i]  = opt_result.f_minimum
    end
    opt_pts, opt_vals
end


text:
push!(LOAD_PATH, "/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/main")
using DHC_2DUtils
using Profile
using BenchmarkTools
using FFTW
using HDF5
using Test
using SparseArrays

function fink_filter_bank_3dizer(hash, cz; nz=256)
    # -------- set parameters
    nx = convert(Int64,hash["npix"])
    n2d = size(hash["filt_value"])[1]
    dz = nz/2-1
    J = size(hash["j_value"])[1]
    L = size(hash["theta_value"])[1]

    im_scale = convert(Int8,log2(nz))
    # -------- number of bins in radial direction (size scales)
    K = (im_scale-2)*cz

    # -------- allocate output array of zeros
    J_L_K     = zeros(Int32, n2d*K,3)
    psi_index = zeros(Int32, J+1, L, K)
    k_value   = zeros(Float64, K)
    logr      = zeros(Float64,nz)
    filt_temp = zeros(Float64,nx,nx,nz)
    F_z       = zeros(Float64,nx,nx,nz)
    #F_p       = zeros(Float64,nx,nx,nz)

    filtind = fill(CartesianIndex{3}[], K*n2d)
    filtval = fill(Float64[], K*n2d)

    for z = 1:convert(Int64,dz+1)
        sz = mod(z+dz,nz)-dz-1
        z2 = sz^2
        logr[z] = 0.5*log2(max(1,z2))
    end

    @inbounds for k_ind = 1:K
        k = k_ind/cz
        k_value[k_ind] = k  # store for later
        krad  = im_scale-k-1
        Δk    = abs.(logr.-krad)
        kmask = (Δk .<= 1/cz)
        k_count = count(kmask)
        k_vals = findall(kmask)

        # -------- radial part
        #I have somehow dropped a factor of sqrt(2) which has been reinserted as a 0.5... fix my brain
        @views F_z[:,:,kmask].= reshape(1/sqrt(2).*cos.(Δk[kmask] .* (cz*π/2)),1,1,k_count)
        @inbounds for index = 1:n2d
            p_ind = hash["filt_index"][index]
            p_filt = hash["filt_value"][index]
            #F_p[p_ind,:] .= p_filt
            f_ind    = k_ind + (index-1)*K
            @views filt_tmp = p_filt.*F_z[p_ind,kmask]
            #temp_ind = findall((F_p .> 1E-13) .& (F_z .> 1E-13))
            #@views filt_tmp = F_p[temp_ind].*F_z[temp_ind]
            ind = findall(filt_tmp .> 1E-13)
            @views filtind[f_ind] = map(x->CartesianIndex(p_ind[x[1]][1],p_ind[x[1]][2],k_vals[x[2]]),ind)
            #filtind[f_ind] = temp_ind[ind]
            @views filtval[f_ind] = filt_tmp[ind]
            @views J_L_K[f_ind,:] = [hash["J_L"][index,1],hash["J_L"][index,2],k_ind-1]
            psi_index[hash["J_L"][index,1]+1,hash["J_L"][index,2]+1,k_ind] = f_ind
        end
    end

    # -------- metadata dictionary
    info=Dict()
    info["2d_npix"]         = hash["npix"]
    info["2d_j_value"]      = hash["j_value"]
    info["2d_theta_value"]  = hash["theta_value"]
    info["2d_psi_index"]    = hash["psi_index"]
    info["2d_phi_index"]    = hash["phi_index"]
    info["2d_J_L"]          = hash["J_L"]
    info["2d_pc"]           = hash["pc"]
    info["2d_wd"]           = hash["wd"]
    info["2d_fs_center_r"]  = hash["fs_center_r"]
    info["2d_filt_index"]   = hash["filt_index"]
    info["2d_filt_value"]   = hash["filt_value"]

    info["nz"]              = nz
    info["cz"]              = cz
    info["J_L_K"]           = psi_index
    info["psi_index"]       = psi_index
    info["k_value"]         = k_value
    info["filt_index"]      = filtind
    info["filt_value"]      = filtval

    return info
end

filter_hash = fink_filter_hash(1, 8, nx=256, pc=1, wd=2)
filt_3d = fink_filter_bank_3dizer(filter_hash, 1, nz=256)

function DHC_compute_3d(image::Array{Float64,3}, filter_hash)
    # Use 2 threads for FFT
    FFTW.set_num_threads(2)

    # array sizes
    (Nx, Ny, Nz)  = size(image)
    (Nf, )    = size(filter_hash["filt_index"])

    # allocate coeff arrays
    out_coeff = []
    S0  = zeros(Float64, 2)
    S1  = zeros(Float64, Nf)
    S20 = zeros(Float64, Nf, Nf)
    S12 = zeros(Float64, Nf, Nf)
    S2  = zeros(Float64, Nf, Nf)  # traditional 2nd order

    # allocate image arrays for internal use
    im_fdf_0_1 = zeros(Float64,           Nx, Ny, Nz, Nf)   # this must be zeroed!
    im_rd_0_1  = Array{Float64, 4}(undef, Nx, Ny, Nz, Nf)

    ## 0th Order
    S0[1]   = mean(image)
    norm_im = image.-S0[1]
    S0[2]   = sum(norm_im .* norm_im)/(Nx*Ny*Nz)
    norm_im ./= sqrt(Nx*Ny*Nz*S0[2])

    append!(out_coeff,S0[:])

    ## 1st Order
    im_fd_0 = fft(norm_im)  # total power=1.0

    # unpack filter_hash
    f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
    f_val   = filter_hash["filt_value"]

    #f_ind2   = filter_hash2["filt_index"]  # (J, L) array of filters represented as index value pairs
    #f_val2   = filter_hash2["filt_value"]

    zarr = zeros(ComplexF64, Nx, Ny, Nz)  # temporary array to fill with zvals

    # make a FFTW "plan" for an array of the given size and type
    P   = plan_ifft(im_fd_0)   # P is an operator, P*im is ifft(im)

    ## Main 1st Order and Precompute 2nd Order
    for f = 1:Nf
        S1tot = 0.0
        f_i = f_ind[f]  # CartesianIndex list for filter
        f_v = f_val[f]  # Values for f_i
        # for (ind, val) in zip(f_i, f_v)   # this is slower!
        if length(f_i) > 0
            for i = 1:length(f_i)
                ind       = f_i[i]
                zval      = f_v[i] * im_fd_0[ind]
                S1tot    += abs2(zval)
                zarr[ind] = zval
                im_fdf_0_1[ind,f] = abs(zval)
            end
            S1[f] = S1tot/(Nx*Ny*Nz)  # image power ###why do we need to normalize again?
            im_rd_0_1[:,:,:,f] .= abs2.(P*zarr)
            zarr[f_i] .= 0
        end
    end
    append!(out_coeff, S1[:])

    # we stored the abs()^2, so take sqrt (this is faster to do all at once)
    im_rd_0_1 .= sqrt.(im_rd_0_1)

    ## 2nd Order
    #Amat = reshape(im_fdf_0_1, Nx*Ny, Nf)
    #S12  = Amat' * Amat
    #Amat = reshape(im_rd_0_1, Nx*Ny, Nf)
    #S20  = Amat' * Amat

    #append!(out_coeff, S20)
    #append!(out_coeff, S12)

    ## Traditional second order
    for f1 = 1:Nf
        thisim = fft(im_rd_0_1[:,:,:,f1])  # Could try rfft here
        # println("  f1",f1,"  sum(fft):",sum(abs2.(thisim))/Nx^2, "  sum(im): ",sum(abs2.(im_rd_0_1[:,:,f1])))
        # Loop over f2 and do second-order convolution
        for f2 = 1:Nf
            f_i = f_ind[f2]  # CartesianIndex list for filter
            f_v = f_val[f2]  # Values for f_i
            # sum im^2 = sum(|fft|^2/npix)
            S2[f1,f2] = sum(abs2.(f_v .* thisim[f_i]))/(Nx*Ny*Nz)
        end
    end
    append!(out_coeff, S2)

    return out_coeff
end

filter_hash = fink_filter_hash(1, 8, nx=64, pc=1, wd=2)
filt_3d = fink_filter_bank_3dizer(filter_hash, 1, nz=64)

test_img = rand(64,64,64)
DHC_compute_3d(test_img,filt_3d)

@time DHC_compute_3d(test_img,filt_3d)

@benchmark DHC_compute_3d(test_img,filt_3d)

filter_hash = fink_filter_hash(1, 16, nx=64, pc=2, wd=2)
filt_3d = fink_filter_bank_3dizer(filter_hash, 1, nz=64)

test_img = rand(64,64,64)
DHC_compute_3d(test_img,filt_3d)

@benchmark DHC_compute_3d(test_img,filt_3d)

fink_filter_hash(1, 8, nx=64, pc=1, wd=1)

A = rand(32,32)
B = rand(32,32)
transpose(B)*A

filter_hash = fink_filter_hash(1, 16, nx=64, pc=2, wd=2)
@benchmark filt_3d = fink_filter_bank_3dizer(filter_hash, 1, nz=64)
test_img = rand(64,64,64)
@benchmark DHC_compute_3d(test_img,filt_3d)

filter_hash = fink_filter_hash(1, 16, nx=128, pc=2, wd=2)
@benchmark filt_3d = fink_filter_bank_3dizer(filter_hash, 1, nz=128)
test_img = rand(128,128,128)
@benchmark DHC_compute_3d(test_img,filt_3d)

filter_hash = fink_filter_hash(1, 16, nx=64, pc=2, wd=2)
@benchmark filt_3d = fink_filter_bank_3dizer(filter_hash, 1, nz=256)
test_img = rand(64,64,256)
@benchmark DHC_compute_3d(test_img,filt_3d)

filter_hash = fink_filter_hash(1, 16, nx=64, pc=2, wd=2)
@benchmark filt_3d = fink_filter_bank_3dizer(filter_hash, 1, nz=512)
test_img = rand(64,64,512)
@benchmark DHC_compute_3d(test_img,filt_3d)
DHC_compute_3d(test_img,filt_3d)

filt_3d = fink_filter_bank_3dizer(filter_hash, 1, nz=512)

filter_hash = fink_filter_hash(1, 16, nx=64, pc=2, wd=2)
@benchmark filt_3d = fink_filter_bank_3dizer(filter_hash, 1, nz=512)
test_img = rand(64,64,513)
@time out = DHC_compute_3d(test_img,filt_3d)

filter_hash = fink_filter_hash(1, 16, nx=64, pc=2, wd=2)
@benchmark filt_3d = fink_filter_bank_3dizer(filter_hash, 1, nz=128)
test_img = rand(64,64,128)
@benchmark DHC_compute_3d(test_img,filt_3d)
DHC_compute_3d(test_img,filt_3d)

filter_hash = fink_filter_hash(1, 16, nx=64, pc=2, wd=2)
@benchmark filt_3d = fink_filter_bank_3dizer(filter_hash, 1, nz=256)
test_img = rand(64,64,256)
@benchmark DHC_compute_3d(test_img,filt_3d)
DHC_compute_3d(test_img,filt_3d)

## Let me try thread testing

filter_hash = fink_filter_hash(1, 16, nx=64, pc=2, wd=2)
filt_3d = fink_filter_bank_3dizer(filter_hash, 1, nz=128)
test_img = rand(64,64,128)
@benchmark DHC_compute_3d(test_img,filt_3d,FFTthreads=1)
@benchmark DHC_compute_3d(test_img,filt_3d,FFTthreads=2)
@benchmark DHC_compute_3d(test_img,filt_3d,FFTthreads=3)
@benchmark DHC_compute_3d(test_img,filt_3d,FFTthreads=4)
@benchmark DHC_compute_3d(test_img,filt_3d,FFTthreads=5)
@benchmark DHC_compute_3d(test_img,filt_3d,FFTthreads=6)
@benchmark DHC_compute_3d(test_img,filt_3d,FFTthreads=7)
@benchmark DHC_compute_3d(test_img,filt_3d,FFTthreads=8)

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
# ref to http://www.genome.ucsc.edu/FAQ/FAQformat.html#format1.8

type Interval
    chrom::ASCIIString
    chromstart::Int
    chromend::Int
end

type BED
    interval::Interval
    score::Int
    strand::ASCIIString
    thickstart::Int
    thickend::Int
    itemrgb::ASCIIString
    blockcount::Int
    blocksizes::Array{Int}
    blockstarts::Array{Int}
end


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

name = "JasPer"
version = v"2.0.28"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/jasper-software/jasper/archive/refs/tags/version-$version.tar.gz", "6b4e5f682be0ab1a5acb0eeb6bf41d6ce17a658bb8e2dbda95de40100939cc88")
]

# Bash recipe for building across all platforms
script = raw"""

cd $WORKSPACE/srcdir/jasper-version-*
mkdir jll_build

cmake -H. \
-B./jll_build \
-DCMAKE_INSTALL_PREFIX=$prefix \
-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
-DJAS_ENABLE_DOC=false \
-DJAS_ENABLE_OPENGL=false \
-DALLOW_IN_SOURCE_BUILD=true \
-DJAS_ENABLE_AUTOMATIC_DEPENDENCIES=false \
-DCMAKE_BUILD_TYPE=Release

cd jll_build/
make clean all
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libjasper", :libjasper),
    ExecutableProduct("imginfo", :imginfo),
    ExecutableProduct("jasper", :jasper),
    ExecutableProduct("imgcmp", :imgcmp)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="JpegTurbo_jll", uuid="aacddb02-875f-59d6-b918-886e6ef4fbf8"))
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
__precompile__()

module TSAnalysis

    # Libraries
    using Dates, Distributed, Logging;
    using LinearAlgebra, Distributions, StableRNGs, Statistics;

    # Custom dependencies
    const local_path = dirname(@__FILE__);
    include("$local_path/types.jl");
    include("$local_path/methods.jl");
    include("$local_path/kalman.jl");
    include("$local_path/subsampling.jl");
    #include("$local_path/uc_models.jl");

    # Export types
    export IntVector, IntMatrix, IntArray, FloatVector, FloatMatrix, FloatArray, JVector, JMatrix, JArray, DiagMatrix, SymMatrix,
           KalmanSettings, KalmanStatus, OnlineKalmanStatus, DynamicKalmanStatus, SizedKalmanStatus, ARIMASettings, VARIMASettings;

    # Export methods
    export check_bounds, nan_to_missing!, error_info, verb_message, interpolate_series, forward_backwards_rw_interpolation, centred_moving_average,
            soft_thresholding, solve_discrete_lyapunov, isconverged, trimmed_mean, sum_skipmissing, mean_skipmissing, std_skipmissing, is_vector_in_matrix, demean, standardise, lag, companion_form;

    # Export functions
    export kfilter!, kfilter_full_sample, kfilter_full_sample!, kforecast, ksmoother, fmin_uc_models, arima, varima, forecast,
           block_jackknife, optimal_d, artificial_jackknife, moving_block_bootstrap, stationary_block_bootstrap;
end


text:
############################################################
# judiVector # #################################################
############################################################

# Authors: Philipp Witte (pwitte@eos.ubc.ca), Henryk Modzelewski (hmodzelewski@eos.ubc.ca)
# Date: January 2017

export judiVector, judiVector_to_SeisBlock, src_to_SeisBlock
export write_shot_record, get_data, convert_to_array, rebuild_jv

############################################################

# structure for seismic data as an abstract vector
mutable struct judiVector{T, AT} <: judiMultiSourceVector{T}
    nsrc::Integer
    geometry::Geometry
    data::Vector{AT}
end

############################################################

## outer constructors

"""
    judiVector
        geometry::Geometry
        data

Abstract vector for seismic data. This vector-like structure contains the geometry and data for either\\
receiver data (shot records) or source data (wavelets).

Constructors
============

Construct vector from `Geometry` structure and cell array of shot records or wavelets. The `data` keyword\\
can also be a single (non-cell) array, in which case the data is the same for all source positions:

    judiVector(geometry, data)

Construct vector for observed data from `SegyIO.SeisBlock`. `segy_depth_key` is the `SegyIO` keyword \\
that contains the receiver depth coordinate:

    judiVector(SegyIO.SeisBlock; segy_depth_key="RecGroupElevation")

Construct vector for observed data from out-of-core data container of type `SegyIO.SeisCon`:

    judiVector(SegyIO.SeisCon; segy_depth_key="RecGroupElevation")

Examples
========

(1) Construct data vector from `Geometry` structure and a cell array of shot records:

    dobs = judiVector(rec_geometry, shot_records)

(2) Construct data vector for a seismic wavelet (can be either cell arrays of individual\\
wavelets or a single wavelet as an array):

    q = judiVector(src_geometry, wavelet)

(3) Construct data vector from `SegyIO.SeisBlock` object:

    using SegyIO
    seis_block = segy_read("test_file.segy")
    dobs = judiVector(seis_block; segy_depth_key="RecGroupElevation")

(4) Construct out-of-core data vector from `SegyIO.SeisCon` object (for large SEG-Y files):

    using SegyIO
    seis_container = segy_scan("/path/to/data/directory","filenames",["GroupX","GroupY","RecGroupElevation","SourceDepth","dt"])
    dobs = judiVector(seis_container; segy_depth_key="RecGroupElevation")

"""
function judiVector(geometry::Geometry, data::Array{T, N}) where {T, N}
    T == Float32 || (data = tof32(data))
    N < 3 || throw(judiMultiSourceException("Only 1D (trace) and 2D (record) input data supported"))
    nsrc = get_nsrc(geometry)
    dataCell = Vector{Array{T, N}}(undef, nsrc)
    for j=1:nsrc
        dataCell[j] = deepcopy(data)
    end
    return judiVector{T, Array{T, N}}(nsrc, geometry, dataCell)
end

# constructor if data is passed as a cell array
function judiVector(geometry::Geometry, data::Vector{Array{T, N}}) where {T, N}
    T == Float32 || (data = tof32.(data))
    nsrcd = length(data)
    nsrcg = get_nsrc(geometry)
    nsrcd == nsrcg || throw(judiMultiSourceException("Number of sources in geometry and data don't match $(nsrcd) != $(nsrcg)"))
    return judiVector{T, Array{T, N}}(nsrcd, geometry, data)
end

# contructor for in-core data container and given geometry
function judiVector(geometry::Geometry, data::SegyIO.SeisBlock)
    # length of data vector
    src = get_header(data,"FieldRecord")
    nsrc = length(unique(src))
    # fill data vector with pointers to data location
    dataCell = Vector{Array{Float32, 2}}(undef, nsrc)
    for j=1:nsrc
        traces = findall(src .== unique(src)[j])
        dataCell[j] = convert(Array{Float32, 2}, data.data[:,traces])
    end
    return judiVector{Float32, Array{Float32, 2}}(nsrc, geometry, dataCell)
end

# contructor for single out-of-core data container and given geometry
function judiVector(geometry::Geometry, data::SegyIO.SeisCon)
    # length of data vector
    nsrc = length(data)
    # fill data vector with pointers to data location
    dataCell = Vector{SegyIO.SeisCon}(undef, nsrc)
    for j=1:nsrc
        dataCell[j] = split(data,j)
    end
    return judiVector{Float32, SegyIO.SeisCon}(nsrc, geometry,dataCell)
end

judiVector(data::SegyIO.SeisBlock; segy_depth_key="RecGroupElevation") = judiVector(Geometry(data; key="receiver", segy_depth_key=segy_depth_key), data)
judiVector(data::SegyIO.SeisCon; segy_depth_key="RecGroupElevation")= judiVector(Geometry(data; key="receiver", segy_depth_key=segy_depth_key), data)
judiVector(data::Vector{SegyIO.SeisCon}; segy_depth_key="RecGroupElevation") = judiVector(Geometry(data; key="receiver", segy_depth_key=segy_depth_key), data)
judiVector(geometry::Geometry, data::Vector{SegyIO.SeisCon}) =  judiVector{Float32, SegyIO.SeisCon}(length(data), geometry, data)

############################################################
## overloaded multi_source functions
time_sampling(jv::judiVector) = jv.geometry.dt

############################################################
# JOLI conversion
jo_convert(::Type{T}, jv::judiVector{T, Array{T, N}}, ::Bool) where {T<:AbstractFloat, N} = jv
jo_convert(::Type{T}, jv::judiVector{vT, Array{vT, N}}, B::Bool) where {T<:AbstractFloat, vT, N} = judiVector{T, Array{T, N}}(jv.nsrc, jv.geometry, jo_convert.(T, jv.data, B))
zero(::Type{T}, v::judiVector{vT, AT}; nsrc::Integer=v.nsrc) where {T, vT, AT} = judiVector{T, AT}(nsrc, deepcopy(v.geometry), T(0) .* v.data[1:nsrc])
function copy!(jv::judiVector, jv2::judiVector)
    jv.geometry = deepcopy(jv2.geometry)
    jv.data .= jv2.data
    jv
end

copyto!(jv::judiVector, jv2::judiVector) = copy!(jv, jv2)
make_input(jv::judiVector{T, Matrix{T}}) where T = jv.data[1]
make_input(jv::judiVector{T, SeisCon}) where T = convert(Matrix{T}, jv.data[1][1].data)

check_compat(ms::Vararg{judiVector, N}) where N = all(y -> compareGeometry(y.geometry, first(ms).geometry), ms)
##########################################################

# Overload needed base function for SegyIO objects
vec(x::SegyIO.SeisCon) = vec(x[1].data)
dot(x::SegyIO.SeisCon, y::SegyIO.SeisCon) = dot(x[1].data, y[1].data)
norm(x::SegyIO.SeisCon, p::Real=2) = norm(x[1].data, p)
abs(x::SegyIO.IBMFloat32) = abs(Float32(x))
*(n::Number, s::SegyIO.SeisCon) = copy(s)

# push!
function push!(a::judiVector{T, mT}, b::judiVector{T, mT}) where {T, mT}
    typeof(a.geometry) == typeof(b.geometry) || throw(judiMultiSourceException("Geometry type mismatch"))
    append!(a.data, b.data)
    a.nsrc += b.nsrc
    push!(a.geometry, b.geometry)
end

# getindex data container
"""
    getindex(x,source_numbers)

getindex seismic data vectors or matrix-free linear operators and extract the entries that correspond\\
to the shot positions defined by `source_numbers`. Works for inputs of type `judiVector`, `judiModeling`, \\
`judiProjection`, `judiJacobian`, `Geometry`, `judiRHS`, `judiPDE`, `judiPDEfull`.

Examples
========

(1) Extract 2 shots from `judiVector` vector:

    dsub = getindex(dobs,[1,2])

(2) Extract geometry for shot location 100:

    geometry_sub = getindex(dobs.geometry,100)

(3) Extract Jacobian for shots 10 and 20:

    Jsub = getindex(J,[10,20])

"""
function getindex(a::judiVector{avDT, AT}, srcnum::RangeOrVec) where {avDT, AT}
    geometry = getindex(a.geometry, srcnum)     # Geometry of getindexd data container
    return judiVector{avDT, AT}(length(srcnum), geometry, a.data[srcnum])
end

# Create SeisBlock from judiVector container to write to file
function judiVector_to_SeisBlock(d::judiVector{avDT, AT}, q::judiVector{avDT, QT};
                                 source_depth_key="SourceSurfaceElevation",
                                 receiver_depth_key="RecGroupElevation") where {avDT, AT, QT}

    typeof(d.geometry) <: GeometryOOC && (d.geometry = Geometry(d.geometry))
    typeof(q.geometry) <: GeometryOOC && (q.geometry = Geometry(q.geometry))

    blocks = Array{Any}(undef, d.nsrc)
    count = 0
    for j=1:d.nsrc

        # create SeisBlock
        blocks[j] = SeisBlock(d.data[j])
        numTraces = size(d.data[j],2)
        traceNumbers = convert(Array{Integer,1},count+1:count+numTraces)

        # set headers
        set_header!(blocks[j], "GroupX", Int.(round.(d.geometry.xloc[j]*1f3)))
        if length(d.geometry.yloc[j]) == 1
            set_header!(blocks[j], "GroupY", Int(round(d.geometry.yloc[j][1]*1f3)))
        else
            set_header!(blocks[j], "GroupY", convert(Array{Integer,1},round.(d.geometry.yloc[j]*1f3)))
        end
        set_header!(blocks[j], receiver_depth_key, Int.(round.(d.geometry.zloc[j]*1f3)))
        set_header!(blocks[j], "SourceX", Int.(round.(q.geometry.xloc[j][1]*1f3)))
        set_header!(blocks[j], "SourceY", Int.(round.(q.geometry.yloc[j][1]*1f3)))
        set_header!(blocks[j], source_depth_key, Int.(round.(q.geometry.zloc[j][1]*1f3)))

        set_header!(blocks[j], "dt", Int(d.geometry.dt[j]*1f3))
        set_header!(blocks[j], "FieldRecord",j)
        set_header!(blocks[j], "TraceNumWithinLine", traceNumbers)
        set_header!(blocks[j], "TraceNumWithinFile", traceNumbers)
        set_header!(blocks[j], "TraceNumber", traceNumbers)
        set_header!(blocks[j], "ElevationScalar", -1000)
        set_header!(blocks[j], "RecSourceScalar", -1000)
        count += numTraces
    end

    # merge into single block
    fullblock = blocks[1]
    for j=2:d.nsrc
        fullblock = merge(fullblock,blocks[j])
        blocks[j] = []
    end
    return fullblock
end

function src_to_SeisBlock(q::judiVector{avDT, QT};
                                 source_depth_key="SourceSurfaceElevation",
                                 receiver_depth_key="RecGroupElevation") where {avDT, QT}

    return judiVector_to_SeisBlock(q, q;
        source_depth_key=source_depth_key,
        receiver_depth_key=receiver_depth_key)
end

function write_shot_record(srcGeometry::GeometryIC, srcData, recGeometry::GeometryIC, recData, options)
    q = judiVector(srcGeometry, srcData)
    d = judiVector(recGeometry, recData)
    pos = [srcGeometry.xloc[1][1], srcGeometry.yloc[1][1],  srcGeometry.zloc[1][1]]
    pos = join(["_"*string(trunc(p; digits=2)) for p in pos])
    file = join([string(options.file_name), pos,".segy"])
    block_out = judiVector_to_SeisBlock(d, q)
    segy_write(join([options.file_path,"/",file]), block_out)
    container = scan_file(join([options.file_path,"/",file]),
                          ["GroupX", "GroupY", "dt", "SourceSurfaceElevation", "RecGroupElevation"];
                          chunksize=256)
    return container
end

####################################################################################################
# Load OOC
function get_data(x::judiVector{T, SeisCon}) where T
    shots = Array{Array{Float32, 2}, 1}(undef, x.nsrc)
    rec_geometry = Geometry(x.geometry)
    for j=1:x.nsrc
        shots[j] = convert(Array{Float32, 2}, x.data[j][1].data)
    end
    return judiVector(rec_geometry, shots)
end

get_data(x::judiVector{T, Array{Float32, 2}}) where T = x
convert_to_array(x::judiVector) = vcat(vec.(x.data)...)

##### Rebuild bad vector
"""
    reuild_jv(v)
rebuild a judiVector from previous version type or JLD2 reconstructed type
"""
rebuild_jv(v::judiVector{T, AT}) where {T, AT} = v
rebuild_jv(v) = judiVector(convgeom(v), convdata(v))

function rebuild_maybe_jld(x::Vector{Any})
    try
        return tof32(x)
    catch e
        if hasproperty(x[1], :offset)
            return [Float32.(StepRangeLen(xi.ref, xi.step, xi.len, xi.offset)) for xi in x]
        end
        return x
    end
end

# Rebuild for backward compatinility
convgeom(x) = GeometryIC{Float32}([getfield(x.geometry, s) for s=fieldnames(GeometryIC)]...)
convdata(x) = convert(Array{Array{Float32, 2}, 1}, x.data)