text:
export MvNormalLH, MvNormalLHSwitches
export means, stds, n_vars
export cov_matrix, check_cov_matrix, draw_from_weights, check_weight_matrix
export conditional_distrib, cond_mean_weights

"""
	$(SIGNATURES)

Switches used to construct Multivariate Normal object.
"""
mutable struct MvNormalLHSwitches{T1} <: AbstractDistributionSwitches{T1}
    meanV :: Vector{T1}
    stdV :: Vector{T1}
end

"""
	$(SIGNATURES)

Multivariation Normal object.

This is a `ModelObject` but currently does not support calibrating parameters.
"""
mutable struct MvNormalLH{T1} <: AbstractDistributionLH{T1}
    switches :: MvNormalLHSwitches{T1}
end

MvNormalLH(meanV :: AbstractVector{T1}, stdV :: AbstractVector{T1}) where T1 =
    MvNormalLH(MvNormalLHSwitches(meanV, stdV));

Base.show(io :: IO, switches :: MvNormalLHSwitches{T1}) where T1 =
    print(io, typeof(switches), " with means $(switches.meanV) and stds $(switches.stdV)");

Base.show(io :: IO, m :: MvNormalLH{T1}) where T1 =
    print(io, typeof(m), " with means $(means(m)) and stds $(stds(m))");

means(m :: MvNormalLH{T1}) where T1 = m.switches.meanV;
stds(m :: MvNormalLH{T1}) where T1 = m.switches.stdV;
n_vars(m :: MvNormalLH{T1}) where T1 = length(means(m));
ModelParams.has_pvector(m :: MvNormalLH{T1}) where T1 = false;

# Make a Distributions.MvNormal that one can sample from
function make_mv_normal(m :: MvNormalLH{T1}, wtM :: AbstractMatrix{T1}) where T1
    return MvNormal(means(m), cov_matrix(m, wtM));
end

"""
	$(SIGNATURES)

Make cov matrix from weight matrix.

The weight matrix is lower triangular with ones on the diagonal.
"""
function cov_matrix(m :: MvNormalLH{T1}, wtM :: AbstractMatrix{T1}) where T1
    @assert check_weight_matrix(m, wtM)
    stdV = stds(m);

    # Scale the weight matrix so that it is the Cholesky decomposition of the Cov matrix
    wt2M = similar(wtM);
    for i1 = 1 : n_vars(m)
       wt2M[i1, :] = wtM[i1, :] ./ sqrt(sum(wtM[i1, :] .^ 2)) .* stdV[i1];
    end

    covMatM = wt2M * (wt2M');
    @assert check_cov_matrix(m, covMatM);
    return covMatM
end

"""
	$(SIGNATURES)

Check that weight matrix is valid.
"""
function check_weight_matrix(m :: MvNormalLH{T1}, wtM :: AbstractMatrix{T1}) where T1
    stdV = stds(m);
    n = length(stdV);

    isValid = true;
    isequal(size(wtM), (n, n))  ||  (isValid = false);
    # Lower triangular?
    isValid = isValid  &&  istril(wtM);
    # Diagonal == 1?
    isValid = isValid  &&  all(isapprox.(diag(wtM), one(T1)));
    # If std deviation is 0, variable must have 0 weights
    for i1 = 1 : n
        if stdV[i1] == zero(T1)
            idxV = trues(n);
            idxV[i1] = false;
            isValid = isValid  &&  all(wtM[i1, idxV] .== zero(T1))  &&
                all(wtM[idxV, i1] .== zero(T1));
        end
    end
    return isValid
end

"""
	$(SIGNATURES)

Check that a covariance matrix is valid for a given `MvNormalLH` object.
"""
function check_cov_matrix(m :: MvNormalLH{T1}, covM :: AbstractMatrix{T1}) where T1
    isValid = true;
    isValid = isValid  && (size(covM) == (n_vars(m), n_vars(m)));
    std2V = sqrt.(diag(covM));
    isValid = isValid  && isapprox(std2V, stds(m), atol = 0.001);
    return isValid
end

"""
	$(SIGNATURES)

Draw random variables given the weight matrix. Returns Matrix by [observation, variable] (unlike Distributions.jl).
"""
function draw_from_weights(m :: MvNormalLH{T1},
    wtM :: AbstractMatrix{T1}, nObs :: Integer, rng :: AbstractRNG) where T1

    mvn = make_mv_normal(m, wtM);
    return rand(rng, mvn, nObs)';
end

## -----------  Conditional distribution

"""
    $(SIGNATURES)

Compute conditional distribution of a set of variables, given the others
For Multiple sets of conditioning observations

# Arguments
- covM
    covariance matrix
- idx2V
    indices of variables on which we condition
- value2M[observation, variable]
    their values

# Outputs
- condMeanM[observation, variable]
- condStdV(variable)
    conditional means and std of each variable, given all others
    for those not in idx2V
- condCovM(variable, variable)
    conditional covariance; for those not in idx2
"""
function conditional_distrib(m :: MvNormalLH{T1},
    idx2V :: AbstractVector{I1}, value2M :: AbstractMatrix{T1},
    covM) where {I1 <: Integer, T1}

    # Input check
    nObs = size(value2M, 1);
    nCond = length(idx2V);

    @assert all(1 .<= idx2V .<= n_vars(m))
    @assert size(value2M) == (nObs, nCond)
    @assert check_cov_matrix(m, covM)

    # Indices of unknowns
    idx1V = [j  for j = 1 : n_vars(m)  if !(j ∈ idx2V)];

    sigma11 = covM[idx1V, idx1V];
    sigma12 = covM[idx1V, idx2V];
    sigma21 = covM[idx2V, idx1V];
    sigma22 = covM[idx2V, idx2V];
    mu1V = means(m)[idx1V];
    mu2V = means(m)[idx2V];

    condCovM  = sigma11 .- (sigma12 / sigma22) * sigma21;
    condStdV  = sqrt.(diag(condCovM));

    # This is equivalent to mu1 + sigma12 / sigma22 * (value2 - mu2)
    condMeanM = mu1V' .+ (value2M .- mu2V') / sigma22 * sigma21;

    # Testing against the original formula; tests matrix expansion
    # This is directly tested by simulation in the test function
    # for i1 = 1 : nObs
    #    condMeanV = (mu1V' .+ ((value2M[i1,:] .- mu2V)' / sigma22) * sigma21)';
    #    @assert (all(abs.(condMeanV .- condMeanM[i1,:]) .< 0.0001));
    # end

    # Output check
    n1 = length(idx1V);
    @assert size(condMeanM) == (nObs, n1)
    @assert size(condStdV) == (n1,)
    @assert size(condCovM) == (n1, n1)

    return condMeanM, condStdV, condCovM
end

"""
    $(SIGNATURES)

Weights of conditioning variables for conditional means: `sigma12/sigma22`.
Such that

    `conditional mean = mu1 + weights * (value2 - mu2)`

# Outputs
- wtM
    n1 x n2 matrix of weights
    each row gives the weights for a "dependent" variable
"""
function cond_mean_weights(m :: MvNormalLH{T1},
    idx2V, covM :: AbstractMatrix{T1}) where T1

    # Input check
    nVars = size(covM, 1);
    nCond = length(idx2V);

    @assert all(1 .<= idx2V .<= nVars);
    @assert check_cov_matrix(m, covM);

    # Indices of unknowns
    idx1V = [j  for j = 1 : n_vars(m)  if !(j ∈ idx2V)];

    sigma12 = covM[idx1V, idx2V];
    sigma22 = covM[idx2V, idx2V];

    # This is equivalent to mu1 + sigma12 / sigma22 * (value2 - mu2)
    wtM = sigma12 / sigma22;
    return wtM
end

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


text:
mutable struct KetState{T<:RLorCX}
    lx :: Int
    d  :: Int
    v  :: Vector{T}
end

function randket(T::Type, rng::AbstractRNG, lx::Int, d::Int)# where {T<:RLorCX}
    v = randn(rng, T, d^lx)
    KetState{T}(lx, d, v)
end

randket(rng::AbstractRNG, lx::Int, d::Int) = randket(Float64, rng, lx, d)

normalize!(v::KetState{T}) where {T<:RLorCX} = normalize!(v.v)

@inline opextend(op, lx::Int, l::Int) = eye(2^(l-1)) ⊗ op ⊗ eye(2^(lx-l))

function apply1site(v::KetState{T}, l::Int, op::SparseMatrixCSC{T,Int}) where {T<:RLorCX}
    size(op) != (v.d, v.d) && error("Operator size not correct!")
    (0 < l <= v.lx) || error("position l is not within the lattice")

    KetState(v.lx, v.d, opextend(op, v.lx, l) * v.v)
end

apply1site(v::KetState{ComplexF64}, l::Int, op::SparseMatrixCSC{Float64, Int}) =
    apply1site(v,l, convert(SparseMatrixCSC{ComplexF64,Int}, op))

function measure1point(v::KetState{T}, op::SparseMatrixCSC{T,Int}) where {T<:RLorCX}
    result = Vector{T}(undef, v.lx)
    for l in 1:v.lx
        result[l] = dot(v.v, opextend(op, v.lx, l) * v.v)
    end
    result
end

measure1point(v::KetState{ComplexF64}, op::SparseMatrixCSC{Float64, Int}) =
    measure1point(v, convert(SparseMatrixCSC{ComplexF64,Int}, op))

ext: jl

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
using CRNT
using Test
using Nemo

@testset "Cones.jl" begin
    @test raysof(cone_positivenullspace([0 0; 0 0])) == [1 0; 0 1]
    @test raysof(cone_positivenullspace([1 -1; 2 -2])) == transpose([1 1])
end

@testset "Base.jl" begin
    a = 0.999
    @test Float16(Nemo.fmpq(a)) == Float16(a)
    @test Nemo.fmpq(a) == 999 // 1000
    a = 0.9999999
    @test Float32(Nemo.fmpq(a)) == Float32(a)
    a = 0.99999999
    @test Float64(Nemo.fmpq(a)) == Float64(a)
    a = 0.11111111111111111111111
    @test_broken BigFloat(Nemo.fmpq(a)) == BigFloat(a)
    @test push!([0.1], Nemo.fmpq(1 // 2)) == [0.1; 0.5]
end


text:
using VectorEngine, Test

function vadd(a, b, c)
    @vectorize for i = 1:(prod(size(a)))
        @inbounds c[i] = a[i] + b[i]
    end
    return
end

dims = (512,)
a = round.(rand(Float64, dims) * 100)
b = round.(rand(Float64, dims) * 100)
c = similar(a)

d_a = VEArray(a)
d_b = VEArray(b)
d_c = VEArray(c)

@synchronize @veda vadd(d_a, d_b, d_c)

c = Array(d_c)
@test a+b ≈ c


text:
using JSON, WordTokenizers, ProgressMeter

function processquestions(dhome, set; word_dic=Dict(), answer_dic=Dict())
    println("original data is loading")
    data = JSON.parsefile(joinpath(dhome,"questions","CLEVR_$(set)_questions.json"))
    println("parsing starts...")
    result = []
    p = Progress(length(data["questions"]))
    for question in data["questions"]
        words   = tokenize(question["question"])
        qtokens = map(x->get!(word_dic,x,length(word_dic)+1),words)
        atoken  = get!(answer_dic,question["answer"],length(answer_dic)+1)
        push!(result,(question["image_filename"],qtokens,atoken,question["question_family_index"]))
	    next!(p)
    end
    open("data/$(set).json", "w") do f
        write(f,json(result))
    end

     return word_dic, answer_dic
end

function preprocess(args)
    root       = args[1]
    adic       = Dict()
    wdic       = Dict()
    if length(args)==2
        dics = JSON.parsefile(args[2])
        wdic = dics["word_dic"]
        adic = dics["answer_dic"]
    end
    println("Parsing training questions...")
    processquestions(root, "train";word_dic=wdic, answer_dic=adic)
    println("Parsing validation questions...")
    processquestions(root, "val"  ;word_dic=wdic, answer_dic=adic)

    open("data/dic.json", "w") do f
        write(f,json(Dict("word_dic"=>wdic,"answer_dic"=>adic)))
    end
end

max_stars_repo_name: feyzaakyurek/MAC-Network

max_issues_repo_name: feyzaakyurek/MAC-Network

max_forks_repo_name: feyzaakyurek/MAC-Network


text:
N = 1000-1
sum = 0
for i = 1:N
	if i % 3 == 0 || i % 5 == 0
#		println(i)
		sum += i
	end
end
println(sum)

max_stars_repo_name: cmey/euler

max_issues_repo_name: cmey/euler

max_forks_repo_name: cmey/euler


text:
module BadelineMovingBlock

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Badeline Boss Moving Block" => Ahorn.EntityPlacement(
        Maple.BadelineMovingBlock,
        "rectangle",
        Dict{String, Any}(),
        function(entity)
            entity.data["nodes"] = [(Int(entity.data["x"]) + Int(entity.data["width"]) + 8, Int(entity.data["y"]))]
        end
    )
)

Ahorn.nodeLimits(entity::Maple.BadelineMovingBlock) = 1, 1
Ahorn.minimumSize(entity::Maple.BadelineMovingBlock) = 8, 8
Ahorn.resizable(entity::Maple.BadelineMovingBlock) = true, true

function Ahorn.selection(entity::Maple.BadelineMovingBlock)
    if entity.name == "finalBossMovingBlock"
        x, y = Ahorn.position(entity)
        nx, ny = Int.(entity.data["nodes"][1])

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

        return [Ahorn.Rectangle(x, y, width, height), Ahorn.Rectangle(nx, ny, width, height)]
    end
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.BadelineMovingBlock, room::Maple.Room)
    Ahorn.drawTileEntity(ctx, room, entity, material='g', blendIn=false)
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.BadelineMovingBlock, room::Maple.Room)
    x, y = Ahorn.position(entity)
    nodes = get(entity.data, "nodes", ())

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

    if !isempty(nodes)
        nx, ny = Int.(nodes[1])
        cox, coy = floor(Int, width / 2), floor(Int, height / 2)

        # Use 'G' instead of 'g', as that is the highlight color of the block (the active color)
        fakeTiles = Ahorn.createFakeTiles(room, nx, ny, width, height, 'G', blendIn=false)
        Ahorn.drawFakeTiles(ctx, room, fakeTiles, room.objTiles, true, nx, ny, clipEdges=true)
        Ahorn.drawArrow(ctx, x + cox, y + coy, nx + cox, ny + coy, Ahorn.colors.selection_selected_fc, headLength=6)
    end
end

end


text:
function hiv_sim!(x, ϵ1, ϵ2, params, dt)
    T1, T2, T1i, T2i, V, E = x
    toti = T1i + T2i
    λ1, λ2, λE, d1, d2, dE, k1, k2, Kb, Kd, δ, δE, ρ1, ρ2, m1, m2, bE, NT, c, f = params

    dT1 = λ1 - (d1 * T1) - (1-ϵ1) * k1 * V * T1
    dT2 = λ2 - (d2 * T2) - (1-f*ϵ1) * k2 * V * T2
    dT1i = (1-ϵ1) * k1 * V * T1 - (δ * T1i) - m1 * E * T1i
    dT2i = (1-f*ϵ1) * k2 * V * T2 - (δ * T2i) - m2 * E * T2i
    dV = (1 - ϵ2) * NT * δ * toti - (c * V) - ((1-ϵ1) * ρ1 * k1 * T1 + (1 - f*ϵ1) * ρ2 * k2 * T2) * V
    dE = λE + bE * toti * E / (toti + Kb) - dE * toti * E / (toti + Kd) - δE * E
    x[1] += dt * dT1
    x[2] += dt * dT2
    x[3] += dt * dT1i
    x[4] += dt * dT2i
    x[5] += dt * dV
    x[6] += dt * dE
    return nothing
end

function hiv_simulate!(x, ϵ, params, dt, sim_steps)
    Δt = dt / sim_steps
    for i in 1:sim_steps
        hiv_sim!(x, ϵ[1], ϵ[2], params, Δt)
    end
    return nothing
end

function hiv_params()
    λ1 = 10000.0
    λ2 = 31.98
    λE = 1.0
    d1 = 0.01
    d2 = 0.01
    dE = 0.25
    k1 = 8.0 * 1e-7
    k2 = 1.0 * 1e-4
    Kb = 100.0
    Kd = 500.0
    δ = 0.7
    δE = 0.1
    ρ1 = 1.0
    ρ2 = 1.0
    m1 = 1.0 * 1e-5
    m2 = 1.0 * 1e-5
    bE = 0.3
    NT = 100.0
    c = 13.0
    f = 0.34
    return λ1, λ2, λE, d1, d2, dE, k1, k2, Kb, Kd, δ, δE, ρ1, ρ2, m1, m2, bE, NT, c, f
end

function hiv_dt(x, p, t)
    T1, T2, T1i, T2i, V, E = x
    toti = T1i + T2i
    ϵ1, ϵ2 = p
    dT1 = 10000.0 - (0.01 * T1) - (1.0 - ϵ1) * (8.0 * 1e-7) * V * T1
    dT2 = 31.98 - (0.01 * T2) - (1.0 - 0.34*ϵ1) * (1.0 * 1e-4) * V * T2
    dT1i = (1-ϵ1) * (8.0 * 1e-7) * V * T1 - (0.7 * T1i) - (1.0 * 1e-5) * E * T1i
    dT2i = (1-0.34*ϵ1) * (1.0 * 1e-4) * V * T2 - (0.7 * T2i) - (1.0 * 1e-5) * E * T2i
    dV = (1 - ϵ2) * 100.0 * 0.7 * (T1i + T2i) - (13.0 * V) - ((1-ϵ1) * 1.0 * (8.0 * 1e-7) * T1 + (1 - 0.34*ϵ1) * 1.0 * (1.0 * 1e-4) * T2) * V
    dE = 1.0 + 0.3 * (T1i + T2i) * E / (T1i + T2i + 100.0) - 0.25 * (T1i + T2i) * E / (T1i + T2i + 500.0) - 0.1 * E
    return @SVector [dT1, dT2, dT1i, dT2i, dV, dE]
end

function hiv_simulate!(prob, x, ϵ, dt)
    tspan = (0.0, dt)  # simulate 5 days at a time
    prob = remake(prob, u0=x, p=ϵ)
    sol = DiffEqBase.solve(prob, BS3(), reltol=1e-1, abstol=1e-1, save_everystep=false)
    return sol.u[2]
end

function hiv_step!(prob::ODEProblem, state, action, dt, actions)
    ϵ = actions[action]
    t, s, x = state
    reward = hiv_reward(s, ϵ)
    s = hiv_simulate!(prob, s, ϵ, dt)
    t += dt
    @. x = log10(s)
    γ = 0.98
    if t ≥ 1000
        γ = 0.0
    end

    return (t, s, x), x, reward, γ

end

function hiv_step!(state, action, dt, actions, params)
    ϵ = actions[action]
    t, s, x = state
    reward = hiv_reward(s, ϵ)
    sim_steps = 5*65
    hiv_simulate!(s, ϵ, params, dt, sim_steps)
    # hiv_simulate!(s, ϵ, dt)
    t += dt
    @. x = log10(s)
    γ = 0.98
    if t ≥ 1000
        γ = 0.0
    end

    return (t, s, x), x, reward, γ

end

function hiv_reward(s, ϵ)
    V, E = s[5], s[6]
    ϵ1, ϵ2 = ϵ[1], ϵ[2]
    return 0.1 * V  + 20000.0 * ϵ1^2 + 2000.0 * ϵ2^2 - 1000.0 * E
end

function hiv_initial()
    t = 0.0
    # s = @SVector [163573.0, 5.0, 11945.0, 46.0, 63919.0, 24.0]
    s = [163573.0, 5.0, 11945.0, 46.0, 63919.0, 24.0]
    x = zeros(6)
    @. x = log10(s)
    return (t, s, x), x
end

function hiv()
    return create_hiv()
end

function create_hiv()
    [(3.0, 8.0), (-1, 4.0), (-2.5, 8.0), (-2.0, 4.0), (-1.0, 8.0), (0.0, 7.0)]
    X = Array{Float64,2}([
         3.0 8
        -1.0 4
        -2.5 8
        -2.0 4
        -1.0 8
         0.0 7
    ])
    S = ([0.0, 1000.0],
        exp10.(X),
        X
    )
    A = 1:4
    # u0 = @SVector [163573.0, 5.0, 11945.0, 46.0, 63919.0, 24.0]
    dt = 5.0
    # tspan = (0.0, dt)
    # ps = (0.0, 0.0)
    # prob = ODEProblem(hiv_dt,u0,tspan,ps)
    params = hiv_params()
    actions = ((0.0, 0.0), (0.7, 0.0), (0.0, 0.3), (0.7, 0.3))
    # p = (s,a)->hiv_step!(prob, s, a, dt, actions)
    p = (s,a)->hiv_step!(s, a, dt, actions, params)
    d0 = hiv_initial

    meta = Dict{Symbol,Any}()
    # meta[:minreward] = -5.0
    # meta[:maxreward] = 5
    # meta[:minreturn] = -5 * ceil(maxT / (dt * 20))
    # meta[:maxreturn] =  5 * ceil(maxT / (dt * 20))
    meta[:stochastic] = false
    meta[:minhorizon] = 200
    meta[:maxhorizon] = 200
    meta[:discounted] = true
    meta[:episodes] = 300

	bplot = HivPlotData()
    render = (state,clearplot=false)->hivplot(bplot,state,clearplot)
    m = SequentialProblem(S,X,A,p,d0,meta,render)
end

struct HivPlotData <: Any
    ts::Vector{Float64}
    xs::Vector{Vector{Float64}}

	function HivPlotData()
		new(Vector{Float64}(),Vector{Vector{Float64}}())
	end
end

@userplot HivPlot
@recipe function f(ap::HivPlot)
	data, state, clearplot = ap.args
	t, state, x = state
	if clearplot
        empty!(data.ts)
        empty!(data.xs)
	end

    push!(data.ts, t)
	push!(data.xs, deepcopy(x))
	legend := false
	xlims := (0, 1000)
	# grid := true
	# ticks := true
	layout := (3,2)
	# foreground_color := :white
	# aspect_ratio := 1.
    xs = hcat(data.xs...)
    names = ["T1 Healthy", "T2 Healthy", "T1 Infected", "T2 Infected", "Free Virus Particles", "Immune Effectors"]
    labels = ["log10 T1", "log10 T2", "log10 T1i", "log10 T2i", "log10 V", "log10 E"]
    ylimits = [(3.0, 8.0), (-1, 4.0), (-2.5, 8.0), (-2.0, 4.0), (-1.0, 8.0), (0.0, 7.0)]
    for i in 1:6
        @series begin
            title --> string(names[i])
            ylims := ylimits[i]
            yguide := labels[i]
            seriestype := :path
            subplot := i
            data.ts, xs[i, :]
        end
    end

end


text:
abstract UniverseModel
abstract MuSigmaModel <: UniverseModel
abstract DistributionModel <: UniverseModel

###################
## MuSigmaModels ##
###################

## SampleMoments
##--------------

type SampleMoments <: MuSigmaModel
    mu::Array{Float64, 1}
    sigma::Array{Float64, 2}
    names::Array{Symbol, 1}

    function SampleMoments(mu::Array{Float64, 1},
                           sigma::Array{Float64, 2},
                           names::Array{Symbol, 1})
        ## TODO: check for equal sizes

        return new(mu, sigma, names)
    end
end

## define empty constructor
##-------------------------

function SampleMoments()
    return SampleMoments(Array(Float64, 0), Array(Float64, 0, 0),
                         Array(Symbol, 0))
end

## Moving window sample moments
##-----------------------------

type MovWinSampleMoments <: MuSigmaModel
    mu::Array{Float64, 1}
    sigma::Array{Float64, 2}
    names::Array{Symbol, 1}

    function MovWinSampleMoments(mu::Array{Float64, 1},
                                 sigma::Array{Float64, 2},
                                 names::Array{Symbol, 1})
        ## TODO: check for equal sizes

        return new(mu, sigma, names)
    end
end

## define empty constructor
##-------------------------

function MovWinSampleMoments()
    return MovWinSampleMoments(Array(Float64, 0),
                               Array(Float64, 0, 0),
                               Array(Symbol, 0))
end

## ExpWeighted
##------------

type ExpWeighted <: MuSigmaModel
    mu::Array{Float64, 1}
    sigma::Array{Float64, 2}
    names::Array{Symbol, 1}
    obsWgts::Array{Float64, 1}

    function ExpWeighted(mu::Array{Float64, 1},
                  sigma::Array{Float64, 2},
                  names::Array{Symbol, 1},
                  obsWgts::Array{Float64, 1})
        ## TODO: check for equal sizes

        return new(mu, sigma, names, obsWgts)
    end
end

## define empty constructor
##-------------------------

function ExpWeighted()
    return ExpWeighted(Array(Float64, 0), Array(Float64, 0, 0),
                Array(Symbol, 0), Array(Float64, 0))
end

####################################
## general MuSigmaModel functions ##
####################################

## check for undefined model
##--------------------------

function isDef(mod::MuSigmaModel)
    if isempty(mod.mu)
        return false
    else
        return true
    end
end

import Base.Multimedia.display
function display(mod::MuSigmaModel)
    ## display type information
    typ = typeof(mod)
    println("\ntype: $typ")

    if !isDef(mod)
        println("Model is still undefined!")
        println("   mu: NA")
        println("sigma: NA")
    else

        ## display mus
        print("mu: ")
        df = DataFrame(name = mod.names, mus = mod.mu)
        display(df)
        print("\n")

        ## display Sigma
        print("sigma: ")
        display(mod.sigma)
        print("\n")
    end
end

##################
## extract info ##
##################

function getVolas(mod::MuSigmaModel)
    nAss = size(mod.sigma, 1)
    variances = Float64[mod.sigma[ii, ii] for ii=1:nAss]
    return sqrt(variances)
end

################
## estimators ##
################

## no estimator may use data of the given date itself!!
## One could think of the estimator as a forecast for the given date
## made in the morning, hence using data only up to there.

## sample moments
##---------------

function fitModel(::Type{SampleMoments}, data::Timematr, dat::Date;
                  minObs::Int = 50)
    ## extract data up to given date
    validDates = idx(data) .< dat

    if sum(validDates) < minObs
        return SampleMoments()
    else
        ## get data
        rets = asArr(data, Float64)[validDates, :]
        mus = mean(rets, 1)[:]
        sigmas = cov(rets)
        return SampleMoments(mus, sigmas, names(data))
    end
end

function fitModel(::Type{MovWinSampleMoments},
                  data::Timematr,
                  dat::Date; windSize::Int = 50)
    ## extract data up to given date
    validDates = idx(data) .< dat

    if sum(validDates) < windSize
        return MovWinSampleMoments()
    else
        ## get data
        rets = asArr(data, Float64)[validDates, :]

        ## only take data equal with length equal to window size
        rets = rets[(end-windSize+1):end, :]

        mus = mean(rets, 1)[:]
        sigmas = cov(rets)
        return MovWinSampleMoments(mus, sigmas, names(data))
    end
end

function fitModel(::Type{ExpWeighted}, data::Timematr, dat::Date;
                  minObs::Int = 50, lambda = 0.96)
    if (lambda < 0) | (1 < lambda)
        error("lambda must be in (0, 1)")
    end

    ## extract data up to given date
    validDates = idx(data) .< dat

    if sum(validDates) < minObs
        return ExpWeighted()
    else
        ## get data
        rets = asArr(data, Float64)[validDates, :]

        nObs, nAss = size(rets)

        ## get weights for observations
        obsWgts = lambda.^((nObs-1):(-1):0)
        wgtsSum = sum(obsWgts)

        ## calculate mus
        mus = sum(rets.*repmat(obsWgts, 1, nAss), 1)/wgtsSum

        ## calculate covariance
        sigmas = rets'*(rets.*repmat(obsWgts, 1, nAss))/wgtsSum

        return ExpWeighted(mus[:], sigmas, names(data), obsWgts)
    end
end

###########################
## get moments over time ##
###########################

macro defineApplyMuSigmaModelEstimators(myType)
    esc(quote
        function applyMuSigmaModelEstimator(modType::Type{$myType},
                                            data::Timematr)
            ## get all days
            allDats = idx(data)

            ## preallocation
            nObs, nAss = size(data)
            musOverTime = Array(Float64, nObs, nAss)
            sigmasOverTime = Array(Float64, nObs, nAss)
            nCovs = int((nAss)*(nAss-1)/2)
            corrOverTime = Array(Float64, nObs, nCovs)
            modelEstimated = falses(nObs)

            for ii=1:length(allDats)
                thisDat = allDats[ii]

                ## estimate moments
                mod = AssetMgmt.fitModel(modType, data, thisDat)

                ## extract mus and sigmas
                if AssetMgmt.isDef(mod)
                    ## set estimation indicator to true
                    modelEstimated[ii] = true

                    musOverTime[ii, :] = mod.mu'
                    sigmasOverTime[ii, :] =
                        (Float64[sqrt(mod.sigma[jj, jj]) for jj=1:nAss])'

                    ## get correlation matrix
                    d = diagm(1./sqrt(diag(mod.sigma)))
                    corrMatr = d*mod.sigma*d

                    ## extract correlations
                    corrs = vcat([corrMatr[(jj+1:end), jj] for jj=1:(nAss-1)]...)
                    corrOverTime[ii, :] = corrs'
                end
            end

            ## remove non-estimatable dates
            musOverTime = musOverTime[modelEstimated, :]
            sigmasOverTime = sigmasOverTime[modelEstimated, :]
            corrOverTime = corrOverTime[modelEstimated, :]

            ## transform to Timematr
            dfMus = DataFrame()
            dfSigmas = DataFrame()
            for ii=1:nAss
                thisNam = names(data)[ii]
                dfMus[thisNam] = musOverTime[:, ii]
                dfSigmas[thisNam] = sigmasOverTime[:, ii]
            end

            dfCorrs = DataFrame()
            for ii=1:size(corrOverTime, 2)
                dfCorrs[ii] = corrOverTime[:, ii]
            end

            dats = idx(data)[modelEstimated]
            musOverTimeTd = Timematr(dfMus, dats)
            sigmasOverTimeTd = Timematr(dfSigmas, dats)
            corrOverTimeTd = Timematr(dfCorrs, dats)
            return (musOverTimeTd, sigmasOverTimeTd, corrOverTimeTd)
        end

    end)
end

for t = (:(SampleMoments), :(MovWinSampleMoments), :(ExpWeighted))
    eval(macroexpand(:(@defineApplyMuSigmaModelEstimators($t))))
end


text:
# Copyright (c) 2013 Jasper den Ouden, under the MIT license,
# see doc/mit.txt from the project directory.

# Parses just a subset, try to get as far as possible with this first.
module JuliaSubsetParse

using Treekenize

# --- end module stuff

export julia_subset_parse, sexpr_print

function julia_subset_parse(input)
    infix_set = ["=", "&&","||", "!","~","==","!=", ">",">=","<=","<",
                 "-","+", "/","*", "::", ":"]

    add_infix(to) = (to[1],to[2],infix_set)
    return treekenize(input,
                      map(add_infix,
                          {("(",")"), ("[","]"), ("{","}"),
                           ("begin","end"),("function","end"),("if","end"),
                           ("type","end"),
                           ("typealias","\n"), ("@","\n"), ("#","\n")}),
                      add_infix(("\n","\n")), 10,1)
end

function sexpr_print(expr::StrExpr)
    if expr.head!="("
        print("($(expr.head)")
    else
        print("(open")
    end
    sexpr_print(expr.body, ",")
    print(")")
end
function sexpr_print(array::Array, between)
    for el in array
        print(between)
        sexpr_print(el)
    end
end
sexpr_print(array::Array) = sexpr_print(array, ",\n")
sexpr_print(other) = print(other)

end # module JuliaSubsetParse


text:
@i function red_reg(reg::ArrayReg{B,T}, Ly::Int, Js, hs, REG_STACK) where {B,T}
    k ← 0
    for j=1:Ly*4
        apply_Gh!(reg, j, hs[j], REG_STACK)
    end
    @invcheckoff for i=1:Ly-1
        for j = 1:4
            apply_Gvb!(reg, (@const (4*(i-1)+j,4*i+j)), Js[k+1], REG_STACK)
            k += 1
        end
    end
    @invcheckoff for i=1:Ly
        apply_G16!(reg, (@const (4i-3,4i-2,4i-1,4i)), Js[k+1:k+16], REG_STACK)
        k += 16
    end
    k → length(Js)
end

@i function isolve_largemem(out!::T, sg::Spinglass{ChimeraLattice, T}, reg::ArrayReg{B,TT}, REG_STACK) where {B,T,TT<:Tropical{T}}
    @invcheckoff begin
    @routine begin
        Lx ← sg.lattice.Nx
        Ly ← sg.lattice.Ny
        Js ← sg.Js
        hs ← sg.hs
        nj_red ← (Ly-1)*4 + 16 * Ly
    end
    k ← 0
    @safe println("Running Chimera, Lx = $Lx, Ly = $Ly, eltype = $T.")
    @safe println("Layer 1/$Lx")
    red_reg(reg, Ly, Js[k+1:k+nj_red], hs[1:4Ly], REG_STACK)
    k += nj_red
    for j=1:Ly*4 # `hs` interated in red->black order
        apply_Gh!(reg, j, hs[j+4Ly], REG_STACK)
    end

    for i=2:Lx
        @routine hk ← (i-1)*Ly*8
        @safe println("Layer $i/$Lx")
        # BLACK
        for j=1:Ly*4
            apply_Ghb!(reg, j, Js[k+1], REG_STACK)
            k += 1
        end
        for j=1:Ly*4 # `hs` interated in red->black order
            apply_Gh!(reg, j, hs[j+hk+4Ly], REG_STACK)
        end

        # Contract with RED
        @routine rr ← _init_reg(Tropical{T}, Ly*4, Val(false))
        red_reg(rr, Ly, Js[k+1:k+nj_red], hs[hk+1:hk+4Ly], REG_STACK)
        reg.state .*= rr.state
        incstack!(REG_STACK)
        swap_state!(REG_STACK, rr.state)
        k += nj_red
        ~@routine
        ~@routine
    end
    summed ← one(TT)
    i_sum(summed, reg.state)
    NiLang.SWAP(summed.n, out!)
    k → length(Js)
    summed → one(TT)
    ~@routine
    end
end

@i function isolve(out!::T, sg::AbstractSpinglass{ChimeraLattice, T}, reg::ArrayReg{B,TT}, A_STACK, B_STACK) where {B,T,TT<:Tropical{T}}
    @invcheckoff begin
    @routine begin
        Lx ← sg.lattice.Nx
        Ly ← sg.lattice.Ny
        Js ← sg.Js
        hs ← sg.hs
        nj_red ← (Ly-1)*4 + 16 * Ly
    end
    k ← 0
    @safe println("Running Chimera, Lx = $Lx, Ly = $Ly, eltype = $T.")
    @safe println("Layer 1/$Lx")

    # initial reg
    @routine begin
        rr ← _init_reg(Tropical{T}, Ly*4, Val(false))
        red_reg(rr, Ly, Js[k+1:k+nj_red], hs[1:4Ly], A_STACK)
    end
    reg.state .*= rr.state
    ~@routine
    k += nj_red

    for j=1:Ly*4 # `hs` interated in red->black order
        apply_Gh!(reg, j, hs[j+4Ly], A_STACK)
    end

    for i=2:Lx
        @safe println("Layer $i/$Lx")
        hk ← (i-1)*Ly*8
        # BLACK
        @routine begin
            for j=1:Ly*4
                apply_Ghb!(reg, j, Js[k+1], A_STACK)
                k += 1
            end
        end
        incstack!(B_STACK)
        store_state!(B_STACK, reg.state)
        ~@routine
        k += 4*Ly
        swap_state!(B_STACK, reg.state)

        for j=1:Ly*4 # `hs` interated in red->black order
            apply_Gh!(reg, j, hs[j+hk+4Ly], A_STACK)
        end

        # RED
        @routine begin
            rr ← _init_reg(Tropical{T}, Ly*4, Val(false))
            red_reg(rr, Ly, Js[k+1:k+nj_red], hs[hk+1:hk+4Ly], A_STACK)
        end
        reg.state .*= rr.state
        ~@routine
        k += nj_red
        hk → (i-1)*Ly*8
    end
    summed ← one(TT)
    i_sum(summed, reg.state)
    NiLang.SWAP(summed.n, out!)
    k → length(Js)
    summed → one(TT)
    ~@routine
    end
end

cachesize_A(lt::ChimeraLattice) = lt.Ny*4
cachesize_B(lt::ChimeraLattice) = lt.Nx-1
cachesize_largemem(lt::ChimeraLattice) = (lt.Nx-1)*lt.Ny*4


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

module APITests
import ..Pkg # ensure we are using the correct Pkg

using Pkg, Test
import Pkg.Types.PkgError, Pkg.Resolve.ResolverError
using UUIDs

using ..Utils

@testset "Pkg.activate" begin
    temp_pkg_dir() do project_path
        cd_tempdir() do tmp
            path = pwd()
            Pkg.activate(".")
            mkdir("Foo")
            cd(mkdir("modules")) do
                Pkg.generate("Foo")
            end
            Pkg.develop(Pkg.PackageSpec(path="modules/Foo")) # to avoid issue #542
            Pkg.activate("Foo") # activate path Foo over deps Foo
            @test Base.active_project() == joinpath(path, "Foo", "Project.toml")
            Pkg.activate(".")
            rm("Foo"; force=true, recursive=true)
            Pkg.activate("Foo") # activate path from developed Foo
            @test Base.active_project() == joinpath(path, "modules", "Foo", "Project.toml")
            Pkg.activate(".")
            Pkg.activate("./Foo") # activate empty directory Foo (sidestep the developed Foo)
            @test Base.active_project() == joinpath(path, "Foo", "Project.toml")
            Pkg.activate(".")
            Pkg.activate("Bar") # activate empty directory Bar
            @test Base.active_project() == joinpath(path, "Bar", "Project.toml")
            Pkg.activate(".")
            Pkg.add("Example") # non-deved deps should not be activated
            Pkg.activate("Example")
            @test Base.active_project() == joinpath(path, "Example", "Project.toml")
            Pkg.activate(".")
            cd(mkdir("tests"))
            Pkg.activate("Foo") # activate developed Foo from another directory
            @test Base.active_project() == joinpath(path, "modules", "Foo", "Project.toml")
            Pkg.activate() # activate home project
            @test Base.ACTIVE_PROJECT[] === nothing
        end
    end
end

@testset "Pkg.status" begin
    temp_pkg_dir() do project_path
        Pkg.add(PackageSpec(name="Example", version="0.5.1"))
        Pkg.add("Random")
        Pkg.status()
        Pkg.status("Example")
        Pkg.status(["Example", "Random"])
        Pkg.status(PackageSpec("Example"))
        Pkg.status(PackageSpec(uuid = "7876af07-990d-54b4-ab0e-23690620f79a"))
        Pkg.status(PackageSpec.(["Example", "Random"]))
        Pkg.status(; mode=PKGMODE_MANIFEST)
        Pkg.status("Example"; mode=PKGMODE_MANIFEST)
        @test_deprecated Pkg.status(PKGMODE_MANIFEST)
        # issue #1183: Test exist in manifest but not in project
        Pkg.status("Test"; mode=PKGMODE_MANIFEST)
        @test_throws PkgError Pkg.status("Test"; mode=Pkg.Types.PKGMODE_COMBINED)
        @test_throws PkgError Pkg.status("Test"; mode=PKGMODE_PROJECT)
        # diff option
        @test_logs (:warn, r"diff option only available") Pkg.status(diff=true)
        git_init_and_commit(project_path)
        @test_logs () Pkg.status(diff=true)
    end
end

@testset "Pkg.test" begin
    temp_pkg_dir() do tmp
        copy_test_package(tmp, "TestArguments")
        Pkg.activate(joinpath(tmp, "TestArguments"))
        # test the old code path (no test/Project.toml)
        Pkg.test("TestArguments"; test_args=`a b`, julia_args=`--quiet --check-bounds=no`)
        Pkg.test("TestArguments"; test_args=["a", "b"], julia_args=["--quiet", "--check-bounds=no"])
        # test new code path
        touch(joinpath(tmp, "TestArguments", "test", "Project.toml"))
        Pkg.test("TestArguments"; test_args=`a b`, julia_args=`--quiet --check-bounds=no`)
        Pkg.test("TestArguments"; test_args=["a", "b"], julia_args=["--quiet", "--check-bounds=no"])
    end
end

end # module APITests


text:
using Distributed
addprocs(2)
@everywhere using DistributedArrays, LinearAlgebra, Random, Test, WaveFD

macro ignore(ex::Expr)
end

@testset "SpaceTime" begin
    @testset "dtmod 2D, T=$(T)" for T in (Float32,Float64)
        for (S,Q,D) in ((:Fornberg,4,1),(:Fornberg,8,1),(:Fornberg,8,2),(:Nihei,8,1))
            dtmod = WaveFD.default_dtmod(S, Q, D, T(10.0), T(10.0), T(1500.0), T(-1.0), T(0.004), T(0.25))
            @test isapprox(0.004/dtmod, round(0.004/dtmod), rtol=1e-6)
        end
    end

    @testset "dtmod 3D, T=$(T)" for T in (Float32,Float64)
        for (S,Q,D) in ((:Fornberg,8,1),(:Fornberg,8,2))
            dtmod = WaveFD.default_dtmod(S, Q, D, T(10.0), T(10.0), T(10.0), T(1500.0), T(-1.0), T(0.004), T(0.25))
            @test isapprox(0.004/dtmod, round(0.004/dtmod), rtol=1e-6)
        end
    end

    @testset "ntmod" begin
        ntrec,dtrec,dtmod=128,0.004,0.001
        ntmod = WaveFD.default_ntmod(dtrec, dtmod, 0.0, ntrec)[2]
        @test (ntrec-1)*dtrec ≈ (ntmod-1)*dtmod
        ntmod = WaveFD.default_ntmod(dtrec, dtmod, ntrec)
        @test (ntrec-1)*dtrec ≈ (ntmod-1)*dtmod

        t0 = Array{Array}(undef, 2)
        t0[1] = [-0.1]
        t0[2] = [0.0]
        it,ntmod = WaveFD.default_ntmod(dtrec, dtmod, t0, ntrec)
        @test -(it-1)*dtmod ≈ t0[1][1]
        t0 = distribute(t0)
        it,ntmod = WaveFD.default_ntmod(dtrec, dtmod, t0, ntrec)
        @test -(it-1)*dtmod ≈ t0[1][1]
    end

    @testset "Time interpolation tests, dot product, T=$(T), mode=$(mode), alg=$(alg), nthreads=$(nthreads), $(length(n)+1)D" for T in (Float32, Float64), mode in (0,1), alg in (WaveFD.LangJulia(),WaveFD.LangC()), n in ((),(4,),(4,5)), nthreads=(1,4)
        dtrec=T(.004)
        dtmod=T(.001)
        m = rand(T,256,n...)
        ms = rand(T,256,n...)
        d = rand(T,64,n...)
        ds = rand(T,64,n...)
        WaveFD.interpforward!(WaveFD.interpfilters(dtmod,dtrec,mode,alg,nthreads), ds, m)
        WaveFD.interpadjoint!(WaveFD.interpfilters(dtmod,dtrec,mode,alg,nthreads), ms, d)
        rhs = dot(ds,d)
        lhs = dot(m,ms)
        err = norm(rhs-lhs)
        write(stdout, "T=$(T), lhs=$(lhs), rhs=$(rhs), mode=$(mode), alg=$(alg), $(length(n)+1)D, err=$(err)\n")
        @test isapprox(err, 0.0, atol=eps(T)*length(m))
    end

    @testset "Time interpolation no-op, mode=$mode, T=$T, n=$n, alg=$alg, nthreads=$nthreads" for mode in (0,1), T in (Float32,Float64), n in ((),(4,),(4,5)), alg in (WaveFD.LangJulia(), WaveFD.LangC()), nthreads=(1,4)
        dtrec=T(.004)
        dtmod=T(.004)
        m = rand(T,64,n...)
        ms = rand(T,64,n...)
        d = rand(T,64,n...)
        ds = rand(T,64,n...)
        WaveFD.interpforward!(WaveFD.interpfilters(dtmod,dtrec,mode,alg,nthreads), ds, m)
        WaveFD.interpadjoint!(WaveFD.interpfilters(dtmod,dtrec,mode,alg,nthreads), ms, d)
        rhs = dot(ds,d)
        lhs = dot(m,ms)
        err = norm(rhs-lhs)
        write(stdout, "T=$(T), lhs=$(lhs), rhs=$(rhs), mode=$(mode), alg=$(alg), $(length(n)+1)D, err=$(err)\n")
        @test isapprox(err, 0.0, atol=eps(T)*length(m))
        @test ds ≈ m
        @test ms ≈ d
    end

    @testset "Time interpolation tests, adjoint, accuracy, T=$(T), alg=$(alg), nthreads=$(nthreads), $(length(n)+1)D" for T in (Float32, Float64), alg=(WaveFD.LangJulia(),WaveFD.LangC()), n in ((),(4,),(4,5)), nthreads=(1,4)
        dtrec=T(.004)
        dtmod=T(.001)

        d = Array{T}(undef,64,n...)
        m = Array{T}(undef,256,n...)
        m_check = Array{T}(undef,256,n...)

        td = .004*collect(0:63)
        tm = .001*collect(0:255)

        dd = sin.(10*td) + .2*sin.(30*td)
        mm = sin.(10*tm) + .2*sin.(30*tm)

        d_expected = Array{T}(undef,64,n...)
        m_expected = Array{T}(undef,256,n...)
        for i in CartesianIndices(n)
            d_expected[:,i] .= dd[:]
            m_expected[:,i] .= mm[:]
        end

        I = (8:256-8,)
        I = length(n) == 1 ? (I...,:) : I
        I = length(n) == 2 ? (I...,:,:) : I

        m = similar(m_expected); m[:] .= 0.0
        WaveFD.interpadjoint!(WaveFD.interpfilters(dtmod,dtrec,0,alg,nthreads),m,d_expected)
        m *= 4
        err = norm((m-m_expected)[I...])/length(m)
        write(stdout, "T=$(T), interpadjoint!, mode=0, nthreads=$(nthreads), $(length(n)+1)D, err=$(err)\n")
        @test isapprox(err, 0.0, atol=.1)

        m = similar(m_expected); m[:] .= 0.0
        WaveFD.interpadjoint!(WaveFD.interpfilters(dtmod,dtrec,1,alg,nthreads),m,d_expected)
        err = norm((m .- m_expected)[I...])/length(m)
        write(stdout, "T=$(T), interpadjoint!, mode=1, nthreads=$(nthreads), $(length(n)+1)D, err=$(err)\n")
        @test isapprox(err, 0.0, atol=.1)
    end

    @testset "Time interpolation tests, forward, accuracy, T=$(T), alg=$(alg), nthreads=$(nthreads), $(length(n)+1)D" for T in (Float32, Float64), alg=(WaveFD.LangJulia(),WaveFD.LangC()), n in ((),(4,),(4,5)), nthreads=(1,4)
        dtrec=T(.004)
        dtmod=T(.001)

        d = Array{T}(undef,64,n...)
        m = Array{T}(undef,256,n...)
        m_check = Array{T}(undef,256,n...)

        td = .004*collect(0:63)
        tm = .001*collect(0:255)

        dd = sin.(10 .* td) .+ .2*sin.(30 .* td)
        mm = sin.(10 .* tm) .+ .2*sin.(30 .* tm)

        d_expected = Array{T}(undef,64,n...)
        m_expected = Array{T}(undef,256,n...)
        for i in CartesianIndices(n)
            d_expected[:,i] .= dd[:]
            m_expected[:,i] .= mm[:]
        end

        I = (8:64-8,)
        I = length(n) == 1 ? (I...,:) : I
        I = length(n) == 2 ? (I...,:,:) : I

        d = similar(d_expected); d[:] .= 0.0
        WaveFD.interpforward!(WaveFD.interpfilters(dtmod,dtrec,0,alg,nthreads),d,m_expected)
        err = norm((d .- d_expected)[I...])/length(d)
        write(stdout, "T=$(T), interpforward!, mode=0, nthreads=$(nthreads), $(length(n)+1)D, err=$(err)\n")
        @test isapprox(err, 0.0, atol=.1)

        d = similar(d_expected); d[:] .= 0.0
        WaveFD.interpforward!(WaveFD.interpfilters(dtmod,dtrec,1,alg,nthreads),d,m_expected)
        d ./= 4
        err = norm((d .- d_expected)[I...])/length(d)
        write(stdout, "T=$(T), interpforward!, mode=1, nthreads=$(nthreads), $(length(n)+1)D, err=$(err)\n")
        @test isapprox(err, 0.0, atol=.1)
    end

    @testset "data injection/extrapolation, 2D off-grid, inner product, T=$(T), F=$(F)" for T in (Float32,Float64), F in (WaveFD.hickscoeffs, WaveFD.linearcoeffs), alg=(WaveFD.LangC(),WaveFD.LangJulia()), nthreads=(1,4)
        nr = 10
        nz, nx = 100, 100
        dz, dx = 10.0, 10.0
        z0, x0 = 0.0, 0.0
        z1 = 5.0*dz; z2 = dx * (nz - 6)
        x1 = 5.0*dx; x2 = dx * (nx - 6)
        z, x = z1 .+ (z2 - z1) .* rand(nr), x1 .+ (x2 - x1) .* rand(nr)
        iz, ix, c = F(T(dz), T(dx), z0, x0, nz, nx, z, x)

        f1 = rand(T, 1, nr)
        f2 = zeros(T, 1, nr)
        g1 = rand(T, nz, nx)
        g2 = zeros(T, nz, nz)

        WaveFD.injectdata!(g2, f1, 1, iz, ix, c, 10)
        WaveFD.extractdata!(f2, g1, 1, iz, ix, c)

        lhs = dot(f1,f2)
        rhs = dot(g1,g2)
        err = abs(lhs - rhs) / abs(lhs + rhs)
        write(stdout, "spacetime, inject/extractdata dot product test (T=$(T), F=$(F)), lhs=$(lhs), rhs=$(rhs), err=$(err)\n")
        @test lhs ≈ rhs
    end

    @testset "data injection/extrapolation, 3D off-grid, inner product, T=$(T), F=$(F)" for T in (Float32,Float64), F in (WaveFD.hickscoeffs, WaveFD.linearcoeffs)
        nr = 10
        nz, ny, nx = 100, 100, 100
        dz, dy, dx = 10.0, 10.0, 10.0
        z0, y0, x0 = 0.0, 0.0, 0.0
        z1 = 5.0*dz; z2 = dx * (nz - 6)
        y1 = 5.0*dy; y2 = dy * (ny - 6)
        x1 = 5.0*dx; x2 = dx * (nx - 6)
        z, y, x = z1 .+ (z2 - z1) .* rand(nr), y1 .+ (y2 - y1) .* rand(nr), x1 .+ (x2 - x1) .* rand(nr)
        iz, iy, ix, c = F(T(dz), T(dy), T(dx), z0, y0, x0, nz, ny, nx, z, y, x)

        f1 = rand(T, 1, nr)
        f2 = zeros(T, 1, nr)
        g1 = rand(T, nz, ny, nx)
        g2 = zeros(T, nz, ny, nz)

        WaveFD.injectdata!(g2, f1, 1, iz, iy, ix, c)
        WaveFD.extractdata!(f2, g1, 1, iz, iy, ix, c)

        lhs = dot(f1,f2)
        rhs = dot(g1,g2)
        err = abs(lhs - rhs) / abs(lhs + rhs)
        write(stdout, "spacetime, inject/extractdata dot product test (T=$(T), F=$(F)), lhs=$(lhs), rhs=$(rhs), err=$(err)\n")
        @test lhs ≈ rhs
    end

    @testset "data injection/extrapolation, 2D on-grid, dot product tests, T=$(T), F=$(F), lang=$(lang), nthreads=$(nthreads)" for T in (Float32,Float64), F in (WaveFD.hickscoeffs, WaveFD.linearcoeffs, WaveFD.ongridcoeffs), lang in (WaveFD.LangC(), WaveFD.LangJulia()), nthreads in (1, 4)
        nz, nx=50, 52
        dz, dx = 10.0, 10.0
        z0, x0 = 0.0, 0.0

        nr = length(5:nx-5)
        rz = 10.0*ones(nr)
        rx = 10.0*collect(5:nx-5)

        m = rand(T,1,nr)    # inject data (source)
        ds = zeros(T,nz,nx) # ... get the field

        d = rand(T,nz,nx)   # given the field
        ms = zeros(T,1,nr)  # ... get the data

        iz, ix, c = F(T(dz), T(dx), z0, x0, nz, nx, rz, rx)

        local optargs
        if F == WaveFD.ongridcoeffs
            optargs = (lang, nthreads)
        else
            optargs = ()
        end
        WaveFD.injectdata!(ds, m, 1, iz, ix, c, optargs...)  # m->ds
        WaveFD.extractdata!(ms, d, 1, iz, ix, c, optargs...) # d->ms
        lhs = dot(m,ms)
        rhs = dot(d,ds)
        err = abs(lhs - rhs) / abs(lhs + rhs)
        write(stdout, "spacetime, inject/extractdata dot product test (T=$(T), F=$(F), lang=$(lang), nthreads=$(nthreads)), lhs=$(lhs), rhs=$(rhs), err=$(err)\n")
        @test lhs ≈ rhs
    end

    @testset "data injection/extrapolation, 3D on-grid, dot product tests, T=$(T), F=$(F)" for T in (Float32,Float64), F in (WaveFD.hickscoeffs, WaveFD.linearcoeffs, WaveFD.ongridcoeffs)
        nz,ny,nx=50,51,52
        dz, dy, dx = 10.0, 10.0, 10.0
        z0, y0, x0 = 0.0, 0.0, 0.0

        rz = 10.0*ones(nx-10,ny-10)
        rx = zeros(nx-10,ny-10)
        ry = zeros(nx-10,ny-10)
        for iy = 1:(ny-10),ix = 1:(nx-10)
            rx[ix,iy] = (ix+5)*dx
            ry[ix,iy] = (iy+5)*dy
        end
        rz = rz[:]
        ry = ry[:]
        rx = rx[:]

        iz, iy, ix, c = F(T(dz), T(dy), T(dx), z0, y0, x0, nz, ny, nx, rz, ry, rx)

        m = rand(T,1,ny*nx)    # inject data (source)
        ds = zeros(T,nz,ny,nx) # ... get the field

        d = rand(T,nz,ny,nx)   # given the field
        ms = zeros(T,1,ny*nx)  # ... get the data

        WaveFD.injectdata!(ds, m, 1, iz, iy, ix, c)  # m->ds
        WaveFD.extractdata!(ms, d, 1, iz, iy, ix, c) # d->ms
        lhs = dot(m,ms)
        rhs = dot(d,ds)
        err = abs(lhs - rhs) / abs(lhs + rhs)
        write(stdout, "spacetime, inject/extractdata dot product test (T=$(T), F=$(F)), lhs=$(lhs), rhs=$(rhs), err=$(err)\n")
        @test lhs ≈ rhs

        if F == WaveFD.ongridcoeffs
            ds_C = copy(ds)
            ms_C = copy(ms)
            ds = zeros(T,nz,ny,nx) # ... get the field
            ms = zeros(T,1,ny*nx)  # ... get the data

            WaveFD.injectdata!(ds, m, 1, iz, iy, ix, c, WaveFD.LangJulia())  # m->ds
            WaveFD.extractdata!(ms, d, 1, iz, iy, ix, c, WaveFD.LangJulia()) # d->ms

            @test ds ≈ ds_C
            @test ms ≈ ms_C

            lhs = dot(m,ms)
            rhs = dot(d,ds)
            err = abs(lhs - rhs) / abs(lhs + rhs)
            write(stdout, "spacetime, inject/extractdata dot product test (T=$(T), F=$(F)), lhs=$(lhs), rhs=$(rhs), err=$(err), Julia code\n")
            @test lhs ≈ rhs
        end
    end

    @testset "data injection, modeling 2D off-grid accuracy tests, T=$(T),F=$(F),freesurface=$(fs)" for
            T in (Float32,), F in (WaveFD.hickscoeffs,WaveFD.linearcoeffs), fs in (true, false)

        function modeling(T,F,ongrid::Bool)
            nthreads = Sys.CPU_THREADS
            z,x,tmax,dt,nbz,nbx,nsponge = 2.0,2.0,0.5,0.001,10,10,10
            dz = ongrid ? 0.01 : 0.02
            dx = ongrid ? 0.01 : 0.02
            fpeak = ongrid ? 1.5/(5*2*dx)/3 : 1.5/(5*dx)/3
            sz,sx = z/2+0.01,x/2+0.01
            nz,nx,nt = round(Int,z/dz)+1, round(Int,x/dx)+1, round(Int,tmax/dt)+1

            prop = Prop2DAcoIsoDenQ_DEO2_FDTD(
                nz=nz,nx=nx,nsponge=nsponge,dz=dx,dx=dx,dt=dt,freesurface=fs)

            WaveFD.V(prop) .= 1.5
            WaveFD.B(prop) .= 1.0
            pcur = WaveFD.PCur(prop)
            pold = WaveFD.POld(prop)

            wavelet = convert(Array{T},reshape(get(WaveFD.WaveletCausalRicker(f=fpeak), dt*collect(0:nt-1)), nt, 1))

            iz, ix, c =  F(T(dz), T(dx), 0.0, 0.0, nz, nx, [sz], [sx])
            blocks = WaveFD.source_blocking(nz, nx, nbz, nbx, iz, ix, c)

            set_zero_subnormals(true)
            for it = 1:nt
                rem(it,10) == 0 && write(stdout, "..$(it/nt*100) percent..\r")
                WaveFD.propagateforward!(prop)
                pcur,pold = pold,pcur
                WaveFD.injectdata!(pcur, blocks, wavelet, it)
            end
            set_zero_subnormals(false)
            pcur./(dz*dx)
        end

        write(stdout, "Computing wavefields, be patient, 2D, T=$(T), F=$(F), freesurface=$(fs) ...\n")
        poff = modeling(T,F,false)
        pon = modeling(T,F,true)[1:2:end,1:2:end]
        r = maximum(vec((poff .- pon) ./ maximum(abs, pon)))
        @test r < .02
    end

    @testset "data injection, modeling 3D off-grid accuracy tests, T=$(T),F=$(F),freesurface=$(fs)" for
            T in (Float32,), F in (WaveFD.hickscoeffs,WaveFD.linearcoeffs), fs in (true, false)

        function modeling(T,F,ongrid::Bool)
            nthreads = Sys.CPU_THREADS
            z,y,x,tmax,dt,nbz,nby,nbx,nsponge = 2.0,2.0,2.0,0.5,0.001,10,10,10,10
            dz = ongrid ? 0.01 : 0.02
            dy = ongrid ? 0.01 : 0.02
            dx = ongrid ? 0.01 : 0.02
            fpeak = ongrid ? 1.5/(5*2*dx)/3 : 1.5/(5*dx)/3
            nz,ny,nx,nt = round(Int,z/dz)+1, round(Int,y/dy)+1, round(Int,x/dx)+1, round(Int,tmax/dt)+1
            sz,sy,sx = z/2+0.01,y/2+0.01,x/2+0.01
            nz,nx,nt = round(Int,z/dz)+1, round(Int,x/dx)+1, round(Int,tmax/dt)+1

            prop = WaveFD.Prop3DAcoIsoDenQ_DEO2_FDTD(
                nz=nz,ny=ny,nx=nx,nsponge=nsponge,dz=dz,dy=dy,dx=dx,dt=dt,freesurface=fs)

            WaveFD.V(prop) .= 1.5
            WaveFD.B(prop) .= 1.0
            pcur = WaveFD.PCur(prop)
            pold = WaveFD.POld(prop)

            wavelet = convert(Array{T},reshape(get(WaveFD.WaveletCausalRicker(f=fpeak), dt*collect(0:nt-1)), nt, 1))

            iz, iy, ix, c = F(T(dz), T(dy), T(dx), 0.0, 0.0, 0.0, nz, ny, nx, [sz], [sy], [sx])
            blocks = WaveFD.source_blocking(nz, ny, nx, nbz, nby, nbx, iz, iy, ix, c)

            set_zero_subnormals(true)
            for it = 1:nt
                rem(it,10) == 0 && write(stdout, "..$(it/nt*100) percent..\r")
                WaveFD.propagateforward!(prop)
                pcur,pold = pold,pcur
                WaveFD.injectdata!(pcur, blocks, wavelet, it)
            end
            set_zero_subnormals(false)
            pcur./(dz*dy*dx)
        end

        write(stdout, "Computing wavefields, be patient, 3D, T=$(T), F=$(F), freesurface=$(fs) ...\n")
        poff = modeling(T,F,false)
        pon = modeling(T,F,true)[1:2:end,1:2:end,1:2:end]
        r = maximum(vec((poff .- pon) ./ maximum(abs, pon)))
        @test r < .1
    end

    @testset "data injection/extrapolation, 2D hicks, on-grid tests" for T in (Float32,Float64)
        nz,nx,dz,dx,z0,x0 = 50,51,10.0,10.0,0.0,0.0
        rng = 6:(nx-6)
        rz = 7 .* dz .* ones(length(rng))
        rx = dx .* collect(rng)
        iz, ix, c = WaveFD.hickscoeffs(T(dz), T(dx), z0, x0, nz, nx, rz, rx)

        for i = 1:length(rz)
            @test size(c[i]) == (1,1)
            @test c[i] ≈ ones(T,1,1)
            @test length(ix[i]) == 1
            @test length(iz[i]) == 1
            @test ix[i][1,1] == round(Int,(rx[i]-x0)/dx) + 1
            @test iz[i][1,1] == round(Int,(rz[i]-z0)/dz) + 1
        end
    end

    @testset "data injection/extrapolation, 3D hicks, on-grid tests" for T in (Float32,Float64)
        nz,ny,nx,dz,dy,dx,z0,y0,x0 = 50,51,52,10.0,10.0,10.0,0.0,0.0,0.0
        rngx = 6:(nx-6)
        rngy = 6:(ny-6)
        nr = length(rngx)*length(rngy)
        rz = 7 .* dz .* ones(nr)
        ry = zeros(nr)
        rx = zeros(nr)
        ir = 1
        for x in rngx, y in rngy
            rx[ir], ry[ir] = dx*x, dy*y
            ir += 1
        end
        iz, iy, ix, c = WaveFD.hickscoeffs(T(dz), T(dy), T(dx), z0, y0, x0, nz, ny, nx, rz, ry, rx)

        c_expected = ones(T, 1, 1)
        for i = 1:length(rz)
            @test c[i] ≈ c_expected
        end
    end

    @testset "data injection/extrapolation, bilinear, on-grid tests" for T in (Float32,Float64)
        nz,nx,dz,dx,z0,x0 = 50,51,10.0,10.0,0.0,0.0
        rng = 2:(nx-1)
        rz = 2 .* dz .* ones(length(rng))
        rx = dx .* collect(rng .- 1)
        iz, iy, c = WaveFD.linearcoeffs(T(dz), T(dx), z0, x0, nz, nx, rz, rx)

        c_expected = zeros(T, 2, 2)
        c_expected[1,1] = 1.0
        for i = 1:length(rz)
            @test c[i] ≈ c_expected
        end
    end

    @testset "data injection/extrapolation, trilinear, on-grid tests" for T in (Float32,Float64)
        nz,ny,nx,dz,dy,dx,z0,y0,x0 = 50,51,52,10.0,10.0,10.0,0.0,0.0,0.0
        rngx = 2:(nx-1)
        rngy = 2:(ny-1)
        nr = length(rngx)*length(rngy)
        rz = 2 .* dz .* ones(nr)
        ry = zeros(nr)
        rx = zeros(nr)
        ir = 1
        for x in rngx, y in rngy
            rx[ir], ry[ir] = dx*(x-1), dy*(y-1)
            ir += 1
        end
        iz, iy, ix, c = WaveFD.linearcoeffs(T(dz), T(dy), T(dx), z0, y0, x0, nz, ny, nx, rz, ry, rx)

        c_expected = zeros(T, 2, 2, 2)
        c_expected[1,1,1] = 1.0
        for i = 1:length(rz)
            @test c[i] ≈ c_expected
        end
    end

    @testset "ongrid detection tests, 2D" for T in (Float32,Float64)
        z0 = 10.0
        x0 = 20.0
        dz = 5.0
        dx = 3.0
        sz = z0 .+ dz.*(randperm(128) .- 1)
        sx = x0 .+ dx.*(randperm(128) .- 1)
        @test WaveFD.allongrid(dz, dx, z0, x0, sz, sx) == true
        sz[10] += .1
        @test WaveFD.allongrid(dz, dx, z0, x0, sz, sx) == false
        sz[10] -= .1
        sx[10] += .1
        @test WaveFD.allongrid(dz, dx, z0, x0, sz, sx) == false
    end

    @testset "ongrid detection tests, 3D" for T in (Float32,Float64)
        z0 = 10.0
        y0 = 10.0
        x0 = 20.0
        dz = 5.0
        dy = 4.0
        dx = 3.0
        sz = z0 .+ dz.*(randperm(128) .- 1)
        sy = y0 .+ dy.*(randperm(128) .- 1)
        sx = x0 .+ dx.*(randperm(128) .- 1)
        @test WaveFD.allongrid(dz, dy, dx, z0, y0, x0, sz, sy, sx) == true
        sz[10] += .1
        @test WaveFD.allongrid(dz, dy, dx, z0, y0, x0, sz, sy, sx) == false
        sz[10] -= .1
        sx[10] += .1
        @test WaveFD.allongrid(dz, dy, dx, z0, y0, x0, sz, sy, sx) == false
        sx[10] -= .1
        sy[10] += .1
        @test WaveFD.allongrid(dz, dy, dx, z0, y0, x0, sz, sy, sx) == false
    end

    @testset "ongrid injection coefficients, 2D" for T in (Float32, Float64)
        z0 = 10.0
        x0 = 20.0
        dz = 5.0
        dx = 3.0
        nz = 128
        nx = 128
        z = z0 .+ dz.*(randperm(nx) .- 1)
        x = x0 .+ dx.*(randperm(nx) .- 1)

        iz, ix, c = WaveFD.ongridcoeffs(dz, dx, z0, x0, nz, nx, z, x)
        @test iz ≈ round.((z .- z0)./dz).+1
        @test ix ≈ round.((x .- x0)./dx).+1
        @test c ≈ ones(size(iz))
    end

    @testset "ongrid injection coefficients, 3D" for T in (Float32, Float64)
        z0 = 10.0
        y0 = 10.0
        x0 = 20.0
        dz = 5.0
        dy = 5.0
        dx = 3.0
        nz = 128
        ny = 128
        nx = 128
        z = z0 .+ dz.*(randperm(nx) .- 1)
        y = y0 .+ dy.*(randperm(nx) .- 1)
        x = x0 .+ dx.*(randperm(nx) .- 1)

        iz, iy, ix, c = WaveFD.ongridcoeffs(dz, dy, dx, z0, y0, x0, nz, ny, nx, z, y, x)
        @test iz ≈ round.((z .- z0)./dz).+1
        @test iy ≈ round.((y .- y0)./dy).+1
        @test ix ≈ round.((x .- x0)./dx).+1
        @test c ≈ ones(size(iz))
    end
end

nothing

max_stars_repo_name: ChevronETC/WaveFD

max_issues_repo_name: ChevronETC/WaveFD

max_forks_repo_name: ChevronETC/WaveFD


text:
"""

# read_optimize

Read optimize output file created by cmdstan.

### Method
```julia
read_optimize(m::Stanmodel)
```

### Required arguments
```julia
* `m::Stanmodel`    : Stanmodel object
```

"""
function read_optimize(model::Stanmodel)
  ## Collect the results in a Dict

  cnames = String[]
  res_type = "optimize"
  tdict = Dict()

  ## tdict contains the arrays of values ##
  tdict = Dict()

  for i in 1:model.nchains
    if isfile("$(model.name)_$(res_type)_$(i).csv")

      # A result type file for chain i is present ##
      instream = open("$(model.name)_$(res_type)_$(i).csv")
      if i == 1
        open("$(model.name)_$(res_type)_$(i).csv") do instream
          str = read(instream, String)
          sstr = split(str)
          tdict[:stan_major_version] = [parse(Int, sstr[4])]
          tdict[:stan_minor_version] = [parse(Int, sstr[8])]
          tdict[:stan_patch_version] = [parse(Int, sstr[12])]
        end
      end

      # After reopening the file, skip all comment lines
      open("$(model.name)_$(res_type)_$(i).csv") do instream
        skipchars(isspace, instream, linecomment='#')
        line = Unicode.normalize(readline(instream), newline2lf=true)

        # Extract samples variable names
        idx = split(strip(line), ",")
        index = [idx[k] for k in 1:length(idx)]
        cnames = convert.(String, idx)

        # Read optimized values
        for i in 1:model.method.iter
          line = Unicode.normalize(readline(instream), newline2lf=true)
          flds = Float64[]
          if eof(instream) && length(line) < 2
            return
          else
            flds = parse.(Float64, split(strip(line), ","))
            for k in 1:length(index)
              if index[k] in keys(tdict)
                # For all subsequent chains the entry should already be in tdict
                append!(tdict[index[k]], flds[k])
              else

                # First chain
                tdict[index[k]] = [flds[k]]
              end
            end
          end
        end
      end
    end
  end

  (tdict, cnames)

end


text:
using SocialNet
using Test

using SocialNet: peek

### First test the required behavior

include("code_test_api.jl")

## We have already tested the interface that translates between persons' names
## and `Int`s representing vertices. Here we test the shortest path algorithm
## without this complication. Furthermore, we test all source and target
## vertices for some basic graphs for which the results are obvious.
## Again, the dot notation broadcasts over all source vertices.
@testset "shortest paths on graph types" begin
    g = Graph(4) # empty graph

    # Map over all possible source vertices, 1:4. For the empty graph, all vertices
    # except the source vertex are unreachable, represented by -1.
    @test shortest_path_length.(g, 1:4) ==
        [[0, -1, -1, -1], [-1, 0, -1, -1], [-1, -1, 0, -1], [-1, -1, -1, 0]]

    g = cycle_graph(5) # This graph is a ring.
    @test shortest_path_length.(g, 1:5) ==
        [[0, 1, 2, 2, 1], [1, 0, 1, 2, 2], [2, 1, 0, 1, 2], [2, 2, 1, 0, 1], [1, 2, 2, 1, 0]]

    g = complete_graph(5)
    @test shortest_path_length.(g, 1:5) ==
        [[0, 1, 1, 1, 1], [1, 0, 1, 1, 1], [1, 1, 0, 1, 1], [1, 1, 1, 0, 1], [1, 1, 1, 1, 0]]

    # The first vertex is the central vertex. It shares an edge with all other vertices.
    # Each of the other vertices shares an edge only with the central vertex.
    g = star_graph(5)
    @test shortest_path_length.(g, 1:5) ==
        [[0, 1, 1, 1, 1], [1, 0, 2, 2, 2], [1, 2, 0, 2, 2], [1, 2, 2, 0, 2], [1, 2, 2, 2, 0]]

    # The star-cycle graph is a cycle graph plus a vertex in the middle with an edge to all vertices.
    # The last vertex is the central vertex.
    g = star_cycle_graph(9)
    @test shortest_path_length.(g, 1:9) ==
        [[0, 1, 2, 2, 2, 2, 2, 1, 1], [1, 0, 1, 2, 2, 2, 2, 2, 1], [2, 1, 0, 1, 2, 2, 2, 2, 1],
         [2, 2, 1, 0, 1, 2, 2, 2, 1], [2, 2, 2, 1, 0, 1, 2, 2, 1], [2, 2, 2, 2, 1, 0, 1, 2, 1],
         [2, 2, 2, 2, 2, 1, 0, 1, 1], [1, 2, 2, 2, 2, 2, 1, 0, 1], [1, 1, 1, 1, 1, 1, 1, 1, 0]] # <-- last vertex

end

### Unit tests / diagnostic tests follow

@testset "basic Graph" begin
    g = Graph()
    @test isempty(g)
    @test length(vertices(g)) == 0
    # Implementation detail, but important
    @test first(vertices(g)) == 1
    @test ne(g) == 0
    @test nv(g) == 0

    add_vertex!(g)
    @test ! isempty(g)
    @test vertices(g) == 1:1
    @test ne(g) == 0
    @test nv(g) == 1

    add_vertex!(g)
    add_edge!(g, 1, 2)
    @test ne(g) == 1
    @test nv(g) == 2
end

@testset "basic PriorityQueue" begin
    q = PriorityQueue()

    @test isempty(q) == true
    insert!(q, 1, 100)
    insert!(q, 10, 99)
    insert!(q, 5, 50)

    # The lowest priority, 1, is for item 100
    @test peek(q) == 100

    # Check if an item is in the queue
    @test (50 in q) == true

    # Remove the lowest priority item
    @test pop!(q) == 100

    # Check that the next lowest priorty item has advanced
    @test peek(q) == 50

end


text:
# TODO: Make the constructor
# IndexSet{N, IndexT, DataT}(::Index...)

# TODO: export this. Maybe make this more
# specific than an alias?
const filterinds = inds

#
# Order singleton type
#

@eval struct Order{x}
  (OrderT::Type{ <: Order})() = $(Expr(:new, :OrderT))
end

Order(x) = Order{x}()

"""
    IndexSet(::Function, ::Order{N})

Construct an IndexSet of length N from a function that accepts
an integer between 1:N and returns an Index.

# Examples
```julia
IndexSet(n -> Index(1, "i\$n"), Order(4))
```
"""
ITensors.IndexSet(f::Function, N::Int) =
  IndexSet(ntuple(f, N))

ITensors.IndexSet(f::Function, ::Order{N}) where {N} =
  IndexSet(ntuple(f, Val(N)))

# TODO: make IndexSet{N} type stable for element type
# and storage type.
Base.filter(O::Order{N},
            f::Function,
            is::IndexSet) where {N} = IndexSet{N}(filter(f,
                                                         Tuple(is)))

Base.filter(O::Order,
            is::IndexSet,
            args...; kwargs...) = filter(O,
                                         ITensors.fmatch(args...;
                                                         kwargs...),
                                         is)

# intersect
ITensors.commoninds(::Order{N},
           A...;
           kwargs...) where {N} =
  IndexSet{N}(intersect(ITensors.itensor2inds.(A)...;
                        kwargs...)...)

# symdiff
ITensors.noncommoninds(::Order{N},
              A...;
              kwargs...) where {N} =
  IndexSet{N}(symdiff(ITensors.itensor2inds.(A)...;
                      kwargs...)...)

# setdiff
ITensors.uniqueinds(::Order{N},
           A...;
           kwargs...) where {N} =
  IndexSet{N}(setdiff(ITensors.itensor2inds.(A)...;
                      kwargs...)...)

# union
ITensors.unioninds(::Order{N},
          A...;
          kwargs...) where {N} =
  IndexSet{N}(union(ITensors.itensor2inds.(A)...;
                    kwargs...)...)


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

function Base.show(io::IO, x::NDArray)
  print(io, "NDArray(")
  Base.show(io, try_get_shared(x, sync = :read))
  print(io, ")")
end

# for REPL
function Base.show(io::IO, ::MIME{Symbol("text/plain")}, x::NDArray{T,N}) where {T,N}
  type_ = split(string(typeof(x)), '.', limit=2)[end]
  n = length(x)
  size_ = N == 1 ? "$n-element" : join(size(x), "×")
  print(io, "$size_ $type_ @ $(context(x))", (n == 0) ? "" : ":\n")
  Base.print_array(io, try_get_shared(x, sync = :read))
end


text:
using BinaryBuilder, Pkg

name = "SuiteSparse"
version = v"5.8.1"

# Collection of sources required to build SuiteSparse
sources = [
    GitSource("https://github.com/DrTimothyAldenDavis/SuiteSparse.git",
              "1869379f464f0f8dac471edb4e6d010b2b0e639d"),
    DirectorySource("./bundled"),
]

# We enable experimental platforms as this is a core Julia dependency
platforms = supported_platforms(;experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libsuitesparseconfig",   :libsuitesparseconfig),
    LibraryProduct("libamd",                 :libamd),
    LibraryProduct("libbtf",                 :libbtf),
    LibraryProduct("libcamd",                :libcamd),
    LibraryProduct("libccolamd",             :libccolamd),
    LibraryProduct("libcolamd",              :libcolamd),
    LibraryProduct("libcholmod",             :libcholmod),
    LibraryProduct("libldl",                 :libldl),
    LibraryProduct("libklu",                 :libklu),
    LibraryProduct("libumfpack",             :libumfpack),
    LibraryProduct("librbio",                :librbio),
    LibraryProduct("libspqr",                :libspqr),
    LibraryProduct("libsuitesparse_wrapper", :libsuitesparse_wrapper),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("libblastrampoline_jll"),
#    Dependency("METIS_jll"),
]

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
# Function for calculating Fermi surface in 2D
function getFermiSurface2D(
            unitcell:: U,
            bondInteractionMatrix:: HB,
            N_points::Int64,
            kappa::Float64;
            fermi_energy::Float64=0.0,
            enforce_hermitian::Bool=false,
            epsilon::Float64=1e-10,
            epsilon_k::Float64=1e-10,
            bounds_lower::Array{Float64,1}=-2*pi.*ones(2),
            bounds_upper::Array{Float64,1}=2*pi.*ones(2),
            refold_to_first_BZ::Bool=true
        ) where {L,NB,S,B<:AbstractBond{L,NB},U<:AbstractUnitcell{S,B},HB<:AbstractBondHamiltonian{L,NS} where {NS}}

    # list of k values that contribute to the BZ
    k_values = zeros(Float64, N_points, 2)

    LTbondInteractionMatrix = newHamiltonian(unitcell,bondInteractionMatrix)

    # funtion of energy
    function energyfunction(kvector)
        eigenvalues = eigvals(getMatrixAtK(LTbondInteractionMatrix,kvector,kappa)) .- fermi_energy
        eigenvalues = eigenvalues .* eigenvalues
        return minimum(eigenvalues)
    end

    # the current search index
    index = 1
    # search until there are enough points
    while index <= N_points

        # find a suitable starting point for the Newton algorithm
        k = Float64[rand(), rand()]
        for j in 1:length(k)
            k[j] = k[j] * bounds_lower[j] + (1-k[j]) * bounds_upper[j]
        end

        # start with the initial energy
        e0 = energyfunction(k)
        # iterate i over 100 newton steps (maximum)
        for i in 1:100
            # check if the energy is already converged
            if e0 < epsilon
                # save the k vector
                k_values[index,:] = k
                # increment the index
                index = index+1
                # break the newton loop
                break
            end
            # the current energy
            H_0 = e0
            # the gradient of the energy
            H_eps =
            [
                energyfunction(k .+ [epsilon_k, 0]),
                energyfunction(k .+ [0, epsilon_k])
            ]
            dH = (H_eps - H_0) ./ epsilon_k
            # absolute value of the gradient
            dHdH = dot(dH, dH)
            # break the Newton loop if the gradient diverges or is flattened too much
            if abs(dHdH) < 1e-20 || abs(dHdH) > 1e20
                break
            end
            # increment k
            dk = dH * (H_0 / dHdH)
            k -= dk
            # calculate a new energy
            e0 = energyfunction(k)
        end
    end

    # if refolding is enabled, refold all points
    if refold_to_first_BZ
        # get the lattice vectors
        a1 = latticeVectors(unitcell)[1]
        a2 = latticeVectors(unitcell)[2]
        # get the reciprocal lattice vectors
        b1 = [a2[2], -a2[1]]
        b2 = [a1[2], -a1[1]]
        # normalize the vectors
        b1 .*= 2*pi/sum(b1.*a1)
        b2 .*= 2*pi/sum(b2.*a2)

        # iterate over all k points
        for index in 1:N_points
            # assumend to be refolded to begin with
            refolded = true
            # iterate while still refolding possible
            while refolded
                # new iteration, not refolded
                refolded = false
                # try all possible refoldings
                for r in [b1, b2, -b1, -b2, b1+b2, -b2-b1, b1-b2, -b2+b1]
                    if sum((k_values[index,:].-r).*(k_values[index,:].-r)) < sum((k_values[index,:]).*(k_values[index,:]))
                        k_values[index,:] .-= r
                        refolded=true
                    end
                end
            end
        end
    end

    # return the list
    return k_values

end

# Function for calculating Fermi surface in 3D
function getFermiSurface3D(
            unitcell:: U,
            bondInteractionMatrix:: HB,
            N_points::Int64,
            kappa::Float64;
            fermi_energy::Float64=0.0,
            enforce_hermitian::Bool=false,
            epsilon::Float64=1e-10,
            epsilon_k::Float64=1e-10,
            bounds_lower::Array{Float64,1}= -2 .*pi .*ones(3),
            bounds_upper::Array{Float64,1}=2 .* pi .*ones(3),
            refold_to_first_BZ::Bool=true
        ) where {L,NB,S,B<:AbstractBond{L,NB},U<:AbstractUnitcell{S,B},HB<:AbstractBondHamiltonian{L,NS} where {NS}}

    # list of k values that contribute to the BZ
    k_values = zeros(Float64, N_points, 3)

    LTbondInteractionMatrix = newHamiltonian(unitcell,bondInteractionMatrix)
    # funtion of energy
    function energyfunction(kvector)
        eigenvalues = eigvals(getMatrixAtK(LTbondInteractionMatrix,kvector,kappa)) .- fermi_energy
        eigenvalues = eigenvalues .* eigenvalues
        return minimum(eigenvalues)
        #return eigmin(getMatrixAtK(LTbondInteractionMatrix,kvector))
    end

    # the current search index
    index = 1
    # search until there are enough points
    while index <= N_points

        # find a suitable starting point for the Newton algorithm
        k = Float64[rand(), rand(), rand()]
        for j in 1:length(k)
            k[j] = k[j] * bounds_lower[j] + (1-k[j]) * bounds_upper[j]
        end

        # start with the initial energy
        e0 = energyfunction(k)
        # iterate i over 100 newton steps (maximum)
        for i in 1:1000
            # check if the energy is already converged
            if e0 < epsilon
                # save the k vector
                k_values[index,:] = k
                # increment the index
                index = index+1
                # break the newton loop
                break
            end
            # the current energy
            H_0 = e0
            # the gradient of the energy
            H_eps =
            [
                energyfunction(k .+ [epsilon_k, 0, 0]),
                energyfunction(k .+ [0, epsilon_k, 0]),
                energyfunction(k .+ [0, 0, epsilon_k])
            ]
            dH = (H_eps .- H_0) ./ epsilon_k
            # absolute value of the gradient
            dHdH = dot(dH, dH)
            # break the Newton loop if the gradient diverges or is flattened too much
            if abs(dHdH) < 1e-20 || abs(dHdH) > 1e20
                break
            end
            # increment k
            dk = dH * (H_0 / dHdH)
            k -= dk
            # calculate a new energy
            e0 = energyfunction(k)
        end
    end

    # if refolding is enabled, refold all points
    if refold_to_first_BZ
        # get the lattice vectors
        a1 = latticeVectors(unitcell)[1]
        a2 = latticeVectors(unitcell)[2]
        a3 = latticeVectors(unitcell)[3]
        # get the reciprocal lattice vectors
        b1 = cross(a2, a3)
        b2 = cross(a3, a1)
        b3 = cross(a1, a2)
        # normalize the vectors
        b1 .*= 2*pi/sum(b1.*a1)
        b2 .*= 2*pi/sum(b2.*a2)
        b3 .*= 2*pi/sum(b3.*a3)
        # put together refolding
        refoldings = Array{Float64,1}[]
        for i in [-1,0,1]
        for j in [-1,0,1]
        for l in [-1,0,1]
            refolding = i*b1 + j*b2 + l*b3
            if sum(abs.(refolding))>1e-8
                push!(refoldings, refolding)
            end
        end
        end
        end
        # iterate over all k points
        for index in 1:N_points
            # assumend to be refolded to begin with
            refolded = true
            # iterate while still refolding possible
            while refolded
                # new iteration, not refolded
                refolded = false
                # try all possible refoldings
                for r in refoldings

                    if sum((k_values[index,:].-r).*(k_values[index,:].-r)) < sum((k_values[index,:]).*(k_values[index,:]))
                        k_values[index,:] .-= r
                        refolded=true
                    end
                end
            end
        end
    end

    # return the list
    return k_values

end

# obtain Fermi surface

"""
    getFermiSurface(
                unitcell::Unitcell,
                N_points::Int64
             [; fermi_energy::Float64=0.0,
                enforce_hermitian::Bool=false,
                epsilon::Float64=1e-10,
                epsilon_k::Float64=1e-10,
                bounds_lower::Array{Float64,1}=-2*pi.*ones(4),
                bounds_upper::Array{Float64,1}=2*pi.*ones(4),
                refold_to_first_BZ::Bool=true ]
            )
calculates `N_points` points which belong to the Fermi surface of the tight-binding model given by
the `Unitcell` object `unitcell`. The Fermi energy can be adjusted by using `fermi_energy=...`.
The points are located using Newton's method with random starting positions. The procedure can be
modified by using some of the optional keywords:
- `enforce_hermitian` determines if the matrix used for energy calculation is made hermitian by construction
- `epsilon` energy threshold for values to be considered close to the Fermi energy
- `epsilon_k` small k increment that is used for calculation of derivatives
- `bounds_upper` / `bounds_lower` the bounds between which the random starting location for the Newton method is searched
After the search has been completed, the values will optionally be refolded into the first brillouin zone by choosing `refold_to_first_BZ=true`.
# Examples
```julia-repl
julia> fermi_surace = getFermiSurface(unitcell, 100)    # 2D unitcell
100x2 Array{Float64}
...
julia> fermi_surace = getFermiSurface(unitcell, 100)    # 3D unitcell
100x3 Array{Float64}
...
julia> fermi_surace = getFermiSurface(unitcell, 100, fermi_energy=1.0)
100x2 Array{Float64}
...
julia> fermi_surace = getFermiSurface(unitcell, 100, refold_to_first_BZ=false)
100x2 Array{Float64}
...
```
"""

function getFermiSurface(
            unitcell:: U,
            bondInteractionMatrix:: HB,
            N_points::Int64,
            kappa::Float64;
            fermi_energy::Float64=0.0,
            enforce_hermitian::Bool=false,
            epsilon::Float64=1e-10,
            epsilon_k::Float64=1e-10,
            bounds_lower::Array{Float64,1}=-2*pi.*ones(3),
            bounds_upper::Array{Float64,1}=2*pi.*ones(3),
            refold_to_first_BZ::Bool=true
        ) where {L,NB,S,B<:AbstractBond{L,NB},U<:AbstractUnitcell{S,B},HB<:AbstractBondHamiltonian{L,NS} where {NS}}

    # check if the unitcell can be put in either method
    if length(sites(unitcell)[1].point) != length(latticeVectors(unitcell))
       println("Unitcell has not the same number of lattice vectors as dimensions")
        return zeros(1,1)
    end

    # check which function to pass to
    if length(latticeVectors(unitcell)) == 2
        # 2D case
        return getFermiSurface2D(
                unitcell,
                bondInteractionMatrix,
                N_points,
                kappa,
                fermi_energy=fermi_energy,
                enforce_hermitian=enforce_hermitian,
                epsilon=epsilon,
                epsilon_k=epsilon_k,
                bounds_lower=bounds_lower,
                bounds_upper=bounds_upper,
                refold_to_first_BZ=refold_to_first_BZ
            )

   elseif length(latticeVectors(unitcell)) == 3
        # 3D case
        return getFermiSurface3D(
                unitcell,
                bondInteractionMatrix,
                N_points,
                kappa,
                fermi_energy=fermi_energy,
                enforce_hermitian=enforce_hermitian,
                epsilon=epsilon,
                epsilon_k=epsilon_k,
                bounds_lower=bounds_lower,
                bounds_upper=bounds_upper,
                refold_to_first_BZ=refold_to_first_BZ
            )
    end
end

export getFermiSurface


text:
import Base: show, IO

"""
    Shunt

Print `s` `n` times before printing `x`.
"""
struct Shunt
    n::Int
    s::Char
    function Shunt(n::Int, s::Char)
        @assert n >= 0
        return new(n, s)
    end
end

function show(io::IO, x::Shunt)
    if x.n > 1
        show(io, x.n == 1 ? x.s : x.s^x.n)
    elseif x.n == 1
        show(io, x.s)
    end
end

extend(shunt::Shunt, n::Int=4) = Shunt(shunt.n + n, shunt.s)

"""
    Shunted{Tx}

Apply `shunt` before printing `x`.
"""
struct Shunted{Tx}
    shunt::Shunt
    x::Tx
end

function show(io::IO, shunted::Shunted)
    show(io, shunted.shunt)
    show(io, shunted.x)
end

dummy_shunt(x) = Shunted(Shunt(0, ' '), x)

function print_shunted_list(io::IO, shunt::Shunt, X::AbstractVector)
    for (j, x) in enumerate(X)
        show(io, Shunted(shunt, x))
        j != length(X) && show('\n')
    end
end
@noinline function print_shunted_list(io::IO, shunt::Shunt, x::Tuple)
    return print_shunted_list(io, shunt, [x...])
end


text:
left(rect::Rect2D) = minimum(rect)[1]
right(rect::Rect2D) = maximum(rect)[1]

bottom(rect::Rect2D) = minimum(rect)[2]
top(rect::Rect2D) = maximum(rect)[2]

Base.getindex(bbox::Rect2D, ::Left) = left(bbox)
Base.getindex(bbox::Rect2D, ::Right) = right(bbox)
Base.getindex(bbox::Rect2D, ::Bottom) = bottom(bbox)
Base.getindex(bbox::Rect2D, ::Top) = top(bbox)

width(rect::Rect2D) = right(rect) - left(rect)
height(rect::Rect2D) = top(rect) - bottom(rect)

bottomleft(bbox::Rect2D{T}) where T = Point2{T}(left(bbox), bottom(bbox))
topleft(bbox::Rect2D{T}) where T = Point2{T}(left(bbox), top(bbox))
bottomright(bbox::Rect2D{T}) where T = Point2{T}(right(bbox), bottom(bbox))
topright(bbox::Rect2D{T}) where T = Point2{T}(right(bbox), top(bbox))

topline(bbox::BBox) = (topleft(bbox), topright(bbox))
bottomline(bbox::BBox) = (bottomleft(bbox), bottomright(bbox))
leftline(bbox::BBox) = (bottomleft(bbox), topleft(bbox))
rightline(bbox::BBox) = (bottomright(bbox), topright(bbox))

function BBox(left::Number, right::Number, bottom::Number, top::Number)
    mini = (left, bottom)
    maxi = (right, top)
    return BBox(mini, maxi .- mini)
end

function IRect2D(bbox::Rect2D)
    return IRect2D(
        round.(Int, minimum(bbox)),
        round.(Int, widths(bbox))
    )
end

function RowCols(ncols::Int, nrows::Int)
    return RowCols(
        zeros(ncols),
        zeros(ncols),
        zeros(nrows),
        zeros(nrows)
    )
end

Base.getindex(rowcols::RowCols, ::Left) = rowcols.lefts
Base.getindex(rowcols::RowCols, ::Right) = rowcols.rights
Base.getindex(rowcols::RowCols, ::Top) = rowcols.tops
Base.getindex(rowcols::RowCols, ::Bottom) = rowcols.bottoms

"""
    eachside(f)
Calls f over all sides (Left, Right, Top, Bottom), and creates a BBox from the result of f(side)
"""
function eachside(f)
    return BBox(map(f, (Left(), Right(), Bottom(), Top()))...)
end

"""
mapsides(
       f, first::Union{Rect2D, RowCols}, rest::Union{Rect2D, RowCols}...
   )::BBox
Maps f over all sides of the rectangle like arguments.
e.g.
```
mapsides(BBox(left, right, bottom, top)) do side::Side, side_val::Number
    return ...
end::BBox
```
"""
function mapsides(
        f, first::Union{Rect2D, RowCols}, rest::Union{Rect2D, RowCols}...
    )
    return eachside() do side
        f(side, getindex.((first, rest...), (side,))...)
    end
end


text:
using Oceananigans.Diagnostics: AbstractDiagnostic

import Oceananigans.Utils: time_to_run
import Oceananigans.Fields: AbstractField, compute!
import Oceananigans.OutputWriters: fetch_output
import Oceananigans.Diagnostics: run_diagnostic

"""
    WindowedTimeAverage{R, OP, FS} <: AbstractDiagnostic

An object for computing 'windowed' time averages, or moving time-averages
of a `operand` over a specified `time_window`, collected on `time_interval`.
"""
mutable struct WindowedTimeAverage{R, OP, FS} <: AbstractDiagnostic
                         result :: R
                        operand :: OP
                    time_window :: Float64
                  time_interval :: Float64
                         stride :: Int
                   field_slicer :: FS
              window_start_time :: Float64
         window_start_iteration :: Int
       previous_collection_time :: Float64
    previous_interval_stop_time :: Float64
                     collecting :: Bool
end

"""
    WindowedTimeAverage(operand; time_window,
                                 time_interval,
                                        stride = 1,
                                  field_slicer = FieldSlicer())

Returns an object for computing running averages of `operand` over `time_window`,
recurring on `time_interval`. During the collection period, averages are computed
every `stride` iteration.

`operand` may be an `Oceananigans.Field`, `Oceananigans.AbstractOperations.Computation,
or `Oceananigans.Diagnostics.Average`.

Calling `wta(model)` for `wta::WindowedTimeAverage` object returns `wta.result`.
"""
function WindowedTimeAverage(operand, model=nothing; time_window, time_interval, stride=1,
                                                     field_slicer = FieldSlicer())

    output = fetch_output(operand, model, field_slicer)
    result = similar(output) # convert views to arrays
    result .= output # initialize `result` with initial output

    return WindowedTimeAverage(result,
                               operand,
                               Float64(time_window),
                               Float64(time_interval),
                               stride,
                               field_slicer,
                               0.0, 0, 0.0, 0.0, false)
end

function time_to_run(clock, wta::WindowedTimeAverage)
    if (wta.collecting ||
        clock.time >= wta.previous_interval_stop_time + wta.time_interval - wta.time_window)

        return true
    else
        return false
    end
end

function accumulate_result!(wta, model)

    # Time increment:
    Δt = model.clock.time - wta.previous_collection_time

    # Time intervals:
    T_current = model.clock.time - wta.window_start_time
    T_previous = wta.previous_collection_time - wta.window_start_time

    # Accumulate left Riemann sum
    integrand = fetch_output(wta.operand, model, wta.field_slicer)
    @. wta.result = (wta.result * T_previous + integrand * Δt) / T_current

    # Save time of integrand collection
    wta.previous_collection_time = model.clock.time

    return nothing
end

function run_diagnostic(model, wta::WindowedTimeAverage)

    if model.clock.iteration == 0 # initialize previous interval stop time
        wta.previous_interval_stop_time = model.clock.time
    end

    # Don't start collecting if we are *only* "initializing" run_diagnostic at the beginning
    # of a Simulation.
    #
    # Note: this can be false at the zeroth iteration if time_interval == time_window (which
    # implies we are always collecting)

    initializing = model.clock.iteration == 0 &&
        model.clock.time < wta.previous_interval_stop_time + wta.time_interval - wta.time_window

    if !(wta.collecting) && !(initializing)
        # run_diagnostic has been called, but we are not currently collecting data.
        # Initialize data collection:

        # Start averaging period
        wta.collecting = true

        # Zero out result
        wta.result .= 0

        # Save averaging start time and the initial data collection time
        wta.window_start_time = model.clock.time
        wta.window_start_iteration = model.clock.iteration
        wta.previous_collection_time = model.clock.time

    elseif model.clock.time >= wta.previous_interval_stop_time + wta.time_interval
        # Output is imminent. Finalize averages and cease data collection.
        accumulate_result!(wta, model)

        # Averaging period is complete.
        wta.collecting = false

        # Reset the "previous" interval time,
        # subtracting a sliver that presents window overshoot from accumulating.
        wta.previous_interval_stop_time = model.clock.time - rem(model.clock.time, wta.time_interval)

    elseif mod(model.clock.iteration - wta.window_start_iteration, wta.stride) == 0
        # Collect data as usual
        accumulate_result!(wta, model)
    end

    return nothing
end

function (wta::WindowedTimeAverage)(model)

    # For the paranoid
    wta.collecting &&
        model.clock.iteration > 0 &&
        @warn "Returning a WindowedTimeAverage before the collection period is complete."

    return wta.result
end


text:
function get_data_paths(config)
    data_paths = []
    for i = 1:nrow(config)
        row = config[i, :]
        miss = replace(row.missingness, "," => ".")
        file = string(
            "n_factors_",
            row.n_factors,
            "_n_items_",
            row.n_items,
            "_missing_",
            miss,
            ".csv"
            )
        push!(data_paths, file)
    end
    return data_paths
end

function read_files(dir, data_paths)
    data = Vector{Any}()
    for i in 1:length(data_paths)
        push!(data, DataFrame(CSV.read(dir*"/"*data_paths[i], DataFrame; missingstring = "NA")))
    end
    return data
end

function gen_SEM_RAM(nfact, nitem)
    nfact = Int64(nfact)
    nitem = Int64(nitem)

    ## Model definition
    nobs = nfact*nitem
    nnod = nfact+nobs
    npar = 3nobs + nfact-1
    Symbolics.@variables x[1:npar]

    #F
    Ind = collect(1:nobs)
    Jnd = collect(1:nobs)
    V = fill(1,nobs)
    F = sparse(Ind, Jnd, V, nobs, nnod)

    #A
    Ind = collect(1:nobs)
    Jnd = vcat([fill(nobs+i, nitem) for i in 1:nfact]...)
    V = [x...][1:nobs]
    A = sparse(Ind, Jnd, V, nnod, nnod)
    xind = nobs+1
    for i in nobs+1:nnod-1
        A[i+1,i] = x[xind]
        xind = xind+1
    end

    #S
    Ind = collect(1:nnod)
    Jnd = collect(1:nnod)
    V = [[x...][nobs+nfact:2nobs+nfact-1]; fill(1.0, nfact)]
    S = sparse(Ind, Jnd, V, nnod, nnod)

	M = [[x...][2nobs+nfact:3nobs+nfact-1]...; fill(0.0, nfact)]

    return RAMMatrices(;A = A, S = S, F = F, M = M, parameters = x)
end

function gen_model(nfact, nitem, data, backend)

    ram_matrices = gen_SEM_RAM(nfact, nitem)

    if backend == "Optim.jl"
        semdiff = SemDiffOptim(
            LBFGS(
                ;linesearch = BackTracking(order=3),
                alphaguess = InitialHagerZhang()
                ),
            Optim.Options(
                ;f_tol = 1e-10,
                x_tol = 1.5e-8)
            )
    elseif backend =="NLopt.jl"
        semdiff = SemDiffNLopt(
            :LD_LBFGS,
            nothing
            )
    end

    model = Sem(
            ram_matrices = ram_matrices,
            data = data,
            imply = RAM,
            diff = semdiff,
            loss = (SemFIML,),
            start_val = start_simple,
            observed = SemObsMissing
        )

    return model
end

function gen_models(config, data_vec)
    models = []
    for i = 1:nrow(config)
        row = config[i, :]
        model = gen_model(row.n_factors, row.n_items, Matrix(data_vec[i]), row.backend)
        push!(models, model)
    end
    return models
end

function benchmark_models(models)
    benchmarks = []
    for model in models
        bm = @benchmark sem_fit($model)
        push!(benchmarks, bm)
    end
    return benchmarks
end

function get_fits(models)
    fits = []
    for model in models
        fit = sem_fit(model)
        push!(fits, fit)
    end
    return fits
end

function compare_estimates(fits, par_vec, config)
    correct =
        [compare_estimate(fit, estimate, n_factors, n_items) for
                (fit, estimate, n_factors, n_items) in zip(fits, par_vec, config.n_factors, config.n_items)]
    return correct
end

function compare_estimate(fit, estimate, n_factors, n_items)

    nfact = Int64(n_factors)
    nitem = Int64(n_items)

    nobs = nfact*nitem

    par_ind = [1:(2nobs+nfact-1)..., (2nobs+2nfact):(3nobs+2nfact-1)...]

    solution = fit.solution

    return StructuralEquationModels.compare_estimates(estimate.est[par_ind], solution, 0.01)

end