text:
"""
    pdf_slope(u, Δ)
    pdf_slope(prim, sw, inK)

Calculate slope of particle distribution function `a = a1 + u * a2 + 0.5 * u^2 * a3`
"""
pdf_slope(u, Δ) = Δ / (u + 1e-7)

function pdf_slope(
    prim::AbstractVector{T1},
    sw::AbstractVector{T2},
    inK,
) where {T1<:Real,T2<:Real}

    sl = similar(sw, axes(prim))

    if length(prim) == 3

        sl[3] =
            4.0 * prim[3]^2 / (inK + 1.0) / prim[1] * (
                2.0 * sw[3] - 2.0 * prim[2] * sw[2] +
                sw[1] * (prim[2]^2 - 0.5 * (inK + 1.0) / prim[3])
            )
        sl[2] = 2.0 * prim[3] / prim[1] * (sw[2] - prim[2] * sw[1]) - prim[2] * sl[3]
        sl[1] =
            sw[1] / prim[1] - prim[2] * sl[2] -
            0.5 * (prim[2]^2 + 0.5 * (inK + 1.0) / prim[3]) * sl[3]

    elseif length(prim) == 4

        sl[4] =
            4.0 * prim[4]^2 / (inK + 2.0) / prim[1] * (
                2.0 * sw[4] - 2.0 * prim[2] * sw[2] - 2.0 * prim[3] * sw[3] +
                sw[1] * (prim[2]^2 + prim[3]^2 - 0.5 * (inK + 2.0) / prim[4])
            )
        sl[3] = 2.0 * prim[4] / prim[1] * (sw[3] - prim[3] * sw[1]) - prim[3] * sl[4]
        sl[2] = 2.0 * prim[4] / prim[1] * (sw[2] - prim[2] * sw[1]) - prim[2] * sl[4]
        sl[1] =
            sw[1] / prim[1] - prim[2] * sl[2] - prim[3] * sl[3] -
            0.5 * (prim[2]^2 + prim[3]^2 + 0.5 * (inK + 2.0) / prim[4]) * sl[4]

    elseif length(prim) == 5

        sl[5] =
            4.0 * prim[5]^2 / (inK + 3.0) / prim[1] * (
                2.0 * sw[5] - 2.0 * prim[2] * sw[2] - 2.0 * prim[3] * sw[3] -
                2.0 * prim[4] * sw[4] +
                sw[1] * (prim[2]^2 + prim[3]^2 + prim[4]^2 - 0.5 * (inK + 3.0) / prim[5])
            )
        sl[4] = 2.0 * prim[5] / prim[1] * (sw[4] - prim[4] * sw[1]) - prim[4] * sl[5]
        sl[3] = 2.0 * prim[5] / prim[1] * (sw[3] - prim[3] * sw[1]) - prim[3] * sl[5]
        sl[2] = 2.0 * prim[5] / prim[1] * (sw[2] - prim[2] * sw[1]) - prim[2] * sl[5]
        sl[1] =
            sw[1] / prim[1] - prim[2] * sl[2] - prim[3] * sl[3] - prim[4] * sl[4] -
            0.5 * (prim[2]^2 + prim[3]^2 + prim[4]^2 + 0.5 * (inK + 3.0) / prim[5]) * sl[5]

    end

    return sl

end

"""
    mixture_pdf_slope(prim, sw, inK)

Calculate slope of multi-component particle distribution function `a = a1 + u * a2 + 0.5 * u^2 * a3`

"""
function mixture_pdf_slope(
    prim::AbstractMatrix{T1},
    sw::AbstractMatrix{T2},
    inK,
) where {T1<:Real,T2<:Real}

    sl = similar(sw, axes(prim))

    for j in axes(sl, 2)
        sl[:, j] .= pdf_slope(prim[:, j], sw[:, j], inK)
    end

    return sl

end

"""
    reduce_distribution(f, weights, dim)
    reduce_distribution(f, v, w, weights, dim)

Reduced distribution function

* @arg : particle distribution function with full velocity space
* @arg : quadrature weights with reduced velocity setting (v & w by default)

"""
function reduce_distribution(
    f::AbstractMatrix{T1},
    weights::AbstractVector{T2},
    dim = 1,
) where {T1<:Real,T2<:Real}

    if dim == 1
        h = similar(f, axes(f, 1))
        for i in eachindex(h)
            h[i] = sum(@. weights * f[i, :])
        end
    elseif dim == 2
        h = similar(f, axes(f, 2))
        for j in eachindex(h)
            h[j] = sum(@. weights * f[:, j])
        end
    else
        throw("dimension dismatch")
    end

    return h

end

function reduce_distribution(
    f::AbstractArray{T1,3},
    weights::AbstractMatrix{T2},
    dim = 1,
) where {T1<:Real,T2<:Real}

    if dim == 1
        h = similar(f, axes(f, 1))
        for i in eachindex(h)
            h[i] = sum(@. weights * f[i, :, :])
        end
    elseif dim == 2
        h = similar(f, axes(f, 2))
        for j in eachindex(h)
            h[j] = sum(@. weights * f[:, j, :])
        end
    elseif dim == 3
        h = similar(f, axes(f, 3))
        for k in eachindex(h)
            h[k] = sum(@. weights * f[:, :, k])
        end
    else
    end

    return h

end

function reduce_distribution(
    f::AbstractArray{X,3},
    v::Y,
    w::Y,
    weights::AbstractMatrix{Z},
    dim = 1,
) where {
    X<:Real,
    Y<:AbstractArray{<:Real,3},
    Z<:Real,
}

    if dim == 1
        h = similar(f, axes(f, 1))
        b = similar(h)
        for i in eachindex(h)
            h[i] = sum(@. weights * f[i, :, :])
            b[i] = sum(@. weights * (v[i, :, :]^2 + w[i, :, :]^2) * f[i, :, :])
        end
    elseif dim == 2
        h = similar(f, axes(f, 2))
        b = similar(h)
        for j in eachindex(h)
            h[j] = sum(@. weights * f[:, j, :])
            b[j] = sum(@. weights * (v[:, j, :]^2 + w[:, j, :]^2) * f[:, j, :])
        end
    elseif dim == 3
        h = similar(f, axes(f, 3))
        b = similar(h)
        for k in eachindex(h)
            h[k] = sum(@. weights * f[:, :, k])
            b[k] = sum(@. weights * (v[:, :, k]^2 + w[:, :, k]^2) * f[:, :, k])
        end
    end

    return h, b

end

"""
    full_distribution(h, b, u, weights, v, w, ρ, γ)
    full_distribution(h, b, u, weights, v, w, prim, γ)

Recover full distribution function from reduced ones

* @arg h & b : reduced particle distribution function with 1D velocity space
* @arg u : quadrature nodes in 1D velocity space
* @arg weights : quadrature weights in 1D velocity space
* @arg v & w : quadrature nodes in the rest velocity space (with 3D setting)
* @return f : particle distribution function with 3D velocity space

"""
function full_distribution(
    h::X,
    b::X,
    u::Y,
    weights::Y,
    v::Z,
    w::Z,
    ρ,
    γ = 5 / 3,
) where {
    X<:AbstractArray{<:AbstractFloat,1},
    Y<:AbstractArray{<:AbstractFloat,1},
    Z<:AbstractArray{<:AbstractFloat,3},
}

    @assert length(h) == size(v, 1) throw(
        DimensionMismatch("reduced and full distribution function mismatch"),
    )

    Ei = 0.5 * discrete_moments(b, u, weights, 0)
    λi = 0.5 * ρ / (γ - 1.0) / Ei / 3.0 * 2.0

    f = similar(v)
    for k in axes(f, 3), j in axes(f, 2), i in axes(f, 1)
        f[i, j, k] = h[i] * (λi / π) * exp(-λi * v[i, j, k]^2) * exp(-λi * w[i, j, k]^2)
    end

    return f

end

full_distribution(
    h::X,
    b::X,
    u::Y,
    weights::Y,
    v::Z,
    w::Z,
    prim::A,
    γ = 5 / 3,
) where {
    X<:AbstractArray{<:AbstractFloat,1},
    Y<:AbstractArray{<:AbstractFloat,1},
    Z<:AbstractArray{<:AbstractFloat,3},
    A<:AbstractArray{<:Real,1},
} = full_distribution(h, b, u, weights, v, w, prim[1], γ)

"""
    shift_pdf!(f, a, du, dt)

Shift distribution function by external force
"""
function shift_pdf!(f::AbstractVector{T}, a, du, dt) where {T<:Real}

    q0 = eachindex(f) |> first # for OffsetArray
    q1 = eachindex(f) |> last

    if a > 0
        shift = Int(floor(a * dt / du)) # only for uniform velocity grid
        for k = q1:-1:q0+shift
            f[k] = f[k-shift]
        end
        for k = q0:shift+q0-1
            f[k] = 0.0
        end

        for k = q0+1:q1
            f[k] += (dt * a - du * shift) * (f[k-1] - f[k]) / du
        end
    else
        shift = Int(floor(-a * dt / du))
        for k = q0:q1-shift
            f[k] = f[k+shift]
        end
        for k = q1-shift+1:q1
            f[k] = 0.0
        end

        for k = q0:q1-1
            f[k] += (dt * a + du * shift) * (f[k] - f[k+1]) / du
        end
    end

    f[q0] = f[q0+1]
    f[q1] = f[q1-1]

    return nothing

end

#--- multi-component gas ---#
function shift_pdf!(
    f::AbstractMatrix{X},
    a::AbstractVector{Y},
    du::AbstractVector{Z},
    dt,
) where {X<:Real,Y<:Real,Z<:Real}

    for j in axes(f, 2)
        _f = @view f[:, j]
        shift_pdf!(_f, a[j], du[j], dt)
    end

    return nothing

end

"""
    chapman_enskog(u, prim, a, A, τ)
    chapman_enskog(u, v, prim, a, b, A, τ)

Recover discrete Chapman-Enskog expansion
"""
function chapman_enskog(
    u::AbstractVector{T1},
    prim::AbstractVector{T2},
    a::AbstractVector{T3},
    A::AbstractVector{T3},
    τ::Real,
) where {T1<:AbstractFloat,T2<:Real,T3<:Real}

    M = maxwellian(u, prim)
    f = @. M * (
        1 -
        τ * (a[1] * u + a[2] * u^2 + 0.5 * a[3] * u^3 + A[1] + A[2] * u + 0.5 * A[3] * u^2)
    )

    return f

end

function chapman_enskog(
    u::AbstractMatrix{T1},
    v::AbstractMatrix{T1},
    prim::AbstractVector{T2},
    a::AbstractVector{T3},
    b::AbstractVector{T3},
    A::AbstractVector{T3},
    τ::Real,
) where {T1<:AbstractFloat,T2<:Real,T3<:Real}

    M = maxwellian(u, v, prim)
    f = @. M * (
        1.0 - τ * (a[1] * u + a[2] * u^2 + a[3] * u * v + 0.5 * a[4] * u * (u^2 + v^2)) -
        τ * (b[1] * v + b[2] * u * v + b[3] * v^2 + 0.5 * b[4] * v * (u^2 + v^2)) -
        τ * (A[1] + A[2] * u + A[3] * v + 0.5 * A[4] * (u^2 + v^2))
    )

    return f

end


text:
using Random
using Distributions
using Knet

Random.seed!(0);

struct LM
  b0
  b1
  sig
end

function f(x, b0, b1)
  return b0 .+ b1 * x
end

function loglike(b0, b1, sig, x, y)
  # return sum(logpdf.(Normal.(f(x, b0, b1), sig), y))
  return -log(sig) - sum(abs2, y - f(x, b0, b1)) / (2 * sig ^2)
end

# DATA
N = 100
x = randn(N)
b0_true = 2.0
b1_true = 3.0
sig_true = .4
y = f(x, b0_true, b1_true) .+ randn(N) * sig_true

# SGD training loop:

function mini(x, y, n)
  idx = rand(1:N, n)
  return (x[idx], y[idx])
end

grad_ll = grad(loglike)

b0 = 0.0
b1 = 0.0
log_sig = 0.0
lr = 1e-2
for i in 1:1000
  g = grad(loglike)(b0, b1, exp(log_sig), x, y)
end

ext: jl


text:
module KeboolaConnectionTransformation
using Pkg
using KeboolaConnectionComponent

function run()::Integer
    config = KeboolaConnectionComponent.config()
    datadir = KeboolaConnectionComponent.datadir()
    parameters = config.parameters

    # prepare the script file
    script = join(get(parameters, "script", []), "\n")
    scriptPath = joinpath(datadir, "script.jl")
    open(scriptPath, "w") do file
        write(file, script)
    end
    print("Script file: " * scriptPath)

    # install packages
    packages = get(parameters, "packages", [])
    if length(packages) > 0
        for package in Iterators.Stateful(packages)
            try
                Pkg.add(package)
            catch e
                throw(ArgumentError("Failed to install package '" * package * "' error: " * e.msg))
            end
        end
        Pkg.resolve()
    end

    # Change current working directory so that relative paths work
    cd(datadir)

    try
        include(scriptPath)
    catch e
        throw(ArgumentError("The transformation script failed: " * sprint(showerror, e)))
    end
    println("Script finished")
    return 0
end

end # module

max_stars_repo_name: keboola/julia-transformation

max_issues_repo_name: keboola/julia-transformation

max_forks_repo_name: keboola/julia-transformation


text:
"""
    GaloisKeys

Stores Galois keys generated by a `KeyGenerator` instance.

See also: [`KeyGenerator`](@ref)
"""
mutable struct GaloisKeys <: SEALObject
  handle::Ptr{Cvoid}

  function GaloisKeys()
    handleref = Ref{Ptr{Cvoid}}(C_NULL)
    # GaloisKeys are created as KSwitchKeys since they share the same data
    retval = ccall((:KSwitchKeys_Create1, libsealc), Clong,
                   (Ref{Ptr{Cvoid}},),
                   handleref)
    @check_return_value retval
    return GaloisKeys(handleref[])
  end

  function GaloisKeys(handle::Ptr{Cvoid})
    object = new(handle)
    finalizer(destroy!, object)
    return object
  end
end

function destroy!(object::GaloisKeys)
  if isallocated(object)
    @check_return_value ccall((:KSwitchKeys_Destroy, libsealc), Clong, (Ptr{Cvoid},), object)
    sethandle!(object, C_NULL)
  end

  return nothing
end

function parms_id(key::GaloisKeys)
  parms_id = zeros(UInt64, 4)
  retval = ccall((:KSwitchKeys_GetParmsId, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt64}),
                 key, parms_id)
  @check_return_value retval
  return parms_id
end


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# Hock-Schittkowski Nonlinear Test Suite  -  HS110
# This file is JuMP implementation of the model described in
#  W. Hock, K. Schittkowski, Test Examples for Nonlinear Programming
#  Codes, Lecture Notes in Economics and Mathematical Systems,
#  Springer, No, 187, 1981
# More information, including original model description, at
# http://www.ai7.uni-bayreuth.de/downloads.htm
#
# This problem has an objective with squared logarithms and a product
# of variables squared.
#############################################################################

using JuMP
using Base.Test

let

m = Model()
@variable(m, -2.001 <= x[1:10] <= 9.999, start = 9)

@NLobjective(m, Min,
    sum{ log(x[j] - 2)^2 + log(10 - x[j])^2, j=1:10} -
    prod{x[i],i=1:10} ^ 0.2
)

solve(m)

@test_approx_eq_eps getobjectivevalue(m) -45.77846971 1e-5

end


text:
import Base: @deprecate

@deprecate by(d::AbstractDataTable, cols, s::Vector{Symbol}) aggregate(d, cols, map(eval, s))
@deprecate by(d::AbstractDataTable, cols, s::Symbol) aggregate(d, cols, eval(s))
@deprecate nullable!(colnames::Array{Symbol,1}, dt::AbstractDataTable) nullable!(dt, colnames)
@deprecate nullable!(colnums::Array{Int,1}, dt::AbstractDataTable) nullable!(dt, colnums)

import Base: keys, values, insert!
@deprecate keys(dt::AbstractDataTable) names(dt)
@deprecate values(dt::AbstractDataTable) DataTables.columns(dt)
@deprecate insert!(dt::DataTable, dt2::AbstractDataTable) merge!(dt, dt2)

@deprecate pool categorical
@deprecate pool! categorical!

@deprecate complete_cases! dropnull!
@deprecate complete_cases completecases

@deprecate sub(dt::AbstractDataTable, rows) view(dt, rows)

@deprecate stackdf stackdt
@deprecate meltdf meltdt


text:
using MortalityTables
using Test
using JSON
using CSV
using Pkg.Artifacts

# This is the path to the Artifacts.toml we will manipulate
MT_artifact_toml = joinpath(pkgdir(MortalityTables), "Artifacts.toml")
tbl_dir = artifact_hash("SOA_Tables", MT_artifact_toml) |> artifact_path

include("CSV.jl")
include("basic.jl")
include("XTbML.jl")
include("parameterized_models.jl")
include("distribution.jl")

# load tables to be used in subsequent tests
tables = MortalityTables.tables()
@test length(tables) > 0

include("projection_scale.jl")
include("MortalityTables.jl")
include("get_SOA_tables.jl")


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

name = "OCCT"
version = v"7.6.2"

# Collection of sources required to build Open CASCADE Technology (OCCT)
sources = [
    GitSource("https://github.com/Open-Cascade-SAS/OCCT.git",
              "bb368e271e24f63078129283148ce83db6b9670a"), # V7_6_2 @ April 29, 2022
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/OCCT
if [[ ${target} == *musl* ]]; then
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/OSD_MemInfo.cxx.patch"
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/OSD_signal.cxx.patch"
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/Standard_StackTrace.cxx.patch"
elif [[ ${target} == *mingw* ]]; then
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/CMakeLists.txt.patch"
elif [[ ${target} == *freebsd* ]]; then
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/Standard_CString.cxx.patch"
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/Standard_StackTrace.cxx.patch"
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/STEPConstruct_AP203Context.cxx.patch"
fi
mkdir build
cd build
cmake -Wno-dev .. \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DBUILD_LIBRARY_TYPE=Shared \
    -DBUILD_MODULE_Draw=0 \
    -DBUILD_MODULE_Visualization=0 \
    -DBUILD_MODULE_ApplicationFramework=0
make -j${nproc}
make install
install_license ../LICENSE_LGPL_21.txt ../OCCT_LGPL_EXCEPTION.txt
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(filter!(p -> arch(p) != "armv6l", supported_platforms()))

# The products that we will ensure are always built
products = [
    LibraryProduct("libTKBO", :libTKBO),
    LibraryProduct("libTKBRep", :libTKBRep),
    LibraryProduct("libTKBin", :libTKBin),
    LibraryProduct("libTKBinL", :libTKBinL),
    LibraryProduct("libTKBinXCAF", :libTKBinXCAF),
    LibraryProduct("libTKBool", :libTKBool),
    LibraryProduct("libTKCAF", :libTKCAF),
    LibraryProduct("libTKCDF", :libTKCDF),
    LibraryProduct("libTKFeat", :libTKFeat),
    LibraryProduct("libTKFillet", :libTKFillet),
    LibraryProduct("libTKG2d", :libTKG2d),
    LibraryProduct("libTKG3d", :libTKG3d),
    LibraryProduct("libTKGeomAlgo", :libTKGeomAlgo),
    LibraryProduct("libTKGeomBase", :libTKGeomBase),
    LibraryProduct("libTKHLR", :libTKHLR),
    LibraryProduct("libTKIGES", :libTKIGES),
    LibraryProduct("libTKLCAF", :libTKLCAF),
    LibraryProduct("libTKMath", :libTKMath),
    LibraryProduct("libTKMesh", :libTKMesh),
    LibraryProduct("libTKOffset", :libTKOffset),
    LibraryProduct("libTKPrim", :libTKPrim),
    LibraryProduct("libTKRWMesh", :libTKRWMesh),
    LibraryProduct("libTKSTEP", :libTKSTEP),
    LibraryProduct("libTKSTEP209", :libTKSTEP209),
    LibraryProduct("libTKSTEPAttr", :libTKSTEPAttr),
    LibraryProduct("libTKSTEPBase", :libTKSTEPBase),
    LibraryProduct("libTKSTL", :libTKSTL),
    LibraryProduct("libTKService", :libTKService),
    LibraryProduct("libTKShHealing", :libTKShHealing),
    LibraryProduct("libTKTopAlgo", :libTKTopAlgo),
    LibraryProduct("libTKV3d", :libTKV3d),
    LibraryProduct("libTKVCAF", :libTKVCAF),
    LibraryProduct("libTKVRML", :libTKVRML),
    LibraryProduct("libTKXCAF", :libTKXCAF),
    LibraryProduct("libTKXDEIGES", :libTKXDEIGES),
    LibraryProduct("libTKXDESTEP", :libTKXDESTEP),
    LibraryProduct("libTKXMesh", :libTKXMesh),
    LibraryProduct("libTKXSBase", :libTKXSBase),
    LibraryProduct("libTKXml", :libTKXml),
    LibraryProduct("libTKXmlL", :libTKXmlL),
    LibraryProduct("libTKXmlXCAF", :libTKXmlXCAF),
    LibraryProduct("libTKernel", :libTKernel)
]

# Some dependencies are needed only on Linux and FreeBSD
x11_platforms = filter(p ->Sys.islinux(p) || Sys.isfreebsd(p), platforms)

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Xorg_xorgproto_jll"; platforms=x11_platforms),
    Dependency("FreeType2_jll"),
    Dependency("Libglvnd_jll"; platforms=x11_platforms),
    Dependency("Xorg_libX11_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXext_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXfixes_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXft_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXinerama_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXrender_jll"; platforms=x11_platforms)
]

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


text:
using MusicXML

exampledir = joinpath(dirname(dirname(@__DIR__)), "examples")

include("$exampledir/parsing.jl")
include("$exampledir/creating.jl")


text:
using StatsBase
using Distributions
using StatPlots,Plots
using LaTeXStrings
include("self_confidence.jl")

# https://github.com/joshday/AverageShiftedHistograms.jl
# might be useful at some point

function hist_diff(h1::Histogram,h2::Histogram)
    @assert length(h1.weights) == length(h2.weights)
    h_diff = h1.weights - h2.weights

    bin_widths = diff(h1.edges[1])
    return h_diff, h1.edges[1][1:end-1]+bin_widths/2, bin_widths[1]
end

function moment_around_reference(val::Array{Float64},loc::Array{Float64},r_star::Float64,
                            s_range::Float64;return_vec::Bool=false)

    r_mom = []
    for i = 1:length(val)
        push!(r_mom,val[i]*loc[i])
    end

    if return_vec
        return r_mom/s_range
    else
        return sum(r_mom/s_range)
    end
end

function residual_moment(c::Array{Float64},t::Array{Float64};
                         r_star::Float64=mean(t),solver_rwd_range::Float64=1.,return_vec::Bool=false,
                        return_hists::Bool=true)

    #  display("c_exp:$(mean(c)), t_exp: $(mean(t))")
    edg_max = maximum([maximum(c), maximum(t)])
    edg_min = minimum([minimum(c), minimum(t)])
    #  display("edg_max: $edg_max, edg_min: $edg_min")
    hist_edges = linspace(edg_min,edg_max,100)
    #  display(hist_edges)

    # based on distribution with widest range, use that as the edges for the second distribution
    c_norm = normalize(fit(Histogram,c,hist_edges,closed=:left),mode=:probability)
    t_norm = normalize(fit(Histogram,t,hist_edges,closed=:left),mode=:probability)
    display(c_norm)
    display("max weight: $(maximum(c_norm.weights))")

    # r for residual
    r_val, r_loc,r_width = hist_diff(c_norm,t_norm)
    r_mom = moment_around_reference(r_val,r_loc,r_star,solver_rwd_range,return_vec=return_vec)
    display("r_mom: $r_mom")
    if return_hists
        return r_val, r_loc, r_width, r_mom, c_norm, t_norm
    else
        return r_val, r_loc, r_width, r_mom
    end
end
function general_logistic(x::Float64;k::Float64=1.,x0::Float64=0.,L::Float64=2.)
    return L/(1+exp(-k*(x-x0)))
end
function general_logistic(x::Array{Float64};k::Float64=1.,x0::Float64=0.,L::Float64=2.)
    return L./(1+exp(-k.*(x-x0)))
end

function X3(c::Array{Float64},t::Array{Float64};
            solver_rwd_range::Float64=1.,scale::Float64=2.,offset::Float64=0.)::Float64
    c = Normal(mean(c),stc(c))
    t = Normal(mean(t),stc(t))
    val = X3(c,t;solver_rwd_range=solver_rwd_range)
    return val
end
function X3(c::Distributions.Distribution,t::Distributions.Distribution;
            solver_rwd_range::Float64=1.,L::Float64=2.,x0::Float64=0.,k::Float64=1.)::Float64
    # using version of 'standardized moment' https://en.wikipedia.org/wiki/Standardized_moment
    # in this version we are using 0. as the absolute mean on which both distributions are being compared
    # also we scale the moment a second time by the total range of known solutions for all solvers
    inv_stdc = 1/c.σ
    inv_stdt = 1/t.σ
    t1 = c.μ*inv_stdc
    t2 = t.μ*inv_stdt

    diff_std_mom = t1-t2 #difference in standard moments about 0

    val = NaN
    if c.μ == 0. && t.μ == 0.
        val = general_logistic(diff_std_mom,k=k,x0=x0,L=L)
    else
        scaled_standardized_moment = diff_std_mom/solver_rwd_range
        val = general_logistic(scaled_standardized_moment,k=k,x0=x0,L=L)
    end
    return val
end

function max_lik(c::Array{Float64},t::Array{Float64};
                         r_star::Float64=mean(t),solver_rwd_range::Float64=1.,return_vec::Bool=false,
                        return_hists::Bool=true)

    #  display("c_exp:$(mean(c)), t_exp: $(mean(t))")
    edg_max = maximum([maximum(c), maximum(t)])
    edg_min = minimum([minimum(c), minimum(t)])
    #  display("edg_max: $edg_max, edg_min: $edg_min")
    hist_edges = linspace(edg_min,edg_max,100)
    #  display(hist_edges)

    # based on distribution with widest range, use that as the edges for the second distribution
    c_norm = normalize(fit(Histogram,c,hist_edges,closed=:left),mode=:probability)
    t_norm = normalize(fit(Histogram,t,hist_edges,closed=:left),mode=:probability)

    ind_max_c = indmax(c_norm.weights)
    c_max = c_norm.weights[ind_max_c]
    c_max_loc = hist_edges[ind_max_c]

    ind_max_t = indmax(t_norm.weights)
    t_max = t_norm.weights[ind_max_t]
    t_max_loc = hist_edges[ind_max_t]

    # r for residual
    r_val, r_loc,r_width = hist_diff(c_norm,t_norm)
    r_mom = moment_around_reference(r_val,r_loc,r_star,solver_rwd_range,return_vec=return_vec)
    display("r_mom: $r_mom")
    if return_hists
        return r_val, r_loc, r_width, r_mom, c_norm, t_norm
    else
        return r_val, r_loc, r_width, r_mom
    end
end

function residual_moment(c::Float64,t::Float64;
                         r_star::Float64=mean(t),solver_rwd_range::Float64=1.,return_vec::Bool=false)

    r_mom = moment_around_reference([1. -1.],[c t],t,solver_rwd_range)
    if return_vec
        return r_mom/solver_rwd_range
    else
        return sum(r_mom/solver_rwd_range)
    end
end

function X3_emp_histogram()
    srand(12345)
    # t for trusted, c for candidate
    t1 = randn(500)*50
    c = t1
    c2 = (randn(1000)+10)*30
    c3 = [randn(1000)-10;randn(1000)+10]*30

    T = [t1,t1,t1]
    C = [c,c2,c3]
    solver_reward_range = [3000.,100.,100.]
    p = []
    c_hists = []
    t_hists = []

    for (t,c,sr) in zip(T,C,solver_reward_range)
        #  NORMALIZE THE VECTORS TO GET A BETTER HISTOGRAM
        r_val, r_loc, r_width, r_mom, c_norm, t_norm = residual_moment(c,t,solver_rwd_range=sr)

        ylimits = (minimum([minimum(r_val),0.]), maximum([maximum(r_val),maximum(t_norm.weights),maximum(c_norm.weights)]))

        #  pyplot()
        b1 = bar(r_loc,r_val,bar_width=r_width,label="residual")
        vline!([mean(t)],lw=3,line=:dot,lc=:green,label=L"$r^*$")
        title!("Residual Moment: $(@sprintf("%0.3f",r_mom))\n sr: $sr")
        yaxis!(ylim=ylimits)

        p1 = plot(t_norm,fillalpha=0.4,label="trusted")
        plot!(c_norm,fillalpha=0.4,label="candidate")
        vline!([mean(t)],lw=3,line=:dot,lc=:green,label=L"$r^*$")
        title!("Trusted vs Candidate Reward Probabilities")
        yaxis!(ylim=ylimits)

        plt = plot(p1,b1,size=(1000,900))
        push!(p,plt)
    end
    plot(p[1],p[2],p[3],layout = (3,1))
end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
module GoldenBlock

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Golden Block" => Ahorn.EntityPlacement(
        Maple.GoldenBlock,
        "rectangle"
    ),
)

frame = "objects/goldblock"
strawberry = "collectables/goldberry/idle00"

Ahorn.minimumSize(entity::Maple.GoldenBlock) = 16, 16
Ahorn.resizable(entity::Maple.GoldenBlock) = true, true

Ahorn.selection(entity::Maple.GoldenBlock) = Ahorn.getEntityRectangle(entity)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.GoldenBlock, room::Maple.Room)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, frame, (i - 1) * 8, 0, 8, 0, 8, 8)
        Ahorn.drawImage(ctx, frame, (i - 1) * 8, height - 8, 8, 16, 8, 8)
    end

    for i in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, 0, (i - 1) * 8, 0, 8, 8, 8)
        Ahorn.drawImage(ctx, frame, width - 8, (i - 1) * 8, 16, 8, 8, 8)
    end

    for i in 2:tilesWidth - 1, j in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, (i - 1) * 8, (j - 1) * 8, 8, 8, 8, 8)
    end

    Ahorn.drawImage(ctx, frame, 0, 0, 0, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, width - 8, 0, 16, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, 0, height - 8, 0, 16, 8, 8)
    Ahorn.drawImage(ctx, frame, width - 8, height - 8, 16, 16, 8, 8)

    Ahorn.drawSprite(ctx, strawberry, width / 2, height / 2)
end

end


text:
@static if VERSION < v"1.6.0-DEV"
    macro load_preference(name, default)
        esc(default)
    end
else
    using Preferences:@load_preference
end

macro pref_log_destruction(expr)
    if @load_preference("LOG_DESTRUCTION", "false") == "true"
        quote
            @ccall jl_safe_printf($(esc(:("Destroying $handle\n")))::Cstring)::Cvoid
            $(esc(expr))
        end
    else
        esc(expr)
    end
end


text:
__precompile__(false)

"""
TODO
"""
module App

using InteractiveUtils
using Distributed
############################
using CodeManager

#names(CodeManager) #shows external links

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

#Base.PROGRAM_FILE

module Game end #init

function getInput()
  print("\x1b[91m>\x1b[0m\x1b[96m ")
  input = readline()
  print("\x1b[0m")
  split(input)
end

const GAME_FILE = joinpath(@__DIR__,"game.jl") #../backup/backup_
const WEBAPP_FILE = joinpath(@__DIR__,"webapp.jl")

load_debug(name,file,app_run,app_clean; clean=true) = print_backtrace(()->load(name,file,app_run,app_clean; clean=clean))

window_is_loaded=false

function load(name,file,app_run,app_clean; clean=true)
  println("=== Load ===")
  loaded = include_module(@__MODULE__, file)[1]

  if loaded
    println("=== Start ===")
    safe_invoke(@__MODULE__, app_run)
    #safe_call(@__MODULE__, app_run)
    if clean
      println("=== CleanUp ===")
      safe_invoke(@__MODULE__, app_clean)
    end
    #safe_call(@__MODULE__, app_clean)
  end

  if clean safe_clean!(@__MODULE__, name) end
end

function run()
  global window_is_loaded
  input = ""
  cmds = Dict([:HELP=>"help"=>"This info",:QUIT=>"quit"=>"Quit program",
  :GAME=>"game"=>"",:APP=>"app"=>"",:CLEAN=>"clean"=>"",:INFO=>"info"=>"",
  :BINDS=>"binds"=>"",:EDIT=>"edit"=>"",:CMD=>"eval"=>"",:RESTART=>"restart"=>"",:EXIT=>"exit"=>""])

  getCMD(cmd::Symbol) = cmds[cmd][1]

  while true
    list = getInput()
    input = list[1]
    value = length(list) > 1 ? string(list[2:end]...) : nothing
    if input == "?" || input == getCMD(:HELP)
        println("[ Command List ]")
        [println("- $cmd") for (_,cmd) in cmds]
    elseif input == "q" || input == getCMD(:QUIT) break
    elseif input == getCMD(:GAME)
      load_debug(:Game, GAME_FILE, :(Game.run), :(Game.cleanUp))
    elseif input == getCMD(:APP)
      if window_is_loaded
        safe_invoke(@__MODULE__, :(WebApp.reload))
      else
        load_debug(:WebApp, WEBAPP_FILE, :(WebApp.run), :(WebApp.cleanUp); clean=false)
        window_is_loaded = true
      end
    elseif input == getCMD(:CLEAN)
        println("<Free Memory>")
        free_memory()
    elseif input == getCMD(:INFO)
        println(varinfo(Game)) #isdefined(@__MODULE__, :Game) ? println(varinfo(Game)) : println("No info for 'Game'")
    elseif input == getCMD(:BINDS)
        println(names(Game,all=true))
    elseif input == getCMD(:EDIT)
        file=joinpath(@__DIR__,value*".jl")
        isfile(file) ? print_backtrace(()->edit(file)) : println("Error: File '"*value*"' does not exist")
    elseif input == getCMD(:CMD)
        result = safe_eval(@__MODULE__, value)
        if result[1]; println(result[2]); end;
    elseif input == getCMD(:RESTART)
        print_backtrace(()->begin this=joinpath(@__DIR__,"../test.bat")
        Base.run(`$this`); exit(); end);
    elseif input == getCMD(:EXIT)
        exit()
    else
        println("Unknown Command")
    end
  end
  println("Close App.")
end #run()

end #App

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
"""
`chem2diff(s::String)`

From chemical equations get differential equations in the form needed for `@diff_eq`

Example:
```
s = "1.0,2.0: 2A -> B"
a = chem2diff(s)
@eval @diff_eq F \$a
```
"""
function chem2diff(s::String)
    #TODO Jacobian
    reactions = split(strip(s),'\n')
    for r in reactions
        @assert sum([i == ':' for i in r]) == 1 "Missing ':' separator for rates in $(r)"
    end
    species = unique(vcat(get_species.(reactions)...))
    diff_eq = ["dy[$(i)] = " for i in 1:length(species) ]
    @info("\n"*join([ "$s <---> dy[$i]"  for (i,s) in enumerate(species) ],'\n'))

    for (i,s) in enumerate(species)
        for r in reactions
            ai, bi = get_stoch(r,s)
            c = bi - ai
            if c != 0.0
                kon, koff = get_rates(r)
                diff_eq[i] *= "$(c)*("
                lhs,rhs = "",""
                for (j,t) in enumerate(species)
                    aj,bj = get_stoch(r,t)
                    if aj == 0.0
                    elseif aj == 1.0
                        lhs *= "*y[$(j)]"
                    else
                        lhs *= "*y[$(j)]^$(aj)"
                    end
                    if bj == 0.0
                    elseif bj == 1.0
                        rhs *= "*y[$(j)]"
                    else
                        rhs *= "*y[$(j)]^$(bj)"
                    end
                end
                if !isempty(lhs)
                    diff_eq[i] *= "$(kon)"*lhs
                else
                    diff_eq[i] *= "0"
                end
                if !isempty(rhs)
                    diff_eq[i] *= "-$(koff)"*rhs
                else
                    diff_eq[i] *= "0"
                end
                diff_eq[i] *= ") + "
            end
        end
        diff_eq[i] = diff_eq[i][1:end-3]*"   #y[$(i)] <--> $(s)"
    end
    a = join("#".*reactions,'\n')*"\n"*join(diff_eq,'\n')
    println(a)
    a = "begin\n"*a*"\nend"
    Meta.parse(a)
end

function get_rates(reaction::AbstractString)
    rates = split(reaction,':')[1]
    a = sum([ i == ',' for i in rates ])
    @assert any(a .== [0,1]) "Incorrect number ($(a)) of rates in $(reaction)"
    if a == 0
        @warn("Reaction $(reaction) has only one rate. Assuming this is the forward rate...")
        return (parse(Float64,rates),0.0)
    else
        rates = split(rates,',')
        return (parse(Float64,rates[1]),parse(Float64,rates[2]))
    end
end

function get_lhs(reaction::AbstractString)
    c = split(reaction,':')[2]
    lhs = strip(split(c,"->")[1])
end

function get_rhs(reaction::AbstractString)
    c = split(reaction,':')[2]
    rhs = strip(split(c,"->")[2])
end

function sep_coefficient(s::AbstractString)
    a = '1'; i = 0; coef = ""
    while isnumeric(a) || a == '.'
        i += 1
        a = s[i]
        coef *= string(a)
    end
    if isempty(coef[1:end-1])
        return 1.0,strip(s)
    end
    coef = coef[1:end-1]
    c = parse(Float64,coef)
    sp = strip(replace(s,coef=>""))
    c,sp
end

function get_species(reaction::AbstractString)
    c = split(reaction,':')[2]
    a = replace(c,"->"=>" + ")
    s = [ String(sep_coefficient(i)[2]) for i in strip.(split(a,'+')) ]
end

function get_stoch(reaction::AbstractString, species::AbstractString)
    a = get_lhs(reaction)
    st_lhs = 0.0
    for i in strip.(split(a,'+'))
        c,sp = sep_coefficient(i)
        if sp == species
            st_lhs += c
        end
    end
    a = get_rhs(reaction)
    st_rhs = 0.0
    for i in strip.(split(a,'+'))
        c,sp = sep_coefficient(i)
        if sp == species
            st_rhs += c
        end
    end
    return st_lhs,st_rhs
end


text:
include("../src/HssMatrices.jl")
using .HssMatrices
using LinearAlgebra
using AbstractTrees
using Plots
using BenchmarkTools

n = 1000
K(x,y) = (x-y) != 0 ? 1/(x-y) : 1000.
A = [ K(x,y) for x=0:1/(n-1):1, y=0:1/(n-1):1];
b = randn(n, 1);

m, n = size(A)
lsz = 64;
rcl = bisection_cluster(1:m, leafsize=lsz)
ccl = bisection_cluster(1:n, leafsize=lsz)
hssA = compress(A, rcl, ccl);

hssA\hssA

n = 20000
K(x,y) = (x-y) != 0 ? 1/(x-y) : 1000.
A = [ K(x,y) for x=0:1/(n-1):1, y=0:1/(n-1):1];
b = randn(n, 1);

# test the simple implementation of cluster trees
m, n = size(A)
lsz = 64;
rcl = bisection_cluster(1:m, leafsize=lsz)
ccl = bisection_cluster(1:n, leafsize=lsz)
hssA = compress(A, rcl, ccl);

@profview hssA\hssA

# # test ULV normally
# x = ulvfactsolve(hssA, b);
# xcor = A\b;
# println(norm(x-xcor)/norm(xcor))

# # test on schewed cluster trees
# lsz = 701;
# rcl = bisection_cluster(1:m, leafsize=lsz)
# ccl = bisection_cluster(1:n, leafsize=lsz)
# rcl.left.left.data = 1:700
# rcl.left.right.data = 701:1001
# #print_tree(rcl)
# hssA = compress(A, rcl, ccl);
# x = ulvfactsolve(hssA, b);
# xcor = A\b;
# println(norm(x-xcor)/norm(xcor))

# # test on schewed cluster trees
# lsz = 701;
# rcl = bisection_cluster(1:m, leafsize=lsz)
# ccl = bisection_cluster(1:n, leafsize=lsz)
# ccl.left.left.data = 1:700
# ccl.left.right.data = 701:1001
# #print_tree(ccl)
# hssA = compress(A, rcl, ccl);
# x = ulvfactsolve(hssA, b);
# xcor = A\b;
# println(norm(x-xcor)/norm(xcor))

plot = plotranks(hssA)


text:
using StochasticSemiDiscretizationMethod

function createSLDOProblem(A,B,ζ,α,β,σ)
    AMx =  ProportionalMX(@SMatrix [0. 1.;-A -2ζ]);
    τ1=2π
    BMx1 = DelayMX(τ1,@SMatrix [0. 0.; B 0.]);
    cVec = Additive(2)
    noiseID = 1
    αMx1 = stCoeffMX(noiseID,ProportionalMX(@SMatrix [0. 0.; α 0.]))
    βMx11 = stCoeffMX(noiseID,DelayMX(τ1,@SMatrix [0. 0.; β 0.]))
    σVec = stAdditive(1,Additive(@SVector [0., σ]))
    LDDEProblem(AMx,[BMx1],[αMx1],[βMx11],cVec,[σVec])
end

SLDOP_lddep=createSLDOProblem(1.,0.1,0.1,0.1,0.1,0.5); # LDDE problem for Hayes equation
method=SemiDiscretization(5,(2π+100eps())/10) # 5th order semi discretization with Δt=2π/10
τmax=2π # the largest τ of the system
# Second Moment mapping
mapping=DiscreteMapping_M2(SLDOP_lddep,method,τmax,n_steps=10,calculate_additive=true); #The discrete mapping of the system

@show spectralRadiusOfMapping(mapping); # spectral radius ρ of the mapping matrix (ρ>1 unstable, ρ<1 stable)
statM2=VecToCovMx(fixPointOfMapping(mapping), length(mapping.M1_Vs[1])); # stationary second moment matrix of the hayes equation (equilibrium position)
@show statM2[1,1] |> sqrt;

# Stability chart for the Delay Oscillator
using MDBM
using Plots
gr();
using LaTeXStrings

method=SemiDiscretization(5,2π/30);
τmax=2π+100eps()
idxs=[1,2,3:2:StochasticSemiDiscretizationMethod.rOfDelay(τmax,method)*2...]

# ζ=0.05, α=0.3*A, β=0.3*B
foo(A,B) = log(spectralRadiusOfMapping(DiscreteMapping_M2(createSLDOProblem(A,B,0.05,0.3*A,0.3*B,0.),method,τmax,idxs,
    n_steps=30),nev=8)); # No additive term calculated

axis=[Axis(-1.0:0.6:5.0,:A),
    Axis(LinRange(-1.5,1.5,12),:B)]

iteration=4;
stab_border_points=getinterpolatedsolution(solve!(MDBM_Problem(foo,axis),iteration));

scatter(stab_border_points...,xlim=(-1.,5.),ylim=(-1.5,1.4),
    label="",title="2nd moment stability border of the Delay Oscillator",xlabel=L"A",ylabel=L"$B$",
    guidefontsize=14,tickfont = font(10),markersize=2,markerstrokewidth=0)


text:
# 3D layer utilities
conv = (stride, in, out) -> Conv((3, 3, 3), in=>out, stride=stride, pad=SamePad())
tran = (stride, in, out) -> ConvTranspose((3, 3, 3), in=>out, stride=stride, pad=SamePad())

conv1 = (in, out) -> Chain(conv(1, in, out), BatchNorm(out), x -> leakyrelu.(x))
conv2 = (in, out) -> Chain(conv(2, in, out), BatchNorm(out), x -> leakyrelu.(x))
tran2 = (in, out) -> Chain(tran(2, in, out), BatchNorm(out), x -> leakyrelu.(x))

# TODO: 2D layer utilities


text:
# Functions for LSRTM
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: December 2016
#

export softThresholding, subsample, Snapshot

mutable struct Snapshot
	path::String
	filename::Array{Any,1}
end

function Snapshot(path::String, filenames...)
	filename_cell = Array{Any}(length(filenames))
	for j=1:length(filenames)
		filename_cell[j] = filenames[j]
	end
	return Snapshot(path, filename_cell)
end

function softThresholding(x, lambda)
# Pointwise soft thresholding function
	return xThresh = sign.(x).*max.(abs.(x)-lambda, 0.f0)
end

max_stars_repo_name: SINBADconsortium/opesciSLIM

max_issues_repo_name: SINBADconsortium/opesciSLIM

max_forks_repo_name: SINBADconsortium/opesciSLIM


text:
import SHA.sha1

function check_sha1(file, hash)
    info("checking SHA1...")
    f = open(file)
    filehash = sha1(f)
    # work around for SHA.sha1() incompatibility
    if isa(filehash, Array{UInt8,1})
        filehash = bytes2hex(filehash)
    end
    if lowercase(filehash) == lowercase(hash)
        info("SHA1 OK")
        return true
    else
        warn("wrong hash, expect $hash, but got $filehash")
        return false
    end
    return false
end

function verify_human_genome_folder(folder, subdir)
    chr_folder = joinpath(folder, subdir)
    for chr in 1:22
        file = joinpath(chr_folder, "chr$chr.fa")
        if !isfile(file)
            return false
        end
    end
    other_chr = ["X", "Y"]
    for chr in other_chr
        file = joinpath(chr_folder, "chr$chr.fa")
        if !isfile(file)
            return false
        end
    end
    return true
end

function download_assembly(assembly)
    if !haskey(human_genomes, assembly)
        error("$assembly is not supported, please use hg19/hg38")
    end

    fileinfo = human_genomes[assembly]

    # check if the assembly already exists
    folder = joinpath(assembly_dir(), assembly)
    if isdir(folder)
        if verify_human_genome_folder(folder, fileinfo["subdir"])
            return true
        end
    end

    # download the reference
    if !isfile(fileinfo["localfile"]) || !check_sha1(fileinfo["localfile"], fileinfo["sha1"])
        info("start to download $assembly from " * fileinfo["source"])
        download(fileinfo["source"], fileinfo["localfile"]*".part")
        mv(fileinfo["localfile"]*".part", fileinfo["localfile"])
        # check sha1
        if !check_sha1(fileinfo["localfile"], fileinfo["sha1"])
            error("Failed to check sha1 integrity for " * fileinfo["localfile"])
        end
    end

    # decompress
    mkpath(folder)
    gz = fileinfo["localfile"]
    info("decompressing...")
    run(`tar -xvzf $gz -C $folder`)

    if verify_human_genome_folder(folder, fileinfo["subdir"])
        return true
    else
        return false
    end
end


text:
using LinearAlgebra

function flag_poles(poles)
    N = length(poles);
    cindex = zeros(Int, N);
    for m=1:N
        if abs(imag(poles[m])) >= eps(Float64)
            if m==1
                cindex[m] = 1;
            else
                if cindex[m-1]==0 || cindex[m-1]==2
                    cindex[m] = 1;
                    cindex[m+1] = 2;
                else
                    cindex[m] = 2;
                end
            end
        end
    end
    return cindex
end
function vectfit3(f, s, poles, weight; relax=true, stable=true, asymp=2,
                  skip_pole=false, skip_res=false)
                  """
                  ===========================================================
                  =   Fast Relaxed Vector Fitting                           =
                  =   Version 1.0                                           =
                  =   Last revised: 08.08.2008                              =
                  =   Written by: Bjorn Gustavsen                           =
                  =   SINTEF Energy Research, N-7465 Trondheim, NORWAY      =
                  =   bjorn.gustavsen@sintef.no                             =
                  =   http://www.energy.sintef.no/Produkt/VECTFIT/index.asp =
                  =   Note: RESTRICTED to NON-COMMERCIAL use                =
                  ===========================================================
              PURPOSE : Approximate f(s) with a state-space model
                      f(s)=C*(s*I-A)^(-1)*B +D +s*E
                        where f(s) is a singe element or a vector of elements.
                        When f(s) is a vector, all elements become fitted with a common
                        pole set.
              INPUT :
              f(s) : function (vector) to be fitted.
                     dimension : (Nc,Ns)
                                  Nc : number of elements in vector
                                  Ns : number of frequency samples
              s : vector of frequency points [rad/sec]
                     dimension : (1,Ns)
              poles : vector of initial poles [rad/sec]
                      dimension : (1,N)
              weight: the rows in the system matrix are weighted using this array. Can be used
                      for achieving higher accuracy at desired frequency samples.
                      If no weighting is desired, use unitary weights: weight=ones(1,Ns).
                      Two dimensions are allowed:
                        dimension : (1,Ns) --> Common weighting for all vector elements.
                        dimension : (Nc,Ns)--> Individual weighting for vector elements.
              relax: optional, default=true
                 true --> Use relaxed nontriviality constraint
                 false --> Use nontriviality constraint of "standard" vector fitting
              stable: optional, default=true
                  false --> unstable poles are kept unchanged
                  true --> unstable poles are made stable by 'flipping' them
                                into the left half-plane
              asymp: optional, default=2
                  1 --> Fitting with D=0,  E=0
                  2 --> Fitting with D!=0, E=0
                  3 --> Fitting with D!=0, E!=0
              skip_pole: optional, default=false
                  true --> The pole identification part is skipped, i.e (C,D,E)
                            are identified using the initial poles (A) as final poles.
              skip_res: optional, default=false
                  true --> The residue identification part is skipped, i.e. only the
                             poles (A) are identified while C,D,E are returned as zero.
              ========
              OUTPUT :
                  fit(s) = C*(s*I-(A)^(-1)*B +D +s.*E
              poles(1,N)    : new poles
              rmserr(1) : root-mean-square error of approximation for f(s).
                                (0 is returned if skip_res==1)
              fit(Nc,Ns): Rational approximation at samples. (0 is returned if
                          skip_res==1).
              APPROACH:
              The identification is done using the pole relocating method known as Vector Fitting [1],
              with relaxed non-triviality constraint for faster convergence and smaller fitting errors [2],
              and utilization of matrix structure for fast solution of the pole identifion step [3].
             ********************************************************************************
              NOTE: The use of this program is limited to NON-COMMERCIAL usage only.
              If the program code (or a modified version) is used in a scientific work,
              then reference should be made to the following:
              [1] B. Gustavsen and A. Semlyen, "Rational approximation of frequency
                  domain responses by Vector Fitting", IEEE Trans. Power Delivery,
                  vol. 14, no. 3, pp. 1052-1061, July 1999.
              [2] B. Gustavsen, "Improving the pole relocating properties of vector
                  fitting", IEEE Trans. Power Delivery, vol. 21, no. 3, pp. 1587-1592,
                  July 2006.
              [3] D. Deschrijver, M. Mrozowski, T. Dhaene, and D. De Zutter,
                  "Macromodeling of Multiport Systems Using a Fast Implementation of
                  the Vector Fitting Method", IEEE Microwave and Wireless Components
                  Letters, vol. 18, no. 6, pp. 383-385, June 2008.
             ********************************************************************************
              This example script is part of the vector fitting package (VFIT3.zip)
              Last revised: 08.08.2008.
              Created by:   Bjorn Gustavsen.
             """
    TOLlow = 1e-18;
    TOLhigh = 1e18;
    a = 0;
    b = 0;
    try
        a, b = size(poles);
    catch
        a = 1;
        b = size(poles)[1];
    end
    if s[1]==0 && a==1
        if poles[1]==0 && poles[2]!=0
            poles[1] = -1;
        elseif poles[2]==0 && poles[1]!=0
            poles[2] = -1;
        elseif poles[1]==0 && poles[2]==0
            poles[1] = -1+10im;
            poles[2] = -1-10im;
        end
    end
    rmserr = [];
    #try
    #    a, b = size(s);
    #catch
    #    a = 1;
    #    b = size(s)[1];
    #end
    #if a < b
        s = transpose(s);
    #end
    #TODO sanity check the arguments
    LAMBD = diagm(0 => poles);
    Ns = length(s);
    N = size(LAMBD)[1]; #length(LAMBD);
    Nc = length(f[:,1]);
    B = ones(N,1);
    SERA = poles;
    SERC = zeros(Nc,N);
    SERD = zeros(Nc,1);
    SERE = zeros(Nc,1);
    roetter = poles;
    fit = zeros(ComplexF64, Nc, Ns);

    weight = transpose(weight);
    if length(weight[1,:])==1
        common_weight = 1;
    elseif length(weight[1,:])==Nc
        common_weight = 0;
    else
        display("ERROR in vectfit3.m: Invalid size of array weight")
        return
    end

    if asymp==1
        offs = 0;
    elseif asymp==2
        offs = 1;
    else
        offs = 2;
    end
    # =======================================================
    ## POLE IDENTIFICATION:
    # =======================================================
    if !(skip_pole)
        Escale = zeros(1,Nc+1);
        # =======================================================
        # Finding out which starting poles are complex :
        # =======================================================
        cindex = flag_poles(poles);
        # =======================================================
        ## Building system - matrix :
        # =======================================================
        if asymp==1 || asymp==2
            Dk = zeros(ComplexF64, Ns,N+1);
            for i=1:Ns
                Dk[i,N+1] = 1;
            end
        elseif asymp==3
            Dk = zeros(ComplexF64, Ns,N+2);
            for i=1:Ns
                Dk[i,N+1] = 1;
                Dk[i,N+2] = s[i];
            end
        else
            Dk = zeros(ComplexF64, Ns,N);
        end
        for m=1:N
            if cindex[m]==0      #real pole
                Dk[:,m] = 1 ./(s .- LAMBD[m,m]);
            elseif cindex[m]==1  #complex pole, 1st part
                Dk[:,m] = 1 ./(s .- LAMBD[m,m]) + 1 ./(s .- conj(LAMBD[m,m]));
                Dk[:,m+1] = 1im ./(s .- LAMBD[m,m]) - 1im ./(s .- conj(LAMBD[m,m]));
            end
        end

        #Scaling for last row of LS-problem (pole identification)
        scale = 0;
        for m=1:Nc
            if length(weight[1,:])==1
                scale = scale + (norm(weight[:,1].*transpose(f[m,:])))^2;
            else
                scale = scale + (norm(weight[:,m].*transpose(f[m,:])))^2;
            end
        end
        scale = sqrt(scale)/Ns;

        if relax
            #Escale=zeros(1,Nc*(N+offs)+N+1);
            #scale=norm(f);#/Ns; #Scaling for sigma in LS problem
            AA = zeros(ComplexF64, Nc*(N+1), N+1);
            bb = zeros(ComplexF64, Nc*(N+1), 1);
            Escale = zeros(ComplexF64, 1, length(AA[1,:]));

            for n=1:Nc
                A = zeros(ComplexF64, Ns, (N+offs)+N+1); #b=zeros(Ns*Nc+1,1);
                if common_weight==1
                    weig = weight;
                else
                    weig = weight[:,n];
                end
                for m=1:N+offs #left block
                    A[1:Ns,m] = weig .* Dk[1:Ns,m];
                end
                inda = N + offs;
                for m=1:N+1 #right block
                    A[1:Ns,inda+m] = -weig .* Dk[1:Ns,m] .* f[n,1:Ns];
                end
                A = [real(A); imag(A); zeros(1, (N+offs)+N+1)];
                #Integral criterion for sigma:
                offset = (N+offs);
                if n==Nc
                    for mm=1:N+1
                        A[2*Ns+1, offset+mm] = real(scale*sum(Dk[:,mm]));
                    end
                end
                F = qr(A);
                Q = Matrix(F.Q);
                R = Matrix(F.R);
                ind1 = N + offs + 1;
                ind2 = N + offs + N + 1;
                R22 = R[ind1:ind2,ind1:ind2];
                AA[(n-1)*(N+1)+1:n*(N+1),:] = R22;
                if n==Nc
                    bb[(n-1)*(N+1)+1:n*(N+1),1] = (Q[end,N+offs+1:end]')*Ns*scale;
                end
            end #for n=1:Nc
            for col=1:length(AA[1,:])
                Escale[col] = 1/norm(AA[:,col]);
                AA[:,col] = Escale[col] .* AA[:,col];
            end
            x = AA\bb;
            #size(x),size(Escale)
            x = x .* transpose(Escale);
        end #if relax

        #Situation: No relaxation, or produced D of sigma extremely small and large. Solve again, without relaxation
        if !relax | (abs(x[end]) < TOLlow) | (abs(x[end]) > TOLhigh)
            AA = zeros(ComplexF64, Nc*(N),N);
            bb = zeros(ComplexF64, Nc*(N),1);
            if !relax
                Dnew = 1;
            else
                if x[end]==0
                    Dnew = 1;
                elseif abs(x[end])<TOLlow
                    Dnew = sign(x[end])*TOLlow;
                elseif abs(x[end])>TOLhigh
                    Dnew = sign(x[end])*TOLhigh;
                end
            end

            for n=1:Nc
                A=zeros(Ns,(N+offs) +N); #b=zeros(Ns*Nc+1,1);
                Escale=zeros(1,N);
                if common_weight==1
                    weig = weight;
                else
                    weig=weight[:,n];
                end

                for m=1:N+offs #left block
                    A[1:Ns,m] = weig.*Dk[1:Ns,m];
                end
                inda = N+offs;
                for m=1:N #right block
                    A[1:Ns,inda+m] = -weig.*Dk[1:Ns,m].*transpose(f[n,1:Ns]);
                end
                b = Dnew*weig.*transpose(f[n,1:Ns]);
                A = [real(A);imag(A)];
                b = [real(b);imag(b)];
                offset = (N+offs);
                Q, R = qr(A,0);
                ind1 = N+offs+1;
                ind2 = N+offs+N;
                R22 = R[ind1:ind2,ind1:ind2];
                AA[(n-1)*N+1:n*N,:] = R22;
                bb[((n-1)*N+1:n*N,1)] = transpose(Q[:,ind1:ind2])*b;
            end #for n=1:Nc
            for col=1:length(AA[1,:])
                Escale[col] = 1 ./ norm(AA[:,col]);
                AA[:,col] = Escale[col] .* AA[:,col];
            end
            #if use_normal==1
            #  x=AA.'*AA\(AA.'*bb);
            #else
            x = AA\bb;
            #end
            x = x .* transpose(Escale);
            x = [x;Dnew];
        end  #if relax==0 | abs(x[end])<TOLlow | abs(x[end])>TOLhigh
        # ************************************
        C = x[1:end-1];
        D = x[end]; #NEW!!
        #We now change back to make C complex :
        # **************
        for m=1:N
            if cindex[m]==1
                for n=1:1 #Nc+1
                    r1 = C[m];
                    r2 = C[m+1];
                    C[m] = r1 + 1im*r2;
                    C[m+1] = r1 - 1im*r2;
                end
            end
        end
        # **************

        # =============================================================================
        # We now calculate the zeros for sigma :
        # =============================================================================
        #oldLAMBD=LAMBD;oldB=B;oldC=C;7
        m = 0;
        for n=1:N
            m = m+1;
            if m<N
                if ( abs(LAMBD[m,m]) > abs(real(LAMBD[m,m])) ) #complex number?
                    LAMBD[m+1,m] = -imag(LAMBD[m,m]);
                    LAMBD[m,m+1] = imag(LAMBD[m,m]);
                    LAMBD[m,m] = real(LAMBD[m,m]);
                    LAMBD[m+1,m+1] = LAMBD[m,m];
                    B[m,1] = 2;
                    B[m+1,1] = 0;
                    koko = C[m];
                    C[m] = real(koko);
                    C[m+1] = imag(koko);
                    m = m+1;
                end
            end
        end
        ZER = LAMBD - B*transpose(C)/D;
        ZER = Array{Float64,2}(ZER); #FIXME that should not be needed, but gives non-conjugate complex pairs...
        roetter = transpose(eigvals(ZER));
        unstables = map(r -> real(r) > 0, roetter);
        if stable
            roetter[unstables] = roetter[unstables] - 2*real(roetter[unstables]); #Forcing unstable poles to be stable...
        end
        roetter = transpose(roetter);
        sort!(roetter, by = x -> (real(x), imag(x)));
        N = length(roetter);
        # =============================================
        #Sorterer polene s.a. de reelle kommer først:
        for n=1:N
            for m=n+1:N
                if (abs(imag(roetter[m])) <= eps(Float64)
                    && abs(imag(roetter[n])) >= eps(Float64))
                    trans = roetter[n];
                    roetter[n] = roetter[m];
                    roetter[m] = trans;
                end
            end
        end
        N1=0;
        for m=1:N
            if imag(roetter[m])==0
                N1 = m;
            end
        end
        if N1<N
            #roetter[N1+1:N] = sort(roetter[N1+1:N]);
            sort!(roetter[N1+1:N], by = x -> (real(x), imag(x)));
        end       # N1: n.o. real poles
        #N2=N-N1; # N2: n.o. imag.poles
        roetter = roetter - 2im*imag(roetter); #10.11.97 !!!
        SERA = transpose(roetter);
    end #if !skip_pole
    # =========================================================================
    #  RESIDUE IDENTIFICATION:
    # =========================================================================
    if !(skip_res)
        # ============================================================================
        # We now calculate SER for f, using the modified zeros of sigma as new poles :
        # ============================================================================
        LAMBD = diagm(0 => roetter);
        #B=ones(N,1);
        # Finding out which poles are complex :
        cindex = flag_poles(roetter);

        # ===============================================================================
        # We now calculate the SER for f (new fitting), using the above calculated
        # zeros as known poles :
        # ===============================================================================
        if asymp==1
            A = zeros(ComplexF64, 2*Ns,N);
            BB = zeros(ComplexF64, 2*Ns,Nc);
        elseif asymp==2
            A = zeros(ComplexF64, 2*Ns,N+1);
            BB = zeros(ComplexF64, 2*Ns,Nc);
        else
            A = zeros(ComplexF64, 2*Ns,N+2);
            BB = zeros(ComplexF64, 2*Ns,Nc);
        end
        for m=1:N
            if cindex[m]==0      #real pole
                Dk[:,m] = 1 ./(s .- LAMBD[m,m]);
            elseif cindex[m]==1  #complex pole, 1st part
                Dk[:,m] = 1 ./(s .- LAMBD[m,m]) + 1 ./(s .- conj(LAMBD[m,m]));
                Dk[:,m+1] = 1im ./(s .- LAMBD[m,m]) - 1im ./(s .- conj(LAMBD[m,m]));
            end
        end
        if common_weight==1
            weight = transpose(weight);
            #I3=diag(ones(1,Nc));I3(:,Nc)=[];
            Dk = zeros(ComplexF64, Ns,N);
            for m=1:N
                if cindex[m]==0      #real pole
                    Dk[:,m] = weight ./(s .- LAMBD[m,m]);
                elseif cindex[m]==1  #complex pole, 1st part
                    Dk[:,m] = weight ./(s .- LAMBD[m,m]) + weight ./(s .- conj(LAMBD[m,m]));
                    Dk[:,m+1] = 1im .*weight ./(s .- LAMBD[m,m]) - 1im .*weight ./(s .- conj(LAMBD[m,m]));
                end
            end
            if asymp==1
                A[1:Ns,1:N] = Dk;
            elseif asymp==2
                A[1:Ns,1:N] = Dk;
                A[1:Ns,N+1] = weight;
            else
                A[1:Ns,1:N] = Dk;
                A[1:Ns,N+1] = weight;
                A[1:Ns,N+2] = weight.*s;
            end
            for m=1:Nc
                BB[1:Ns,m] = weight.*transpose(f[m,:]);
            end
            A[Ns+1:2*Ns,:] = imag(A[1:Ns,:]);
            A[1:Ns,:] = real(A[1:Ns,:]);
            BB[Ns+1:2*Ns,:] = imag(BB[1:Ns,:]);
            BB[1:Ns,:] = real(BB[1:Ns,:]);

            if asymp==2
                A[1:Ns,N+1] = A[1:Ns,N+1];
            elseif asymp==3
                A[1:Ns,N+1] = A[1:Ns,N+1];
                A[Ns+1:2*Ns,N+2] = A[Ns+1:2*Ns,N+2];
            end
            Escale = zeros(ComplexF64, 1,length(A[1,:]));
            for col=1:length(A[1,:]);
                Escale[col] = norm(A[:,col],2);
                A[:,col] = A[:,col]./Escale[col];
            end
            X = A\BB;
            for n=1:Nc
                X[:,n] = X[:,n]./transpose(Escale);
            end
            X = transpose(X);
            C = X[:,1:N];

            if asymp==2
                SERD = X[:,N+1];
            elseif asymp==3
                SERE = X[:,N+2];
                SERD = X[:,N+1];
            end

        else #if common_weight==1
            SERD = zeros(ComplexF64, Nc,1);
            SERE = zeros(ComplexF64, Nc,1);
            C = zeros(ComplexF64, Nc,N);
            for n=1:Nc
                if asymp==1
                    A[1:Ns,1:N] = Dk;
                elseif asymp==2
                    A[1:Ns,1:N] = Dk;
                    A[1:Ns,N+1] = 1;
                else
                    A[1:Ns,1:N] = Dk;
                    A[1:Ns,N+1] = 1;
                    A[1:Ns,N+2] = s;
                end
                for m=1:length(A[1,:])
                    A[1:Ns,m] = weight[:,n].*A[1:Ns,m];
                end
                BB = weight[:,n].*transpose(f(n,:));
                A[Ns+1:2*Ns,:] = imag(A[1:Ns,:]);
                A[1:Ns,:] = real(A[1:Ns,:]);
                BB[Ns+1:2*Ns] = imag(BB[1:Ns]);
                BB[1:Ns] = real(BB[1:Ns]);

                if asymp==2
                    A[1:Ns,N+1] = A[1:Ns,N+1];
                elseif asymp==3
                    A[1:Ns,N+1] = A[1:Ns,N+1];
                    A[Ns+1:2*Ns,N+2] = A[Ns+1:2*Ns,N+2];
                end
                Escale = zeros(1,length(A[1,:]));
                for col=1:length(A[1,:]);
                    Escale[col] = norm(A[:,col],2);
                    A[:,col] = A[:,col]./Escale[col];
                end
                x = A\BB;
                x = x./transpose(Escale);
                C[n,1:N] = transpose(x[1:N]);
                if asymp==2
                    SERD[n] = x[N+1];
                elseif asymp==3
                    SERE[n] = x[N+2];
                    SERD[n] = x[N+1];
                end
            end #for n=1:Nc
        end #if common_weight==1
        # =========================================================================
        #We now change back to make C complex.
        for m=1:N
            if cindex[m]==1
                for n=1:Nc
                    r1 = C[n,m];
                    r2 = C[n,m+1];
                    C[n,m] = r1 + 1im*r2;
                    C[n,m+1] = r1 - 1im*r2;
                end
            end
        end
        # **************
        B = ones(N,1);
        # ====================================================
        SERA = LAMBD;
        SERB = B;
        SERC = C;
        Dk = zeros(ComplexF64, Ns,N);
        for m=1:N
            Dk[:,m] = 1 ./(s .- SERA[m,m]);
        end
        for n=1:Nc
            fit[n,:] = transpose(Dk*SERC[n,:]);
            if asymp==2
                fit[n,:] = fit[n,:] .+ SERD[n];
            elseif asymp==3
                fit[n,:] = fit[n,:] .+ SERD[n] .+ transpose(s).*SERE[n];
            end
        end
        #fit = transpose(fit);
        #f = transpose(f);
        dif1 = fit - f;
        rmserr = sqrt( sum( sum(map(abs, dif1.^2)) ) )/sqrt(Nc*Ns);
        #fit = transpose(fit);
    end # if !skip_res
    A = SERA;
    poles = A;
    if !(skip_res)
        B = SERB;
        C = SERC;
        D = SERD;
        E = SERE;
    else
        B = ones(N,1);
        C = zeros(Nc,N);
        D = zeros(Nc,Nc);
        E = zeros(Nc,Nc);
        rmserr = 0;
    end;
    #A = sparse(diag(A));    % A is complex, make it diagonal
    #SER.A=A; SER.B=B; SER.C=C; SER.D=D; SER.E=E;
    return SERD[1] , diag(poles)[:], SERC[:]
end


text:
# Types and functions for passing around sets of physical values

# Types

"Holds a set of planetary physical values (mass, radius, etc)"
abstract type ValueSet{mc<:ModelComplexity} end

"Holds physical values of mass, radius, pressure, and temperature."
mutable struct PhysicalValues{Ma<:Mass, Ra<:Radius, Pr<:Pressure, Te<:Temperature} <: ValueSet{WithTemp}
    m::Ma
    r::Ra
    P::Pr
    T::Te
end
ValueSet(m, r, P, T) = PhysicalValues(m, r, P, T)

"Holds physical values of mass, radius, and pressure"
mutable struct MassRadiusPressure{Ma<:Mass,Ra<:Radius,Pr<:Pressure} <: ValueSet{NoTemp}
    m::Ma
    r::Ra
    P::Pr
end
ValueSet(m, r, P) = MassRadiusPressure(m, r, P)

# Query and work with sets of physical values

"Number of physical variables (mass, radius, pressure...) used"
function nvars end
nvars(::Type{NoTemp}) = 3  # M, R, P
nvars(::Type{WithTemp}) = 4  # M, R, P, T

"Number of independent physical variables (radius, pressure...) used"
ndeps{mc<:ModelComplexity}(::Type{mc}) = nvars(mc) - 1

# Print them in terms of Earth masses and radii
function Base.show(io::IO, pv::PhysicalValues)
    m, r, P, T = pv.m, pv.r, pv.P, pv.T
    show(io, "$(m/M_earth) M⊕, $(r/R_earth) R⊕, $(P/Pa) Pa, $(T/K) K")
end
function Base.show(io::IO, mrp::MassRadiusPressure)
    m, r, P = mrp.m, mrp.r, mrp.P
    show(io, "$(m/M_earth) M⊕, $(r/R_earth) R⊕, $(P/Pa) Pa")
end

# Get individual values from the ValueSets
mass(vs::ValueSet) = vs.m::Mass
radius(vs::ValueSet) = vs.r::Radius
pressure(vs::ValueSet) = vs.P::Pressure
temperature(vs::PhysicalValues) = vs.T::Temperature
function gravity(vs::ValueSet)
    m = mass(vs)
    r = radius(vs)
    (G * m ./ (r.^2))::Gravity
end
"Get dependent physical values (radius, pressure, [temperature])"
function nonmass end
nonmass(pv::PhysicalValues) = [radius(pv), pressure(pv), temperature(pv)]
nonmass(mrp::MassRadiusPressure) = [radius(mrp), pressure(mrp)]

"Is a given `ValueSet` physical (all positive)?"
function isphysical end
function isphysical(vs::PhysicalValues)
    mass(vs) > 0 && radius(vs) > 0 && pressure(vs) > 0 && temperature(vs) > 0
end
function isphysical(vs::MassRadiusPressure)
    mass(vs) > 0 && radius(vs) > 0 && pressure(vs) > 0
end

# Zero values
Base.zero(::Type{MassRadiusPressure}) = MassRadiusPressure(0, 0, 0)
Base.zero(::Type{PhysicalValues}) = PhysicalValues(0, 0, 0, 0)
Base.zero(::Type{ValueSet{NoTemp}}) = zero(MassRadiusPressure)
Base.zero(::Type{ValueSet{WithTemp}}) = zero(PhysicalValues)

# Set individual values from the ValueSets
setmass!(vs::ValueSet, m) = (vs.m = m)
setradius!(vs::ValueSet, r) = (vs.r = r)
setpressure!(vs::ValueSet, P) = (vs.P = P)
settemperature!(vs::PhysicalValues, T) = (vs.T = T)


text:
Base.:*(x::Union{Num,Symbolic},y::Unitful.AbstractQuantity) = x * y
Base.:/(x::Union{Num,Symbolic},y::Unitful.AbstractQuantity) = x / y

struct ValidationError <: Exception
    message::String
end

function screen_unit(result)
    result isa Unitful.Unitlike || throw(ValidationError("Unit must be a subtype of Unitful.Unitlike, not $(typeof(result))."))
    result isa Unitful.ScalarUnits || throw(ValidationError("Non-scalar units such as $result are not supported. Use a scalar unit instead."))
    result == u"°" && throw(ValidationError("Degrees are not supported. Use radians instead."))
end
"Find the unit of a symbolic item."
get_unit(x::Real) = unitless
function get_unit(x::Unitful.Quantity)
    result = Unitful.unit(x)
    screen_unit(result)
    return result
end
equivalent(x,y) = isequal(1*x,1*y)
unitless = Unitful.unit(1)

get_unit(x::Num) = get_unit(value(x))
function get_unit(x::Symbolic)
    if x isa Sym || operation(x) isa Sym || (operation(x) isa Term && operation(x).f == getindex) || x isa Symbolics.ArrayOp
        if x.metadata !== nothing
            symunits = get(x.metadata, VariableUnit, unitless)
            screen_unit(symunits)
        else
            symunits = unitless
        end
        return symunits
    elseif operation(x) isa Differential
        return get_unit(arguments(x)[1]) / get_unit(operation(x).x)
    elseif  operation(x) isa Difference
        return get_unit(arguments(x)[1]) / get_unit(operation(x).t) #TODO: make this same as Differential
    elseif x isa Pow
        pargs = arguments(x)
        base,expon = get_unit.(pargs)
        @assert expon isa Unitful.DimensionlessUnits
        if base == unitless
            unitless
        else
            pargs[2] isa Number ? operation(x)(base, pargs[2]) : operation(x)(1*base, pargs[2])
        end
    elseif x isa Add # Cannot simply add the units b/c they may differ in magnitude (eg, kg vs g)
        terms = get_unit.(arguments(x))
        firstunit = terms[1]
        for other in terms[2:end]
            termlist = join(map(repr,terms),", ")
            equivalent(other,firstunit) || throw(ValidationError(", in sum $x, units [$termlist] do not match."))
        end
        return firstunit
    elseif operation(x) == Symbolics._mapreduce
        if x.arguments[2] == +
            get_unit(x.arguments[3])
        else
            throw(ValidationError("Unsupported array operation $x"))
        end
    else
        return get_unit(operation(x)(1 .* get_unit.(arguments(x))...))
    end
end

"Get unit of term, returning nothing & showing warning instead of throwing errors."
function safe_get_unit(term, info)
    side = nothing
    try
        side = get_unit(term)
    catch err
        if err isa Unitful.DimensionError
            @warn("$info: $(err.x) and $(err.y) are not dimensionally compatible.")
        elseif err isa ValidationError
            @warn(info*err.message)
        elseif err isa MethodError
            @warn("$info: no method matching $(err.f) for arguments $(typeof.(err.args)).")
        else
            rethrow()
        end
    end
    side
end

function _validate(terms::Vector, labels::Vector{String}; info::String = "")
    equnits = safe_get_unit.(terms, info*" ".*labels)
    allthere = all(map(x -> x!==nothing, equnits))
    allmatching = true
    first_unit = nothing
    if allthere
        for idx in 1:length(equnits)
            if !isequal(terms[idx],0)
                if first_unit === nothing
                    first_unit = equnits[idx]
                elseif !equivalent(first_unit, equnits[idx])
                    allmatching = false
                    @warn("$info: units [$(equnits[1])] for $(labels[1]) and [$(equnits[idx])] for $(labels[idx]) do not match.")

                end
            end
        end
    end
    allthere && allmatching
end

function validate(jump::Union{ModelingToolkit.VariableRateJump, ModelingToolkit.ConstantRateJump}, t::Symbolic; info::String = "")
    newinfo = replace(info,"eq."=>"jump")
    _validate([jump.rate, 1/t], ["rate", "1/t"], info = newinfo) && # Assuming the rate is per time units
    validate(jump.affect!,info = newinfo)
end

function validate(jump::ModelingToolkit.MassActionJump, t::Symbolic; info::String = "")
    left_symbols = [x[1] for x in jump.reactant_stoch] #vector of pairs of symbol,int -> vector symbols
    net_symbols = [x[1] for x in jump.net_stoch]
    all_symbols = vcat(left_symbols,net_symbols)
    allgood = _validate(all_symbols, string.(all_symbols), info = info)
    n = sum(x->x[2],jump.reactant_stoch,init = 0)
    base_unitful = all_symbols[1] #all same, get first
    allgood && _validate([jump.scaled_rates, 1/(t*base_unitful^n)], ["scaled_rates", "1/(t*reactants^$n))"], info = info)
end

function validate(jumps::ArrayPartition{<:Union{Any, Vector{<:JumpType}}}, t::Symbolic)
    labels = ["in Mass Action Jumps,", "in Constant Rate Jumps,", "in Variable Rate Jumps,"]
    all([validate(jumps.x[idx], t, info = labels[idx]) for idx in 1:3])
end

validate(eq::ModelingToolkit.Equation; info::String = "") = _validate([eq.lhs, eq.rhs], ["left", "right"], info = info)
validate(eq::ModelingToolkit.Equation, term::Union{Symbolic,Unitful.Quantity,Num}; info::String = "") = _validate([eq.lhs, eq.rhs, term], ["left","right","noise"], info = info)
validate(eq::ModelingToolkit.Equation, terms::Vector; info::String = "") = _validate(vcat([eq.lhs, eq.rhs], terms), vcat(["left", "right"], "noise  #".*string.(1:length(terms))), info = info)

"Returns true iff units of equations are valid."
validate(eqs::Vector; info::String = "") = all([validate(eqs[idx], info = info*" in eq. #$idx") for idx in 1:length(eqs)])
validate(eqs::Vector, noise::Vector; info::String = "") = all([validate(eqs[idx], noise[idx], info = info*" in eq. #$idx") for idx in 1:length(eqs)])
validate(eqs::Vector, noise::Matrix; info::String = "") = all([validate(eqs[idx], noise[idx, :], info = info*" in eq. #$idx") for idx in 1:length(eqs)])
validate(eqs::Vector, term::Symbolic; info::String = "") = all([validate(eqs[idx], term, info = info*" in eq. #$idx") for idx in 1:length(eqs)])
validate(term::Symbolics.SymbolicUtils.Symbolic) = safe_get_unit(term,"") !== nothing

"Throws error if units of equations are invalid."
check_units(eqs...) = validate(eqs...) || throw(ValidationError("Some equations had invalid units. See warnings for details."))


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

# %% [markdown]
# https://twitter.com/kado_judo0312/status/1412968953450622977

# %%
using Plots
using LsqFit

# %%
model(t, p) = p[1] * exp.(-p[2] * t)
tdata = range(0, 10, length=21)

noise = [0.075383981064656, -0.8126158128658965, 1.5758753705225542, 1.4656391027715079, -0.5883715498465719, -0.680892588706559, 1.2067535098436946, -1.9517689516680197, 0.2532440785082103, 0.7943300540482112, -1.0897646372738818, -2.176470865841504, 0.9492433626731639, 1.3235403906519267, 0.13798106753050526, 1.0267333913533003, 0.03511342272679487, -1.5026374968513874, 0.5459381122556283, 0.6502966880952431, 1.455355259859152];

ydata = model(tdata, [1.0, 2.0]) .* exp.(noise);

# %%
scatter(tdata, ydata)

# %%
scatter(tdata, ydata; yscale=:log10)

# %%
p0 = [0.5, 0.5]
fit = curve_fit(model, tdata, ydata, p0)
param = fit.param
yFit(t) = model(t, param)

# %%
scatter(tdata, ydata)
plot!(tdata , yFit)

# %%
scatter(tdata, ydata; yscale=:log10)
plot!(tdata , yFit)

# %%
logmodel(t, a) = a[1] .- a[2] * t
logfit = curve_fit(logmodel, tdata, log.(ydata), [0.0, 0.0])
@show logfit.param
logyFit(t) = logmodel(t, logfit.param)

# %%
scatter(tdata, ydata; yscale=:log10)
plot!(tdata , exp∘logyFit)

# %%
scatter(tdata, ydata)
plot!(tdata , exp∘logyFit)

# %%

max_stars_repo_path: 0009/LsqFit example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0009/LsqFit example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0009/LsqFit example.jl

max_forks_repo_name: genkuroki/public


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

JLLWrappers.@generate_main_file_header("AlgRemez")
JLLWrappers.@generate_main_file("AlgRemez", UUID("acb6dc63-88f0-54c7-a126-ccdc963b8b3f"))
end  # module AlgRemez_jll


text:
# testing environment: Julia 1.0.1 on Intel i7-4710HQ, 12GB DDR3, Samsung SSD 860 EVO 250GB, GTX 850M

using Revise  # for debugging
push!(LOAD_PATH, "./")  # loading working directory to search paths

import StochInt  # import our library of: stochastic process & integral
# -------------
# DEMO: generate (discretized) daily stock returns with Brownian motion
# NOTE: usually, returns are estimated in percentages rather than digits
# TIME: 1 time: 0.000044 seconds (24 allocations: 9.094 KiB)
@time StkRet = StochInt.simu_BM(
    250,   # 250 trading days
    N = 1,  # periods to discretize in one trading day, here equals one to simulate daily returns
    mu = 0.1,  # drift, estimated by percentage daily returns data
    sigma = 1.0,  # diffusion, estimated by percentage daily returns data
) .+ 2.5;  # add an initial level/returns (%) at time 0
# -------------
# DEMO: generate (discretized) daily stock prices with Geometric Brownian motion
# NOTE: in price simulation, GBM requires returns in digits but not in percentages
# TIME: 1 time: 0.000018 seconds (9 allocations: 2.359 KiB)
@time StkPri = exp.( StkRet ./ 100 );
# -------------
# DEMO: generate (discretized) long-term interest rates (e.g. 10-year fed bond returns) with Ornstein-Uhlenbeck process
# NOTE: interest rates are usually estimated in percecntages rather than in digits
# TIME: 1 time: 0.000046 seconds (24 allocations: 9.094 KiB)
@time IntRat = StochInt.simu_OU(
    250,  # 250 trading days
    N = 1,  # periods to discretize in one trading day
    mu = 0.1,  # drift
    sigma = 1.0,  # diffusion
) .+ 7.1;  # add benchmark interest rate level
# ------------
# # DEMO: solve the steady state distribution of a given Markov Chain (denoted by one-step transition probability matrix)
# # NOTE: element p_{i,j} in one-step matrix P means the probability from state i to state j
# # TIME:
# Pmat = [ 0.25 0.75; 0.2 0.8 ]
# @time SState = StochInt.solve_SS_MarkovChain( Pmat )
# ------------
# DEMO: solve the steady state distribution of a given Markov Chain through a fucking savage method
# NOTE: yes, just do P*P*P... until it converges
# TIME: 1 time (50 * 50 P): 0.0077 seconds (1.27 k allocations: 6.002 MiBs)
Pmat = rand(50,50); for x in 1:50; Pmat[x,:] ./= sum(Pmat[x,:]); end;  # construct sample one-step transition matrix
@time SState = StochInt.solve_SS_MarkovChain_FUCK( Pmat, atol = 1E-6, maxiter = 100 );
# ------------
# DEMO: simulate a Markov chain when given a one-step transition matrix, an initial state and time periods to simulate
# TIME: 1 time (50 * 50 P): 0.00018 seconds (319 allocations: 199.313 KiB)
Pmat = rand(50,50); for x in 1:50; Pmat[x,:] ./= sum(Pmat[x,:]); end;  # construct sample one-step transition matrix
@time RndMC = StochInt.simu_MarkovChain(
    Pmat, # one-step transition matrix, an Array{Float64,2}
    2,  # initial state to start simulation, from 1 ~ 50
    100, # periods to generate
)

import EasyIO  # import our library: easy I/O process
# -----------------
# DEMO: create a time tagging to name files:
# NOTE: it is useful when running the same model for many times
# TIME: 10000 times: 0.030428 seconds (460.00 k allocations: 22.430 MiB, 8.55% gc time)
# NOTE:
@time for x in 1:10000; EasyIO.timeTag(); end;
print(string( "FileSample_", EasyIO.timeTag(), ".txt" ))
# -----------------

#

max_stars_repo_name: Clpr/NumToolLibs

max_issues_repo_name: Clpr/NumToolLibs

max_forks_repo_name: Clpr/NumToolLibs


text:
using Traits
# Check Cmp-trait (comparison) which is implemented in Traits.jl/src/commontraits.jl
@assert istrait(Cmp{Int,Float64})        # Int and Float64 can be compared
@assert istrait(Cmp{Int,String})==false  # Int and String cannot be compared

# make a new trait and add a type to it:
@traitdef MyTr{X,Y} begin
    foobar(X,Y) -> Bool # All type-tuples for which there is a method foo
                        # with that signature belong to MyTr
end
type A
    a::Int
end
@assert istrait(MyTr{A,A})==false  # foobar not implement yet
foobar(a::A, b::A) = a.a==b.a      # implement it
@assert istrait(MyTr{A,A})         # voila!
@assert istrait(MyTr{Int,Int})==false

# make a function which dispatches on traits:
@traitfn ft1{X,Y; Cmp{X,Y}}(x::X,y::Y)  = x>y ? 5 : 6
@traitfn ft1{X,Y; MyTr{X,Y}}(x::X,y::Y) = foobar(x,y) ? -99 : -999

ft1(4,5)        # ==6    i.e. dispatches to first definition
ft1(A(5), A(6)) # ==-999 i.e. dispatches to second definition

try
    ft1("asdf", 5)
catch err
    println(err)
end
foobar(a::String, b::Int) = length(a)==b
ft1("asdf", 5)

## dispatch using subtraits
@traitdef MyTr2{X,Y} <: MyTr{X,Y} begin
    bar(X,Y) -> Bool
end

@traitfn gt1{X,Y; MyTr2{X,Y}}(x::X,y::Y)  = "MyTr2"
@traitfn gt1{X,Y; MyTr{X,Y}}(x::X,y::Y)   = "MyTr"

type B1{T}
    a::T
end
foobar(a::B1, b::B1) = a.a==b.a
type B2{T}
    a::T
end
foobar(a::B2, b::B2) = a.a==b.a
bar(a::B2, b::B2) = a.a==b.a
@assert istrait(MyTr{B1{Int},B1{Int}})  # true
@assert istrait(MyTr2{B1,B1})==false
@assert istrait(MyTr2{B2{Int},B2{Int}}) # need the type parameters
                                        # here for Julia to infer the
                                        # return type

@assert gt1(B1(1), B1(1))=="MyTr"
@assert gt1(B2(1), B2(1))=="MyTr2"


text:
immutable Monetary{C}
    val::Int64

    Monetary(minor::Integer) = new(minor)
    Monetary(major::Integer, minor::Integer) = new(major*10^scale(Monetary{C})+minor)
end

@generated show(io::IO, m::Monetary) =
    let s = scale(m)
        if s == 0
            quote
                print(io, m.val)
            end
        else
            fmt = "%.0$(s)f"
            f = 10^s
            quote
                @printf(io, $fmt, m.val/$f)
            end
        end
    end

currency{C}(::Union{Type{Monetary{C}},Monetary{C}}) = C
scale{C}(::Type{Monetary{C}}) = 2
scale(m::Monetary) = scale(typeof(m))

dollars(m::Monetary) = div(m.val, 10^scale(m))
cents(m::Monetary) = rem(m.val, 10^scale(m))

Base.zero{M<:Monetary}(::Union{Type{M},M}) = M(0)
Base.one{M<:Monetary}(::Union{Type{M},M}) = M(1, 0)

Base.isfinite{M<:Monetary}(::Union{Type{M},M}) = true
Base.eps{M<:Monetary}(::M) = M(1)
Base.typemax{M<:Monetary}(::Union{Type{M},M}) = M(typemax(Int64))
Base.typemin{M<:Monetary}(::Union{Type{M},M}) = M(typemin(Int64))
Base.isless{M<:Monetary}(x::M, y::M) = isless(x.val, y.val)
(==){M<:Monetary}(x::M, y::M) = x.val == y.val

Base.round{M<:Monetary}(::Type{M}, x::AbstractFloat) =
    M(round(Int64, x*10^scale(M)))
Base.round{M<:Monetary}(::Type{M}, x::AbstractFloat, r::RoundingMode) =
    M(round(Int64, x*10^scale(M), r))

Base.convert{M<:Monetary}(::Type{M}, x::Integer) = M(x, 0)
Base.convert{M<:Monetary}(::Type{M}, x::Real) = M(round(Int64, x*10^scale(M)))
Base.convert{T<:Number}(::Type{T}, x::Monetary) = convert(T, x.val/10^scale(x))

(+){M<:Monetary}(m::M) = m
(-){M<:Monetary}(m::M) = M(-m.val)

(+){M<:Monetary}(m1::M, m2::M) = M(m1.val+m2.val)
(-){M<:Monetary}(m1::M, m2::M) = M(m1.val-m2.val)

(*){M<:Monetary}(m::M, k::Integer) = M(k*m.val)
(*){M<:Monetary}(m::M, k::Real) = M(round(Int64, k*m.val))
(*)(k::Number, m::Monetary) = m*k
(/)(m::Monetary, k::Number) = (one(k)/k)*m
(/){M<:Monetary}(m1::M, m2::M) = m1.val/m2.val

const USD = convert(Monetary{:USD}, 1)
const CAD = convert(Monetary{:CAD}, 1)
const GBR = convert(Monetary{:GBR}, 1)
const EUR = convert(Monetary{:EUR}, 1)

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
"""
	Options <: AbstractOptions

Type that wraps a dictionary used for options

### Fields

"""

struct Options <: AbstractOptions
	dict::Dict{Symbol, Any}
	Options(args::Pair{Symbol,Any}...) = new(Dict{Symbol,Any}(args))
	Options(args::Dict{Symbol,<:Any}) = new(args)
end

"""
	keys

"""

keys(opt::Options) = keys(opt.dict)

"""
	values

"""
values(opt::Options) = values(opt.dict)


text:
"""
    prsabrs(RR::Vector{<:Real}, SBP::Vector{<:Real}; L::Int = 15)

Calculate a measure for BRS based on phase-rectified signal averaging as defined by Bauer et al. 2010.

# Keywords

* L: defines the segment length as 2L+1
"""
function prsabrs(RR::Vector{<:Real}, SBP::Vector{<:Real}; L::Int = 15)
    L > 0 || throw(DomainError("L needs to be greater 0, 15 is recommended"))
    # find all anchor points (SBP[i] >= SBP[i-1])
    isAnchor = [false; diff(SBP) .> 0]
    meanΔSBP = mean(filter(x -> x > 0, diff(SBP)))
    anchorInds =  collect(1:length(SBP))[isAnchor]

    # Select segments
    len = length(RR)
    segments = zeros(Float64, 2L+1, length(anchorInds))
    for (i, ind) in enumerate(anchorInds)
        inds = max(ind-L, 1):min(ind+L, len)
        segments[inds .+ (L+1 -ind), i] = RR[max(ind-L, 1):min(ind+L, len)] .- RR[ind]
    end

    meanSegment = vec(mean(segments, dims = 2))
    prsaBRSv = 0.25(meanSegment[L+1] + meanSegment[L+2] - meanSegment[L] - meanSegment[L-1])
    prsaBRSNormv = prsaBRSv / meanΔSBP

    return prsaBRS(value = prsaBRSv, valueNorm = prsaBRSNormv, meanSegment = meanSegment, segments = segments, L = L, anchorInds = anchorInds)
end

"""
Struct that stores all information regarding the prsaBRS estimation. The final result is stored in `value` or `valueNorm` for a normalized value. It can be plotted for visual inspection.
"""
@with_kw mutable struct prsaBRS
    L::Int = 0
    value::Real = 0
    valueNorm::Real = 0
    segments::Matrix{<:Real} = zeros(Float64, 1,1)
    meanSegment::Vector{<:Real} = Float64[]
    anchorInds::Vector{<:Real} = Float64[]
end

@recipe function f(res::prsaBRS)

    layout := (2,1)
    L = res.L

    @series begin
        subplot := 2
        label := ""
        seriescolor := :black

        -L:L, res.meanSegment
    end

    for i in 1:size(res.segments, 2)
        @series begin
            subplot := 1
            label := ""
            seriescolor := :grey
            seriesalpha := 0.2

            -L:L, res.segments[:, i]
        end
    end

    @series begin

        subplot := 1
        yguide := "ΔRR [ms]"
        label := ""
        seriescolor := :black
        linewidth := 2

        -L:L, res.meanSegment
    end

    @series begin
        subplot := 2
        label := ""
        marker := :d
        yguide := "ΔRR [ms]"
        seriescolor := :red

        (-2:1), res.meanSegment[(-2:1) .+ (L+1)]
    end

    @series begin
        subplot := 2
        seriestype := :scatter
        label := "prsaBRS = $(round(res.value, digits = 2)) ms"
        marker := :circ
        yguide := "ΔRR [ms]"
        seriescolor := :black
        ylims := (minimum(res.meanSegment), maximum(res.meanSegment))

        [0,0], [minimum(res.meanSegment)-10, minimum(res.meanSegment)-10]
    end

    @series begin
        subplot := 2
        seriestype := :scatter
        label := "prsaBRSNorm = $(round(res.valueNorm, digits = 2)) ms/mmHg"
        marker := :circ
        yguide := "ΔRR [ms]"
        seriescolor := :black
        ylims := (minimum(res.meanSegment), maximum(res.meanSegment))

        [0,0], [minimum(res.meanSegment)-10, minimum(res.meanSegment)-10]
    end
end


text:
const DECIMALS = 4 # the maximum number of decimals to show
const PADDING = 2  # number of spaces minimum between fields

str_width(s::AbstractString) = length(s)
str_width(s::Symbol) = length(string(s))
str_width(n::Number) = length(string(n))
str_width(f::Function) = length(string(f))
str_width(b::Bool) = b ? 4 : 5
str_width(c::Char) = 1
str_width(::Missing) = 7

_round(x, n=4) = ismissing(x) ? missing : round(x, digits=n)

function getshowrows(io::IO, x::TS)
    nrow = size(x,1)
    display_rows, _ = displaysize(io)
    display_rows -= 3
    if nrow > display_rows
        row_chunk_size = fld(display_rows-4, 2) - 1
        return (collect(1:row_chunk_size), collect(nrow-row_chunk_size:nrow))
    else
        return (collect(1:nrow), Int[])
    end
end

function getwidths(io::IO, X::TS{V,T}, padding::Int) where {V,T}
    width_index = T==Date ? 10 : 19
    width_values = str_width.(X.fields)
    if V <: Bool
        width_values .= 5
    elseif V <: Number
        @inbounds for j in 1:size(X,2)
            width_values[j] = max(width_values[j], maximum(str_width.(_round.(X.values[:,j], DECIMALS))))
        end
    else
        @inbounds for j in 1:size(X,2)
            width_values[j] = max(width_values[j], maximum(str_width.(X.values[:,j])))
        end
    end
    return width_index + padding, width_values .+ padding
end

hasnegs(x::Vector)::Bool = eltype(x)<:Number ? any(x.<zero(eltype(x))) : false
hasnegs(X::Matrix)::Vector{Bool} = [hasnegs(X[:,j]) for j in 1:size(X,2)]

function show(io::IO, X::TS, padding::Int=PADDING, digits::Int=DECIMALS)::Nothing
    # print summary of data structure
    if (isempty(X))
        print(io, "Empty $(typeof(X))\n")
        return nothing
    else
        print(io, "$(size(X,1))x$(size(X,2)) $(typeof(X)): $(X.index[1]) to $(X.index[end])\n\n")
    end
    # partition rows if data too large
    toprows, bottomrows = getshowrows(io, X)
    negatives = hasnegs(X.values)
    width_index, width_values = getwidths(io, X, padding)
    headerline = [rpad("Index", width_index); [rpad(string(X.fields[j]), width_values[j]) for j in 1:size(X,2)]]
    print(io, join(headerline), '\n')
    if isempty(bottomrows)
        @inbounds for (t, x) in X
            datarow = [
                       rpad(t, width_index);
                       [rpad(string(round(x[j], digits=digits)), width_values[j]) for j in 1:size(X,2)]
                      ]
            print(io, join(datarow), '\n')
        end
    else
        @inbounds for (t, x) in X[toprows]
            datarow = [
                       rpad(t, width_index);
                       [rpad(string(round(x[j], digits=digits)), width_values[j]) for j in 1:size(X,2)]
                      ]
            print(io, join(datarow), '\n')
        end
        print(io, "⋮\n")
        @inbounds for (t, x) in X[bottomrows]
            datarow = [
                       rpad(t, width_index);
                       [rpad(string(round(x[j], digits=digits)), width_values[j]) for j in 1:size(X,2)]
                      ]
            print(io, join(datarow), '\n')
        end
    end
    return nothing
end


text:
abstract type AbstractNode end

# COMPONENT
struct Component{N} <: AbstractNode
    type::Symbol
    name::Symbol
    freeports::MVector{N,Bool}
    vertex::RefValue{Int}
    parameters::Vector{Num}
    states::Vector{Num}
    equations::Vector{Equation}
    defaults::AbstractDict
    function Component{N}(t, n, v::Int, p::Vector, x::Vector, eq::Vector, def::AbstractDict) where {N}
        new(Symbol(t), Symbol(n), ones(MVector{N,Bool}), Ref(v), p, x, eq, def)
    end
end

function Component(type, name = type; vertex::Int = 0, library = BondGraphs.DEFAULT_LIBRARY,
    comp_dict = haskey(library, type) ? library[type] : Dict(),
    numports::Int = haskey(comp_dict, :numports) ? comp_dict[:numports] : 1,
    parameters = haskey(comp_dict, :parameters) ? collect(keys(comp_dict[:parameters])) : Num[],
    states = haskey(comp_dict, :states) ? collect(keys(comp_dict[:states])) : Num[],
    equations = haskey(comp_dict, :equations) ? comp_dict[:equations] : Equation[],
    defaults = haskey(comp_dict, :defaults) ? copy(comp_dict[:defaults]) : Dict())

    Component{numports}(type, name, vertex, parameters, states, equations, defaults)
end

# JUNCTION
abstract type Junction <: AbstractNode end

struct EqualEffort <: Junction
    name::Symbol
    freeports::Vector{Bool}
    weights::Vector{Int}
    vertex::RefValue{Int}
    EqualEffort(; name = Symbol("0"), v::Int = 0) = new(Symbol(name), [true], [0], Ref(v))
end

struct EqualFlow <: Junction
    name::Symbol
    freeports::Vector{Bool}
    weights::Vector{Int}
    vertex::RefValue{Int}
    EqualFlow(; name = Symbol("1"), v::Int = 0) = new(Symbol(name), [true], [0], Ref(v))
end

# Source-sensor
struct SourceSensor <: AbstractNode
    name::Symbol
    freeports::MVector{1,Bool}
    vertex::RefValue{Int}
    SourceSensor(; name=:SS, v::Int = 0) = new(Symbol(name), ones(MVector{1,Bool}) , Ref(v))
end

# PROPERTIES
# Type
type(n::AbstractNode) = n.type
type(j::Junction) = typeof(j)
type(ss::SourceSensor) = :SS

# Name
name(n::AbstractNode) = n.name

# Ports
freeports(n::AbstractNode) = n.freeports
numports(n::AbstractNode) = length(n.freeports)
updateport!(n::AbstractNode, idx::Int) = freeports(n)[idx] = !freeports(n)[idx]

# Weights
weights(j::Junction) = j.weights
set_weight!(j::Junction, idx::Int, w::Int) = j.weights[idx] = w

nextfreeport(n::AbstractNode) = findfirst(freeports(n))
function nextfreeport(j::Junction)
    freeport = findfirst(freeports(j))
    if isnothing(freeport)
        push!(j.freeports, true)
        push!(j.weights, 0)
        return length(j.freeports)
    else
        return freeport
    end
end

# Vertex
vertex(n::AbstractNode) = n.vertex[]
set_vertex!(n::AbstractNode, v::Int) = n.vertex[] = v

# Parameters
parameters(n::AbstractNode) = n.parameters
parameters(::Junction) = Num[]
parameters(::SourceSensor) = Num[]

# State variables
states(n::AbstractNode) = n.states
states(::Junction) = Num[]
states(::SourceSensor) = Num[]

# Equations
equations(n::AbstractNode) = n.equations
equations(::Junction) = Equation[]
equations(::SourceSensor) = Equation[]

# Defaults
defaults(n::AbstractNode) = n.defaults
defaults(::Junction) = Dict{Num,Any}()
defaults(::SourceSensor) = Dict{Num,Any}()

# Set and get default parameter values
function get_parameter(n::AbstractNode, var)
    p, = @parameters $var
    p in keys(defaults(n)) || error("Component does not have parameter $var")
    defaults(n)[p]
end
function set_parameter!(n::AbstractNode, var, val)
    p, = @parameters $var
    p in keys(defaults(n)) || error("Component does not have parameter $var")
    defaults(n)[p] = val
end

# Set and getinitial conditions
function get_initial_value(n::AbstractNode, var)
    _, x = @variables t, $var(t)
    x in keys(defaults(n)) || error("Component does not have state variable $var")
    defaults(n)[x]
end
function set_initial_value!(n::AbstractNode, var, val)
    _, x = @variables t, $var(t)
    x in keys(defaults(n)) || error("Component does not have state variable $var")
    defaults(n)[x] = val
end

# BASE FUNCTIONS
# This definition will need to expand when equations etc. are added
==(n1::AbstractNode, n2::AbstractNode) = type(n1) == type(n2) && n1.name == n2.name

show(io::IO, node::AbstractNode) = print(io, "$(type(node)):$(node.name)")
show(io::IO, node::Junction) = print(io, "$(node.name)")

max_stars_repo_name: mic-pan/BondGraphs

max_issues_repo_name: mic-pan/BondGraphs

max_forks_repo_name: mic-pan/BondGraphs


text:
using .DualNumbers

function blkimpv_dual(num1,num2,num3,num4,num5)
	@eval function blkimpv(S0::$num1,K::$num2,r::$num3,T::$num4,Price::$num5,FlagIsCall::Bool=true,xtol::Real=1e-14,ytol::Real=1e-15)

		blscheck(S0,K,r,T,0.1,r);

		return blsimpv(S0,K,r,T,Price,r,FlagIsCall);

	end
end

type_blkimpv_dual_1=[Dual,Union{Real,Dual},Union{Real,Dual},Union{Real,Dual},Union{Real,Dual}]
type_blkimpv_dual_2=[Real,Dual,Union{Real,Dual},Union{Real,Dual},Union{Real,Dual}]
type_blkimpv_dual_3=[Real,Real,Union{Real,Dual},Union{Real,Dual},Union{Real,Dual}]
blkimpv_dual(type_blkimpv_dual_1[1],type_blkimpv_dual_1[2],type_blkimpv_dual_1[3],type_blkimpv_dual_1[4],type_blkimpv_dual_1[5])
blkimpv_dual(type_blkimpv_dual_2[1],type_blkimpv_dual_2[2],type_blkimpv_dual_2[3],type_blkimpv_dual_2[4],type_blkimpv_dual_2[5])
blkimpv_dual(type_blkimpv_dual_3[1],type_blkimpv_dual_3[2],type_blkimpv_dual_3[3],type_blkimpv_dual_3[4],type_blkimpv_dual_3[5])

function blsimpv_dual(num1,num2,num3,num4,num5,num6)
	@eval function blsimpv(S0::$num1,K::$num2,r::$num3,T::$num4,Price::$num5,d::$num6=0.0,FlagIsCall::Bool=true,xtol::Real=1e-14,ytol::Real=1e-15)
	if (Price< $num5(0))
		throw(ErrorException("Option Price Cannot Be Negative"));
	end
	FinancialToolbox.blscheck(S0,K,r,T,0.1,d);
	value__(x)=x.value;
	value__(x::Real)=x;
	f(x)=(blsprice(value__(S0),value__(K),value__(r),value__(T),x,value__(d),FlagIsCall)-value__(Price));
	σ=FinancialToolbox.brentMethod(f,0.001,1.2,xtol,ytol);
	out=dual(0.0,0.0);

	if(!isdual(Price))
		der_=-(blsprice(S0,K,r,T,σ,d,FlagIsCall)/blsvega(value__(S0),value__(K),value__(r),value__(T),σ,value__(d),FlagIsCall)).epsilon
		out=dual(σ,der_);
	else
		der_=1/blsvega(value__(S0),value__(K),value__(r),value__(T),σ,value__(d),FlagIsCall)
		out=dual(σ,der_);
	end

	return out;

	end
end

type_blsimpv_dual_1=[Dual,Union{Real,Dual},Union{Real,Dual},Union{Real,Dual},Union{Real,Dual},Union{Real,Dual}]
type_blsimpv_dual_2=[Real,Dual,Union{Real,Dual},Union{Real,Dual},Union{Real,Dual},Union{Real,Dual}]
type_blsimpv_dual_3=[Real,Real,Union{Real,Dual},Union{Real,Dual},Union{Real,Dual},Union{Real,Dual}]
blsimpv_dual(type_blsimpv_dual_1[1],type_blsimpv_dual_1[2],type_blsimpv_dual_1[3],type_blsimpv_dual_1[4],type_blsimpv_dual_1[5],type_blsimpv_dual_1[6])
blsimpv_dual(type_blsimpv_dual_2[1],type_blsimpv_dual_2[2],type_blsimpv_dual_2[3],type_blsimpv_dual_2[4],type_blsimpv_dual_2[5],type_blsimpv_dual_2[6])
blsimpv_dual(type_blsimpv_dual_3[1],type_blsimpv_dual_3[2],type_blsimpv_dual_3[3],type_blsimpv_dual_3[4],type_blsimpv_dual_3[5],type_blsimpv_dual_3[6])