text:
# Benchmark the calculation of a distance in a PDB file
# The distance is the closest distance between any atoms of residues 50 and 60
#   of chain A in 1AKE

using MIToS.PDB

pdb_filepath = "data/1AKE.pdb"
struc = read(pdb_filepath, PDBFile, model="1", chain="A", group="ATOM")

# Run to JIT compile
distance(struc[50], struc[60])

elapsed = @elapsed distance(struc[50], struc[60])

println(elapsed)

max_stars_repo_name: wojdyr/pdb-benchmarks

max_issues_repo_name: wojdyr/pdb-benchmarks

max_forks_repo_name: wojdyr/pdb-benchmarks


text:
# Autogenerated wrapper script for Singular_jll for x86_64-unknown-freebsd11.1-cxx03
export Singular, libfactory, libomalloc, libparse, libpolys, libsingular, libsingular_resources

using cddlib_jll
using FLINT_jll
using GMP_jll
using MPFR_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"
LIBPATH_default = ""

# Relative path to `Singular`
const Singular_splitpath = ["bin", "Singular"]

# This will be filled out by __init__() for all products, as it must be done at runtime
Singular_path = ""

# Singular-specific global declaration
function Singular(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(Singular_path)
    end
end

# Relative path to `libfactory`
const libfactory_splitpath = ["lib", "libfactory.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libfactory_path = ""

# libfactory-specific global declaration
# This will be filled out by __init__()
libfactory_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libfactory = "libfactory-4.1.3.so"

# Relative path to `libomalloc`
const libomalloc_splitpath = ["lib", "libomalloc.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libomalloc_path = ""

# libomalloc-specific global declaration
# This will be filled out by __init__()
libomalloc_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libomalloc = "libomalloc-0.9.6.so"

# Relative path to `libparse`
const libparse_splitpath = ["bin", "libparse"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libparse_path = ""

# libparse-specific global declaration
function libparse(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(libparse_path)
    end
end

# Relative path to `libpolys`
const libpolys_splitpath = ["lib", "libpolys.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libpolys_path = ""

# libpolys-specific global declaration
# This will be filled out by __init__()
libpolys_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libpolys = "libpolys-4.1.3.so"

# Relative path to `libsingular`
const libsingular_splitpath = ["lib", "libSingular.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libsingular_path = ""

# libsingular-specific global declaration
# This will be filled out by __init__()
libsingular_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libsingular = "libSingular-4.1.3.so"

# Relative path to `libsingular_resources`
const libsingular_resources_splitpath = ["lib", "libsingular_resources.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libsingular_resources_path = ""

# libsingular_resources-specific global declaration
# This will be filled out by __init__()
libsingular_resources_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libsingular_resources = "libsingular_resources-4.1.3.so"

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"Singular")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (cddlib_jll.PATH_list, FLINT_jll.PATH_list, GMP_jll.PATH_list, MPFR_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (cddlib_jll.LIBPATH_list, FLINT_jll.LIBPATH_list, GMP_jll.LIBPATH_list, MPFR_jll.LIBPATH_list,))

    global Singular_path = normpath(joinpath(artifact_dir, Singular_splitpath...))

    push!(PATH_list, dirname(Singular_path))
    global libfactory_path = normpath(joinpath(artifact_dir, libfactory_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libfactory_handle = dlopen(libfactory_path)
    push!(LIBPATH_list, dirname(libfactory_path))

    global libomalloc_path = normpath(joinpath(artifact_dir, libomalloc_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libomalloc_handle = dlopen(libomalloc_path)
    push!(LIBPATH_list, dirname(libomalloc_path))

    global libparse_path = normpath(joinpath(artifact_dir, libparse_splitpath...))

    push!(PATH_list, dirname(libparse_path))
    global libpolys_path = normpath(joinpath(artifact_dir, libpolys_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libpolys_handle = dlopen(libpolys_path)
    push!(LIBPATH_list, dirname(libpolys_path))

    global libsingular_path = normpath(joinpath(artifact_dir, libsingular_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libsingular_handle = dlopen(libsingular_path)
    push!(LIBPATH_list, dirname(libsingular_path))

    global libsingular_resources_path = normpath(joinpath(artifact_dir, libsingular_resources_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libsingular_resources_handle = dlopen(libsingular_resources_path)
    push!(LIBPATH_list, dirname(libsingular_resources_path))

    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(vcat(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ':')
end  # __init__()

ext: jl


text:
"""
    EB{TI<:Integer,TR<:Real} <: VagueAlternative

    EB( # default values
      samplespace::SampleSpace,
      p0::Real, pmcrv::Real, prior,
      gamma::Real, lambda::Real;
      a::Real = 1, b::Real = 1, targetpower::Real = .8, k::Real = 1,
      alpha::Real = 0.05,
      minconditionalpower::Real = 0.0, MONOTONECONDITIONALPOWER::Bool = true,
      npriorpivots::Integer = 50, ngpivots::Integer = 15,
      label::String = ""
    )

This type represents a set of parameters for finding optimal two-stage designs
maximizing the expected benefit (EB).

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| samplespace  | a sample space object |
| p0           | upper boundary of the null hypothesis |
| pmcrv        | point alternative to power on |
| prior        | prior distribution on true response rate, expected score under the prior conditional on a relevant effect is minimized |
| gamma        | per patient costs |
| lambda       | benefit of correctly rejecting null if relevant effect exists |
| a            | scaling parameter of nonlinear power transfrom (set to infinity for thresholding) |
| b            | scaling parameter of nonlinear power transfrom |
| targetpower  | target power for thresholding |
| k            | factor for non-responder costs |
| alpha        | maximal type one error rate |
| minconditionalpower | minimal conditional power upon continuation to stage two |
| MONOTONECONDITIONALPOWER | if true, the conditional power must be monotonously increasing, this constraint is only relevant if nmax is set very restrictively |
| npriorpivots | number of pivot points for integration of prior |
| npivots      | number of pivot points piecewise linear approximation of local score |
| label        | descriptive string for parameter set |
"""
mutable struct EB{TI<:Integer,TR<:Real} <: VagueAlternative

  label::String
  samplespace::SampleSpace
  p0::TR
  pmcrv::TR
  prior
  gamma::TR
  lambda::TR
  a::TR
  b::TR
  targetpower::TR
  k::TR
  mtoer::TR
  minconditionalpower::TR
  MONOTONECONDITIONALPOWER::Bool
  npriorpivots::TI # number of pivots for prior evaluation
  ngpivots::TI # number of pivots for approximation of g

  function EB{TI,TR}(
    label::String,
    samplespace::SampleSpace,
    p0::TR, pmcrv::TR, prior,
    gamma::TR, lambda::TR,
    a::TR, b::TR, targetpower::TR, k::TR,
    alpha::TR,
    minconditionalpower::TR, MONOTONECONDITIONALPOWER::Bool,
    npriorpivots::TI, ngpivots::TI
  ) where {TI<:Integer,TR<:Real}

    @checkprob p0 pmcrv alpha minconditionalpower targetpower
    z = quadgk(prior, 0, 1, abstol = .0001)[1]
    abs(z - 1) > .001 ?
      error(@sprintf("prior must integrate to one, is %.6f", z)) : nothing
    gamma < 0 ?
      error(@sprintf("gamma must be positive, is %.6f", gamma)) : nothing
    lambda < 0 ?
      error(@sprintf("lambda must be positive, is %.6f", gamma)) : nothing
    lambda / gamma > 500  ?
      warn(@sprintf("lambda/gamma is large (%.3f), potentially infeasible", lambda/gamma)) : nothing
    a <= 0 ?
      error(@sprintf("a must be positive, is %.3f", a)) : nothing
    b <= 0 ?
      error(@sprintf("b must be positive, is %.3f", b)) : nothing
    npriorpivots < 5 ?
      error(@sprintf("npriorpivots must be at least 5, is %i", npriorpivots)) : nothing
    ngpivots < 5 ?
      error(@sprintf("npivots must be at least 5, is %i", ngpivots)) : nothing
    new(
      label,
      samplespace,
      p0, pmcrv, prior,
      gamma, lambda,
      a, b, targetpower, k,
      alpha,
      minconditionalpower, MONOTONECONDITIONALPOWER,
      npriorpivots, ngpivots
    )

    end # inner constructor

end # EB

function EB( # default values
    samplespace::SampleSpace,
    p0::Real, pmcrv::Real, prior,
    gamma::Real, lambda::Real;
    a::Real = 1, b::Real = 1, targetpower::Real = .8, k::Real = 1,
    alpha::Real = 0.05,
    minconditionalpower::Real = 0.0, MONOTONECONDITIONALPOWER::Bool = true,
    npriorpivots::Integer = 50, ngpivots::Integer = 15,
    label::String = ""
)

  p0, pmcrv, gamma, lambda, alpha, targetpower, a, b, k, minconditionalpower =
    promote(p0, pmcrv, gamma, lambda, alpha, targetpower, a, b, k, minconditionalpower)
  npriorpivots, ngpivots = promote(npriorpivots, ngpivots)

  EB{typeof(ngpivots),typeof(p0)}(
    label,
    samplespace,
    p0, pmcrv, prior,
    gamma, lambda,
    a, b, targetpower, k,
    alpha,
    minconditionalpower, MONOTONECONDITIONALPOWER,
    npriorpivots, ngpivots
  )

end # EB

maxsamplesize(params::EB) = maxsamplesize(params.samplespace)

isgroupsequential(params::EB) = isgroupsequential(params.ss)

hasmonotoneconditionalpower(params::EB) = params.MONOTONECONDITIONALPOWER

minconditionalpower(params::EB) = params.minconditionalpower

Base.show(io::IO, params::EB) = print("EB")

function expectedtransformedpower(
  design::Design, params::Union{EB, MBESS}, x1::TI
) where {TI<:Integer}

  x1 < 0 ?
    error(@sprintf("x1 must be positive, is %i", x1)) : nothing
  n1       = interimsamplesize(design)
  x1 > n1 ?
    error(@sprintf("x1 must be smaller than n1, is %i", x1)) : nothing
  pmcrv    = mcrv(params)
  phi(p)   = prior(params, p)
  z        = quadgk(p -> phi(p) * Distributions.pdf(Distributions.Binomial(n1, p), x1), pmcrv, 1)[1]
  omega(p) = p < pmcrv ? 0 : phi(p) * Distributions.pdf(Distributions.Binomial(n1, p), x1) / z # conditional prior for stage 2
  return quadgk(p -> g(params, power(design, x1, p)) * omega(p), pmcrv, 1)[1]

end # expectedtransformedpower

function expectedtransformedpower(design::Design, params::Union{EB, MBESS})

  pmcrv  = mcrv(params)
  phi(p) = prior(params, p)
  n1     = interimsamplesize(design)
  z      = quadgk(phi, pmcrv, 1)[1]
  omega(p) = p < pmcrv ? 0 : phi(p) / z # conditional prior given effect
  return quadgk(p -> g(params, power(design, p)) * omega(p), pmcrv, 1)[1]

end # expectedtransformedpower

function expectedcost(design::Design, params::EB, p::T) where {T<:Real}

  @checkprob p
  ess = SampleSize(design, p) |> mean
  return params.gamma * ess * (p + params.k * (1 - p))

end # expectedcost

function g(params::EB, power::T) where {T<:Real}

  @checkprob power
  if params.a == params.b == Inf
    if power >= params.targetpower
      return 1.0
    else
      return 0.0
    end
  else
    return Distributions.cdf(Distributions.Beta(params.a, params.b), power)
  end

end # g

function expectedbenefit(design::Design, params::EB, p::T) where {T<:Real}

  @checkprob p
  if p < mcrv(params)
    return 0.0
  else
     return params.lambda * g(params, power(design, p))
  end

end # expectedbenefit

function score(design::Design, params::EB, p::T) where {T<:Real}

  @checkprob p
  return -(expectedbenefit(design, params, p) - expectedcost(design, params, p))

end # score

function score(design::Design, params::EB)

  return quadgk(p -> params.prior(p)*score(design, params, p), 0, 1, reltol = .001)[1]

end # score

function completemodel(ipm::IPModel, params::EB, n1::T) where {T<:Integer}

    ss = samplespace(params)
    !possible(n1, ss) ? warn("completemodel(EB): n1 and sample space incompatible") : nothing
    # extract ip model
    m             = ipm.m
    y             = ipm.y
    nvals         = ipm.nvals
    cvals         = ipm.cvals
    cvalsfinite   = ipm.cvalsfinite
    cvalsinfinite = [-Inf; Inf]
    # extract other parameters
    nmax          = maxsamplesize(ss, n1)
    p0            = null(params)
    pmcrv         = mcrv(params)
    prior         = params.prior
    npriorpivots  = params.npriorpivots
    ngpivots      = params.ngpivots
    # get grid for conditional prior
    cpriorpivots, dccdf   = findgrid(prior, pmcrv, 1, npriorpivots)
    # add type one error rate constraint
    JuMP.@constraint(m,
        sum(dbinom(x1, n1, p0)*_cpr(x1, n1, n, c, p0)*y[x1, n, c] for
            x1 in 0:n1, n in nvals, c in cvals
        ) <= mtoer(params)
    )
    z = zeros(n1 + 1) # normalizing constants for prior conditional on p >= mcrv, X1=x1
    for x1 in 0:n1
        z[x1 + 1] = quadgk(p -> prior(p)*dbinom(x1, n1, p), pmcrv, 1, abstol = .001)[1]
    end
    JuMP.@expression(m, cep[x1 in 0:n1], # define expressions for conditional power given X1=x1
        sum(quadgk(p -> prior(p) * dbinom(x1, n1, p) * _cpr.(x1, n1, n, c, p), pmcrv, 1, abstol = 0.001)[1] / z[x1 + 1] * y[x1, n, c]
            for n in nvals, c in cvals
        )
    )
    for x1 in 0:n1
        JuMP.@constraint(m, # add conditional type two error rate constraint (power)
            cep[x1] >= minconditionalpower(params) * (1 - sum(y[x1, n1, c] for c in cvalsinfinite)) # must be conditional on continuation!
        )
        if (x1 >= 1) & hasmonotoneconditionalpower(params)
            JuMP.@constraint(m, # ensure monotonicity if required
                cep[x1] - cep[x1 - 1] >= 0
            )
        end
    end
    z = zeros(n1 + 1) # normalizing constant for prior conditional on X1=x1
    for x1 in 0:n1
        z[x1 + 1] = quadgk(p -> dbinom(x1, n1, p) * prior(p), 0, 1, abstol = .001)[1]
    end
    JuMP.@expression(m, ec[x1 in 0:n1], # weighted expected sample size conditional on X1=x1
        sum(
            (params.gamma*(x1 + params.k*(n1 - x1)) + # stage one cost
             quadgk(p -> params.gamma*(p + params.k*(1 - p))*(n - n1) * dbinom(x1, n1, p)*prior(p)/z[x1 + 1], 0, 1, abstol = .001)[1]  # stage two expected
            ) *
            y[x1, n, c]
            for n in nvals, c in cvals
        )
    )
    # construct expressions for power
    JuMP.@expression(m, designpower[p in cpriorpivots],
        sum(
            dbinom(x1, n1, p)*_cpr(x1, n1, n, c, p) * y[x1, n, c]
            for x1 in 0:n1, n in nvals, c in cvals
        )
    )
    prob_crv = quadgk(p -> prior(p), mcrv(params), 1)[1]
    if params.a == params.b == Inf
        JuMP.@variable(m, exceedstargetpower[1:npriorpivots], Bin)#pivots = [0; params.targetpower - .001; params.targetpower + .001; 1] # lambda formulaion requires edges! model step function as linear approximation
        for i in 1:npriorpivots
            JuMP.@constraint(m,
                designpower[cpriorpivots[i]] - params.targetpower + (1 - exceedstargetpower[i]) >= .00001
            )
        end
        JuMP.@expression(m, obj,
            sum( ec[x1]*z[x1 + 1] for x1 in 0:n1 ) - prob_crv*params.lambda*sum( exceedstargetpower[i]*dccdf[i] for i in 1:npriorpivots ) # negative score!
        )
    else
        lbound = quantile(Distributions.Beta(params.a, params.b), .025) # TODO: implement a = b = inf
        ubound = quantile(Distributions.Beta(params.a, params.b), .975)
        pivots = [0; collect(linspace(lbound, ubound, max(1, ngpivots - 2))); 1] # lambda formulaion requires edges! exploitn fact that function is locally linear!
        JuMP.@variable(m, 0 <= lambdaSOS2[cpriorpivots, pivots] <= 1)
        JuMP.@variable(m, wpwr[cpriorpivots]) # weighted power
        for p in cpriorpivots
            JuMP.addSOS2(m, [lambdaSOS2[p, piv] for piv in pivots])
            JuMP.@constraint(m, sum(lambdaSOS2[p, piv] for piv in pivots) == 1)
            JuMP.@constraint(m, sum(lambdaSOS2[p, piv]*piv for piv in pivots) == designpower[p]) # defines lambdas!
            if p >= params.pmcrv
                JuMP.@constraint(m, sum(lambdaSOS2[p, piv]*g(params, piv) for piv in pivots) == wpwr[p])
            else
                JuMP.@constraint(m, wpwr[p] == 0)
            end
        end
        JuMP.@expression(m, obj,
            sum( ec[x1]*z[x1 + 1] for x1 in 0:n1 ) - prob_crv*params.lambda*sum( wpwr[cpriorpivots[i]]*dccdf[i] for i in 1:npriorpivots ) # negative score!
        )
    end
    JuMP.@objective(m, Min,
        obj
    )
    return true

end # completemodel

# ToDo: implement!
_isfeasible(design::Design, params::EB) = true

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
using Pkg
Pkg.activate(".") # Create new environment in this folder

# first time you need to install dependencies
Pkg.add("Distributions")
Pkg.add("StatsBase")
Pkg.add(["DataFrames","DataFramesMeta","Chain"])
Pkg.add("Plots")
Pkg.add("CategoricalArrays")
Pkg.add("LightGraphs")
Pkg.add("FixedEffectModels")

# past the first time, you only need to instanciate the current folder
Pkg.instantiate() # Updates packages given .toml file

#----------------------------------------------------

using Distributions
using LinearAlgebra
using StatsBase
using DataFrames
using Plots
using CategoricalArrays
using FixedEffectModels
using Chain
using DataFramesMeta
using LightGraphs
using TikzGraphs

#----------------------------------------------------

α_sd = 1
ψ_sd = 1

csort = 0.5 # Sorting effect
csig  = 0.5 # Cross-sectional standard deviation
w_sigma = 0.2
λ = 0.1 # Moving Probability is fixed
nt = 10

function data_generating_process(α_sd,
                                ψ_sd,
                                csort,
                                csig,
                                w_sigma,
                                λ,
                                nt)

    cnetw = 0.2 # Network effect

    nk = 30
    nl = 10

    # approximate each distribution with some points of support
    ψ = quantile.(Normal(), (1:nk) / (nk + 1)) * α_sd
    α = quantile.(Normal(), (1:nl) / (nl + 1)) * ψ_sd

    # Let's create type-specific transition matrices
    # We are going to use joint normals centered on different values
    G = zeros(nl, nk, nk)
    for l in 1:nl, k in 1:nk
        G[l, k, :] = pdf( Normal(0, csig), ψ .- cnetw * ψ[k] .- csort * α[l])
        G[l, k, :] = G[l, k, :] ./ sum(G[l, k, :])
    end

    # We then solve for the stationary distribution over psis for each alpha value
    # We apply a crude fixed point approach
    H = ones(nl, nk) ./ nk
    for l in 1:nl
        M = transpose(G[l, :, :])
        for i in 1:100
            H[l, :] = M * H[l, :]
        end
    end

    # p1 = plot(G[1, :, :], xlabel="Previous Firm", ylabel="Next Firm", zlabel="G[1, :, :]", st=:wireframe)
    # p2 = plot(G[nl, :, :], xlabel="Previous Firm", ylabel="Next Firm", zlabel="G[nl, :, :]", st=:wireframe, right_margin = 10Plots.mm) # right_margin makes sure the figure isn't cut off on the right
    # plot(p1, p2, layout = (1, 2), size=[600,300])

    #--------------------------------------------------------------

    nt = 10
    ni = 10000

    # We simulate a balanced panel
    ll = zeros(Int64, ni, nt) # Worker type
    kk = zeros(Int64, ni, nt) # Firm type
    spellcount = zeros(Int64, ni, nt) # Employment spell

    for i in 1:ni

        # We draw the worker type
        l = rand(1:nl)
        ll[i,:] .= l

        # At time 1, we draw from H
        kk[i,1] = sample(1:nk, Weights(H[l, :]))

        for t in 2:nt
            if rand() < λ
                kk[i,t] = sample(1:nk, Weights(G[l, kk[i,t-1], :]))
                spellcount[i,t] = spellcount[i,t-1] + 1
            else
                kk[i,t] = kk[i,t-1]
                spellcount[i,t] = spellcount[i,t-1]
            end
        end

    end

    #------------------------------------------------------------

    firms_per_type = 15
    jj = zeros(Int64, ni, nt) # Firm identifiers

    draw_firm_from_type(k) = sample(1:firms_per_type) + (k - 1) * firms_per_type

    for i in 1:ni

        # extract firm type
        k = kk[i,1]

        # We draw the firm (one of firms_per_type in given group)
        jj[i,1] = draw_firm_from_type(k)

        for t in 2:nt
            if spellcount[i,t] == spellcount[i,t-1]
                # We keep the firm the same
                jj[i,t] = jj[i,t-1]
            else
                # We draw a new firm
                k = kk[i,t]

                new_j = draw_firm_from_type(k)
                # Make sure the new firm is actually new
                while new_j == jj[i,t-1]
                    new_j = draw_firm_from_type(k)
                end

                jj[i,t] = new_j
            end
        end
    end
    # Make sure firm ids are contiguous
    contiguous_ids = Dict( unique(jj) .=> 1:length(unique(jj))  )
    jj .= getindex.(Ref(contiguous_ids),jj);

    #----------------------------------------------------------------------------------------

    ii = repeat(1:ni,1,nt)
    tt = repeat((1:nt)',ni,1)
    df = DataFrame(i=ii[:], j=jj[:], l=ll[:], k=kk[:], α=α[ll[:]], ψ=ψ[kk[:]], t=tt[:], spell=spellcount[:]);

    return df, G, H

end

#Number of workers in the firm:

df,G,H = data_generating_process(α_sd,ψ_sd,csort,csig, w_sigma, λ, nt);

crossSectionMean = @chain df begin
    groupby([:j, :t])
    combine(nrow => :count)
    groupby(:j)
    combine(:count => mean)
    @aside println("On average each firm has $(round(mean(_.count_mean))) observations, averaging across time periods.")
    # first(_, 5)
end

#Number of movers at the firm, defined as people that moved into the firm:

println("Maximum number of moves by a single worker: ", maximum(df.spell))

moversDataFrame = @chain df begin
    groupby([:j, :i, :spell])
    combine(:t => mean)
    @transform!(:spell_true = :spell .> 0)
    groupby(:j)
    combine(:spell_true => sum)
    @aside println("On average each firm has $(round(mean(_.spell_true_sum))) movers over the 10 periods.")
end

#---------------------------------------------------------------------------------------------------

df[!, :lw] = df.α + df.ψ + w_sigma * rand(Normal(), size(df)[1]);
df[!, :dummy] .= 1;

#1. Mean wage within firm:
eventStudyPanel =  @chain df begin
    groupby([:j, :t]) # Group by firm and time
    transform(:lw => mean) # Obtain wage at that group level
    @transform!(:wage_percentile = cut(:lw_mean,4)) # Get wage_percentile at same level of agregation
    groupby([:i, :spell])
    transform(nrow => :years_at_firm) # Get number of years individual spend at a single firm
    sort([:i,:t], rev = true) # Sort over time by individual (this helps to get first 2 periods of last firm)
    groupby([:i, :spell])  # by group and spell
    transform(:dummy .=>  cumsum => :years_at_firm_sofar) # Get the number of years spent at a firm so far
    sort([:i,:t], rev = false) # Sort over time by individual (this helps to get first 2 periods of last firm)
    groupby([:i, :spell])  # by group and spell
    transform(:dummy .=>  cumsum => :years_at_firm_sofar_inverse) # Get the number of years spent at a firm so far, but backwards.
    @aside initialFirmDataFrame = @chain _ begin
        subset(:spell => ByRow(==(1)), :years_at_firm_sofar_inverse => ByRow(<=(2))) # Generate dataframe for initial firm, keep last two years (call it initialFirmDataFrame)
    end
    subset(:spell => ByRow(==(0)), :years_at_firm_sofar => ByRow(<=(2))) # Generate dataframe for subsequent firm, keep first two years
    append!(initialFirmDataFrame) # Append both dataframes
    groupby(:i) # group by person
    transform(nrow => :nreps) # and get number of repetitions
    subset(:nreps => ByRow(==(4))) # and get number of repetitions
    sort([:i,:t]) # sort by i and t to check if it worked out
    # Generating event time variable:
    sort([:i,:t], rev = false) # Sort over time by individual (this helps to get first 2 periods of last firm)
    groupby([:i])  # by worker and time
    transform(:dummy .=>  cumsum => :event_time) # Get the number of years spent at a firm so far, but backwards.
end

# Generate dictionary to index cut label:
percentile_cut = Dict(1:4 .=> unique(sort(eventStudyPanel.wage_percentile)))

# Get Card Heining and Kling event Figure:
# There may be a smart way to do this, here I'm brute forcing it ...
initial = 1
final = 4
function generateEventStudy(eventStudyPanel, initial, final)
    return  eventStudy =  @chain eventStudyPanel begin
                            @aside finalJob = @chain _ begin
                                subset(:wage_percentile => ByRow(==(percentile_cut[final])), :spell => ByRow(==(1)))
                            end
                            subset(:wage_percentile => ByRow(==(percentile_cut[initial])), :spell => ByRow(==(0)))
                            append!(finalJob) # Append both dataframes
                            groupby(:i)
                            transform(nrow => :nreps)
                            subset(:nreps => ByRow(==(4)))
                            groupby([:wage_percentile, :event_time])
                            combine(:lw => mean)
                            sort(:event_time)
    end
end

eventStudySwitchers = [generateEventStudy(eventStudyPanel, 1, 1).lw_mean,
                        generateEventStudy(eventStudyPanel, 1, 2).lw_mean,
                        generateEventStudy(eventStudyPanel, 1, 3).lw_mean,
                        generateEventStudy(eventStudyPanel, 1, 4).lw_mean]
eventStudySwitchers2 = [generateEventStudy(eventStudyPanel, 4, 1).lw_mean,
                        generateEventStudy(eventStudyPanel, 4, 2).lw_mean,
                        generateEventStudy(eventStudyPanel, 4, 3).lw_mean,
                        generateEventStudy(eventStudyPanel, 4, 4).lw_mean];

plot(1:4,eventStudySwitchers, label=["1 to 1" "1 to 2" "1 to 3" "1 to 4"],markershape = :square)
plot!(1:4,eventStudySwitchers2,label =  ["4 to 1 " "4 to 2" "4 to 3" "4 to 4"],markershape = :circle)
plot!(legend=:outertopright)

# Calibrating the parameters:
#----------------------------------------------------------------------------------------

# Generate function with variance decomposition:

function variance_decomposition(df)
    varianceDecomposition =  @chain df begin
        groupby([:k]) # Group by firm type
        combine(:α => mean, :ψ => mean)
    end

    std_α = std(varianceDecomposition.α_mean)
    std_ψ = std(varianceDecomposition.ψ_mean)
    std_αψ = cov(varianceDecomposition.α_mean,varianceDecomposition.ψ_mean)

    return [std_α std_ψ std_αψ]
end

α_sd = 1
ψ_sd = 1
csort = 0.5 # Sorting effect
csig  = 0.5 # Cross-sectional standard deviation
w_sigma = 0.2

gap = 0.20
α_sd_list = 0.01:gap:1
ψ_sd_list = 0.01:gap:1
csort_list = 0.01:gap:1 # Sorting effect
csig_list  = 0.01:gap:1 # Cross-sectional standard deviation
w_sigma_list = 0.01:gap:1

flat_gridpoints(grids) = vec(collect(Iterators.product(grids...)))
grid_points = flat_gridpoints((α_sd_list, ψ_sd_list, csort_list, csig_list,w_sigma_list));

jj = 1

true_var_decomp = [0.186 0.101 0.110]
difference_outcome = []

#Apply grid search, this can be optimized for sure...
for jj in 1:length(grid_points)
    α_sd,ψ_sd,csort,csig, w_sigma = grid_points[jj]
    df,_,_ = data_generating_process(α_sd,ψ_sd,csort,csig, w_sigma, λ, nt);
    var_decomp = variance_decomposition(df)
    euclidean_dist = sum((true_var_decomp-var_decomp).^2)
    push!(difference_outcome, euclidean_dist)
end

# Get index of minimum euclidean distance estimation:
_, index_calibration = findmin(difference_outcome)
α_sd,ψ_sd,csort,csig, w_sigma = grid_points[index_calibration]
df, G, H  = data_generating_process(α_sd,ψ_sd,csort,csig, w_sigma, λ, nt);

p1 = plot(G[1, :, :], xlabel="Previous Firm", ylabel="Next Firm", zlabel="G[1, :, :]", st=:wireframe)
p2 = plot(G[nl, :, :], xlabel="Previous Firm", ylabel="Next Firm", zlabel="G[nl, :, :]", st=:wireframe, right_margin = 10Plots.mm) # right_margin makes sure the figure isn't cut off on the right
plot(p1, p2, layout = (1, 2), size=[600,300])

plot(H, xlabel="Worker", ylabel="Firm", zlabel="H", st=:wireframe)

# The output reduces cross sectional variance too much, check what happens when ↑ csig ...
df, G, H  = data_generating_process(α_sd,ψ_sd,csort, .1, w_sigma, λ, nt);
var_decomp = variance_decomposition(df)
euclidean_dist = sum((true_var_decomp-var_decomp).^2)
println("Euclidean distance is: $euclidean_dist")

p1 = plot(G[1, :, :], xlabel="Previous Firm", ylabel="Next Firm", zlabel="G[1, :, :]", st=:wireframe)
p2 = plot(G[nl, :, :], xlabel="Previous Firm", ylabel="Next Firm", zlabel="G[nl, :, :]", st=:wireframe, right_margin = 10Plots.mm) # right_margin makes sure the figure isn't cut off on the right
plot(p1, p2, layout = (1, 2), size=[600,300])

# It seems csig contributes importantly to covariance! Check this again and get feedback from prof.
# I'll leave it at 0.05 which is 5 times it optimal* value.

# Delete method in case of mistake
# m = @which data_generating_process(α_sd,ψ_sd,csort,csig, w_sigma, λ, nt)
# Base.delete_method(m)

#--------------------------------------------------------------------------------------

# Getting connected set:

function individualDeterministicTransitionMatrix(df,ii)

    nfirms = length(unique(df.j))
    shift_i = @chain df begin
                        subset(:i => ByRow(.==(ii)))
                        sort(:t)
                        combine(first, groupby(_,:spell))
                    end

    nshifts = size(shift_i)[1]
    transitionMatrix_i = zeros(Int32, nfirms, nfirms);

    for ii in 1:nshifts
        if ii != nshifts
            current = shift_i.j[ii]
            next = shift_i.j[ii+1]
            transitionMatrix_i[current,next] = 1
        end
    end

    return transitionMatrix_i
end

function getConnectedDataSet(df)

    # Add all shifts happening in the economy during the whole time
    nfirms = length(unique(df.j))

    totalDeterministicShifts = zeros(nfirms, nfirms);

    for ii in unique(df.i)
        totalDeterministicShifts = totalDeterministicShifts + individualDeterministicTransitionMatrix(df,ii)
    end

    adjacencyMatrix = (UpperTriangular(totalDeterministicShifts) + transpose(LowerTriangular(totalDeterministicShifts))).>1

    adjacencyMatrix = adjacencyMatrix + transpose(adjacencyMatrix)

    println("There are $(sum(totalDeterministicShifts)) job shifts across time ...")

    println("There are $(sum(adjacencyMatrix)) edges between nodes ...")

    simpleGraph = SimpleGraph(adjacencyMatrix);

    connectedNetwork = connected_components(simpleGraph)

    connectedSet = connectedNetwork[1]

    println("We have only $(length(connectedSet)) firms fully connected $(length(connectedSet)/nfirms) of the market...") # Double check we might be wrong...

    return df_connected = df[in(connectedSet).(df.j),:]

    println("Due to unconnectedness we eliminated $(size(df)[1]-size(df_connected)[1]) observations, not much")

end

# Estimate AKM model:
"""
This part of the problem set is for you to implement the AKM estimator.
    As discussed in class, this can be done simply by updating, in turn,
    the worker FE and the firm FE.

    Start by appending 2 new columns `alpha_hat` and `psi_hat` to your data.
    Then loop over the following:

    1. Update `alpha_hat` by taking the mean within `i` net of firm FE
    2. Update `psi_hat` by taking the mean within `fid` net of worker FE

    Question 7:

     - Run the previous steps in a loop, and at each step evaluate how much
     the total mean square error has changed.
     Check that is goes down with every step.
     Stop when the MSE decreases by less than 1e-9.

    Note that you can increase speed by focusing on movers only first.
"""
#--------------------------------------------------------------------------------------

function akm_estimation(df_connected)

    df_connected[:,:alpha_hat] .= .0;
    df_connected[:,:psi_hat] .= .0;

    # Compute firm type fixed effects ols model:
    delta = Inf
    tol = 0.00001
    nIter = 1
    msePast = 0
    while delta>tol

        if nIter == 1
            # Regress controling for industry fixed effects...
            model_j = reg(df_connected, term(:lw) ~ fe(:j), save=true);
            # Obtain residuals (which are controled by industry fe):
            df_connected[:,:alpha_hat] = residuals(model_j);
        else
            # Just obtain the alpha parameters (not averaged) by netting out psi_hat from prev iter.
            df_connected[:,:alpha_hat] = df_connected[:,:lw] - df_connected[:,:psi_hat]
        end

        # Average fixed effects by individual:
        df_connected = transform(groupby(df_connected, :i), :alpha_hat => mean => :alpha_hat)

        # Net out individual fixed effects
        df_connected[:,:psi_hat] = df_connected[:,:lw] - df_connected[:,:alpha_hat]

        # Compute average industry fixed effects
        df_connected = transform(groupby(df_connected, :j), :psi_hat => mean => :psi_hat)

        # Model verification
        model_verify = reg(df_connected, @formula(lw ~ alpha_hat + psi_hat), save=true);

        # Compute residuals and mse
        delta = abs(msePast - sum(residuals(model_verify, df_connected).^2))

        msePast = sum(residuals(model_verify, df_connected).^2)

        # Add count
        nIter = nIter + 1

        println("At iteration number $(nIter), the MSE is $(delta)")

    end

    return df_connected

end

df_connected = getConnectedDataSet(df)

df_connected_results = akm_estimation(df_connected)

# %%

# Limited mobility bias:
#---------------------------------------------------------------------

λ_list = [0.1 0.2 0.3 0.4 0.5 0.6]

nt_list = [5 6 8 10 15]

df,_,_ = data_generating_process(α_sd, ψ_sd, csort, csig, w_sigma, λ, nt);

df_connected = getConnectedDataSet(df)

df_connected_results = akm_estimation(df_connected)


text:
import MathProgBase
export can_solve_mip, can_solve_socp, can_solve_sdp, can_solve_exp

# NOTE: Convex has been tested with the following solvers:
#   ECOS:   ECOSSolver
#   SCS:    SCSSolver
#   Gurobi: GurobiSolver
#   Mosek:  MosekSolver
#   GLPKMathProgInterface: GLPKSolverMIP
# It has not been tested with other solvers such a CPLEX.

function can_solve_mip(solver)
    name = typeof(solver).name.name
    if name == :GurobiSolver || name == :MosekSolver || name == :GLPKSolverMIP || name == :CPLEXSolver || name == :CbcSolver
        return true
    else
        @info "$name cannot solve mixed integer programs. Consider using Gurobi, Mosek, or GLPK."
        return false
    end
end

function can_solve_socp(solver)
    if :SOC in MathProgBase.supportedcones(solver)
        return true
    else
        name = typeof(solver).name.name
        @info "$name cannot solve second order cone programs. Consider using SCS, ECOS, Mosek, or Gurobi."
        return false
    end
end

function can_solve_exp(solver)
    if :ExpPrimal in MathProgBase.supportedcones(solver)
        return true
    else
        name = typeof(solver).name.name
        @info "$name cannot solve exponential programs. Consider using SCS or ECOS."
        return false
    end
end

function can_solve_sdp(solver)
    if :SDP in MathProgBase.supportedcones(solver)
        return true
    else
        name = typeof(solver).name.name
        @info "$name cannot solve semidefinite programs. Consider using SCS or Mosek."
        return false
    end
end


text:
# ---
# title: 900. RLE Iterator
# id: problem900
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/rle-iterator/description/>
# hidden: true
# ---
#
# Write an iterator that iterates through a run-length encoded sequence.
#
# The iterator is initialized by `RLEIterator(int[] A)`, where `A` is a run-
# length encoding of some sequence.  More specifically, for all even `i`, `A[i]`
# tells us the number of times that the non-negative integer value `A[i+1]` is
# repeated in the sequence.
#
# The iterator supports one function: `next(int n)`, which exhausts the next `n`
# elements (`n >= 1`) and returns the last element exhausted in this way.  If
# there is no element left to exhaust, `next` returns `-1` instead.
#
# For example, we start with `A = [3,8,0,9,2,5]`, which is a run-length encoding
# of the sequence `[8,8,8,5,5]`.  This is because the sequence can be read as
# "three eights, zero nines, two fives".
#
#
#
# **Example 1:**
#
#
#
#     Input: ["RLEIterator","next","next","next","next"], [[[3,8,0,9,2,5]],[2],[1],[1],[2]]
#     Output: [null,8,8,5,-1]
#     Explanation:
#     RLEIterator is initialized with RLEIterator([3,8,0,9,2,5]).
#     This maps to the sequence [8,8,8,5,5].
#     RLEIterator.next is then called 4 times:
#
#     .next(2) exhausts 2 terms of the sequence, returning 8.  The remaining sequence is now [8, 5, 5].
#
#     .next(1) exhausts 1 term of the sequence, returning 8.  The remaining sequence is now [5, 5].
#
#     .next(1) exhausts 1 term of the sequence, returning 5.  The remaining sequence is now [5].
#
#     .next(2) exhausts 2 terms, returning -1.  This is because the first term exhausted was 5,
#     but the second term did not exist.  Since the last term exhausted does not exist, we return -1.
#
#
#
# **Note:**
#
#   1. `0 <= A.length <= 1000`
#   2. `A.length` is an even integer.
#   3. `0 <= A[i] <= 10^9`
#   4. There are at most `1000` calls to `RLEIterator.next(int n)` per test case.
#   5. Each call to `RLEIterator.next(int n)` will have `1 <= n <= 10^9`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using LowRankModels
# tests basic functionality of glrm.jl
srand(1);
m,n,k,s = 1000,1000,5,1000*1000;
# matrix to encode
X_real, Y_real = randn(m,k), randn(k,n);
A = X_real*Y_real;
losses = fill(QuadLoss(),n)
rx, ry = ZeroReg(), ZeroReg();
glrm = GLRM(A,losses,rx,ry,5, scale=false, offset=false, X=randn(k,m), Y=randn(k,n));

p = Params(1, max_iter=200, abs_tol=0.0000001, min_stepsize=0.001)
@time X,Y,ch = fit!(glrm, params=p);
Ah = X'*Y;
p.abs_tol > abs(vecnorm(A-Ah)^2 - ch.objective[end])

function validate_folds(trf,tre,tsf,tse)
	for i=1:length(trf)
		if length(intersect(Set(trf[i]), Set(tsf[i]))) > 0
			println("Error on example $i: train and test sets overlap")
		end
	end
	for i=1:length(tre)
		if length(intersect(Set(tre[i]), Set(tse[i]))) > 0
			println("Error on feature $i: train and test sets overlap")
		end
	end
	true
end

obs = LowRankModels.flatten_observations(glrm.observed_features)
folds = LowRankModels.getfolds(obs, 5, size(glrm.A)..., do_check = false)
for i in 1:length(folds)
	@assert validate_folds(folds[i]...)
end


text:
# "smart" summary and REPL printing

function Base.summary(sys::System)
    return "System (base power $(sys.basepower))"
end

function Base.show(io::IO, sys::System)
    println(io, "basepower=$(sys.basepower)")
    show(io, sys.data)
end

function Base.show(io::IO, ::MIME"text/plain", sys::System)
    println(io, "System")
    println(io, "======")
    println(io, "Base Power: $(sys.basepower)\n")
    show(io, MIME"text/plain"(), sys.data)
end

function Base.show(io::IO, ::MIME"text/html", sys::System)
    println(io, "
System
")
    println(io, "
Base Power: $(sys.basepower)
")
    show(io, MIME"text/html"(), sys.data)
end

function Base.summary(tech::DeviceParameter)
    return "$(typeof(tech))"
end

function Base.summary(arc::Arc)
    return "$(get_name(get_from(arc))) -> $(get_name(get_to(arc))): ($(typeof(arc)))"
end

function Base.show(io::IO, ::MIME"text/plain", data::PowerSystemTableData)
    println(io, "$(typeof(data)):")
    println(io, "  directory:  $(data.directory)")
    if !isnothing(data.timeseries_metadata_file)
        println(io, "  timeseries_metadata_file:  $(data.timeseries_metadata_file)")
    end
    println(io, "  basepower:  $(data.basepower)")
    for (field, df) in data.category_to_df
        print(io, "  $field:  ")
        println(io, "$(summary(df))")
    end
end


text:
### Parsing utilities

@generated function tryparse_internal{T<:TimeType, S, F}(::Type{T}, str::AbstractString, df::DateFormat{S, F}, raise::Bool=false)
    token_types = Type[dp <: DatePart ? SLOT_RULE[first(dp.parameters)] : Void for dp in F.parameters]
    N = length(F.parameters)

    types = slot_order(T)
    num_types = length(types)
    order = Vector{Int}(num_types)
    for i = 1:num_types
        order[i] = findfirst(token_types, types[i])
    end

    field_defaults = slot_defaults(T)
    field_order = tuple(order...)
    tuple_type = slot_types(T)

    # `slot_order`, `slot_defaults`, and `slot_types` return tuples of the same length
    assert(num_types == length(field_order) == length(field_defaults))

    quote
        R = Nullable{$tuple_type}
        t = df.tokens
        l = df.locale
        pos, len = start(str), endof(str)

        err_idx = 1
        Base.@nexprs $N i->val_i = 0
        Base.@nexprs $N i->(begin
            pos > len && @goto done
            nv, next_pos = tryparsenext(t[i], str, pos, len, l)
            isnull(nv) && @goto error
            val_i, pos = unsafe_get(nv), next_pos
            err_idx += 1
        end)
        pos <= len && @goto error

        @label done
        parts = Base.@ntuple $N val
        return R(reorder_args(parts, $field_order, $field_defaults, err_idx)::$tuple_type)

        @label error
        # Note: Keeping exception generation in separate function helps with performance
        raise && throw(gen_exception(t, err_idx, pos))
        return R()
    end
end

function gen_exception(tokens, err_idx, pos)
    if err_idx > length(tokens)
        ArgumentError("Found extra characters at the end of date time string")
    else
        ArgumentError("Unable to parse date time. Expected token $(tokens[err_idx]) at char $pos")
    end
end

#    reorder_args(val, idx, default, default_from)
#
# reorder elements of `val` tuple according to `idx` tuple. Use `default[i]`
# when `idx[i] == 0` or i >= default_from
#
# returns a tuple `xs` of the same length as `idx` where `xs[i]` is
# `val[idx[i]]` if `idx[i]` is non zero, `default[i]` if `idx[i]` is zero.
#
# `xs[i]` is `default[i]` for all i >= `default_from`.
#
#
function reorder_args{N}(val::Tuple, idx::NTuple{N}, default::Tuple, default_from::Integer)
    ntuple(Val{N}) do i
        if idx[i] == 0 || idx[i] >= default_from
            default[i]
        else
            val[idx[i]]
        end
    end
end

function Base.tryparse{T<:TimeType}(::Type{T}, str::AbstractString, df::DateFormat)
    nt = tryparse_internal(T, str, df, false)
    if isnull(nt)
        return Nullable{T}()
    else
        return Nullable{T}(T(unsafe_get(nt)...))
    end
end

default_format(::Type{Date}) = ISODateFormat
default_format(::Type{DateTime}) = ISODateTimeFormat

function Base.parse{T<:TimeType}(::Type{T},
                                 str::AbstractString,
                                 df::DateFormat=default_format(T))
    nt = tryparse_internal(T, str, df, true)
    T(unsafe_get(nt)...)
end

@inline function tryparsenext_base10(str::AbstractString, i::Int, len::Int, min_width::Int=1, max_width::Int=0)
    i > len && (return Nullable{Int64}(), i)
    min_pos = min_width <= 0 ? i : i + min_width - 1
    max_pos = max_width <= 0 ? len : min(i + max_width - 1, len)
    d::Int64 = 0
    @inbounds while i <= max_pos
        c, ii = next(str, i)
        if '0' <= c <= '9'
            d = d * 10 + (c - '0')
        else
            break
        end
        i = ii
    end
    if i <= min_pos
        return Nullable{Int64}(), i
    else
        return Nullable{Int64}(d), i
    end
end

@inline function tryparsenext_word(str::AbstractString, i, len, locale, maxchars=0)
    word_start, word_end = i, 0
    max_pos = maxchars <= 0 ? len : min(chr2ind(str, ind2chr(str,i) + maxchars - 1), len)
    @inbounds while i <= max_pos
        c, ii = next(str, i)
        if isalpha(c)
            word_end = i
        else
            break
        end
        i = ii
    end
    if word_end == 0
        return Nullable{SubString}(), i
    else
        return Nullable{SubString}(SubString(str, word_start, word_end)), i
    end
end

function Base.parse(::Type{DateTime}, s::AbstractString, df::typeof(ISODateTimeFormat))
    i, end_pos = start(s), endof(s)

    dm = dd = Int64(1)
    th = tm = ts = tms = Int64(0)

    nv, i = tryparsenext_base10(s, i, end_pos, 1)
    dy = isnull(nv) ? (@goto error) : unsafe_get(nv)
    i > end_pos && @goto error

    c, i = next(s, i)
    c != '-' && @goto error
    i > end_pos && @goto done

    nv, i = tryparsenext_base10(s, i, end_pos, 1, 2)
    dm = isnull(nv) ? (@goto error) : unsafe_get(nv)
    i > end_pos && @goto done

    c, i = next(s, i)
    c != '-' && @goto error
    i > end_pos && @goto done

    nv, i = tryparsenext_base10(s, i, end_pos, 1, 2)
    dd = isnull(nv) ? (@goto error) : unsafe_get(nv)
    i > end_pos && @goto done

    c, i = next(s, i)
    c != 'T' && @goto error
    i > end_pos && @goto done

    nv, i = tryparsenext_base10(s, i, end_pos, 1, 2)
    th = isnull(nv) ? (@goto error) : unsafe_get(nv)
    i > end_pos && @goto done

    c, i = next(s, i)
    c != ':' && @goto error
    i > end_pos && @goto done

    nv, i = tryparsenext_base10(s, i, end_pos, 1, 2)
    tm = isnull(nv) ? (@goto error) : unsafe_get(nv)
    i > end_pos && @goto done

    c, i = next(s, i)
    c != ':' && @goto error
    i > end_pos && @goto done

    nv, i = tryparsenext_base10(s, i, end_pos, 1, 2)
    ts = isnull(nv) ? (@goto error) : unsafe_get(nv)
    i > end_pos && @goto done

    c, i = next(s, i)
    c != '.' && @goto error
    i > end_pos && @goto done

    nv, j = tryparsenext_base10(s, i, end_pos, 1, 3)
    tms = isnull(nv) ? (@goto error) : unsafe_get(nv)
    tms *= 10 ^ (3 - (j - i))

    j > end_pos || @goto error

    @label done
    return DateTime(dy, dm, dd, th, tm, ts, tms)

    @label error
    throw(ArgumentError("Invalid DateTime string"))
end


text:
# Main file for solving the naive SAA problem

using Distributions, StatsFuns, StatsBase, JuMP, Gurobi

const gurobi_env = Gurobi.Env()

# read these as inputs from the corresponding "instance file"
modRepNum = parse(ARGS[1])
degNum = parse(ARGS[2])
dataRepNum = parse(ARGS[3])
sampSizeNum = parse(ARGS[4])

#*========= MODEL INFORMATION ==========
const caseNum = 1
const paramsFile = "params_case" * string(caseNum) * "_nsaa.jl"
#*======================================

#*========= INCLUDE FILES ====================
include("maxParameters.jl")
include("randomSeeds.jl")

include(paramsFile)
# the next two data files are defined in paramsFile
include(modelFile)
include(dataFile)

include("solveModel.jl")
include("estimateSolnCost.jl")
include("evaluateSoln.jl")
#*======================================

# set seed for reproducibility
srand(startingSeed)

# directory name for storing results
const baseDirName = "case" * string(caseNum) * "_nsaa/" * "mod_" * string(modRepNum) * "/" * "deg_" * string(degNum) * "/"
const subDirName = baseDirName * "rep_" * string(dataRepNum) * "/"
const subDirName2 = subDirName * "samp_" * string(sampSizeNum) * "/"
mkpath(subDirName2)

#*========= PRINT OPTIONS ====================
const storeResults = true

const infoFile = "ddsp.txt"
const modelDataFile = "model_data.txt"
const numSampFile = "num_samples.txt"

const nsaaObjFile = "nsaa_obj.txt"
const nsaaDDObjFile = "nsaa_ddobj.txt"
const nsaaTimeFile = "nsaa_time.txt"
#*======================================

#*========= GENERATE MODEL PARAMETERS ====================

srand(randomSeeds_models[modRepNum])
numCovariates = 3

covariate_mean, covariate_covMat, coeff_true = generateBaseDemandModel(numCovariates)

#*========= STORE RESULTS ====================
if(storeResults)

	# write details to text file, including some key details about the test instance
	details_file = baseDirName * infoFile
	open(details_file, "w") do f
		write(f,"case number: $caseNum \n")
		write(f,"numResources: $numResources \n")
		write(f,"numCustomers: $numCustomers \n")
		write(f,"model replicate number: $modRepNum \n")
		write(f,"degree: $(degree[degNum]) \n")
		write(f,"sample sizes: $numDataSamples \n")
		write(f,"demand_errors_scaling: $demand_errors_scaling \n\n")

		write(f,"randomSeeds_MC: $randomSeeds_MC \n")
		write(f,"randomSeeds_models: $randomSeeds_models \n")
		write(f,"randomSeeds_data: $randomSeeds_data \n")
	end

	mod_data_file = baseDirName * modelDataFile
	open(mod_data_file, "w") do f
		write(f,"covariate_mean = $covariate_mean \n")
		write(f,"covariate_covMat = $covariate_covMat \n")
		write(f,"coeff_true = $coeff_true \n")
	end

	samp_size_file = baseDirName * numSampFile
	open(samp_size_file, "w") do f
		for i = 1:length(numDataSamples)
			write(f,"$(numDataSamples[i]) \n")
		end
	end

end
#*============================================

#*========= GENERATE DATA REPLICATE ====================

srand(randomSeeds_data[dataRepNum])

covariate_obs = generateCovariateReal(numCovariates,covariate_mean,covariate_covMat)

demand_data, ~ = generateDemandData(numCovariates,numDataSamples[sampSizeNum],degree[degNum],covariate_mean,covariate_covMat,coeff_true)

#*========= SOLVE NSAA MODEL ====================

tic()

z_soln_NSAA, objDDNSAA = solveSAAModel(demand_data)

nsaaObjEstimates = estimateSolnQuality(z_soln_NSAA,covariate_obs,degree[degNum],numMCScenarios,numMCReplicates,coeff_true)

NSAATime = toq()

#*========= STORE RESULTS ====================
if(storeResults)

	# write data to text file
	obj_est_file = subDirName2 * nsaaObjFile
	open(obj_est_file, "w") do f
		for i = 1:length(nsaaObjEstimates)
			write(f,"$(nsaaObjEstimates[i]) \n")
		end
	end

	ddobj_est_file = subDirName2 * nsaaDDObjFile
	open(ddobj_est_file, "w") do f
		write(f,"$objDDNSAA")
	end

	time_obj_file = subDirName2 * nsaaTimeFile
	open(time_obj_file, "w") do f
		write(f,"$NSAATime")
	end

end
#*============================================

max_stars_repo_name: rohitkannan/DD-SAA

max_issues_repo_name: rohitkannan/DD-SAA

max_forks_repo_name: rohitkannan/DD-SAA


text:
"""
    load_graphml(filename)

Read a graphML format file with file name `filename`.
Note: only support DiGraph{Node{String}} for now.
"""
function load_graphml(filename)
    open(EzXML.StreamReader, filename) do reader
        serial = 0
        g = DiGraph{Node{String}}()
        # scan nodes
        for typ in reader
            if typ == EzXML.READER_ELEMENT
                elname = EzXML.nodename(reader)
                if elname == "node"
                    add_node!(g, reader["id"])
                elseif elname == "edge"
                    add_edge!(g, reader["source"], reader["target"])
                end
            end
        end
    return g
    end
end


text:
using Benchmark
################################################################################
# While the C++ implementation is only slightly faster than the julia version,
# my laptop is actually faster than the server......
#
# On my laptop (OS X):
#  Row | Function         | Average   | Relative | Replications |
# |-----|------------------|-----------|----------|--------------|
# | 1   | "pooling_julia"  | 0.100337  | 2.23666  | 100          |
# | 2   | "pooling_native" | 0.0448601 | 1.0      | 100          |
#
# On a linux server
# | Row | Function         | Average   | Relative | Replications |
# |-----|------------------|-----------|----------|--------------|
# | 1   | "pooling_julia"  | 0.129456  | 1.30713  | 100          |
# | 2   | "pooling_native" | 0.0990383 | 1.0      | 100          |
################################################################################

################################################################################
# OpenMP significantly improved the result.
#
# On a linux server with OMP_NUM_THREADS=16
# | Row | Function         | Average   | Relative | Replications |
# |-----|------------------|-----------|----------|--------------|
# | 1   | "pooling_julia"  | 0.12026   | 9.82832  | 100          |
# | 2   | "pooling_native" | 0.0122361 | 1.0      | 100          |
#
# On the same linux server without setting OMP_NUM_THREADS
# | Row | Function         | Average   | Relative | Replications |
# |-----|------------------|-----------|----------|--------------|
# | 1   | "pooling_julia"  | 0.130325  | 7.58476  | 100          |
# | 2   | "pooling_native" | 0.0171824 | 1.0      | 100          |
################################################################################

function pooling(input::Array, output::Array, mask::Array, kernel, pad, stride)
  width, height, channels, num = size(input)
  pooled_width  = int(ceil(float(width +2*pad[1]-kernel[1]) / stride[1]))+1
  pooled_height = int(ceil(float(height+2*pad[2]-kernel[2]) / stride[2]))+1

  for n = 1:num
    for c = 1:channels
      for ph = 1:pooled_height
        for pw = 1:pooled_width
          hstart = (ph-1)*stride[2] - pad[2] + 1
          wstart = (pw-1)*stride[1] - pad[1] + 1
          hend = min(hstart + kernel[2] - 1, height)
          wend = min(wstart + kernel[1] - 1, width)
          hstart = max(1, hstart)
          wstart = max(1, wstart)

          maxval = -Inf
          maxw = 0
          maxh = 0
          for w = wstart:wend
            for h = hstart:hend
              @inbounds val = input[w,h,c,n]
              if val > maxval
                maxval = val
                maxw = w
                maxh = h
              end
            end
          end
          @inbounds output[pw,ph,c,n] = maxval
          @inbounds mask[pw,ph,c,n] = (maxh-1) * width + maxw-1
        end
      end
    end
  end
end

library = dlopen("./libpooling.so")
func_handle = dlsym(library, :max_pooling_impl_double)
function pooling_native(input::Array, output::Array, mask::Array, kernel, pad, stride)
  width, height, channels, num = size(input)
  pooled_width  = int(ceil(float(width +2*pad[1]-kernel[1]) / stride[1]))+1
  pooled_height = int(ceil(float(height+2*pad[2]-kernel[2]) / stride[2]))+1

  ccall(func_handle, Void, (Ptr{Float64}, Ptr{Float64}, Ptr{Csize_t}, Cint, Cint, Cint, Cint,
        Cint, Cint, # pooled_width, pooled_height
        Cint, Cint, Cint, Cint, Cint, Cint, # kernel, pad, stride
      ), input, output, mask, width, height, channels, num, pooled_width, pooled_height,
      kernel[1], kernel[2], pad[1], pad[2], stride[1], stride[2])
end

input = rand(28, 28, 50, 128)
kernel = (5, 5)
pad = (2, 2)
stride = (2, 2)
width, height, channels, num = size(input)
pooled_width  = int(ceil(float(width +2*pad[1]-kernel[1]) / stride[1]))+1
pooled_height = int(ceil(float(height+2*pad[2]-kernel[2]) / stride[2]))+1

output1 = zeros(pooled_width, pooled_height, channels, num)
output2 = zeros(pooled_width, pooled_height, channels, num)

mask1 = zeros(Csize_t, size(output1))
mask2 = zeros(Csize_t, size(output2))

pooling_julia() = pooling(input, output1, mask1, kernel, pad, stride)
pooling_native() = pooling_native(input, output2, mask2, kernel, pad, stride)

# make sure results are correct
pooling_julia()
pooling_native()
@assert all(abs.(output1-output2) .< 1e-10)
@assert mask1 == mask2

# compare performance
df = compare([pooling_julia, pooling_native], 100)
println("$df")


text:
############################################################################################
struct SMC2Constructor{W<:AbstractConstructor,J<:AbstractConstructor} <:AbstractConstructor
    "Constructor for (online) state propagation over time."
    propagation::W
    "Constructor for jittering step."
    jitter::J
    function SMC2Constructor(propagation::W, jitter::J) where {W<:AbstractConstructor,J<:AbstractConstructor}
        return new{W,J}(propagation, jitter)
    end
end
function get_sym(constructor::SMC2Constructor)
#    return get_sym(constructor.jitter)
    sym1 = get_sym(constructor.jitter)
    sym2 = get_sym(constructor.propagation)
    #Return unique symbols as tuple
    sym = unique((sym1..., sym2...))
    return Tuple(sym)
end

function SMC2(
    propagation::J, jitter::W; kwargs...
) where {J<:AbstractConstructor, W<:AbstractConstructor}
    return SMC(SMC2Constructor(propagation, jitter); kwargs...)
end

############################################################################################
"""
$(TYPEDEF)

SMC2 Kernel, consisting of a kernel for the particle trajectory and kernel for all model parameter.

# Fields
$(TYPEDFIELDS)
"""
struct SMC2Kernel{P<:ParticleFilter,M<:PMCMC} <: BaytesCore.AbstractAlgorithm
    "Kernel for (online) state propagation over time."
    pf::P   #θ dimension
    "Kernel for jittering step."
    pmcmc::M   #particle dimension
    function SMC2Kernel(pf::P, pmcmc::M) where {P<:ParticleFilter,M<:PMCMC}
        @argcheck isa(pf.tune.referencing, Marginal) "Cannot use Conditional Filter in θ dimension, use Marginal referencing instead"
        return new{P,M}(pf, pmcmc)
    end
end

function result!(kernel::SMC2Kernel, result::L) where {L<:ℓObjectiveResult}
    result!(kernel.pmcmc, result)
    return nothing
end

function get_result(kernel::SMC2Kernel)
    return get_result(kernel.pmcmc)
end
function predict(_rng::Random.AbstractRNG, kernel::SMC2Kernel, objective::Objective)
    return predict(_rng, kernel.pf, objective)
end

function infer(
    _rng::Random.AbstractRNG, kernel::SMC2Kernel, model::ModelWrapper, data::D
) where {D}
    return BaytesCore.infer(_rng, kernel.pf, model, data)
end
function infer(
    _rng::Random.AbstractRNG,
    diagnostics::Type{AbstractDiagnostics},
    kernel::SMC2Kernel,
    model::ModelWrapper,
    data::D,
) where {D}
    return BaytesCore.infer(_rng, diagnostics, kernel.pmcmc, model, data)
end

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

function jitterkernel(particles::SMCParticles{<:SMC2Kernel}, iter::Int64)
    return particles.kernel[iter].pmcmc
end

function SMCweight(_rng::Random.AbstractRNG, objective::Objective, algorithm::SMC2Kernel, cumweightsₜ₋₁)
    #!NOTE: Use cumulative weight from propagated pf, taking into account new data
    cumweightsₜ = objective.temperature * algorithm.pf.particles.ℓobjective.cumulative #get_ℓweight(algorithm.pf)
    #!NOTE: cumweightsₜ₋₁ already has correct temperature at t-1
    return cumweightsₜ, cumweightsₜ - cumweightsₜ₋₁
end
function SMCreweight(_rng::Random.AbstractRNG, objective::Objective, algorithm::SMC2Kernel, cumweightsₜ₋₁)
    #!NOTE: Run Particle Filter wrt to jittered model parameter, and compute new weight
    propose!(
        _rng,
        algorithm.pf,
        objective.model,
        objective.data,
        objective.temperature,
        BaytesCore.UpdateTrue()
    )
    cumweightsₜ = objective.temperature * algorithm.pf.particles.ℓobjective.cumulative #get_ℓweight(algorithm.pf)
    #!NOTE: cumweightsₜ₋₁ already has correct temperature at t-1
    return cumweightsₜ, cumweightsₜ - cumweightsₜ₋₁
end

############################################################################################
"""
$(SIGNATURES)
Initialize SMC² sampler.

# Examples
```julia
```

"""
function SMCParticles(
    _rng::Random.AbstractRNG,
    kernel::SMC2Constructor,
    objective::Objective,  #Initial Model with prior (Parameter values are overwritten by prior samples)
    info::BaytesCore.SampleDefault,
    tune::SMCTune
)
    @unpack propagation, jitter = kernel
    ## Initialize individual PF and PMCMC algorithm
    @unpack model, data, temperature = objective
    @unpack Ntuning = tune
    Nchains = tune.chains.Nchains
    modelᵗᵉᵐᵖ = deepcopy(model)
    tagged = Tagged(model, propagation.sym)
    ModelWrappers.fill!(
        modelᵗᵉᵐᵖ,
        tagged,
        BaytesCore.to_NamedTuple(
            keys(tagged.parameter),
            getfield(modelᵗᵉᵐᵖ.val, keys(tagged.parameter)[1])[1:maximum(size(data))],
        ),
    )
    ## Initialize individual Models ~ Assign copies for  Model and Algorithms ~ necessary in initiation so no pointer error
    modelᵛ = [deepcopy(modelᵗᵉᵐᵖ) for _ in Base.OneTo(Nchains)]
    algorithmᵛ = [
        SMC2Kernel(
            propagation(_rng, modelᵛ[iter], data, temperature, info),
            jitter(_rng, modelᵛ[iter], data, temperature, info),
        ) for iter in Base.OneTo(Nchains)
    ]
    ## Assign weights
    weights = BaytesCore.ParameterWeights(Nchains)
    ## Assign buffer for inplace resampling
    buffer = SMCBuffer(
        _rng, algorithmᵛ[1].pmcmc, modelᵗᵉᵐᵖ, data, Nchains
    )
    ## Loop through all models
    #Polyester.@batch per=thread minbatch=tune.batchsize for iter in eachindex(algorithmᵛ)
    Base.Threads.@threads for iter in eachindex(algorithmᵛ)
        ## Tune PMCMC algorithm
        propose!(_rng, algorithmᵛ[iter].pmcmc, modelᵛ[iter], data, temperature, BaytesCore.UpdateTrue())
        for _ in 2:Ntuning
            propose!(_rng, algorithmᵛ[iter].pmcmc, modelᵛ[iter], data, temperature, tune.capture)
        end
        ## Assign a Particle Filter for tuned theta
        propose!(_rng, algorithmᵛ[iter].pf, modelᵛ[iter], data, temperature, BaytesCore.UpdateTrue())
        ## Calculate initial weights
        _objective = Objective(modelᵛ[iter], objective.data, tune.tagged, objective.temperature)
        buffer.cumweights[iter], weights.buffer[iter] = SMCweight(
            _rng,
            _objective,
            algorithmᵛ[iter],
            _objective(_objective.model.val)
        )
    end
    ## Normalize weights
    weights(weights.buffer)
    ## Return container
    return SMCParticles(modelᵛ, algorithmᵛ, weights, buffer)
end

############################################################################################
"""
$(SIGNATURES)
Propagate data forward over time.

# Examples
```julia
```

"""
function propagate!(_rng::Random.AbstractRNG, particles::SMCParticles{<:SMC2Kernel}, tune::SMCTune, data::D, temperature::F) where {D, F<:AbstractFloat}
    #!NOTE: Here, smc.buffer.jitterdiagnostics will not be updated, and for SMCDiagnostics, last available jitterdiagnostics are provided.
    ## Propagate series forward with recent particle
    #Polyester.@batch per=thread minbatch=tune.batchsize for iter in eachindex(particles.model)
    Base.Threads.@threads for iter in eachindex(particles.model)
        #_, diagnostics = propagate!(
        propagate!(
            _rng,
            particles.kernel[iter].pf,
            particles.model[iter],
            data,
            temperature
        )
        #particles.buffer.predictions[iter] = diagnostics.base.prediction
    end
    return nothing
end

############################################################################################
#export
export
    SMC2,
    SMC2Constructor,
    SMC2Kernel,
    SMCweight,
    resample!,
    propagate!,
    weight!


text:
module GraphBasedSystems

using SparseArrays
using StaticArrays
using LightGraphs

export System,
    full_matrix,
    full_vector,
    children,
    connections,
    parents,

    ldu_solve!

include(joinpath("util", "custom_static.jl"))

include(joinpath("system", "entry.jl"))
include(joinpath("system", "system.jl"))
include(joinpath("system", "setup_functions.jl"))

include(joinpath("solvers", "lu.jl"))
include(joinpath("solvers", "llt.jl"))
include(joinpath("solvers", "ldlt.jl"))
include(joinpath("solvers", "ldu.jl"))

end


text:
function aca_compression(
    matrix::Function,
    rowindices,
    colindices;
    tol=1e-14,
    T=ComplexF64,
    maxrank=40,
    svdrecompress=true
)
    function smartmaxlocal(roworcolumn, acausedindices)
        maxval = -1
        index = -1
        for i=1:length(roworcolumn)
            if !acausedindices[i]
                if abs(roworcolumn[i]) > maxval
                    maxval = abs(roworcolumn[i])
                    index = i
                end
            end
        end
        return index, maxval
    end

    acausedrowindices = zeros(Bool,length(rowindices))
    acausedcolumnindices = zeros(Bool,length(colindices))
    acarowindicescounter = 1
    acacolumnindicescounter = 1

    nextrowindex = 1
    #acarowindices = [nextrowindex]
    acausedrowindices[nextrowindex] = true

    U = zeros(T, length(rowindices), maxrank)
    V = zeros(T, maxrank, length(colindices))

    i = 1
    #acacolumnindices = Integer[]

    @views matrix(V[acarowindicescounter:acarowindicescounter, :],
                    rowindices[nextrowindex:nextrowindex],
                    colindices[:])

    @views nextcolumnindex, maxval = smartmaxlocal(V[acarowindicescounter, :], acausedcolumnindices)
    acausedcolumnindices[nextcolumnindex] = true

    #push!(acacolumnindices, nextcolumnindex)

    @views V[acarowindicescounter:acarowindicescounter, :] /= V[acarowindicescounter, nextcolumnindex]

    @views matrix(
        U[:, acacolumnindicescounter:acacolumnindicescounter],
        rowindices,
        colindices[nextcolumnindex:nextcolumnindex]
    )

    @views normUVlastupdate = norm(U[:, 1])*norm(V[1, :])
    normUVsqared = normUVlastupdate^2

    while normUVlastupdate > sqrt(normUVsqared)*tol &&
        i <= length(rowindices)-1 &&
        i <= length(colindices)-1 &&
        acacolumnindicescounter < maxrank

        i += 1

        @views nextrowindex, maxval = smartmaxlocal(U[:,acacolumnindicescounter], acausedrowindices)

        #if nextrowindex == -1
        #    error("Failed to find new row index: ", nextrowindex)
        #end

        #if isapprox(maxval, 0.0)
        #    println("Future V entry is close to zero. Abort.")
        #    return U[:,1:acacolumnindicescounter], V[1:acarowindicescounter, :]
        #end

        acausedrowindices[nextrowindex] = true

        acarowindicescounter += 1
        @views matrix(
            V[acarowindicescounter:acarowindicescounter, :],
            rowindices[nextrowindex:nextrowindex],
            colindices[:]
        )

        @views V[acarowindicescounter:acarowindicescounter, :] -= U[nextrowindex:nextrowindex, 1:acacolumnindicescounter]*V[1:(acarowindicescounter-1), :]
        @views nextcolumnindex, maxval = smartmaxlocal(V[acarowindicescounter, :], acausedcolumnindices)

        if (isapprox(V[acarowindicescounter, nextcolumnindex],0.0))
            normUVlastupdate = 0.0
            V[acarowindicescounter:acarowindicescounter, :] .= 0.0
            acarowindicescounter -= 1
            println("Matrix seems to have exact rank: ", acarowindicescounter)
        else
            #push!(acarowindices, nextrowindex)

            @views V[acarowindicescounter:acarowindicescounter, :] /= V[acarowindicescounter, nextcolumnindex]

            #if nextcolumnindex == -1
            #    error("Failed to find new column index: ", nextcolumnindex)
            #end

            #if isapprox(maxval, 0.0)
            #    println("Future U entry is close to zero. Abort.")
            #    return U[:,1:acacolumnindicescounter], V[1:acarowindicescounter-1,:]
            #end

            acausedcolumnindices[nextcolumnindex] = true

            #push!(acacolumnindices, nextcolumnindex)

            acacolumnindicescounter += 1
            @views matrix(U[:, acacolumnindicescounter:acacolumnindicescounter],
                            rowindices,
                            colindices[nextcolumnindex:nextcolumnindex])

            @views U[:, acacolumnindicescounter] -= U[:, 1:(acacolumnindicescounter-1)]*V[1:(acarowindicescounter-1),nextcolumnindex:nextcolumnindex]

            @views normUVlastupdate = norm(U[:,acacolumnindicescounter])*norm(V[acarowindicescounter,:])

            normUVsqared += normUVlastupdate^2
            for j = 1:(acacolumnindicescounter-1)
                @views normUVsqared += 2*abs(dot(U[:,acacolumnindicescounter], U[:,j])*dot(V[acarowindicescounter,:], V[j,:]))
            end
        end
    end

    if acacolumnindicescounter == maxrank
        println("WARNING: aborted ACA after maximum allowed rank")
    end

    if svdrecompress && acacolumnindicescounter > 1
        @views Q,R = qr(U[:,1:acacolumnindicescounter])
        @views U,s,V = svd(R*V[1:acarowindicescounter,:])

        opt_r = length(s)
        for i in eachindex(s)
            if s[i] < tol*s[1]
                opt_r = i
                break
            end
        end

        A = (Q*U)[:,1:opt_r]
        B = (diagm(s)*V')[1:opt_r,:]
        return A, B
    else
        return U[:,1:acacolumnindicescounter], V[1:acarowindicescounter,:]
    end
end

function aca_compression(matrix::Function, testnode::BoxTreeNode, sourcenode::BoxTreeNode;
                        tol=1e-14, T=ComplexF64, maxrank=40, svdrecompress=true)
    U, V = aca_compression(matrix, testnode.data, sourcenode.data,
                             tol=tol, T=T, maxrank=maxrank, svdrecompress=svdrecompress)

    if false == true && size(U,2) >= 30
        println("Rank of compressed ACA is too large")
        println("Size of U: ", size(U))
        println("   Testnode box")
        println("    - Center: ", testnode.boundingbox.center)
        println("    - Halflength: ", testnode.boundingbox.halflength)
        println("    - Nodes: ", length(testnode.data))
        println("    - Level: ", testnode.level)
        println("   Sourcenode box")
        println("    - Center: ", sourcenode.boundingbox.center)
        println("    - Halflength: ", sourcenode.boundingbox.halflength)
        println("    - Nodes: ", length(sourcenode.data))
        println("    - Level: ", sourcenode.level)
        #UU, SS, VV = svd(U*V)
        #k = 1
        #while k < length(SS) && SS[k] > SS[1]*tol
        #    k += 1
        #end
        #println("   Actual rank: ", k)
        #error("emergency stop")
    end

    return U, V
end

max_stars_repo_name: sbadrian/FastBEAST

max_issues_repo_name: sbadrian/FastBEAST

max_forks_repo_name: sbadrian/FastBEAST


text:
function get_current_state(obj::AbstractODEStateInterpolator)
    return jcall(obj, "getCurrentState", ODEStateAndDerivative, ())
end

function get_global_current_state(obj::AbstractODEStateInterpolator)
    return jcall(obj, "getGlobalCurrentState", ODEStateAndDerivative, ())
end

function get_global_previous_state(obj::AbstractODEStateInterpolator)
    return jcall(obj, "getGlobalPreviousState", ODEStateAndDerivative, ())
end

function get_interpolated_state(obj::AbstractODEStateInterpolator, arg0::jdouble)
    return jcall(obj, "getInterpolatedState", ODEStateAndDerivative, (jdouble,), arg0)
end

function get_previous_state(obj::AbstractODEStateInterpolator)
    return jcall(obj, "getPreviousState", ODEStateAndDerivative, ())
end

function is_current_state_interpolated(obj::AbstractODEStateInterpolator)
    return jcall(obj, "isCurrentStateInterpolated", jboolean, ())
end

function is_forward(obj::AbstractODEStateInterpolator)
    return jcall(obj, "isForward", jboolean, ())
end

function is_previous_state_interpolated(obj::AbstractODEStateInterpolator)
    return jcall(obj, "isPreviousStateInterpolated", jboolean, ())
end

function restrict_step(obj::AbstractODEStateInterpolator, arg0::ODEStateAndDerivative, arg1::ODEStateAndDerivative)
    return jcall(obj, "restrictStep", AbstractODEStateInterpolator, (ODEStateAndDerivative, ODEStateAndDerivative), arg0, arg1)
end

function taylor(::Type{AdamsStateInterpolator}, arg0::EquationsMapper, arg1::ODEStateAndDerivative, arg2::jdouble, arg3::jdouble, arg4::Vector{jdouble}, arg5::Array2DRowRealMatrix)
    return jcall(AdamsStateInterpolator, "taylor", ODEStateAndDerivative, (EquationsMapper, ODEStateAndDerivative, jdouble, jdouble, Vector{jdouble}, Array2DRowRealMatrix), arg0, arg1, arg2, arg3, arg4, arg5)
end


text:
"""
    struct GridPortion{Dc,Dp,G} <: Grid{Dc,Dp}
      oldgrid::G
      cell_to_oldcell::Vector{Int}
      node_to_oldnode::Vector{Int}
    end
"""
struct GridPortion{Dc,Dp,G} <: Grid{Dc,Dp}
  oldgrid::G
  cell_to_oldcell::Vector{Int}
  node_to_oldnode::Vector{Int}
  cell_to_nodes::Table{Int,Vector{Int},Vector{Int32}}
  @doc """
      GridPortion(oldgrid::Grid{Dc,Dp},cell_to_oldcell::Vector{Int}) where {Dc,Dp}
  """
  function GridPortion(oldgrid::Grid{Dc,Dp},cell_to_oldcell::Vector{Int}) where {Dc,Dp}

    oldcell_to_oldnodes = get_cell_nodes(oldgrid)
    noldnodes = num_nodes(oldgrid)
    oldnode_to_coords = get_node_coordinates(oldgrid)

    node_to_oldnode, oldnode_to_node = _find_active_nodes(oldcell_to_oldnodes,cell_to_oldcell,noldnodes)
    cell_to_nodes = _renumber_cell_nodes(oldcell_to_oldnodes,oldnode_to_node,cell_to_oldcell)

    new{Dc,Dp,typeof(oldgrid)}(oldgrid,cell_to_oldcell,node_to_oldnode,cell_to_nodes)
  end
end

function OrientationStyle(::Type{GridPortion{Dc,Dp,G}}) where {Dc,Dp,G}
  OrientationStyle(G)
end

function RegularityStyle(::Type{GridPortion{Dc,Dp,G}}) where {Dc,Dp,G}
  RegularityStyle(G)
end

function get_node_coordinates(grid::GridPortion)
  oldnode_to_coords = get_node_coordinates(grid.oldgrid)
  reindex(oldnode_to_coords,grid.node_to_oldnode)
end

function get_cell_nodes(grid::GridPortion)
  grid.cell_to_nodes
end

function get_reffes(grid::GridPortion)
  get_reffes(grid.oldgrid)
end

function get_cell_type(grid::GridPortion)
  reindex(get_cell_type(grid.oldgrid),grid.cell_to_oldcell)
end

# Helpers

function _find_active_nodes(oldcell_to_oldnodes,cell_to_oldcell,noldnodes)
  oldnode_is_active = fill(false,noldnodes)
  cache = array_cache(oldcell_to_oldnodes)
  for oldcell in cell_to_oldcell
    oldnodes = getindex!(cache,oldcell_to_oldnodes,oldcell)
    for oldnode in oldnodes
      oldnode_is_active[oldnode] = true
    end
  end
  node_to_oldnode = findall(oldnode_is_active)
  oldnode_to_node = fill(UNSET,noldnodes)
  oldnode_to_node[node_to_oldnode] = 1:length(node_to_oldnode)
  (node_to_oldnode, oldnode_to_node)
end

function _renumber_cell_nodes(oldcell_to_oldnodes,oldnode_to_node,cell_to_oldcell)
  ncells = length(cell_to_oldcell)
  cell_to_nodes_ptrs = zeros(Int32,ncells+1)
  cache = array_cache(oldcell_to_oldnodes)
  for (cell,oldcell) in enumerate(cell_to_oldcell)
    oldnodes = getindex!(cache,oldcell_to_oldnodes,oldcell)
    cell_to_nodes_ptrs[cell+1] = length(oldnodes)
  end
  length_to_ptrs!(cell_to_nodes_ptrs)
  ndata = cell_to_nodes_ptrs[end]-1
  cell_to_nodes_data = zeros(Int,ndata)
  for (cell,oldcell) in enumerate(cell_to_oldcell)
    oldnodes = getindex!(cache,oldcell_to_oldnodes,oldcell)
    a = cell_to_nodes_ptrs[cell]-1
    for (lnode,oldnode) in enumerate(oldnodes)
      node = oldnode_to_node[oldnode]
      @assert node > 0
      cell_to_nodes_data[a+lnode] = node
    end
  end
  Table(cell_to_nodes_data,cell_to_nodes_ptrs)
end