text:
# Mean Field Aligned (MFA) coordinates transformation.

using EmpiricalModeDecomposition
using EmpiricalModeDecomposition: zerocrossing!
using Statistics: mean

export mfa, emd_average

"""
    mfa(x, y, z, Bx, By, Bz; method="MAVG", nw=15, fmax=-1.0, fsample=-1.0, threshold=0.1,
       verbose=false) -> Bμ, Bϕ, Bν

Obtain the time series of rotation matrix from geocentric solar ecliptic (GSE) to Mean
field-aligned (MFA) reference frame, given the positions of the spacecraft `x`, `y`, `z` and
magnetic fields `Bx`, `By`, `Bz` in GSE.

# Optional Arguments
- `method::String="MAVG"`: approach to obtain background ambient field, from moving average
`"MAVG"` with `nw` points, to empirical mode decomposition `"EMD"` with cutoff frequency
`fmax`, sampling frequency `fsample` and power selection `threshold` ∈ (0,1).
"""
function mfa(x::AbstractVector, y::AbstractVector, z::AbstractVector,
   Bx::AbstractVector, By::AbstractVector, Bz::AbstractVector;
   method="MAVG", nw=15, fsample=-1.0, fmax=fsample/50 , threshold=0.1, verbose=false)

   if !(length(x) == length(y) == length(z) ==
      length(Bx) == length(By) == length(Bz))
      throw(DimensionMismatch("requires same vector length for inputs!"))
   end

   μx, μy, μz = similar(x), similar(x), similar(x)
   ϕx, ϕy, ϕz = similar(x), similar(x), similar(x)
   νx, νy, νz = similar(x), similar(x), similar(x)
   Bμ, Bϕ, Bν = similar(x), similar(x), similar(x)

   if method == "MAVG"
     B̄x = sma(Bx, nw)
     B̄y = sma(By, nw)
     B̄z = sma(Bz, nw)
   elseif method == "EMD"
     B̄x = emd_average(Bx, fmax, fsample, threshold; verbose)
     B̄y = emd_average(By, fmax, fsample, threshold; verbose)
     B̄z = emd_average(Bz, fmax, fsample, threshold; verbose)
   end

   for i in axes(x, 1)
      B̄mag = hypot(B̄x[i], B̄y[i], B̄z[i])
      μx[i] = B̄x[i] / B̄mag
      μy[i] = B̄y[i] / B̄mag
      μz[i] = B̄z[i] / B̄mag
      ϕx[i] = y[i]*μz[i] - z[i]*μy[i]
      ϕy[i] = z[i]*μx[i] - x[i]*μz[i]
      ϕz[i] = x[i]*μy[i] - y[i]*μx[i]
      rBmag = hypot(ϕx[i], ϕy[i], ϕz[i])
      ϕx[i] /= rBmag
      ϕy[i] /= rBmag
      ϕz[i] /= rBmag
      νx[i] = μy[i]*ϕz[i] - μz[i]*ϕy[i]
      νy[i] = μz[i]*ϕx[i] - μx[i]*ϕz[i]
      νz[i] = μx[i]*ϕy[i] - μy[i]*ϕx[i]

      Bμ[i] = μx[i]*Bx[i] + μy[i]*By[i] + μz[i]*Bz[i]
      Bϕ[i] = ϕx[i]*Bx[i] + ϕy[i]*By[i] + ϕz[i]*Bz[i]
      Bν[i] = νx[i]*Bx[i] + νy[i]*By[i] + νz[i]*Bz[i]
   end

   if verbose
      @info "method = $method"
   end

   Bμ, Bϕ, Bν
end

"""
    emd_average(y, fmax, fsample, threshold=0.1; verbose=false)

Extracts low frequency component from a signal using Empirical Mode Decomposition (EMD).
The signal is reconstructed by using Intrinsic Mode Functions (IMF) which have an average
frequency lower than a specified cut-off frequency `fmax` and a percentage of the power of
the signal, attributable to frequencies lower than `fmax` and greater than `thresold`.

Reference:
The use of the empirical mode decomposition for the identification of mean field aligned
reference frames, M Regi, Alfredo Del Corpo, Marcello De Lauretis, 2016
https://doi.org/10.4401/ag-7067
"""
function emd_average(y, fmax, fsample, threshold=0.1; verbose=false)

   len = length(y)
   x = range(1, length=len, step=1/fsample)

   imfs = emd(y, x)  # intrinsic mode functions

   if length(imfs) == 1
      verbose && @warn "No oscillation detected!"
      return imfs[1]
   end

   res = imfs[end]
   imf = imfs[1:end-1]

   ## Zero-crossing
   # maximum number of zero-crossing is the number of zero_crossing of a sin wave of
   # frequency fmax in a sampling time interval T0 = len/fsample
   nmax = floor(Int, length(y)*(2*fmax)/fsample+1)

   # number of zeros in each IMF
   crossings = [Int[] for _ = 1:length(imfs)-1] # initialization

   for i in 1:length(crossings)
      zerocrossing!(imf[i], crossings[i])
   end

   nzeros = [length(crossings[i]) for i in 1:length(crossings)]

   ## Average period and signal power
   tmin = floor(Int, fsample/fmax)    # miminum period [# of samples]
   tminhalf = floor(Int, fsample/(2*fmax)) # miminum half period [# of samples]
   #T = zeros(length(crossings))
   ind1_ = Vector{Bool}(undef, length(crossings))

   for i = 1:length(imf)
      # indices of zero crossings plus first and last entries
      ind_ = [1, crossings[i]..., len]
      # lengths of consecutive zero-crossing [# of samples]
      tp = diff(ind_)
      #T[i] = mean(2*tp) # average period of the mode [# of samples]

      p = zeros(length(tp))
      # Evaluating the IMF power in each segment
      for j = 1:length(tp)
         p[j] = sum( abs.( imf[i][ind_[j]:ind_[j+1]-1] ) ) / fsample
      end

      ps = p[tp .> tminhalf] # select powers linked to periods greater than tmin

      # select IMF if the percentage of the selected power is greater than the thresold
      ind1_[i] = sum(ps) / sum(p) > threshold

      verbose && @info "IMF $i, threshold = $threshold, power = $(sum(ps) / sum(p))"
   end

   ## Reconstructing ambient signal
   ind2_ = nzeros .< nmax
   imf_bk = imf[ind1_ .| ind2_] # number of IMFs to add to residue to obtain ambient signal

   y_bk = res + sum(imf_bk)

   return y_bk
end


text:
# ----------------------
# -- read             --
# ----------------------

function read_uns(skipnum,uns)
    # skipnum = 9
    ite = Int64(0)

    fff=[]
    open(uns, "r") do f
        fff=read(f,String)
    end
    fff=split(fff,"\n",keepempty=false)
    num_nodes=length(fff)-skipnum

    for i in 1+skipnum:length(fff)
        fff[i]=replace(fff[i]," \r" => "")
        fff[i]=replace(fff[i],"\r" => "")
    end

    ite += skipnum

    # boundary
    ite += 1
    temp = split(fff[ite]," ")
    bdnum = parse(Int64,temp[1])
    for i in 1:bdnum
        # 境界条件処理
        ite += 1
    end

    # skip
    ite += 4

    # Nodes
    ite += 1
    temp = split(fff[ite]," ")
    nodes_num = parse(Int64,temp[1])
    nodes = zeros(nodes_num,3)

    for i in 1:nodes_num
        ite +=1
        temp = split(fff[ite]," ")

        k = 1
        for j in 1:length(temp)
            if temp[j] != ""
                nodes[i,k] = parse(Float64,temp[j])
                k += 1
            end
        end
    end

    println("read nodes ")

    #=
    k = 1
    for j in 1:length(temp)
        if temp[j] != ""
            readQ[i,k] = parse(Float64,temp[j])
            k += 1
        end
    end
    =#

    # boundary face
    ite += 2
    temp = split(fff[ite]," ")
    bdface_num = parse(Int64,temp[1])
    for i in 1:bdface_num
        # 境界条件処理
        ite += 1
    end

    # elements
    ite += 1
    temp = split(fff[ite]," ")

    k=0
    ite_temp = ite
    while k == 0
        if occursin("Variables",fff[ite_temp]) == true
            k=1
        end

        ite_temp += 1
    end

    println("read elements ")

    ele_num = ite_temp - ite -2
    element = zeros(Int64,ele_num,8)
    element_shape = zeros(Int64,ele_num)

    for i in 1:ele_num
        ite +=1
        temp = split(fff[ite]," ")

        nnn = parse(Int64,temp[1])
        if nnn == 1
            element_shape[i] = 4
        elseif nnn == 2
            element_shape[i] = 8
        elseif nnn == 3
            element_shape[i] = 6
        elseif nnn == 4
            element_shape[i] = 5
        else
            " check element in .uns !! \n"
            throw(UndefVarError(:x))
        end

        for j in 1:element_shape[i]
            element[i,j] = parse(Int64,temp[2+j])
        end
    end

    return nodes, nodes_num, element, element_shape
end

function read_result(skipnum)
    fff=[]
    open("test333.dat", "r") do f
        fff=read(f,String)
    end
    fff=split(fff,"\n",keepempty=false)
    num_point=length(fff)-skipnum

    for i in 1+skipnum:length(fff)
        fff[i]=replace(fff[i]," \r" => "")
    end

    readQ = zeros(num_point,5)
    for i in 1:num_point
        temp=split(fff[i+skipnum]," ")

        k = 1
        for j in 1:length(temp)
            if temp[j] != ""
                readQ[i,k] = parse(Float64,temp[j])
                k += 1
            end
        end

    end
    return readQ
end

function read_allgrid(uns)
    skip=9
    nodes, nodes_num, element, element_shape = read_uns(skip,uns)

    #readQ     = read_result(skip)
    println("fin read grid")
    return nodes, nodes_num, element, element_shape
end

ext: jl


text:
## U updates

"""
    updateU!(M::HePPCATModel,Y,method)

Update the factor eigenvectors `M.U` with respect to data `Y` using `method`.
"""
function updateU! end

# Update method: Minorize Maximize (with linear minorizer)
"""
    updateU!(M::HePPCATModel,Y,::MinorizeMaximize)

Minorize maximize update of `M.U` using a linear minorizer.
"""
function updateU!(M::HePPCATModel,Y,::MinorizeMaximize)
    M.U .= polar(gradF(M.U,M.λ,M.v,Y))
    return M
end

# Update method: Projected Gradient Ascent
"""
    updateU!(M::HePPCATModel,Y,pga::ProjectedGradientAscent)

Projected gradient ascent update of `M.U`.
Supported step size types:
`Number` (i.e., constant step size),
[`InverseLipschitz`](@ref),
and
[`ArmijoSearch`](@ref) (experimental).
"""
function updateU!(M::HePPCATModel,Y,pga::ProjectedGradientAscent{<:Number})
    if pga.stepsize == Inf
        M.U .= polar(gradF(M.U,M.λ,M.v,Y))
    else
        M.U .= polar(M.U + pga.stepsize*gradF(M.U,M.λ,M.v,Y))
    end
    return M
end
function updateU!(M::HePPCATModel,Y,pga::ProjectedGradientAscent{<:InverseLipschitz})
    return updateU!(M,Y,ProjectedGradientAscent(inv(pga.stepsize.bound(M,Y))))
end
function updateU!(M::HePPCATModel,Y,pga::ProjectedGradientAscent{<:ArmijoSearch})
    params = pga.stepsize

    dFdU = gradF(M.U,M.λ,M.v,Y)
    F0, FΔ = F(M.U,M.λ,M.v,Y), params.tol * norm(dFdU)^2
    m = findfirst(IdentityUnitRange(0:params.maxsearches-1)) do m
        Δ = params.contraction^m * params.stepsize
        F(polar(M.U + Δ*dFdU),M.λ,M.v,Y) >= F0 + Δ * FΔ
    end
    if isnothing(m)
        @warn "Exceeded maximum line search iterations. Accuracy not guaranteed."
        m = params.maxsearches
    end

    Δ = params.contraction^m * params.stepsize
    M.U .= polar(M.U + Δ*dFdU)
    return M
end

# Update method: Stiefel Gradient Ascent
"""
    updateU!(M::HePPCATModel,Y,sga::StiefelGradientAscent)

Stiefel gradient ascent update of `M.U`.
Supported step size types:
[`ArmijoSearch`](@ref).
"""
function updateU!(M::HePPCATModel,Y,sga::StiefelGradientAscent{<:ArmijoSearch})
    params = sga.stepsize

    dFdU = gradF(M.U,M.λ,M.v,Y)
    ∇F = dFdU - M.U*(dFdU'M.U)

    F0, FΔ = F(M.U,M.λ,M.v,Y), params.tol * norm(∇F)^2
    m = findfirst(IdentityUnitRange(0:params.maxsearches-1)) do m
        Δ = params.contraction^m * params.stepsize
        F(geodesic(M.U,∇F,Δ),M.λ,M.v,Y) >= F0 + Δ * FΔ
    end
    if isnothing(m)
        @warn "Exceeded maximum line search iterations. Accuracy not guaranteed."
        m = params.maxsearches
    end

    Δ = params.contraction^m * params.stepsize
    M.U .= geodesic(M.U,∇F,Δ)
    return M
end

# Objective with respect to U
"""
    F(U,λ,v,Y)

Objective for the optimization problem w.r.t `U`.
"""
F(U,λ,v,Y) = 1/2*sum(norm(sqrt(Diagonal(λ./vl./(λ.+vl)))*U'*Yl)^2 for (Yl,vl) in zip(Y,v))

"""
    gradF(U,λ,v,Y)

Euclidean gradient of objective for the optimization problem w.r.t `U`.
"""
gradF(U,λ,v,Y) = sum(Yl * Yl' * U * Diagonal(λ./vl./(λ.+vl)) for (Yl,vl) in zip(Y,v))

"""
    LipBoundU1(M::HePPCATModel,Y)

Lipschitz bound of objective w.r.t. `U` at `M` with data `Y`.
"""
function LipBoundU1(M::HePPCATModel,Y)
    L, λmax = length(M.v), maximum(M.λ)
    return sum(norm(Y[l])^2*λmax/M.v[l]/(λmax+M.v[l]) for l in 1:L)
end

"""
    LipBoundU2(M::HePPCATModel,Y)

Lipschitz bound of objective w.r.t. `U` at `M` with data `Y`.
"""
function LipBoundU2(M::HePPCATModel,Y)
    L, λmax = length(M.v), maximum(M.λ)
    return sum(opnorm(Y[l])^2*λmax/M.v[l]/(λmax+M.v[l]) for l in 1:L)
end

# Utilities
"""
    polar(A)

Polar factor of `A`.
"""
function polar(A)
    F = svd(A)
    return F.U*F.Vt
end

"""
    skew(A)

Skew-symmetrize `A` by computing `(A-A')/2`.
"""
skew(A) = (A-A')/2

"""
    geodesic(U,X,t)

Geodesic step of size `t` from `U` in direction of `X` along the Stiefel manifold.
"""
function geodesic(U,X,t)
    k = size(U,2)

    A = skew(U'X)
    Q,R = qr(X - U*(U'X))

    MN = exp(t*[A -R'; R zeros(k,k)])[:,1:k]
    M, N = MN[1:k,:], MN[k+1:end,:]

    return U*M + Matrix(Q)*N
end


text:
using HDF5
using FastGaussQuadrature
using LaTeXStrings
using LinearAlgebra
using Particles
using Plots
using Random
using ReducedBasisMethods
using Statistics

# HDF5 file containing training data
fpath = "../runs/BoT_Np5e4_k_010_050_np_10_T25.h5"

# read sampling parameters
params = read_sampling_parameters(fpath)

# read training parameters
μₜ = h5read(fpath, "parameters/mu_train")

# read integrator parameters
IP = IntegratorParameters(fpath)

# create spline Poisson solver
poisson = PoissonSolverPBSplines(fpath)

# read snaptshot data
X = h5read(fpath, "snapshots/X")
V = h5read(fpath, "snapshots/V")
E = h5read(fpath, "snapshots/E")
# D = h5read(fpath, "snapshots/D")
# Φ = h5read(fpath, "snapshots/Phi")

# Reference draw
P₀ = ParticleList(X[:,1], V[:,1], ones(IP.nₚ) .* poisson.L ./ IP.nₚ)

# EVD

# XV = hcat(X, V);
XV = copy(X)
for p in 1:IP.nparam
    # XV = hcat(XV, V[:,1+(p-1)*IP.nₛ])
    global XV = hcat(XV, V[:,1+(p-1)*IP.nₛ]) #, E[:,p*IP.nₛ])
    #print(1+(p-1)*IP.nₛ, " ", p*IP.nₛ, "\n" )
end

# clear
X, V = 0, 0
GC.gc()

# @time F = eigen(Xₑₓₜ' * Xₑₓₜ)
# Λ, Ω = sorteigen(F.values, F.vectors)
@time F = eigen(XV' * XV)
@time Λ, Ω = sorteigen(F.values, F.vectors)

@time Fₑ = eigen(E' * E)
@time Λₑ, Ωₑ = sorteigen(Fₑ.values, Fₑ.vectors)

# clear
F, Fₑ = 0, 0
GC.gc()

# plot
plot(xlabel = L"$i$", ylabel = L"$\lambda_i$", yscale = :log10,
     grid = true, gridalpha = 0.5)
plot!(abs.(Λ )[1:1000], linewidth = 2, alpha = 0.25, label = L"$X$")
plot!(abs.(Λₑ)[1:1000], linewidth = 2, alpha = 0.5,  label = L"$F$")
savefig("../runs/BoT_Np5e4_k_010_050_np_10_T25_SVDs_BoT_1.pdf")

# plot
plot(xlabel = L"$i$", ylabel = L"$\lambda_i$", yscale = :log10,
     grid = true, gridalpha = 0.5, legend = :none)
plot!(abs.(Λ ), linewidth = 2, alpha = 0.25, label = L"$X_v$")
plot!(abs.(Λₑ), linewidth = 2, alpha = 0.5,  label = L"$E$")
savefig("../runs/BoT_Np5e4_k_010_050_np_10_T25_SVDs_BoT_2.pdf")

# Projection Matrices
@time k, Ψ = get_Ψ(XV, Λ, Ω, 1e-9, 0)
@time kₑ, Ψₑ = get_Ψ(E, Λₑ, Ωₑ, 1e-5, 0)

# clear
XV, E = 0, 0
Λ, Λₑ, Ω, Ωₑ = 0, 0, 0, 0
GC.gc()

# DEIM
@time Πₑ = deim_get_Π(Ψₑ)

# save to HDF5
h5save("../runs/BoT_Np5e4_k_010_050_np_10_T25_projections.h5", IP, poisson, params, μₜ, k, kₑ, Ψ, Ψₑ, Πₑ, P₀)


text:
# using Base.Test
# using MIToS.MSA

print("""
Tests for Mappings on Annotated Multiple Sequence Alignments
------------------------------------------------------------
""")

let pfam = read(joinpath(pwd(), "data", "PF09645_full.stockholm"), Stockholm,
                generatemapping=true, useidcoordinates=true),
    F112_SSV1 = collect(string(".....QTLNSYKMAEIMYKILEKKGELTLEDILAQFEISVPSAYNIQRALKAICERHPDECEVQY",
                               "KNRKTTFKWIKQEQKEEQKQEQTQDNIAKIFDAQPANFEQTDQGFIKAKQ....."))

  @test minimum(getsequencemapping(pfam, 4)) == 3
  @test maximum(getsequencemapping(pfam, "F112_SSV1/3-112")) == 112
  @test getcolumnmapping(pfam) == getindex([ i for i in 1:length(F112_SSV1) ], F112_SSV1 .!= '.')
  @test length(pfam.annotations.sequences) == 9
end

let fasta = read(joinpath(pwd(), "data", "Gaoetal2011.fasta"), FASTA, generatemapping=true)

  @test getsequencemapping(fasta, 1) == [1, 2, 3, 4, 5, 6]
  @test getsequencemapping(fasta, "SEQ2") == [1, 2, 3, 4, 5, 6]
  @test getcolumnmapping(fasta) == [1, 2, 3, 4, 5, 6]
end

print("""

Tests for annotation of modifications
-------------------------------------
""")

let pfam = read(joinpath(pwd(), "data", "PF09645_full.stockholm"), Stockholm,
                generatemapping=true, useidcoordinates=true), io = IOBuffer()

  print("""
  Test printing annotations and modifications
  """)
  print(io, pfam.annotations)
  printed = split(takebuf_string(io), '\n')
  @test length(printed) == 17
  @test printed[1] == string("#=GF NCol   120")
  @test printed[2] == string("#=GF ColMap   6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22",
                             ",23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,",
                             "43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,",
                             "63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,",
                             "84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,",
                             "103,104,105,106,107,108,109,110,111,112,113,114,115")
  @test ismatch(r"MIToS_", printed[3])
  @test ismatch(r"MIToS_", printed[4])

  print("""
  Test delete MIToS annotations
  """)
  delete_annotated_modifications!(pfam)
  print(io, pfam.annotations)
  printed = split(takebuf_string(io), '\n')
  @test length(printed) == 17 - 2 # 2 MIToS annotations
  @test !ismatch(r"MIToS_", printed[3])
  @test !ismatch(r"MIToS_", printed[4])
end


text:
module BEASTDataPrep

using CSV, DataFrames, Gadfly, Statistics, PhyloNetworks
using HypothesisTests

export parse_continuous_data,
       plot_transformed_data,
       log_transform_data!,
       standardize_data!,
       conform_tree_and_data,
       conform_traitdata,
       merge_traitdata

include("parse_data.jl")
include("transform_and_scale.jl")
include("trim_taxa.jl")

end


text:
function init_observable_mappings!(m::RealBondMkup)

    observables = OrderedDict{Symbol,Observable}()
    population_mnemonic = get(get_setting(m, :population_mnemonic))

    ############################################################################
    ## 1. GDP
    ############################################################################
    gdp_fwd_transform = function (levels)
        # FROM: Level of nominal GDP (FRED :GDP series)
        # TO:   The cyclical component of log real GDP

        levels[:rGDP] = nominal_to_real(:GDP, levels)
        gdp = oneqtrpctchange(levels[:rGDP])
        _, levels[:GDP_cyc_comp] = hpfilter(gdp, 1600)
        levels[:GDP_cyc_comp]
        #=levels[:lGDP] = log.(levels[:rGDP])
        _, levels[:lGDP_cyc_comp] = hpfilter(levels[:lGDP], 1600)

        levels[:lGDP_cyc_comp]=#
    end

    gdp_rev_transform = logleveltopct_annualized

    observables[:obs_gdp] = Observable(:obs_gdp, [:GDP__FRED, population_mnemonic, :GDPDEF__FRED],
                                       gdp_fwd_transform, gdp_rev_transform,
                                       "Log Real GDP", "")

    ############################################################################
    ## 2. Core PCE Inflation
    ############################################################################

    pce_fwd_transform = function (levels)
        # FROM: Core PCE index
        # INTO: Approximate quarter-to-quarter percent change of Core PCE,
        # i.e. quarterly core pce inflation

        oneqtrpctchange(levels[:PCEPILFE])
    end

    pce_rev_transform = loggrowthtopct_annualized

    observables[:obs_corepce] = Observable(:obs_corepce, [:PCEPILFE__FRED],
                                           pce_fwd_transform, pce_rev_transform,
                                           "Core PCE Inflation",
                                           "Core PCE Inflation")

    ############################################################################
    ## 3. Nominal short-term interest rate (3 months)
    ############################################################################

    nominalrate_fwd_transform = function (levels)
        # FROM: Nominal effective federal funds rate (aggregate daily data at a
        #       quarterly frequency at an annual rate)
        # TO:   Nominal effective fed funds rate, at a quarterly rate

        annualtoquarter(levels[:DFF])
    end

    nominalrate_rev_transform = quartertoannual

    observables[:obs_nominalrate] = Observable(:obs_nominalrate, [:DFF__FRED],
                                               nominalrate_fwd_transform, nominalrate_rev_transform,
                                               "Nominal FFR",
                                               "Nominal Effective Fed Funds Rate")

    m.observable_mappings = observables
end


text:
using LibBSON

bsonObject = BSONObject(Dict(
    "null" => nothing,
    "bool" => true,
    "int" => 42,
    "double" => 3.141,
    "string" => "Hello, Jérôme",
    "oid" => BSONOID(),
    "minkey" => :minkey,
    "maxkey" => :maxkey,
    "array" => Any[5.41, false]
    ))
println(bsonObject)
println(bsonObject["string"])
for (k, v) in bsonObject
    println("$k => $v")
end
bsonArray = BSONArray(Any[
    "one",
    Dict("key" => 6.7)
    ])
for e in bsonArray
    println(e)
end


text:
using Crystallography
using Test

@testset "Crystallography.jl" begin
end


text:
#InspectDR: Defaults
#-------------------------------------------------------------------------------

#==Constants: Initial default values
===============================================================================#
#Default values for data area dimensions (used to save single plot):
const DEFAULT_DATA_WIDTH = 500.0
const DEFAULT_DATA_HEIGHT = DEFAULT_DATA_WIDTH / φ #Use golden ratio

#Default values for plot dimensions (used to save multi-plot):
const DEFAULT_PLOT_WIDTH = 600.0
const DEFAULT_PLOT_HEIGHT = DEFAULT_PLOT_WIDTH / φ #Use golden ratio

#Default font:
const DEFAULT_FONTNAME = @windows? "Cambria": "Serif"
#Cairo "built-in": Serif, Sans, Serif, Fantasy, Monospace

#==Types
===============================================================================#
type Defaults
	rendersvg::Bool #Might want to dissalow SVG renderings for performance reasons
	showtimestamp::Bool
	fontname::AbstractString
	fontscale::Float64

	#Default values for data-area dimensions (saving single plot):
	wdata::Float64
	hdata::Float64

	#Default values for plot dimensions (saving multi-plot):
	wplot::Float64
	hplot::Float64
end
Defaults() =
	Defaults(true, false, DEFAULT_FONTNAME, 1,
		DEFAULT_DATA_WIDTH, DEFAULT_DATA_HEIGHT,
		DEFAULT_PLOT_WIDTH, DEFAULT_PLOT_HEIGHT,
	)

#==Data
===============================================================================#
function __init__()
	global const defaults = Defaults()
	local userdefaults

	try; userdefaults = Main.DEFAULTS_INSPECTDR
	catch; return;	end

	for (k, v) in userdefaults
		try
			setfield!(defaults, k, v)
		catch
			warn("Cannot set InspectDR.defaults.$k = $v.")
		end
	end
end

#==Functions
===============================================================================#
#Forces a function call to read *current* settings during parameter initialization:
#=Not needed:
function getdefaults(field::Symbol)
	global defaults
	return getfield(defaults, field)
end
=#

#Last line


text:
"""
    `getmatfail(sigma1::Real, sigma2::Real, tau12::Real, xt::Real,
        xc::Real, yt::Real, yc::Real, s::Real, method::String)`

    `getmatfail(sigma1::Real, sigma2::Real, tau12::Real, mat::Material,
        method::String)`

    `getmatfail(stress::AbstractArray{<:Real,1}, mat::Material, method::String)`

    `getmatfail(stress::AbstractArray{<:Real,1},xt::Real,xc::Real,yt::Real,
        yc::Real,s::Real, method::String)`

    `getmatfail(stress::AbstractArray{<:AbstractArray{<:Real,1},1},
        mat::AbstractArray{Material}, lam::Laminate, method::String)`

Determines ply failure according to specified theory. Values greater than one
signify ply failure. Choose between "maxstress", "tsaiwu", and "hashinrotem".
Safety factors are stored in the second return argument.

# Theory Specific Output
Maximum Ply Stress:
1. Longitudinal tension
2. Longitudinal compression
3. Transverse tension
4. Transverse compression
5. Positive shear
6. Negative shear

Tsai-Wu
1. Failure Criterion

Hashin-Rotem
1. Fiber failure in tension
2. Fiber failure in compression
3. Matrix failure in tension
4. Matrix failure in compression
"""
function getmatfail(sigma1::Real, sigma2::Real, tau12::Real, xt::Real,
    xc::Real, yt::Real, yc::Real, s::Real, method::AbstractString)

    if method == "maxstress"
        return maxstress(sigma1, sigma2, tau12, xt, xc, yt, yc, s)
    elseif method == "tsaiwu"
        return tsaiwu(sigma1, sigma2, tau12, xt, xc, yt, yc, s)
    elseif method == "hashinrotem"
        return hashinrotem(sigma1, sigma2, tau12, xt, xc, yt, yc, s)
    else
        error("Chosen method ($method) not implemented. Choose between `maxstress`,
            `tsaiwu`, and `hashinrotem`")
    end
end

getmatfail(stress::AbstractArray{<:Real,1}, mat::Material, method::AbstractString) =
    getmatfail(stress[1], stress[2], stress[3], mat.xt, mat.xc, mat.yt, mat.yc,
    mat.s, method)
getmatfail(stress::AbstractArray{<:Real,1}, xt::Real, xc::Real, yt::Real,
    yc::Real, s::Real, method::AbstractString) = getmatfail(stress[1], stress[2],
    stress[3], xt, xc, yt, yc, s, method)
getmatfail(sigma1::Real, sigma2::Real, tau12::Real, mat::Material, method::AbstractString) =
    getmatfail(sigma1, sigma2, tau12, mat.xt, mat.xc, mat.yt, mat.yc, mat.s, method)

function getmatfail(stress::AbstractArray{<:AbstractArray{<:Real,1},1},
    mat::AbstractArray{<:Material,1}, lam::Laminate, method::AbstractString)

    result = [getmatfail(stress[i][1], stress[i][2], stress[i][3],
        mat[lam.matid[i]].xt, mat[lam.matid[i]].xc,
        mat[lam.matid[i]].yt, mat[lam.matid[i]].yc,
        mat[lam.matid[i]].s, method) for i in 1:length(lam.matid)]

    matfail = [result[i][1] for i in 1:length(lam.matid)]
    sf = [result[i][2] for i in 1:length(lam.matid)]

    return matfail, sf
end

function maxstress(sigma1::Real, sigma2::Real, tau12::Real, xt::Real,
    xc::Real, yt::Real, yc::Real, s::Real)

    fail = [sigma1/xt, -sigma1/xc, sigma2/yt, -sigma2/yc, tau12/s, -tau12/s]

    safetyfactor = 1.0./fail

    return fail, safetyfactor
end
maxstress(stress::AbstractArray{<:Real,1}, mat::Material) =
    maxstress(stress[1], stress[2], stress[3], mat.xt, mat.xc, mat.yt, mat.yc, mat.s)
maxstress(stress::AbstractArray{<:Real,1}, xt::Real, xc::Real, yt::Real,
    yc::Real, s::Real) = maxstress(stress[1], stress[2], stress[3], xt, xc, yt, yc, s)
maxstress(sigma1::Real, sigma2::Real, tau12::Real, mat::Material) =
    maxstress(sigma1, sigma2, tau12, mat.xt, mat.xc, mat.yt, mat.yc, mat.s)

function maxstress(stress::AbstractArray{<:AbstractArray{<:Real,1},1},
    mat::AbstractArray{<:Material,1}, lam::Laminate)

    result = [maxstress(stress[i][1], stress[i][2], stress[i][3],
        mat[lam.matid[i]].xt, mat[lam.matid[i]].xc,
        mat[lam.matid[i]].yt, mat[lam.matid[i]].yc,
        mat[lam.matid[i]].s, method) for i in 1:length(lam.matid)]

    matfail = [result[i][1] for i in 1:length(lam.matid)]
    sf = [result[i][2] for i in 1:length(lam.matid)]

    return matfail, sf
end

function tsaiwu(sigma1::Real, sigma2::Real, tau12::Real, xt::Real,
    xc::Real, yt::Real, yc::Real, s::Real)

    a = (sigma1^2.0/(xt*xc))+(sigma2^2.0/(yt*yc))-
         sqrt(1.0/(xt*xc)*1.0/(yt*yc))*sigma1*sigma2+tau12^2.0/s^2.0
    b = (1.0/xt-1.0/xc)*sigma1+(1.0/yt-1.0/yc)*sigma2
    c = -1
    fail = a+b
    safetyfactor = abs(fail) > 1e-20 ? (-b+sqrt(b^2-4*a*c))/(2a) : Inf

    return fail, safetyfactor
end
tsaiwu(stress::AbstractArray{<:Real,1}, mat::Material) =
    tsaiwu(stress[1], stress[2], stress[3], mat.xt, mat.xc, mat.yt, mat.yc, mat.s)
tsaiwu(stress::AbstractArray{<:Real,1}, xt::Real, xc::Real, yt::Real,
    yc::Real, s::Real) = tsaiwu(stress[1], stress[2], stress[3], xt, xc, yt, yc, s)
tsaiwu(sigma1::Real, sigma2::Real, tau12::Real, mat::Material) =
    tsaiwu(sigma1, sigma2, tau12, mat.xt, mat.xc, mat.yt, mat.yc, mat.s)

function tsaiwu(stress::AbstractArray{<:AbstractArray{<:Real,1},1},
    mat::AbstractArray{<:Material,1}, lam::Laminate)

    result = [tsaiwu(stress[i][1], stress[i][2], stress[i][3],
        mat[lam.matid[i]].xt, mat[lam.matid[i]].xc,
        mat[lam.matid[i]].yt, mat[lam.matid[i]].yc,
        mat[lam.matid[i]].s, method) for i in 1:length(lam.matid)]

    matfail = [result[i][1] for i in 1:length(lam.matid)]
    sf = [result[i][2] for i in 1:length(lam.matid)]

    return matfail, sf
end

function hashinrotem(sigma1::Real, sigma2::Real, tau12::Real, xt::Real,
    xc::Real, yt::Real, yc::Real, s::Real)

    tensionfail = (sigma2^2.0/yt^2.0+tau12^2.0/s^2.0)
    compressionfail = (sigma2^2.0/yc^2.0+tau12^2.0/s^2.0)
    fail = [sigma1/xt, -sigma1/xc, sign(sigma2)*tensionfail, sign(sigma2)*-compressionfail]
    safetyfactor = [xt/sigma1, -xc/sigma1, 1/sqrt(tensionfail), sign(sigma2)*-1/sqrt(compressionfail)]

    return fail, safetyfactor
end
hashinrotem(stress::AbstractArray{<:Real,1}, mat::Material) =
    hashinrotem(stress[1], stress[2], stress[3], mat.xt, mat.xc, mat.yt, mat.yc, mat.s)
hashinrotem(stress::AbstractArray{<:Real,1}, xt::Real, xc::Real, yt::Real,
    yc::Real, s::Real) = hashinrotem(stress[1], stress[2], stress[3], xt, xc, yt, yc, s)
hashinrotem(sigma1::Real, sigma2::Real, tau12::Real, mat::Material) =
    hashinrotem(sigma1, sigma2, tau12, mat.xt, mat.xc, mat.yt, mat.yc, mat.s)

function hashinrotem(stress::AbstractArray{<:AbstractArray{<:Real,1},1},
    mat::AbstractArray{<:Material,1}, lam::Laminate)

    result = [hashinrotem(stress[i][1], stress[i][2], stress[i][3],
        mat[lam.matid[i]].xt, mat[lam.matid[i]].xc,
        mat[lam.matid[i]].yt, mat[lam.matid[i]].yc,
        mat[lam.matid[i]].s, method) for i in 1:length(lam.matid)]

    matfail = [result[i][1] for i in 1:length(lam.matid)]
    sf = [result[i][2] for i in 1:length(lam.matid)]

    return matfail, sf
end


text:
### Visualize the whole search structure of a problem

function search_structure(st)

end

"""
    search_structure(v0, trace)

"""
function search_structure(v0, trace, draw_vertex, draw_edge, draw_active; cell_width = 20, cell_height = 20, fname = "sstructure.svg", dir = :min, rev = false)
    X, Y, Qa, Qd = generate_coordinates(v0, trace, dir = dir, rev = rev)
    m = last(X)
    n = maximum(Y)

    @show n, m
    @show X
    @show Y
    @show Qa
    @show Qd
    d = Drawing((m+3)*cell_width, (n+3)*cell_height, fname)

    origin()
    sethue("black")
    setfont("STIX Two Math", 12)
    background("white")

    t = Table(n, m, cell_width, cell_height)
    drawgrid(t)

    # Drawing
    foreach((x,y) -> draw_vertex(t[x,y]), Y,X)
    foreach(e -> draw_edge(t[e[1]...], t[e[2]...]), Qd)
    foreach(e -> draw_active(t[e[1]...], t[e[2]...]), Qa)

    finish()
    preview()
end

function drawgrid(t)
    n,m = size(t)
    o = t[n,1]
    x = t[n,m] + Point(20,0)
    y = t[1,1] - Point(0,20)
    sethue("lightgray")
    setline(1)
    hline = y -> line(t[1,y], t[n,y], :stroke)
    vline = x -> line(t[x,1], t[x,m], :stroke)
    foreach(idx -> hline(idx), 1:m)
    foreach(idx -> vline(idx), 1:n)
    setline(4)
    sethue("black")
    tickdrawer = p -> circle(p, 1, :fill)
    axis(o, x, "𝓷"; ticks = t[n,:], drawticks = tickdrawer)
    axis(o, y, "𝓡"; ticks = t[:,1], drawticks = tickdrawer)

end

function axis(t0, t1, l1; ticks = nothing, drawticks = x -> x, valign = "center", halign = "center")
    arrow(t0, t1)
    foreach( d -> drawticks(d), ticks)
    settext(l1, t1)
end

"""
    generate_searchcoordinates(v0, trace::TraceLogger)

Given an initial solution and a search trace, compute the coordinates required to display the search process as the search for an improving path through the ranked search space.
"""
function generate_coordinates(v0, trace; dir = :min, rev = true)
    A = trace.accepted
    V = compute_objectives(v0, trace.deltas, A)
    R = rank_objectives(V, dir = dir, rev = rev)
    curr = (1,R[v0])
    it = 1
    X = Vector{Int}(undef,length(V))
    Y = Vector{Int}(undef,length(V))
    Qa = Vector{Tuple{Tuple{Int,Int},Tuple{Int,Int}}}()
    Qd = Vector{Tuple{Tuple{Int,Int},Tuple{Int,Int}}}()
    for (i,v) in enumerate(V)
        X[i] = it
        Y[i] = R[v]
        node = (R[v], it)
        if i != 1
            if A[i]
                it += 1
                push!(Qa, (curr,node))
                curr = node
            else
                push!(Qd, (curr,node))
            end
        end
    end
    return X,Y,Qa,Qd
end

"""

"""
function rank_objectives(V; dir = :min, rev = false)
    if dir == :min
        Vsorted = unique(sort(V, rev = rev))
    else
        Vsorted = unique(sort(V, rev = !rev))
    end
    return Dict(Vsorted .=> 1:length(Vsorted))
end

##
"""
    compute_objectives(v0, Δv, acceptance)

Compute the objective value of each solution considered in a search trace, baced on the v0, move-Δs and acceptance.
"""
function compute_objectives(v0, ΔV, A)
    vc = v0
    V = Array{eltype(ΔV), 1}(undef, length(ΔV)+1)
    V[1] = vc
    pushfirst!(A, true)
    for (i,Δv,a) in Iterators.zip(1:length(ΔV), ΔV, A)
        V[i+1] = vc + Δv
        if a
            vc += Δv
        end
    end
    return V
end

###

max_stars_repo_name: Michiel-VL/PolyaViz

max_issues_repo_name: Michiel-VL/PolyaViz

max_forks_repo_name: Michiel-VL/PolyaViz


text:
@everywhere using POMDPToolbox
@everywhere using POMDPModels

pomdp = BabyPOMDP()
fwc = FeedWhenCrying()
rnd = solve(RandomSolver(MersenneTwister(7)), pomdp)

q = [] # vector of the simulations to be run
push!(q, Sim(pomdp, fwc, max_steps=32, rng=MersenneTwister(4), metadata=Dict(:policy=>"feed when crying")))
push!(q, Sim(pomdp, rnd, max_steps=32, rng=MersenneTwister(4), metadata=Dict(:policy=>"random")))

# this creates two simulations, one with the feed-when-crying policy and one with a random policy

data = run_parallel(q)

# by default, the dataframe output contains the reward and the contents of `metadata`
@show data
# data = 2×2 DataFrames.DataFrame
# │ Row │ policy             │ reward   │
# ├─────┼────────────────────┼──────────┤
# │ 1   │ "feed when crying" │ -4.5874  │
# │ 2   │ "random"           │ -27.4139 │

# to perform additional analysis on each of the simulations one can define a processing function with the `do` syntax:
data2 = run_parallel(q,progress=false) do sim, hist
    println("finished a simulation - final state was $(last(state_hist(hist)))")
    return [:steps=>n_steps(hist), :reward=>discounted_reward(hist)]
end

@show data2

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
#=
SPUNG = State Research Program for Utilization of Natural Gas

based on:
Unification of the two-parameter equation of state and the principle ofcorresponding states
Jørgen Mollerup

the idea is to use mollerup shape factors and PropaneRef as the EoS
but any function that provides a_scaled and shape_factors will suffice
=#

struct SPUNG{S<:EoSModel,E<:EoSModel} <: EoSModel
    shape_model::S
    shape_ref::S
    model_ref::E
end

function SPUNG(components::Vector{String},refmodel=PropaneRef(),shapemodel::SHAPE=SRK(components)) where SHAPE<:EoSModel
    refname = component_names(refmodel)
    shape_ref = SHAPE.name.wrapper(refname)
    model = SPUNG(shapemodel,shape_ref,refmodel)
    return model
end

function eos(model::SPUNG,V,T,z=SA[1.0],phase="unknown")
    f,h = shape_factors(model,V,T,z)
    T0 = T/f
    V0 = V/h
    return eos(model.model_ref,V0,T0)
end

function shape_factors(model::SPUNG{<:ABCubicModel},V,T,z=SA[1.0])
    n = sum(z)
    x = z * (1/n)
    a,b = cubic_ab(model.shape_model,T,x)
    a0,b0 = cubic_ab(model.shape_ref,T,SA[1.0])
    h = b/b0
    fh = n*a/a0
    f = fh/h
    return f,h
end

mw(model::SPUNG) = mw(model.shape_model)

function Base.show(io::IO,mime::MIME"text/plain",model::SPUNG)
    println(io,"Extended Corresponding States model")
    println(io," reference model: ",model.model_ref)
    print(io," shape model: ",model.shape_model)
end

function Base.show(io::IO,model::SPUNG)
    print(io,"SPUNG(",model.shape_ref,",",model.model_ref,")")
end

function lb_volume(model::SPUNG,z=SA[1.0];phase=:l)
    lb_v0 = lb_volume(model.model_ref)
    T0 = T_scale(model.model_ref)
    f,h = shape_factors(model,lb_v0,T0,z) #h normaly should be independent of temperature
    return lb_v0*h
end

function T_scale(model::SPUNG,z=SA[1.0])
    lb_v0 = lb_volume(model.model_ref)
    T0 = T_scale(model.model_ref)
    f,h = shape_factors(model,lb_v0,T0,z) #h normaly should be independent of temperature
    return T0*f
end

function p_scale(model::SPUNG,z=SA[1.0])
    lb_v0 = lb_volume(model.model_ref)
    T0 = T_scale(model.model_ref)
    p0 = p_scale(model.model_ref)
    f,h = shape_factors(model,lb_v0,T0,z) #h normaly should be independent of temperature
    ps = p0*f/h
    return ps
end

#=
ideally we could perform SPUNG only providing x0, but i cant find the error here
function x0_sat_pure(model::SPUNG,T,z=SA[1.0])
    lb_v = lb_volume(model,z)
    f,h = shape_factors(model,lb_v,T,z)
    T0 = T/f
    @show T0
    vl0,vv0 = exp10.(x0_sat_pure(model.model_ref,T0,SA[1.0]))
    @show vl = vl0*h
    @show vv = vv0*h
    return [log10(vl),log10(vv)]
end
=#

#overloading sat_pure for SPUNG directly seems to be the way
function sat_pure(model::SPUNG,T::Real;v0=nothing)
    lb_v = lb_volume(model,SA[1.0])
    f,h = shape_factors(model,lb_v,T,SA[1.0])
    T0 = T/f
    psat0,vl0,vv0 = sat_pure(model.model_ref,T0;v0)
    p = pressure(model,vv0*h,T)
     return (p,vl0*h,vv0*h)
end

function sat_pure_p(model::SPUNG,p::Real)
    return naive_sat_pure_p(model,p)
end

function general_shape_factors(model::SPUNG,V,T,z=SA[1.0])
    n = sum(z)
    x = z * (1/n)
    RT = R̄*T
    b = lb_volume(model.shape_model,x)
    b0 = lb_volume(model.shape_ref,x)
    B = second_virial_coefficient(model.shape_model,T,x)
    B0 = second_virial_coefficient(model.shape_ref,T)
    #B = b-a/RT
    #a/RT = b-B
    #a = RT(b-B)
    a = RT*(b-B)
    a0 = RT*(b0-B0)
    tau = 1/(1-4(B/b))
    tau0 = 1(1-4(B0/b0))
    Tc = T_scale(model.shape_model,x)
    Tc0 = T_scale(model.shape_ref,x)

    f0 = (tau*Tc)/(tau0*Tc0)
    #f0 = tau/tau0
    #@show T/f0
     #T0 = Roots.find_zero(f0_f,T/f0)
    #@show T0
     #B0 = second_virial_coefficient(model.shape_ref,T0)
    #tau0 = 1-(B0/b0)
    #tau0 = 1-(B0/b0)
    #f = tau/tau0

    #a,b = cubic_ab(model.shape_model,T,x)
    #a0,b0 = cubic_ab(model.shape_ref,T,SA[1.0])
    h = b/b0
    #fh = n*a/a0
    #f = fh/h
    f = f0
    return f,h
end

#=
Tc = 0.26+2.1R
R = λ-1

=#

export SPUNG


text:
!isempty(ARGS) || error("No config supplied.")
isfile(ARGS[1]) || error("Cannot read '$(ARGS[1])'")
isabspath(ARGS[1]) || error("Please use an absolute path for the config.")
println("Config supplied: '$(ARGS[1])'")
config_file = ARGS[1]
include(config_file)

using MLDataUtils
using Random
using DelimitedFiles

function downsample_class(data, labels, target_class, other_class, target_percentage)
    target_indicies = shuffle!(findall(x -> x .== target_class, labels))
    mask = falses(length(labels))
    mask[labels .== other_class] .= true
    mask[target_indicies[1:ceil(Int, target_percentage * sum(labels .== other_class) / (1 - target_percentage))]] .= true
    return data[mask, :], labels[mask]
end

function process_file(input_file, output_file; num_versions=1)
    raw = readdlm(input_file, ',')
    num_attributes = length(findall(x -> occursin("@ATTRIBUTE", string(x)), raw[:, 1])) - 2
    id_column = findfirst(x -> occursin("@ATTRIBUTE 'id'", string(x)), raw[:, 1]) - 1
    label_column = findfirst(x -> occursin("@ATTRIBUTE 'outlier'", string(x)), raw[:, 1]) - 1
    data_start_row = findlast(x -> x == "@DATA", raw[:, 1]) + 1
    raw[:, label_column] = map(x -> x == "'yes'" ? :outlier : :inlier, raw[:, label_column])
    data, labels = raw[data_start_row:end, [i for i in 1:size(raw, 2) if i != id_column && i != label_column]], raw[data_start_row:end, label_column]
    data = hcat(data, labels)
    @assert size(data, 2) - 1 == num_attributes
    @assert size(data, 1) == length(labels)

    for i in 1:num_versions
        Random.seed!(i)
        outlier_percentage = sum(labels .== :outlier) / length(labels)
        resampling = outlier_percentage != TARGET_OUTLIER_PERCENTAGE || length(labels) > MAX_VALUES
        target_output_file = resampling ? "$(output_file[1:end-4])_r0$i.csv" : output_file
        @info "Generating '$target_output_file'."
        res_data, res_labels = copy(data), copy(labels)
        if outlier_percentage > TARGET_OUTLIER_PERCENTAGE
            @info "Downsampling outlier class (outlier_percentage = $(outlier_percentage))."
            res_data, res_labels = downsample_class(res_data, res_labels, :outlier, :inlier, TARGET_OUTLIER_PERCENTAGE)
        elseif outlier_percentage < TARGET_OUTLIER_PERCENTAGE
            @info "Downsampling inlier class (outlier_percentage = $(outlier_percentage))."
            res_data, res_labels = downsample_class(res_data, res_labels, :inlier, :outlier, 1 - TARGET_OUTLIER_PERCENTAGE)
        end
        if length(res_labels) > MAX_VALUES
            @info "Downsampling from $(length(res_labels)) to $MAX_VALUES observations."
            p = MAX_VALUES / length(res_labels)
            (res_data, res_labels), _ = stratifiedobs((res_data, res_labels), p=p, obsdim=1)
        end
        outlier_percentage = sum(res_labels .== :outlier) / length(res_labels)
        @info "Final outlier_percentage = $(outlier_percentage))."
        @assert size(res_data, 1) == length(res_labels)
        @assert abs(outlier_percentage - TARGET_OUTLIER_PERCENTAGE) < 0.01
        writedlm(target_output_file, res_data, ',')
    end
end

Random.seed!(0)
MAX_VALUES = 1000
TARGET_OUTLIER_PERCENTAGE = 0.05
target_versions_semantic = r"withoutdupl_norm_05_v0[1-3]"
target_versions_literature = r"withoutdupl_norm"
dataset_dir = normpath(joinpath(data_root, "input", "raw"))
output_path = normpath(joinpath(data_root, "input", "processed"))

mkpath(output_path)
@info "Saving processed files to $output_path."

for dataset_class in ["semantic", "literature"]
    for d in data_dirs[dataset_class]
        @info d
        outdir = joinpath(output_path, d)
        isdir(outdir) || mkpath(outdir)
        if dataset_class == "semantic"
            target_files = filter(x -> occursin(target_versions_semantic, x), readdir(joinpath(dataset_dir, dataset_class, d)))
            @assert length(target_files) == 3
            @info "[$(d)] Found $(length(target_files)) files."
            for f in target_files
               process_file(joinpath(dataset_dir, "semantic", d, f), joinpath(outdir, f[1:end-5] * ".csv"))
            end
        else
            target_files = filter(x -> occursin(target_versions_literature, x), readdir(joinpath(dataset_dir, dataset_class, d)))
            if (i = findfirst(x -> occursin("catremoved", x), target_files)) !== nothing
                target_file = target_files[i]
            else
                target_file = first(target_files)
            end
            @info "[$(d)] Resampling '$target_file'."
            process_file(joinpath(dataset_dir, dataset_class, d, target_file), joinpath(outdir, target_file[1:end-5] * ".csv"), num_versions=3)
        end
    end
end

max_stars_repo_name: kit-dbis/ocal-evaluation

max_issues_repo_name: kit-dbis/ocal-evaluation

max_forks_repo_name: kit-dbis/ocal-evaluation


text:
module TestDigitalCurrency
using AlphaVantage
using Test

TEST_SLEEP_TIME =  parse(Float64, get(ENV, "TEST_SLEEP_TIME", "15"))
MAX_TESTS = parse(Int64, get(ENV, "MAX_TESTS", "1"))

@testset "Digital Currencies" begin
    for f in (:digital_currency_daily, :digital_currency_weekly, :digital_currency_monthly, :digital_currency_intraday)
        @eval begin
            testname = string($f)
            @testset "$testname" begin
                @testset "AlphaVantageResponse" begin
                    data = $f("MSFT")
                    @test isa(data, AlphaVantageResponse)
                end
                sleep(TEST_SLEEP_TIME + 2*rand())
                @testset "JSON" begin
                    data = $f("BTC", datatype = "json")
                    @test typeof(data) === Dict{String,Any}
                    @test length(data) === 2
                end
                sleep(TEST_SLEEP_TIME + 2*rand())
                @testset "CSV" begin
                    data = $f("BTC", datatype = "csv")
                    @test typeof(data) === Tuple{Array{Any, 2}, Array{AbstractString, 2}}
                    @test length(data) === 2
                end
            end
        end
        sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
    end

    @testset "crypto_rating" begin
        data = crypto_rating("BTC")
        @test typeof(data) === Dict{String, Any}
        @test length(data) === 1
        @test data["Crypto Rating (FCAS)"]["1. symbol"] === "BTC"
        sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
    end

end

end # module


text:
@testset "code generation (relying on a device)" begin

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

@testset "LLVM" begin

@testset "stripping invariant.load" begin
    function kernel(ptr, x)
        i = CUDAnative.threadIdx_x()
        @inbounds unsafe_store!(ptr, x[i], 1)
        return
    end

    arr = CuArray(zeros(Float64))
    ptr = pointer(arr)

    @cuda kernel(ptr, (1., 2., ))
    @test Array(arr)[] == 1.
end

@testset "stripping const TBAA" begin
    # this one is particularly nasty because it occurs in a nested function

    _a = rand(Int, 2, 1)
    b = ((1,9999),(1,9999))

    out = CuArray(zeros(Int, 2,1))
    a = Tuple(_a)

    function kernel(out, a, b)
        i = threadIdx().x
        blockIdx().x
        @inbounds out[i,1] = a[i] + b[i][1]
        return
    end

    @cuda threads=2 kernel(out, a, b)
    @test Array(out) == (_a .+ 1)
end

@testset "ptxas-compatible control flow" begin
    @noinline function throw_some()
        throw(42)
        return
    end

    @inbounds function kernel(input, output, n)
        i = threadIdx().x

        temp = @cuStaticSharedMem(Int, 1)
        if i == 1
            1 <= n || throw_some()
            temp[1] = input
        end
        sync_threads()

        1 <= n || throw_some()
        unsafe_store!(output, temp[1], i)

        return
    end

    function gpu(input)
        output = CuArray(zeros(eltype(input), 2))
        ptr = pointer(output)
        ptr = reinterpret(Ptr{eltype(input)}, ptr)

        @cuda threads=2 kernel(input, ptr, 99)

        return Array(output)
    end

    function cpu(input)
        output = zeros(eltype(input), 2)

        for j in 1:2
            @inbounds output[j] = input
        end

        return output
    end

    input = rand(1:100)
    @test cpu(input) == gpu(input)
end

end

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

@testset "SASS" begin

@testset "basic reflection" begin
    valid_kernel() = return
    invalid_kernel() = 1

    @test CUDAnative.code_sass(devnull, valid_kernel, Tuple{}) == nothing
    @test_throws KernelError CUDAnative.code_sass(devnull, invalid_kernel, Tuple{})
end

@testset "function name mangling" begin
    @eval @noinline $(Symbol("dummy_^"))(x) = x

    @eval kernel_341(ptr) = (@inbounds unsafe_store!(ptr, $(Symbol("dummy_^"))(unsafe_load(ptr))); nothing)

    CUDAnative.code_sass(devnull, kernel_341, Tuple{Ptr{Int}})
end

end

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

end


text:
module ResultFrames

   include( "result_frames.jl" )

end


text:
using Test
using Liblet

test_transition_totalorder() = Transition("a", "b", "c") > Transition("a", "b", "b")

test_transition_hash() = begin
    S = Dict(
            Transition("a", "b", "c") => 1,
            Transition("a", "b", "c") => 2
    )
    1 == length(S)
end

test_transition_str() = "frm-label->to" == string(Transition("frm", "label", "to"))

test_transition_lt() = "{frm}-label->{to}" == string(Transition(Set(["frm"]), "label", Set(["to"])))

test_transition_set() = "{frm}-label->{to}" == string(Transition(Set(["frm"]), "label", Set(["to"])))

test_transition_tupleofitems() = "(A -> •B)-label->(C -> •D)" == string(Transition([Item("A", ["B"])], "label", [Item("C", ["D"])]))

test_transition_wrong_label1() = Transition("frm", 1, "to")

test_transition_wrong_label2() = Transition("frm", "", "to")

test_transition_wrong_from1() = Transition(1, "label", "to")

test_transition_wrong_from2() = Transition("", "label", "to")

test_transition_wrong_from3() = Transition(Set(), "label", "to")

test_transition_wrong_from4() = Transition(Set([""]), "label", "to")

test_transition_wrong_to1() = Transition("frm", "label", 1)

test_transition_wrong_to2() = Transition("frm", "label", "")

test_transition_wrong_to3() = Transition("frm", "label", Set())

test_transition_wrong_to4() = Transition("frm", "label", Set([""]))

# Tests

function runtransitiontests()
    @testset "transition_tests" begin
        @test test_transition_totalorder()
        #@test test_transition_hash()
        #@test test_transition_str()
        #@test test_transition_lt()
        #@test test_transition_set()
        #@test_throws ArgumentError test_transition_tupleofitems()
        @test_throws MethodError test_transition_wrong_label1()
        @test_throws ArgumentError test_transition_wrong_label2()
        @test_throws MethodError test_transition_wrong_from1()
        @test_throws ArgumentError test_transition_wrong_from2()
        @test_throws ArgumentError test_transition_wrong_from3()
        @test_throws ArgumentError test_transition_wrong_from4()
        @test_throws MethodError test_transition_wrong_to1()
        @test_throws ArgumentError test_transition_wrong_to2()
        @test_throws ArgumentError test_transition_wrong_to3()
        @test_throws ArgumentError test_transition_wrong_to4()
    end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Libdl, GMP_jll

@testset "GMP_jll" begin
    vn = VersionNumber(unsafe_string(unsafe_load(cglobal((:__gmp_version, libgmp), Ptr{Cchar}))))
    @test vn == v"6.2.0"
end

max_stars_repo_name: udesou/julia

max_issues_repo_name: udesou/julia

max_forks_repo_name: udesou/julia


text:
using ArgParse
using JSON
using DataStructures
using NCDatasets

function runOneCmd(cmd)
    println(">> ", string(cmd))
    run(cmd)
end

function pleaseRun(cmd)
    if isa(cmd, Array)
        for i = 1:length(cmd)
            runOneCmd(cmd[i])
        end
    else
        runOneCmd(cmd)
    end
end

println("""
This file produce an netcdf file with only `time` and `time_bound` variables
for users to append variables to other files.

Notice: This program uses no-leap calendar (365 days fixed)

""")

s = ArgParseSettings()
@add_arg_table s begin

    "--output"
        help = "The file you want to produce."
        arg_type = String
        required = true

    "--years"
        help = "How many years?"
        arg_type = Int64
        default = 1

end

parsed = DataStructures.OrderedDict(parse_args(ARGS, s))

JSON.print(parsed, 4)

dom = [31.0, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
sum_dom = sum(dom)
sum_dom == 365 || throw(ErrorException("Sum of dom is $(sum_dom) rather than 365."))

_t    = zeros(Float64, length(dom))
_bnds = zeros(Float64, 2, length(dom))

for m=1:length(dom)
    #bnds[m, 1] = beg of month  m
    #bnds[m, 2] = end of month  m
    if m==1
        _bnds[1, m] = 0.0
    else
        _bnds[1, m] = _bnds[2, m-1]
    end

    _bnds[2, m] = _bnds[1, m] + dom[m]

    _t[m] = (_bnds[1, m] + _bnds[2, m]) / 2.0
end

t    = zeros(Float64, 12*parsed["years"])
bnds = zeros(Float64, 2, 12*parsed["years"])

for y = 1:parsed["years"]
    i_offset = (y-1)*12
    t_offset = (y-1)*sum_dom
    t[i_offset+1:i_offset+12]       .+= _t    .+ t_offset
    bnds[:, i_offset+1:i_offset+12] .+= _bnds .+ t_offset
end

Dataset(parsed["output"], "c") do ds

    defDim(ds, "time", Inf)
    defDim(ds, "d2", 2)

    defVar(ds, "time", t, ("time", ), ; attrib = Dict(
        "long_name" => "time",
        "bounds"    => "time_bound",
        "calendar"  =>  "noleap",
        "units"     => "days since 0001-01-01 00:00:00",
    ))

    defVar(ds, "time_bound", bnds, ("d2", "time"), ; attrib = Dict(
        "long_name" => "boundaries for time-averaging interval",
        "units"     => "days since 0001-01-01 00:00:00",
    ))

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using Term.Consoles
import Term: tprint, tprintln

println("\nTesting tprint, stdout temporarily disabled")

@suppress_out begin
    @testset "\e[34mTPRINT" begin
        @test_nothrow tprint("string")
        @test_nothrow tprint("{red}adasd{/red}")
        @test_nothrow tprint("{blue on_green}adasda")
        @test_nothrow tprint("{red}dadas{green}insdai{/green}outssdrse{blue}fsfsf{/blue}")
        @test_nothrow tprint(Panel("test"))
        @test_nothrow tprint(TextBox("test"))
        @test_nothrow tprint(1)
        @test_nothrow tprint(:x)
        @test_nothrow tprint(1, Panel("test"), "test")
        @test_nothrow tprintln("test")
        @test_nothrow tprintln(Panel("test"))
        @test_nothrow tprintln(1, Panel("test"), "test")

        @test_nothrow tprint(stdout, "string")
        @test_nothrow tprint(stdout, "{red}adasd{/red}")
        @test_nothrow tprint(stdout, "{blue on_green}adasda")
        @test_nothrow tprint(
            stdout,
            "{red}dadas{green}insdai{/green}outssdrse{blue}fsfsf{/blue}",
        )
        @test_nothrow tprint(stdout, Panel("test"))
        @test_nothrow tprint(stdout, TextBox("test"))
        @test_nothrow tprint(stdout, 1)
        @test_nothrow tprint(stdout, :x)
        @test_nothrow tprint(stdout, 1, Panel("test"), "test")
        @test_nothrow tprintln(stdout, "test")
        @test_nothrow tprintln(stdout, Panel("test"))
        @test_nothrow tprintln(stdout, 1, Panel("test"), "test")
    end
end

@testset "\e[34mCONSOLE" begin
    @test console_height() == displaysize(stdout)[1]

    @test console_width() == displaysize(stdout)[2]

    @test sprint(cursor_position) == "\e[6n"
    @test sprint(up) == "\e[1A"
    @test sprint(beginning_previous_line) == "\e[F"
    @test sprint(down) == "\e[1B"
    @test sprint(clear) == "\e[2J"
    @test sprint(hide_cursor) == "\e[?25l"
    @test sprint(show_cursor) == "\e[?25h"
    @test sprint(line) == "\n"
    @test sprint(prev_line) == "\e[1F"
    @test sprint(next_line) == "\e[1E"
    @test sprint(erase_line) == "\e[2K"
    @test sprint(cleartoend) == "\e[0J"
    @test sprint(move_to_line) == "\e[1;1H"
    @test sprint(change_scroll_region) == "\e[1;1r\e[1B"
    @test sprint(savecursor) == "\e[s"
    @test sprint(restorecursor) == "\e[u"
end


text:
# dynamics tests
tf = 10.0
t = LinRange(0.,tf,Nt)
sim2 = Sim(sim;γ = 0.0,tf=tf,t=t)

#make initial state
ϕi = sol[end]
@pack! sim2 = ϕi

sol2,err = testsim(sim2)
@test err == false

psitest = xspace(sol2[end],sim2)
npeak = abs2.(psitest[256])
@test isapprox(g*npeak,μ,rtol=1e-2)


text:
export randnCPR, randnP, randnC

function randnC(args...;kwargs...)
    return randn(args...;kwargs...) .+ im * randn(args...;kwargs...)
end

function randnCPR(args...;kwargs...)
    return randnP(args...;kwargs...) .+ im * randn(args...;kwargs...)
end

function randnP(args...;kwargs...)
    return abs.(randn(args...;kwargs...))
end


text:
#in module Units

meter = Physical(1, m=1)
kilogram = Physical(1, kg=1)
second = Physical(1, s=1)
ampere = Physical(1, A=1)
kelvin = Physical(1, K=1)
mole = Physical(1, mol=1)
candela = Physical(1, cd=1)
hertz = 1/second
newton = kilogram*meter/second^2
pascal = newton/meter^2
joule = newton*meter
watt = joule/second
coulomb = ampere * second
volt = watt/ampere
farad = coulomb/volt
ohm = volt/ampere
weber = volt*second
tesla = weber/meter^2
henry = weber/ampere
lumen = candela
lux = lumen/meter^2

NA = 6.02e23/mole
planck = 6.62607004081e-34joule/second


text:
module tools

export unfolding, folding, mode_n_mult, tten, hosvd
export get_xi, get_rhs_sin, get_rhs_norm, fullten
export rhs_sin, rhs_norm

export aca, cur, cur_fun
export unfolding_fun, aca_fun, approx_aca

using LinearAlgebra

#### New interface for right-hand side
# This can be used more like functions.
# For any other function on the right hand side,
# just make a non-abstract subtype of rhs and
# implement the desired function as a callable.

# The functions (r::T)(i,j,k) where T <: rhs is some subtype of rhs
# make the objects callable, i.e. one can write
# > r = rhs_sin(100)
# > r(1,35,83)
# this makes the code very convenient.  The Base.getindex
# implementation for rhs is for a more natural indexign via
# > r[1,35,83]
# which returns the same value as the call above.

abstract type rhs end

struct rhs_sin <: rhs
    n::Int
    den::Float64
    rhs_sin(n::Int) = new(n, 1/(n+1))
end

function (r::rhs_sin)(i::Int,j::Int,k::Int)
    # multiplication with reciprocal denominator is much faster than
    # division
    sin((i+j+k)*r.den)
end

struct rhs_norm <: rhs
    n::Int
    den::Float64
    rhs_norm(n::Int) = new(n, 1/(n+1))
end

function (r::rhs_norm)(i::Int,j::Int,k::Int)
    sqrt( (i*r.den)^2 + (j*r.den)^2 + (k*r.den)^2)
end

function Base.getindex(r::rhs, i::Int, j::Int, k::Int)
    # cheap boundscheck; note this does not use the Julia convetnions
    # of bounds checking with several nested function calls
    # check may be removed for better performance
    correct = 1<=i<=r.n && 1<=j<=r.n && 1<=k<=r.n
    if ~correct
        throw(BoundsError(r,(i,j,k)))
    end
    return r(i,j,k)
end

# alternative to getindex for callables This makes arrayindexing with
# Fortran syntax possible, i.e. a(i,j,k)
#### use with care!!
#### uncomment only if you know what you do
# function (a::Array)(i::Int...)
#     a[i...]
# end

"""
    unfolding_fun(a, n, dims) -> a_u(i,j)

Returns a function for evaluating the unfolding (matricisation,
flattening) of a tensor that is given only via a function.

"""
function unfolding_fun_expr(a, n, dims)
    # this is just for getting the first index after permutation, it
    # is assumed that the rank of the tensors is not huge
    # This is a bit easy here, because I assume tensors of fixed order
    # 3 only. Then products just vanish.
    temp = dims[1:end .!= n]
    ind = Tuple(temp)

    # here one-based indexing is a PITA. Maybe this can be done
    # easier, but at least this works. It has been tested with the
    # array based version.
    if n == 1
        ret = :( (i,j) ->
                 # $a(i, ind2sub($ind,j)...)
                 $a(i, CartesianIndices($ind)[j][1], CartesianIndices($ind)[j][2])
                 # if rem(j,$ind)==0
                 # a(i,$ind,div(j,$ind))
                 # else
                 # a(i,rem(j,$ind),1+div(j,$ind))
                 # end
                 )
    elseif n == 2
        ret = :( (i,j) ->
                 $a(CartesianIndices($ind)[j][1], i, CartesianIndices($ind)[j][2])
                 # if rem(j,$ind)==0
                 # a($ind,i,div(j,$ind))
                 # else
                 # a(rem(j,$ind),i,1+div(j,$ind))
                 # end
                 )
    elseif n==3
        ret = :( (i,j) ->
                 $a(CartesianIndices($ind)[j][1], CartesianIndices($ind)[j][2], i)
                 # if rem(j,$ind)==0
                 # a($ind,div(j,$ind),i)
                 # else
                 # a(rem(j,$ind),1+div(j,$ind),i)
                 # end
                 )
    end
    # show(ret)                   # for seeing what substitutions are made
    return ret
end

function unfolding_fun(a, n, dims)
    return eval(unfolding_fun_expr(a, n, dims))
end

"""
    get_xi(n::Int)
returns a vector xi of length n with entries
```julia
xi[i] = 1/(n+1)
```
"""
function get_xi(n::Int)
    collect(1:n)/(n+1)
end

@inline function ξ(i, n)
    i/(n+1)
end

function get_rhs_sin(xi::Array{Float64,1})
    n = length(xi)
    rhs = Array{Float64,3}(undef, (n,n,n))
    @inbounds for i in 1:n
        for j in 1:n
            @simd for k in 1:n
                rhs[k,j,i] = sin(xi[k]+xi[j]+xi[i])
            end
        end
    end
    return rhs
end

function get_rhs_sin(n::Int)
    rhs = Array{Float64,3}(undef, (n,n,n))
    @inbounds for i in 1:n
        for j in 1:n
            @simd for k in 1:n
                rhs[k,j,i] = sin(ξ(k,n)+ξ(j,n)+ξ(i,n))
            end
        end
    end
    return rhs
end

function get_rhs_sin_opt(n::Int)
    # this function retruns different values
    # this is due to machine arithmetic and
    # the precomputation of the denominator
    rhs = Array{Float64,3}(undef, (n,n,n))
    den::Float64 = 1.0/(n+1)
    @inbounds for i in 1:n
        for j in 1:n
            @simd for k in 1:n
                rhs[k,j,i] = sin((k+j+i)*den)
            end
        end
    end
    return rhs
end

function get_rhs_norm(n::Int)
    rhs = Array{Float64,3}(undef, (n,n,n))
    @inbounds for i in 1:n
        for j in 1:n
            @simd for k in 1:n
                rhs[k, j, i] = sqrt(ξ(k,n)*ξ(k,n)+ξ(j,n)*ξ(j,n)+ξ(i,n)*ξ(i,n))
            end
        end
    end
    return rhs
end

function get_rhs_norm(xi::Array{Float64,1})
    n = length(xi)
    rhs = Array{Float64,3}(undef, (n,n,n))
    @inbounds for i in 1:n
        for j in 1:n
            @simd for k in 1:n
                rhs[k,j,i] = sqrt(xi[k]*xi[k] + xi[j]*xi[j] + xi[i]*xi[i])
            end
        end
    end
    return rhs
end

function approx_sin(n)
    rhs = get_rhs_sin(n)
    hosvd(rhs, 1e-4)
end

function approx_norm(n)
    rhs = get_rhs_norm(n)
    hosvd(rhs, 1e-4)
end

"""
n-mode unfolding or matricisation of a tensor a
"""
function unfolding(a, n)
    dims = size(a)
    indices = setdiff(1:ndims(a),n)
    reshape(permutedims(a,[n;indices]),dims[n],prod(dims[indices]))
end

"""
fold matrix into tensor by mode n
reverse of unfolding
"""
function folding(m, n, dims)
    indices = setdiff(1:length(dims), n)
    a = reshape(m, (dims[n],dims[indices]...))
    permutedims(a,invperm([n;indices]))
end

"""
mode-n multiplication of tensor a with matrix m

# inputs
- a :  Array
- n :  mode for multiplication
- m :  matrix that gets multiplied with a in mode n
"""
function mode_n_mult(a, n, m)
    dims = size(a)
    # new_dims = [dims...]
    new_dims = collect(dims)
    new_dims[n] = size(m,1)
    b = m*unfolding(a, n)
    folding(b,n,new_dims)
end

"""
Type for Tucker tensor
core:   Array, this is the core tensor of size (r₁, r₂, r₃, ...)
frames: Array of Arrays, these are the mode frames
"""
struct tten{T<:AbstractFloat}
    core::Array{T}
    frames::Array{Array{T,2},1}

end

"""
    hosvd(a, ϵ)

Computes the Higher Order Singular Value Decomposition

The computation is done via sequentially truncated HOSVD, this is much faster,
as the tensor gets compressed early.
Returns a `tten` object. The tten is a tucker tensor consisting of a
core tensor and an array of mode frames.

input:
- a: full tensor
- ϵ: relative error bound

output:
- `tten` object representing the Tucker approximation
"""
function hosvd(a, ϵ::T=1e-4, sv::Bool=false) where {T<:AbstractFloat}
    d = ndims(a)
    # dims = [size(a)...]
    dims = collect(size(a))
    ϵ = ϵ/sqrt(d) * norm(a[:],2) # divide relative error equal to all dimensions
    ϵ² = ϵ*ϵ                     # use the sqare for cheaper comparison

    frames = Array{Array{Float64,2},1}(undef, d)
    if sv
        svs = Array{Array{Float64,1},1}(undef, d)
    end
    for i in 1:d
        U,S,V = svd(unfolding(a,i))
        rᵢ = 1
        temp = 0.0
        for j in length(S):-1:1
            temp += S[j]*S[j]
            if temp > ϵ²
                # take the one satisfying the error bound, but stay
                # within the bounds => check whether on bound or not
                rᵢ = j==length(S) ? j : (j)
                break
            end
        end
        frames[i] = U[:, 1:rᵢ]
        #a = mode_n_mult(a,i,frames[i]')
        # use a more compact representation of the remaining matrix
        # the implementation of mode_n_mult copys the data around more than once
        # this works as the remainder S*V is just the same as
        # U'*A = U'*U*A*V', as U is orthogonal
        dims[i] = rᵢ
        a = folding(diagm(0=>S[1:rᵢ])*V[:,1:rᵢ]',i,dims);
        if sv
            svs[i] = S
        end
        # garbage collection makes life easier and code faster
        GC.gc()
    end
    if sv
        return tten(a,frames), svs
    else
        return tten(a,frames)
    end
end

# the hosvd2 functions are worse because they reduce the core at the
# end when done during the first loop, simultaneous to svd's, it
# reduces the size and thus cost much
function hosvd2(a, ranks::Vector{T}) where {T<:Int}
    d = ndims(a)
    @assert(d==length(ranks), "dimension mismatch")

    frames = Array{Array{Float64,2},1}(undef, d)
    for i in 1:d
        U = svd(unfolding(a,i))
        frames[i] = U[:, 1:ranks[i]]
    end

    for i in 1:d
        a = mode_n_mult(a,i,frames[i]')
    end
    return tten(a, frames)

end

function hosvd2(a, ranks::Tuple{Vararg{T}}) where {T<:Int}
    # rank_arr = [ranks...]
    rank_arr = collect(ranks)
    hosvd(a, rank_arr)
end

function hosvd2(a, tol::T) where {T<:AbstractFloat}
    d = ndims(a)
    tol = tol/sqrt(d) * norm(a[:],2)

    frames = Array{Array{Float64,2},1}(undef, d)
    for i in 1:d
        U,S = svd(unfolding(a,i))
        r_i = length(S)
        for j in 1:length(S)
            # rework this
            # this is not exactly what we want
            # better is ||S[j:end]|| < tol or even
            # sum(S[j:end].^2) < tol^2
            if S[j] < tol
                r_i = j
                break
            end
        end
        frames[i] = U[:, 1:r_i]
        gc()
    end

    for i in 1:d
        a = mode_n_mult(a,i,frames[i]')
        gc()
    end
    return tten(a,frames)
end

function fullten(a::tten{T}) where {T<:AbstractFloat}
    res = a.core
    for i in 1:length(a.frames)
        res = mode_n_mult(res, i, a.frames[i])
        GC.gc()
    end
    return res
end

# """
#     approx_aca(a, dims, ϵ)

# Compute an approximation of a in Tucker format. dim is an array
# holding the size of a.
# """
# function approx_aca(a, dims, ϵ::Float64=1e-4)
#     d = length(dims)
#     frames =  Array{Array{Float64,2},1}(d)

#     a_u = unfolding_fun(a, 1, dims)
#     dims_u = [dims[1], prod(dims[2:end]) ]
#     I,J = aca_fun(a_u, dims_u, ϵ)
#     C,U,R = cur_fun(a_u, dims_u, I, J)
#     frames[1] = C
#     dims[1] = length(I)
#     core = folding(U*R, 1, dims)

#     for i = 2:3
#         core_u = unfolding(core, i)
#         I,J = aca(core_u, ϵ)
#         C, U, R = cur(core_u, I, J)
#         frames[i] = C
#         dims[i] = length(I)
#         core = folding(U*R, i, dims)
#         gc()
#     end

#     return tten(core, frames)
# end

"""
    aca(a, ϵ)

Computes adaptive cross approximation of matrix a with relative accuracy ϵ.

This implementation uses partial pivoting, this implies the array has
to be evaluated along a row or column. For the computation not the
whole array is needed.

"""
function aca(a::Array{Float64,2},ϵ::Float64=1e-4)
    Rk = a
    maxrk = min(size(a)...)
    I = Int[]
    J = Int[]
    I_good = Int[]
    J_good = Int[]
    k = 1
    i = 1
    # i = 70
    normk2 = 0.0
    us = []
    vs = []
    ϵ = ϵ*ϵ
    # use ↓ this crappy construct to simulate a do-while-loop
    while true
        # maxval is not really used here but comes for free in mymax
        # function
        # maxval,j = mymax(abs.(Rk[i,:]),J)
        # use better pivoting compared to pseudo-alg from lecture
        i,j = find_pivot(Rk, i, I, J)
        # reorder execution compared to algorith from lecture
        # now this is not really important anymore, but it does not
        # matter where this is done. For an early (not working)
        # version this was important
        append!(I, i)
        append!(J, j)

        δ = Rk[i,j]
        # improve the error here. ϵ is maybe not optimal, as it is
        # just the sqared error bound for the overall
        # approximation. This just works for the tested cases.
        if abs(δ) < ϵ
            if length(I) == maxrk - 1
                # if this happens, no proper LRA has been calculated
                @warn("Return with maximal rank => no Low Rank Approximation")
                return I_good,J_good
            end

            # try again, if still bad exit note: this is a trade off
            # between accuracy and speed.
            # It may happen, that there are still good candidates, but
            # those cannot be reached with the search used for the
            # pivot elements. E.g. for block diagonal matrices this
            # easily happens.
            i,j = find_pivot(Rk, i, I, J)
            δ = Rk[i,j]
            if abs(δ) < ϵ
                @warn("Could not find good pivot anymore.")
                return I_good,J_good
            end
        else
            uk = Rk[:,j]
            vk = Rk[i,:]'/δ
            Rk = Rk - uk*vk
            k = k + 1
            # use only the indices that contribute
            append!(I_good, i)
            append!(J_good, j)
            push!(us, uk)
            push!(vs, vk)

            # stopping criterion
            # calculate the parts seperately
            # first interactions of uk, vk with old u and v then the
            # sqared norm, at the end add everything up
            temp = 0
            # k-2 due to early incrementation of k
            for l in 1:k-2
                temp += uk'*us[l] * vs[l]*vk'
            end
            #                         ↓ sqared norms
            normk2 = normk2 + 2temp + (uk'*uk)*(vk*vk')

            # check for stopping criterion
            # also if it is not fulfilled check for high rank, if rank
            # gets too high just return silently and don't complain
            ### maybe it is a good idea to add some special warning or
            ### throw an error if rank exceeds the bounds, after all
            ### this is for _low rank_ approximation
            if uk'*uk * vk*vk' <= ϵ*normk2 || k==maxrk
                return I_good,J_good
            end
        end
        # if abs(δ)<2eps() does not make a
        # difference, if uk and vk are not updated, the error estimate
        # remains the same. So this can be checked in the else
        # statement

        # needed when directly implementing pseudo-algo from lecture
        # that is not the best choice
        # maxval, i = mymax(abs.(uk), I)
    end
end

function find_pivot(a, i, I, J, n::Int=3)
    j = 1
    for l = 1:n
        maxval, j = mymax(abs.(a[i,:]), J)
        maxval, i = mymax(abs.(a[:,j]), I)
    end
    return i,j
end

"""
    mymax(v, I) -> max,i

Computes the maximal value and the index of v, where i ∉ I. I is some
indexset, every element of I is a valid index for v. (The
implementation assumes that at least the minimal index of I is valid.)
"""
function mymax(v, I)
    indices = setdiff(1:length(v),I)
    ind = min(indices...)
    maxval = v[ind]
    for i ∈ indices
        if v[i] > maxval
            maxval = v[i]
            ind = i
        end
    end
    return maxval, ind
end

"""
    update_rk(a, I, J) -> r

takes a function a(i,j), arrays if indices I, J and returns function
calculating the residual a - u*s*v'. Where u are columns J of a, v are
rows I of a, and s is inverse of a(I,J).
"""
function update_rk(a, I, J)
    # check this again! Not sure this works
    #
    # Idea: compute s = inv(a(I,J))
    # then compute a - u*s*v'
    # for the last computation use only the rows and columns that are
    # required
    # this function returns a closure, if not used to it looks strange
    # at first

    if length(I) != length(J)
        error("I and J have to be of equal size")
        return
    end

    s = zeros(length(I),length(J))
    for (j,jj) ∈ enumerate(J)
        for (i,ii) ∈ enumerate(I)
            s[i,j] = Base.invokelatest(a, ii,jj)
        end
    end
    s = inv(s)

    function u(ii)
        # check sizes!
        ret = zeros(length(J))
        for (j,jj) ∈ enumerate(J)
            ret[j] = Base.invokelatest(a, ii,jj)
        end
        return ret
    end

    function v(jj)
        ret = zeros(length(I))
        for (i,ii) ∈ enumerate(I)
            ret[i] = Base.invokelatest(a, ii,jj)
        end
        return ret
    end

    # define the closure
    function rk(ii,jj)
        return Base.invokelatest(a, ii,jj) - u(ii)'*s*v(jj)
    end

    return rk
end

"""
   find_pivot_fun(a, dims, i, I, J)

retunr indices i,j that maximise a(i,j) using a heuristic (search
along axes of a). dims array contains the size of a, i.e.
a ∈ ℝ^(dims[1],dims[2]).

"""
function find_pivot_fun(a, dims, i, I, J)
    j = 1
    for n = 1:3
        a1 = jj -> abs(a(i,jj))
        # maxval, j = mymax_fun( jj-> abs(a(i,jj)), dims[2], J)
        maxval, j = mymax_fun(a1, dims[2], J)
        a2 = ii -> abs(a(ii,j))
        # maxval, i = mymax_fun( ii-> abs(a(ii,j)), dims[1], I)
        maxval, i = mymax_fun(a2, dims[1], I)
    end
    return i,j
end

"""
    mymax_fun(v, n, I)

return maximum value of function v and index ∈ {1...n}∖I
"""
function mymax_fun(v, n, I)
    indices = setdiff(1:n, I)
    ind = min(indices...)
    # maxval = v(ind)
    maxval = Base.invokelatest(v, ind)
    for i ∈ indices
        # if v(i) > maxval
        if Base.invokelatest(v, i) > maxval
            # maxval = v(i)
            maxval = Base.invokelatest(v, i)
            ind = i
        end
    end
    return maxval, ind
end

function aca_fun(a, dims, ϵ::Float64=1e-4)
    #R = Array{Any,1}([a])
    # R = []
    # push!(R, a)
    # Rk = a
    maxrk = min(dims...)
    I = Int[]
    J = Int[]
    I_good = Int[]
    J_good = Int[]
    k = 1
    i = 1
    normk2 = 0.0
    us = []
    vs = []
    Rk = []
    push!(Rk, a)
    ϵ = ϵ*ϵ

    while true
        i,j = find_pivot_fun(Rk[k], dims, i, I, J)
        append!(I, i)
        append!(J, j)

        # δ = R[k](i,j)
        δ = Base.invokelatest(Rk[k], i,j)
        if abs(δ) < ϵ
            # if abs(δ) == 0
            #     return sort!(I),sort!(J)
            # end
            if length(I) == maxrk-1
                @warn("Return with maximal rank")
                return I_good, J_good
            end

            i,j = find_pivot_fun(Rk[k], dims, i, I, J)
            δ = Rk[k](i,j)
            if abs(δ) < ϵ
                @warn("Could not find good pivot aymore.")
                return I_good, J_good
            end
        else
            # uk = (ii)->R[k](ii,j)
            # vk = (jj)->R[k](i,jj)/δ
            # uk = (ii)->Rk[k](ii,j)
            # vk = (jj)->Rk[k](i,jj)/δ
            uk = (ii)->Base.invokelatest(Rk[k], ii,j)
            vk = (jj)->Base.invokelatest(Rk[k],i,jj)/δ
            push!(us, uk)
            push!(vs, vk)

            # FIX: calculate the whole function in every step
            # this costs a loop of length k in each step
            # push!(R, (ii,jj)->R[k](ii,jj)-uk(ii)*vk(jj) )
            # Rk = update_rk(Rk, uk, vk)
            # Rk = update_rk(a, us, vs)

            # update later, otherwise uk and vk get zero!!!!!
            # Rk = update_rk(a, I, J)
            # k = k+1

            temp = 0.0
            for l in 1:k-2
                # temp += prod(map( ii->uk(ii)*us[l](ii), 1:dims[1])) *
                #     prod(map( jj->ns[l](jj)*vk(jj), 1:dims[2]))
                # use the handy |> piping syntax
                temp += sum(1:dims[1] .|> ii->uk(ii)*us[l](ii)) *
                    sum(1:dims[2] .|> jj->vs[l](jj)*vk(jj))
            end

            # calculate product of sqared norms
            # sqnormprod = prod(map( ii->uk(ii)*uk(ii), 1:dims[1])) * prod(map( jj->vk(jj)*vk(jj), 1:dims[2]))
            sqnormprod = sum(1:dims[1] .|> ii->uk(ii)*uk(ii)) * sum(1:dims[2] .|> jj->vk(jj)*vk(jj))
            normk2 = normk2 + 2temp + sqnormprod

            if sqnormprod <= ϵ*normk2 || k==maxrk
                return I_good, J_good
            end

            # use only the indices that contribute to the solution
            append!(I_good, i)
            append!(J_good, j)

            # update here, otherwise uk and vk would be zero
            # k already increased
            k = k+1
            push!(Rk, update_rk(a, I_good, J_good))
        end

    end
end

"""
    cur_fun(a, dims, I, J) -> c, u, r

returns CUR decomposition of `a` given as function a(i,j). dims array
defines sizes of a, i.e. a ∈ ℝ^(dims[1],dims[2]). I, J are the row and
columns indices as returned from aca_fun.

When using aca with functions for huge sizes be careful. This function creates **full matrices**!

"""
function cur_fun(a, dims, I, J)
    m = dims[1]
    n = dims[2]
    c = zeros(m, length(J))
    u = zeros(length(I),length(J))
    r = zeros(length(I), n)

    # set c to coloumns of a
    for (j,jj) ∈ enumerate(J)
        for i in 1:m
            c[i,j] = Base.invokelatest(a, i,jj)
        end
    end

    # set r to rows of a
    for j in 1:n
        for (i,ii) ∈ enumerate(I)
            r[i,j] = Base.invokelatest(a, ii,j)
        end
    end

    # first get a[I,J], then invert it
    for (j,jj) ∈ enumerate(J)
        for (i,ii) ∈ enumerate(I)
            u[i,j] = Base.invokelatest(a, ii,jj)
        end
    end
    u = inv(u)

    return c, u, r
end

"""
   cur(a, i, j) -> C,U,R

Compute the CUR approximation given the original matrix a and the
indiex sets i and j.
"""
function cur(a, i, j)
    return a[:,j], inv(a[i,j]), a[i,:]
end

"""
    approx_aca(a, dims, ϵ)

Compute an approximation of a in Tucker format. dim is an array
holding the size of a.
"""
function approx_aca(a, dims, ϵ::Float64=1e-4)
    d = length(dims)
    frames =  Array{Array{Float64,2},1}(undef, d)

    a_u = unfolding_fun(a, 1, dims)
    dims_u = [dims[1], prod(dims[2:end]) ]
    I,J = aca_fun(a_u, dims_u, ϵ)
    C,U,R = cur_fun(a_u, dims_u, I, J)
    frames[1] = C
    dims[1] = length(I)
    core = folding(U*R, 1, dims)

    for i = 2:3
        core_u = unfolding(core, i)
        I,J = aca(core_u, ϵ)
        C, U, R = cur(core_u, I, J)
        frames[i] = C
        dims[i] = length(I)
        core = folding(U*R, i, dims)
        gc()
    end

    return tten(core, frames)
end

end


text:
using {PKGNAME}
using Test

# write tests here


text:
################################################################################
#  Copyright 2020, Marta Vanin, Tom Van Acker                                  #
################################################################################
# PowerModelsDistributionStateEstimation.jl                                    #
# An extention package of PowerModelsDistribution.jl for Static Power System   #
# State Estimation.                                                            #
################################################################################

## Exponential
# functions
heslogpdf(d::_DST.Exponential{T}, x::Real) where T<:Real = 0

## Weibull
# functions
function heslogpdf(d::_DST.Weibull{T}, x::Real) where T<:Real
    if _DST.insupport(_DST.Weibull, x)
        θ, α = _DST.params(d)
        - (α - 1) / x^2 - (α - 1) * α * x^(α - 2) / θ^(α)
    else
        zero(T)
    end
end

## Normal
# functions
heslogpdf(d::_DST.Normal{T}, x::Real) where T<:Real = -1/_DST.std(d)^2

## LogNormal
# functions
function heslogpdf(d::_DST.LogNormal{T}, x::Real) where T<:Real
    if _DST.insupport(_DST.LogNormal, x)
        μ, σ = _DST.params(d)
        ( log(x - μ) - 1 ) / ( 2σ^2 * ( x- μ )^2 ) + 1 / x^2
    else
        zero(T)
    end
end

## Beta
# functions
function heslogpdf(d::_DST.Beta{T}, x::Real) where T<:Real
    if _DST.insupport(_DST.Beta, x)
        α, β = _DST.params(d)
        - (α - 1) / x^2 - (β - 1) / (1 - x)^2
    else
        zero(T)
    end
end

## Extended Beta
# struct
"""
    ExtendedBeta

The [**extended beta distribution**](https://www.vosesoftware.com/riskwiki/Beta4distribution.php)
with shape parameters `α` and `β`, and optional support parameters `min` and
`max` has a probability density function
```math
f(x, α, β, \\text{min}, \\text{max}) =
    \\begin{cases}
        0,                                                                                &\\text{if:}~x < \\text{min},               \\\\
        \\frac{(x-\\text{min})^{α-1} (\\text{max}-x)^{β-1}}{B(α,β) (\\text{max}-\\text{min})^{α+β-1}},  &\\text{if:}~\\text{min} ≤ x ≤ \\text{max}, \\\\
        0,                                                                                &\\text{if:}~x > max,
    \\end{cases}
```
where B(α,β) is a Beta function.
"""
struct ExtendedBeta{T<:Real} <: _DST.ContinuousUnivariateDistribution
    α::T
    β::T
    min::T
    max::T
end

# additional constructors
ExtendedBeta(α::Real, β::Real) = ExtendedBeta(α, β, 0.0, 1.0)

# functions
scale(dst::ExtendedBeta) = (dst.α, dst.β)
params(dst::ExtendedBeta) = (dst.α, dst.β, dst.min, dst.max)
insupport(dst::ExtendedBeta, x::Real) = dst.min ≤ x ≤ dst.max

mean(dst::ExtendedBeta) = dst.min + dst.α * (dst.max - dst.min) / (dst.α + dst.β)
function mode(dst::ExtendedBeta)
    α, β, min, max = params(dst)
    if α > 1 && β > 1       return min + (α - 1) * (max - min) / (α + β - 2)
    elseif α < 1  && β ≥ 1  return min
    elseif α == 1 && β > 1  return min
    elseif α ≥ 1  && β < 1  return max
    elseif α > 1  && β == 1 return max
    elseif α ≤ 1  && β ≤ 1  error("mode is defined only when α > 1 and/or β > 1")
    end
end
skewness(dst::ExtendedBeta) =
    2 * (dst.β - dst.α) / (dst.α + dst.β + 2) * sqrt((dst.α + dst.β + 1) / (dst.α * dst.β))

function pdf(dst::ExtendedBeta{T}, x::Real) where T<:Real
    α, β, min, max = params(dst)
    if insupport(dst, x)
        (x - min)^(α - 1) * (max - x)^(β - 1) / _SF.beta(α,β) / (max - min)^(α + β - 1)
    else
        zero(T)
    end
end
function logpdf(dst::ExtendedBeta{T}, x::Real) where T<:Real
    α, β, min, max = params(dst)
    if insupport(dst, x)
        (α - 1) * log(x - min) + (β - 1) * log(max - x) - _SF.logbeta(α,β) - (α + β - 1) * log(max - min)
    else
        -T(Inf)
    end
end
function gradlogpdf(dst::ExtendedBeta{T}, x::Real) where T<:Real
    α, β, min, max = params(dst)
    if insupport(dst, x)
        (α - 1) / (x - min) - (β - 1) / (max - x)
    else
        zero(T)
    end
end
function heslogpdf(dst::ExtendedBeta{T}, x::Real) where T<:Real
    α, β, min, max = params(dst)
    if insupport(dst, x)
        - (α - 1) / (x - min)^2 - (β - 1) / (max - x)^2
    else
        zero(T)
    end
end
function maximum(dst::ExtendedBeta{T}) where T<:Real
    return params(dst)[4]
end
function minimum(dst::ExtendedBeta{T}) where T<:Real
    return params(dst)[3]
end

rand(dst::ExtendedBeta, N::Int) =
    (dst.max - dst.min) .* rand(_DST.Beta(dst.α, dst.β), N) .+ dst.min

rand(r::_RAN.AbstractRNG, dst::ExtendedBeta, N::Int) =
    (dst.max - dst.min) .* rand(r, _DST.Beta(dst.α, dst.β), N) .+ dst.min

## Gamma heslogpdf
function heslogpdf(d::_DST.Gamma{T}, x::Real) where T<:Real
    if _DST.insupport(_DST.Gamma, x)
        α, θ = _DST.params(d)
        - (α - 1) / x^2
    else
        zero(T)
    end
end
##
# gradlogpdf and heslogpdf for Gaussian mixtures done directly
# with Distributions.jl
function gradlogpdf(d::_DST.MixtureModel, x::Real)
    @assert all([c isa _DST.Normal for c in d.components]) "Only Gaussian mixture models are supported!"
    σ = [c.σ for c in d.components]
    μ = [c.μ for c in d.components]
    γ = d.prior.p./(sqrt(2*π)*σ)
    sum([-γ[i]*(x-μ[i])*exp(-0.5*(x-μ[i])^2/σ[i]^2)*σ[i]^(-2) for i in 1:length(d.components)])/sum([γ[i]*exp(-0.5*(x-μ[i])^2/σ[i]^2) for i in 1:length(d.components)])
end

function heslogpdf(d::_DST.MixtureModel, x::Real)
    σ = [c.σ for c in d.components]
    μ = [c.μ for c in d.components]
    γ = d.prior.p./(sqrt(2*π)*σ)
    p1 = sum([γ[i]*(x-μ[i])^2*exp(-0.5*(x-μ[i])^2/σ[i]^2)*σ[i]^(-4) - γ[i]*exp(-0.5*(x-μ[i])^2/σ[i]^2)*σ[i]^(-2) for i in 1:length(d.components)])
    p2 = sum([γ[i]*exp(-0.5*(x-μ[i])^2/σ[i]^2) for i in 1:length(d.components)])
    p3 = sum([-γ[i]*(x-μ[i])*exp(-0.5*(x-μ[i])^2/σ[i]^2)*σ[i]^(-2) for i in 1:length(d.components)])^2
    p4 = sum([γ[i]*exp(-0.5*(x-μ[i])^2/σ[i]^2) for i in 1:length(d.components)])^2
    return p1/p2-p3/p4
end


text:
import Base: eltype

# Gregorian calendar-ic globals
const second = 1.0
const minute = 60second
const hour = 60minute
const day = 24hour
const year = 365day
const stellaryear = 23hour + 56minute + 4.098903691
const Ω = 2π/stellaryear

@inline minuszero(args...) = -0
@inline maxsqrt(ϕ::T) where T = sqrt(max(zero(T), ϕ))
@inline maxsqrt(ϕ, i) = @inbounds maxsqrt(ϕ[i])
@inline nan2inf(a::T) where T = ifelse(isnan(a), T(Inf), a)
@inline inf2zero(a::T) where T = ifelse(isinf(a), zero(T), a)

function pressenter()
    println("\nPress enter to continue.")
    chomp(readline())
end

macro zeros(T, dims, vars...)
    expr = Expr(:block)
    append!(expr.args, [:($(esc(var)) = zeros($(esc(T)), $(esc(dims))); ) for var in vars])
    expr
end

macro def(name, definition)
    return quote
        macro $(esc(name))()
            esc($(Expr(:quote, definition)))
        end
    end
end

function build_homogeneous_solution(names, typeprefix=:CellField; prefix=Symbol(""))
    nfields = length(names)
    fieldtype = Expr(:curly, typeprefix, :A, :G, :T)
    exprs = [ :( $(names[i]) :: $fieldtype ) for i = 1:nfields ]

    fullname = Symbol(prefix, :Solution)
    signature = Expr(:curly, fullname, :A, :G, :T)

    return quote
        #import StaticArrays: FieldVector

        struct $signature <: AbstractSolution{$(nfields), $fieldtype}
            $(exprs...)
        end
    end
end

function build_accessories(names; prefix=Symbol(""))
    nfields = length(names)
    bcsfields = [ :( $(names[i]) :: FieldBoundaryConditions ) for i = 1:nfields ]
    accfields = [ :( $(names[i]) :: T ) for i = 1:nfields ]
    lhsfields = [ :( $(names[i]) :: Tridiagonal{T, A} ) for i = 1:nfields ]

    bcs_signature = Symbol(prefix, :BoundaryConditions)
    acc_signature = Expr(:curly, Symbol(prefix, :Accessory), :T)
    lhs_signature = Expr(:curly, Symbol(prefix, :LeftHandSide), :T, :A)

    return quote
        import StaticArrays: FieldVector
        import LinearAlgebra: Tridiagonal
        import OceanTurb: build_lhs

        struct $bcs_signature <: FieldVector{$(nfields), FieldBoundaryConditions}
            $(bcsfields...)
        end

        struct $acc_signature <: FieldVector{$(nfields), T}
            $(accfields...)
        end

        struct $lhs_signature <: FieldVector{$(nfields), Tridiagonal{T, A}}
            $(lhsfields...)
            function LeftHandSide(solution::AbstractSolution{1, AbstractField})
                lhs = build_lhs(solution)[1]
                new{eltype(solution[1]), arraytype(solution[1])}(lhs)
            end
        end

        function LeftHandSide(solution::AbstractSolution{N, AbstractField}) where N
            lhs = build_lhs(solution)
            LeftHandSide{eltype(solution[1]), arraytype(solution[1])}(lhs...)
        end
    end
end

macro solution(names...)
    esc(quote
        $(build_homogeneous_solution(names))
        $(build_accessories(names))
    end
    )
end

macro prefixed_solution(prefix, names...)
    esc(quote
        $(build_homogeneous_solution(names; prefix=prefix))
        $(build_accessories(names; prefix=prefix))
    end
    )
end

@def add_standard_model_fields begin
    clock       :: Clock{T}
    grid        :: G
    timestepper :: TS
    solution    :: Solution
    bcs         :: BoundaryConditions
end

@def add_clock_grid_timestepper begin
    clock       :: Clock{T}
    grid        :: G
    timestepper :: TS
end

"""
    run_until!(model, dt, tfinal)

Run `model` until `tfinal` with time-step `dt`.
"""
function run_until!(model, dt, tfinal)
    nt = floor(Int, (tfinal - time(model))/dt)
    time_step!(model, dt, nt)

    last_dt = tfinal - time(model)
    time_step!(model, last_dt)

    return nothing
end

"""
    diffusive_flux!(flux, fieldname, model)

Calculates the diffusive flux associated with `fieldname` in `model.solution.fieldname`
and stores the result in the `FaceField` `flux`.
"""
function diffusive_flux!(flux, fieldname, model)
    field = getproperty(model.solution, fieldname)
    K = getproperty(model.timestepper.eqn.K, fieldname)

    for i in interiorindices(flux)
        @inbounds flux[i] = - K(model, i) * ∂z(field, i)
    end

    return nothing
end

"""
    diffusive_flux(fieldname, model)

Returns `flux::FaceField` with the diffusive flux associated with `fieldname`
in `model.solution.fieldname`.
"""
function diffusive_flux(fieldname, model)
    flux = FaceField(model.grid)
    diffusive_flux!(flux, fieldname, model)
    return flux
end

"""
    prettytime(t)

Convert a floating point value `t` representing an amount of time in seconds to a more
human-friendly formatted string with three decimal places. Depending on the value of `t`
the string will be formatted to show `t` in nanoseconds (ns), microseconds (μs),
milliseconds (ms), seconds (s), minutes (min), hours (hr), or days (day).
"""
function prettytime(t)
    # Modified from: https://github.com/JuliaCI/BenchmarkTools.jl/blob/master/src/trials.jl
    iszero(t) && return "0.000 s"

    if t < 1e-6
        value, units = t * 1e9, "ns"
    elseif t < 1e-3
        value, units = t * 1e6, "μs"
    elseif t < 1
        value, units = t * 1e3, "ms"
    elseif t < minute
        value, units = t, "s"
    elseif t < hour
        value, units = t / minute, "min"
    elseif t < day
        value, units = t / hour, "hr"
    else
        value, units = t / day, "day"
    end

    return @sprintf("%.3f", value) * " " * units
end


text:
@doc """
    findid(id::String, S::GphysData)
    findid(S::GphysData, id::String)

Get the index of the first channel in S where `id.==S.id` is true. Returns 0
for failure.

    findid(S::GphysData, T::GphysData)

Get index corresponding to the first channel in T that matches each ID in S;
equivalent to [findid(id,T) for id in S.id].

    findid(C::SeisChannel, S::SeisData)
    findid(S::SeisData, C::SeisChannel)

Get the index to the first channel `c` in S where `S.id[c]==C.id`.
""" findid
function findid(id::String, ID::Array{String,1})
  c = 0
  for i = 1:length(ID)
    if ID[i] == id
      c = i
      break
    end
  end
  return c
end

function findid(id::DenseArray{UInt8,1}, ID::Array{String,1})
  c = 0
  for i = 1:length(ID)
    if codeunits(ID[i]) == id
      c = i
      break
    end
  end
  return c
end

function findid(id::Union{Regex,String}, S::T) where {T<:GphysData}
  j=0
  for i=1:length(S.id)
    if S.id[i] == id
      j=i
      break
    end
  end
  return j
end
findid(S::T, id::Union{String,Regex})  where {T<:GphysData} = findid(id, S)
# DND ...why the fuck is findfirst so fucking slow?!

findid(C::TC, S::TS) where {TC<:GphysChannel, TS<:GphysData} = findid(getfield(C, :id), S)
findid(S::TS, C::TC) where {TC<:GphysChannel, TS<:GphysData} = findid(getfield(C, :id), S)


text:
using SparseArrays, LinearAlgebra, ProgressMeter, Serialization

function spcT()
    @info "Reading the data and genotyped ID"
    Ai = deserialize("Ainv.ser")
    jx = begin
        tmp = Int[]
        for i in eachline("genotyped_renumid_sort.dat")
            push!(tmp, parse(Int, i))
        end
        tmp
    end
    N = size(Ai, 1)
    n2 = length(jx)
    ix = sort(collect(setdiff(Set(1:N), Set(jx))))
    A11 = Ai[ix, ix]
    A12 = Ai[ix, jx]
    A22 = Matrix(Ai[jx, jx])
    nr = n2 ÷ 10
    stops = collect(nr:nr:n2)
    stops[end] = n2
    start = 1
    @info "Calculating by blocks"
    @showprogress for stop in stops
        rhs = Matrix(A12[:, start:stop])
        sol = A11 \ rhs
        blk = A12'sol
        A22[:, start:stop] -= blk
        start = stop + 1
    end
    serialize("a22.txt", A22)
end


text:
@doc """
### AdvectionEquationMod.applyGLS2

  This function updates the residual with the GLS term, using the definition
  of the parameter tau from Hughes Part III.

  Inputs:
    mesh
    sbp
    eqn
    opts
    src_func:  The functor that calculates the source term at a node

  Aliasing restrictions: none
"""->
function applyGLS2(mesh::AbstractMesh{Tmsh}, sbp::AbstractOperator, eqn::AdvectionData{Tsol, Tres, 2}, opts, src_func::SRCType) where {Tmsh, Tsol, Tres}

#  println("----- entered applyGLS2 -----")
  Dx = zeros(Float64, mesh.numNodesPerElement, mesh.numNodesPerElement)
  Dy = zeros(Dx)

  # working vectors that do the reduction
  red_vec1 = zeros(Tsol, mesh.numNodesPerElement)
  red_vec2 = zeros(Tsol, mesh.numNodesPerElement)
  red_vec3 = zeros(red_vec1)
  red_vec4 = zeros(red_vec2)
  u = zeros(Tsol, mesh.numNodesPerElement)
  dxidx = zeros(Tmsh, 2,2, mesh.numNodesPerElement)
  tau_vec = zeros(Tres, mesh.numNodesPerElement)
  # calculate Dxi and Deta
  Dxi = diagm(1./sbp.w)*sbp.Q[:, :, 1]
  Deta = diagm(1./sbp.w)*sbp.Q[:, :, 2]

  for i=1:mesh.numEl
    # calculate the quantities needed for this element
    dxidx_hat = sview(mesh.dxidx[:, :, :, i])
    jac = sview(mesh.jac, :, i)
    res = sview(eqn.res, :, :, i)
    # constant coefficient advection only!
    alpha_x = eqn.params.alpha_x
    alpha_y = eqn.params.alpha_y
    # calclate the true dxidx (not scaled by the mapping jacobian)
    for j=1:mesh.numNodesPerElement
      for k=1:2
        for p=1:2
          dxidx[p, k, j] = dxidx_hat[p, k, j]*jac[j]
        end
      end
    end

    for j=1:mesh.numNodesPerElement
      for k=1:mesh.numNodesPerElement
        Dx[k, j] = dxidx[1,1]*Dxi[k, j] + dxidx[2, 1]*Deta[k, j]
        Dy[k, j] = dxidx[1,2]*Dxi[k, j] + dxidx[2, 2]*Deta[k, j]
      end
    end

   for j=1:mesh.numNodesPerElement
#     alpha_xj = eqn.params.alpha_x[1, j, i]
#     alpha_yj = eqn.params.alpha_y[1, j, i]
     dxidx_j = sview(dxidx, :, :, j)
     p = 2
#     tau_vec[j] = getTau(alpha_xj, alpha_yj, dxidx_j, p)

     tau_vec[j] = getTau(alpha_x, alpha_y, jac[j], mesh.min_node_dist)
#     println("tau_vec[$j] = ", tau_vec[j])
   end

   for j=1:mesh.numNodesPerElement
     u[j] = eqn.q[1, j, i]
   end

   # now compute the GLS term
   # trial space
   smallmatvec!(Dx, u, red_vec1)
   smallmatvec!(Dy, u, red_vec2)

   for j=1:mesh.numNodesPerElement
     red_vec1[j] = alpha_x*red_vec1[j] + alpha_y*red_vec2[j]
   end

   # add source term
   for j=1:mesh.numNodesPerElement
     coords_j = sview(mesh.coords, :, j, i)
     red_vec1[j] -= src_func(eqn.params, coords_j, eqn.t)
   end

   # middle terms
   # also copy pinto red_vec2 at same time
   for j=1:mesh.numNodesPerElement
     red_vec1[j] *= tau_vec[j]*sbp.w[j]/jac[j]
     red_vec2[j] = red_vec1[j]
   end

   smallmatTvec!(Dx, red_vec1, red_vec3)
   smallmatTvec!(Dy, red_vec2, red_vec4)

   for j=1:mesh.numNodesPerElement
     red_vec3[j] = alpha_x*red_vec3[j] + alpha_y*red_vec4[j]
   end

   # update res
   for j=1:mesh.numNodesPerElement
     res[j] -= red_vec3[j]
   end

  end  # end loop over elements

  #  println("----- Finished applyGLS2 -----")
  return nothing

end  # end function

@doc """
### AdvectionEquationMod.getTau

  This function calculates the stabilization parameter tau for GLS at a node,
  using the definition given in Hughes Part III

  Inputs:
    alpha_x: advection coefficient in the x direction at the node
    alpha_y: advection coefficient in the y direction at the node
    dxidx: a 2x2 matrix containing the mapping jacobian at the node
    p : which p norm to use (typically either 1 or 2)

  Outputs:
    tau: the value of tau at the node

"""->
function getTau(alpha_x, alpha_y, dxidx::AbstractMatrix, p)
  fac = 1.0
  b1 = dxidx[1,1]*alpha_x + dxidx[1,2]*alpha_y
  b2 = dxidx[2,1]*alpha_x + dxidx[2,2]*alpha_y
  bp = absvalue(b1)^p + absvalue(b2)^p
  return fac*(bp^(-1/p))
end

@doc """
### AdvectionEquationMod.getTau

  This function calculate the stabilization parameter tau at a node using
  the (heuristic) definition 0.5*h/sqrt(alpha_x^2 + alpha_y^2)

  Inputs:
    alpha_x: advection coefficient in the x direction at the node
    alpha_y: advection coefficient in the y direction at the node
    jac: determinant of the mapping jacobian at the node

  Outputs:
    tau: the value of tau at the node
"""->
function getTau(alpha_x, alpha_y, jac::Tmsh, min_node_dist) where Tmsh

  fac = 1
  h = (1/sqrt(jac))/2  # /2 because reference element is -1 to 1
  alpha_nrm = sqrt(alpha_x*alpha_x + alpha_y*alpha_y)
  return fac*0.5*h/alpha_nrm
end


text:
"""
An n × n almost-banded matrix is defined by

A = B + S,

where B is an n × n banded matrix and S is semi-separable matrix given by

S = triu(uvᵀ, bu+1) + tril(pqᵀ, -bl-1).

Here, bu, bl ⩾ 0. u and v are n × ru matrices and p and q are n × rl matrices.
See [1] for the definition.

[1] Chandrasekaran and Gu, Fast and stable algorithms for banded plus
Semiseparable systems of linear equations, SIMAX, 25 (2003), pp. 373-384.
"""

struct AlmostBandedMatrix{T} <: AbstractMatrix{T}
    B::BandedMatrix{T}
    S::SemiSeparableMatrix{T}

    function AlmostBandedMatrix(B::BandedMatrix{T}, S::SemiSeparableMatrix{T}) where T
        @assert size(B) == size(S)
        new{T}(B, S)
    end
end


text:
using DataFrames, JuMP, CPLEX, Dates
include("utils/extract_data.jl")
include("utils/display_results.jl")
include("utils/misc.jl")
include("single_depot/solve_single_depot.jl")
include("multi_depot/solve_multi_depot.jl")

function solve_problem(file, nb_vehicles, approx, verbose)

    print("\n\n")
    println("===========================================")
    println("GROOMER OPTIMIZATION PROBLEM")
    println("-----------")

    println("Extracting data from file ...")
    (data, graph) = extract_data(file)
    data[1, :nb_vehicles] = nb_vehicles

    println("-----------")

    print_info_dataset(data)

    println("-----------")

    if length(data[1, :depots])>1
        (tours, solving_time) = solve_multi_depot(data, graph, approx, verbose)
    elseif length(data[1, :depots])==1
        (tours, solving_time) = solve_single_depot(data, graph)
    else
        error("Problem with depot list. Check dataset's validity.")
    end

    println("-----------")

    display_results(tours, graph, data, solving_time)

    println("===========================================")
end

max_stars_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_issues_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_forks_repo_name: TheoGuyard/Groomer-Optimization-Problem


text:
import gmsh
using JLD2

gmsh.initialize()
gmsh.open("square.msh")
nodes = gmsh.model.mesh.getNodes()
nodes_label = convert(Array{Int64}, nodes[1])
nodes_N = length(nodes_label)
nodes_coordinate = reshape(nodes[2], (3, nodes_N))
nodes_boundary = convert(Array{Int64}, unique(gmsh.model.mesh.getNodesByElementType(1)[1]))
nodes_boundary_N = length(nodes_boundary)

elements = gmsh.model.mesh.getElements()

elements_bound_label = convert(Array{Int64}, elements[2][1])
elements_bound_N = length(elements_bound_label)

elements_int_label = convert(Array{Int64}, elements[2][2])
elements_int_N = length(elements_int_label)

elements_bound = reshape(convert(Array{Int64}, elements[3][1]), (3, elements_bound_N))
elements_int = reshape(convert(Array{Int64}, elements[3][2]), (9, elements_int_N))
#gmsh.fltk.run()
gmsh.finalize()
@save "square.jld2"