text:
function solve_simplex(input::Input)

    midterm = MidTerm()
    if check_phase_1(input)
        midterm.base =  collect((input.n - input.m + 1):input.n)
        midterm.nbase = collect(1:(input.n-input.m))
        val, i = findmin(input.b)

        aux = midterm.base[i]
        push!(midterm.nbase, aux)
        midterm.base[i] = input.n + 1
        c_mem = copy(input.c)
        input.A = [input.A -ones(input.m)]
        input.c = zeros(input.n + 1)
        input.c[end] = -1
        input.n = input.n + 1

        init_log_simplex1(input)
        while midterm.termination_status == 0 && midterm.iter < input.max_iter
            iterate!(input, midterm)
        end
        update_midterm!(input, midterm)
        if midterm.z < -input.tol
            midterm.termination_status = 3 #infeasible
            output = write_output(input, midterm)
            return output
        end
        midterm.iter = 0
        input.n = input.n - 1
        cache = findfirst(x->x==input.n + 1,midterm.nbase)
        input.c = c_mem
        midterm.termination_status = 0

        if cache !== nothing
            deleteat!(midterm.nbase, cache)
        else
            l = findfirst(x->x==input.n + 1,midterm.base)
            B = view(input.A,:,midterm.base)
            for (k,w) in enumerate(midterm.nbase)
                sol = B \ input.A[:,w]
                if abs.(sol[l]) > input.tol
                    push!(midterm.base, w)
                    deleteat!(midterm.nbase, k)
                    deleteat!(midterm.base, l)
                    break
                end
                if k == length(midterm.nbase)
                    input.A = input.A[1:end .!= l, :]
                    input.b = input.b[1:end .!= l]
                    deleteat!(midterm.base, l)
                    input.m -= 1
                    break
                end
            end
        end
        input.A = input.A[:,1:input.n]
    else
        midterm.base =  collect((input.n - input.m + 1):input.n)
        midterm.nbase = collect(1:(input.n-input.m))
    end

    init_log_simplex2(input)
    while midterm.termination_status == 0 && midterm.iter < input.max_iter
        iterate!(input, midterm)
    end
    update_midterm!(input, midterm)
    output = write_output(input, midterm)
    return output
end

function check_phase_1(input::Input)
    b_min = minimum(input.b)
    return b_min < 0
end

function iterate!(input::Input, midterm::MidTerm)

    midterm.iter += 1
    A = input.A
    b = input.b
    c = input.c
    base = midterm.base
    nbase = midterm.nbase
    tol = input.tol
    B = view(A,:,base)
    N = view(A,:,nbase)
    xB = B \ b
    y = B' \ c[base]
    midterm.red_cost = c[nbase] - N'*y
    val = maximum(midterm.red_cost)

    if val <= tol
        midterm.termination_status = 1
        return #optimal
    end
    midterm.j = findfirst(x->x>tol,midterm.red_cost)

    d = zeros(length(c))
    d_base = B \ N[:,midterm.j]
    d[base] = - d_base
    d[nbase[midterm.j]] = 1
    midterm.d = d
    d_base = max.(d_base, 0)
    r = max.(xB, tol) ./ d_base
    val, midterm.i = findmin(r)

    if val == Inf
        midterm.termination_status = 2
        return #unbounded
    end
    midterm.z = c[base]'xB
    midterm.x = zeros(input.n)
    midterm.x[base] = xB
    iteration_log(input, midterm)
    base[midterm.i], nbase[midterm.j] = nbase[midterm.j], base[midterm.i]
    return #max iteration
end


text:
# Calculator.jl
module Calculator

# include sub-modules
include("Mortgage.jl")
using .Mortgage: payment

# functions for the main module
include("funcs.jl")

end # module

ext: jl


text:
using StatsBase
using StatsBase: PValue, TestStat
using Test, Random

v1 = [1.45666, -23.14, 1.56734e-13]
v2 = ["Good", "Great", "Bad"]
v3 = [1, 56, 2]
v4 = [-12.56, 0.1326, 2.68e-16]
v5 = [0.12, 0.3467, 1.345e-16]
ct = CoefTable(Any[v1, v2, v3, v4, v5],
               ["Estimate", "Comments", "df", "t", "p"],
               ["x1", "x2", "x3"], 5, 4)
ct_noname = CoefTable(Any[v1, v2, v3, v4, v5],
                      ["Estimate", "Comments", "df", "t", "p"],
                      [], 5, 4)

@test sprint(show, ct) == """
───────────────────────────────────────────────
         Estimate  Comments  df       t       p
───────────────────────────────────────────────
x1    1.45666         Good    1  -12.56  0.1200
x2  -23.14            Great  56    0.13  0.3467
x3    1.56734e-13     Bad     2    0.00  <1e-15
───────────────────────────────────────────────"""

@test sprint(show, ct_noname) == """
────────────────────────────────────────────────
          Estimate  Comments  df       t       p
────────────────────────────────────────────────
[1]    1.45666         Good    1  -12.56  0.1200
[2]  -23.14            Great  56    0.13  0.3467
[3]    1.56734e-13     Bad     2    0.00  <1e-15
────────────────────────────────────────────────"""

@test sprint(show, MIME"text/markdown"(), ct) == """
|    |      Estimate | Comments | df |      t |      p |
|:---|--------------:|---------:|---:|-------:|:-------|
| x1 |   1.45666     |    Good  |  1 | -12.56 | 0.1200 |
| x2 | -23.14        |    Great | 56 |   0.13 | 0.3467 |
| x3 |   1.56734e-13 |    Bad   |  2 |   0.00 | <1e-15 |"""

@test sprint(show, MIME"text/markdown"(), ct_noname) == """
|     |      Estimate | Comments | df |      t |      p |
|:----|--------------:|---------:|---:|-------:|:-------|
| [1] |   1.45666     |    Good  |  1 | -12.56 | 0.1200 |
| [2] | -23.14        |    Great | 56 |   0.13 | 0.3467 |
| [3] |   1.56734e-13 |    Bad   |  2 |   0.00 | <1e-15 |"""

@test length(ct) === 3
@test eltype(ct) ==
    NamedTuple{(:Name, :Estimate, :Comments, :df, :t, :p),
               Tuple{String,Float64,String,Int,Float64,Float64}}
@test collect(ct) == [
    (Name = "x1", Estimate = 1.45666, Comments = "Good", df = 1, t = -12.56, p = 0.12)
    (Name = "x2", Estimate = -23.14, Comments = "Great", df = 56, t = 0.1326, p = 0.3467)
    (Name = "x3", Estimate = 1.56734e-13, Comments = "Bad", df = 2, t = 2.68e-16, p = 1.345e-16)
]

m = [0.11258244478647295  0.05664544616214151  0.38181274408522614  0.8197779704008801
     0.36831406658084287  0.12078054506961555  0.8151038332483567   0.6699313951612162
     0.3444540231363058   0.17957407667101322  0.2422083248151139   0.4530583319523316]
ct = CoefTable(m, ["Estimate", "Stderror", "df", "p"], [], 4)
@test sprint(show, ct) == """
──────────────────────────────────────────
     Estimate   Stderror        df       p
──────────────────────────────────────────
[1]  0.112582  0.0566454  0.381813  0.8198
[2]  0.368314  0.120781   0.815104  0.6699
[3]  0.344454  0.179574   0.242208  0.4531
──────────────────────────────────────────"""
@test length(ct) === 3
@test eltype(ct) ==
    NamedTuple{(:Estimate, :Stderror, :df, :p),
               Tuple{Float64,Float64,Float64,Float64}}
@test collect(ct) == [
    (Estimate = 0.11258244478647295, Stderror = 0.05664544616214151,
     df = 0.38181274408522614, p = 0.8197779704008801)
    (Estimate = 0.36831406658084287, Stderror = 0.12078054506961555,
     df = 0.8151038332483567, p = 0.6699313951612162)
    (Estimate = 0.3444540231363058, Stderror = 0.17957407667101322,
     df = 0.2422083248151139, p = 0.4530583319523316)
]

@test sprint(show, PValue(1.0)) == "1.0000"
@test sprint(show, PValue(1e-1)) == "0.1000"
if VERSION > v"1.6.0-DEV"
    @test sprint(show, PValue(1e-5)) == "<1e-04"
else
    @test sprint(show, PValue(1e-5)) == "<1e-4"
end
@test sprint(show, PValue(NaN)) == "NaN"
@test_throws ErrorException PValue(-0.1)
@test_throws ErrorException PValue(1.1)

@test sprint(show, TestStat(NaN)) == "NaN"
@test sprint(show, TestStat(1e-1)) == "0.10"
@test sprint(show, TestStat(1e-5)) == "0.00"
@test sprint(show, TestStat(π)) == "3.14"

@testset "Union{PValue, TestStat} is Real" begin
    vals = [0.0, Rational(1,3), NaN]
    for T in [PValue, TestStat],
        f in (==, <, ≤, >, ≥, isless, isequal),
        lhs in vals, rhs in vals
        # make sure that T behaves like a Real,
        # regardless of whether it's on the LHS, RHS or both
        @test f(T(lhs), T(rhs)) == f(lhs, rhs)
        @test f(lhs, T(rhs)) == f(lhs, rhs)
        @test f(T(lhs), rhs) == f(lhs, rhs)
    end

    # the (approximate) equality operators get a bit more attention
    for T in [PValue, TestStat]
        @test T(Rational(1,3)) ≈ T(1/3)
        @test Rational(1,3) ≈ T(1/3) atol=0.01
        @test T(Rational(1,3)) isa Real
        @test T(T(0.05)) === T(0.05)
        @test hash(T(0.05)) == hash(0.05)
        @test hash(T(0.05), UInt(42)) == hash(0.05, UInt(42))
    end
end

@test sprint(showerror, ConvergenceException(10)) == "failure to converge after 10 iterations."

@test sprint(showerror, ConvergenceException(10, 0.2, 0.1)) ==
    "failure to converge after 10 iterations. Last change (0.2) was greater than tolerance (0.1)."

@test sprint(showerror, ConvergenceException(10, 0.2, 0.1, "Try changing maxIter.")) ==
    "failure to converge after 10 iterations. Last change (0.2) was greater than tolerance (0.1). Try changing maxIter."

err = @test_throws ArgumentError ConvergenceException(10,.1,.2)
@test err.value.msg == "Change must be greater than tol."


text:
@discord_object struct User
    id::Snowflake
    username::String
    discriminator::String
    avatar::String
    bot::Bool
    system::Bool
    mfa_enabled::Bool
    locale::String
    verified::Bool
    email::String
    flags::Int
    premium_type::UserPremiumType.UserPremiumTypeEnum
    public_flags::Int
    # TODO: Add a PartialGuildMember struct or something?
    # https://discord.com/developers/docs/resources/channel#message-object-message-structure
end


text:
"""
    ConstantPoissonInput{N, L}

Iterate over a Poisson input for `N` neurons with rate `λ[n]` for each neuron `n` `L` times.
"""
struct ConstantPoissonInput{N, L} <: AbstractPoissonInput{N, L}
    λs::Vector{Float64}
    buffer::Vector{Int64}
    start::Int64
    stop::Int64

    function ConstantPoissonInput{N, L}(λs, buffer, start, stop) where {N, L}
        if !(typeof(N) <: Integer)
            error("Number of neurons N must be an integer")
        elseif !(typeof(L) <: Integer)
            error("Input length L must be an integer")
        elseif length(λs) != N
            error("λs must have length N.")
        elseif length(buffer) != N
            error("buffer must have length N.")
        end

        new(λs, buffer, start, stop)
    end
end

ConstantPoissonInput(λs::Vector{<:Real}, L::Integer; start::Integer=one(Int64), stop::Integer=L) = ConstantPoissonInput{length(λs), L}(λs, zeros(Int64, N), start, stop)
ConstantPoissonInput(λ::Real, N::Integer, L::Integer; start::Integer=one(Int64), stop::Integer=L) = ConstantPoissonInput{N, L}(fill(λ, N), zeros(Int64, N), start, stop)

"""
    rates(iter::ConstantPoissonInput{N, L}, i::Integer)::Vector{Float64} where {N, L}

Obtain the firing rate for each neuron in `iter` for iteration `i`.
"""
function rates(iter::ConstantPoissonInput{N, L}, i::Integer)::Vector{Float64} where {N, L}
    return iter.λs
end


text:
using SparseArrays, MatrixNetworks, Graphs, SimpleWeightedGraphs
using SimpleValueGraphs
# using IndexedGraphs
include("../src/IndexedGraphs.jl")
using .IndexedGraphs
using Random
using LinearAlgebra

# build directed graph
N = 10_000
W = sprand(N, N, 5/N)
# remove self loops
for i in 1:size(W,2); W[i,i] = 0; end
dropzeros!(W)
W = W + W'
w = sparse(UpperTriangular(W)).nzval    # unique wights for Graph
wD = nonzeros(permutedims(W))           # duplicate weights for DiGraph

# pick source at random
rng = MersenneTwister(0)
s = rand(rng, 1:size(W,2))

g = IndexedGraph(W)
gD = IndexedDiGraph(W)
g_MN = MatrixNetwork(W)
g_Graphs = SimpleGraph(W)
g_SWG = SimpleWeightedGraph(W)
g_SVG = ValGraph(g_Graphs,
    edgeval_types=(Float64,), edgeval_init=(s, d) -> (W[s,d],))

ds_IG = dijkstra_shortest_paths(g, s, w)
d_IG = ds_IG.dists; p_IG  = ds_IG.parents

ds_IDG = dijkstra_shortest_paths(gD, s, wD)
d_IDG = ds_IDG.dists; p_IDG  = ds_IDG.parents

d_MN, p_MN = MatrixNetworks.dijkstra(g_MN, s)

ds_Graphs = dijkstra_shortest_paths(g_Graphs, s, W)
d_Graphs = ds_Graphs.dists; p_Graphs = ds_Graphs.parents

ds_SWG = dijkstra_shortest_paths(g_SWG, s)
d_SWG = ds_SWG.dists; p_SWG = ds_SWG.parents

ds_SVG = SimpleValueGraphs.Experimental.dijkstra_shortest_paths(g_SVG, s)
d_SVG = ds_SVG.dists; p_SVG = ds_SVG.parents

# check that results are correct
@assert d_IG == d_MN == d_Graphs == d_SWG == d_SVG == d_IDG
@assert p_IG == p_MN == p_Graphs == p_SWG == p_SVG == p_IDG

### BENCHMARK
using BenchmarkTools
println("IndexedDiGraph:")
@btime dijkstra_shortest_paths($gD, $s, $wD)
println("IndexedGraph:")
@btime dijkstra_shortest_paths($g, $s, $w)
println("MatrixNetwork:")
@btime MatrixNetworks.dijkstra($g_MN, $s)
println("SimpleGraph")
@btime dijkstra_shortest_paths($g_Graphs, $s, $W)
println("SimpleWeightedGraph:")
@btime dijkstra_shortest_paths($g_SWG, $s);
println("ValGraph (SimpleValueGraphs.Experimental):")
@btime SimpleValueGraphs.Experimental.dijkstra_shortest_paths($g_SVG, $s);


text:
using MathOptInterface

const MOI = MathOptInterface
const CI = MOI.ConstraintIndex
const VI = MOI.VariableIndex

const MOIU = MOI.Utilities

mutable struct MOISolution
    status::Int
    primal::Vector{Float64} # primal of variables
    slack::Vector{Float64}
    dual::Vector{Float64} # dual of constraints
    obj_val::Float64
    dual_objval::Float64
    solve_time::Float64
    iter::Int
end
MOISolution() = MOISolution(0, Float64[], Float64[], Float64[], NaN, NaN, NaN, 0)

mutable struct ModelData
    m::Int # Numbero of constraints
    n::Int # Number of variables
    I::Vector{Int} # List of rows
    J::Vector{Int} # List of cols
    V::Vector{Float64} # List of coefficients
    b::Vector{Float64} # constants
    obj_constant::Float64 # The objective is min c'x + obj_constant
    c::Vector{Float64} # Objective coefficients
end

mutable struct ConeData
    f::Int # length of the zero cone (equality constraints)
    l::Int # length of the nonnegatives cone (<= constraints)
    nrows::Dict{Int, Int} # The number of rows of each vector sets, this is used by `constrrows` to recover the number of rows used by a constraint when getting `ConstraintPrimal` or `ConstraintDual`
    function ConeData()
        return new(0, 0, Dict{Int, Int}())
    end
end

mutable struct Optimizer <: MOI.AbstractOptimizer
    cone::ConeData
    data::Union{Nothing, ModelData} # only non-Void between MOI.copy_to and MOI.optimize!
    maxsense::Bool
    sol::MOISolution
    params::Params
    function Optimizer(kwargs...)
        return new(ConeData(), nothing, false, MOISolution(), Params(kwargs...))
    end
end
function Optimizer(;args...)
    return Optimizer(args)
end

function MOI.is_empty(optimizer::Optimizer)
    return !optimizer.maxsense && optimizer.data === nothing
end
function MOI.empty!(optimizer::Optimizer)
    optimizer.maxsense = false
    optimizer.data = nothing # It should already be nothing except if an error is thrown inside copy_to
    optimizer.sol.status = 0
    return
end

MOI.supports(::Optimizer, ::MOI.Silent) = true
function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool)
    optimizer.params.verbose = value
end
MOI.get(optimizer::Optimizer, ::MOI.Silent) = optimizer.params.verbose

MOIU.supports_allocate_load(model::Optimizer, copy_names::Bool) = true

function MOI.supports(::Optimizer,
                      ::Union{MOI.ObjectiveSense,
                              MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}})
    return true
end

const SF = MOI.VectorAffineFunction{Float64}
const SS = Union{MOI.Zeros, MOI.Nonnegatives}
function MOI.supports_constraint(optimizer::Optimizer, ::Type{<:SF}, ::Type{<:SS})
    return true
end
function MOI.supports_constraint(::Optimizer, ::Type{MOI.AbstractFunction}, ::Type{MOI.AbstractSet})
    return false
end

function MOI.supports(::Optimizer, ::MOI.VariablePrimalStart,
                      ::Type{MOI.VariableIndex})
    return false
end

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; kws...)
    return MOIU.automatic_copy_to(dest, src; kws...)
end

# Computes cone dimensions
constroffset(cone::ConeData, ci::CI{<:MOI.AbstractFunction, MOI.Zeros}) = ci.value
function _allocate_constraint(cone::ConeData, f::MOI.AbstractFunction, s::MOI.Zeros)
    ci = cone.f
    cone.f += MOI.dimension(s)
    return ci
end
constroffset(cone::ConeData, ci::CI{<:MOI.AbstractFunction, MOI.Nonnegatives}) = cone.f + ci.value
function _allocate_constraint(cone::ConeData, f::MOI.AbstractFunction, s::MOI.Nonnegatives)
    ci = cone.l
    cone.l += MOI.dimension(s)
    return ci
end
constroffset(optimizer::Optimizer, ci::CI) = constroffset(optimizer.cone, ci::CI)
function MOIU.allocate_constraint(optimizer::Optimizer, f::F, s::S) where {F <: MOI.AbstractFunction, S <: MOI.AbstractSet}
    return CI{F, S}(_allocate_constraint(optimizer.cone, f, s))
end
constrrows(s::MOI.AbstractVectorSet) = 1:MOI.dimension(s)
# When only the index is available, use the `optimizer.ncone.nrows` field
constrrows(optimizer::Optimizer, ci::CI{<:MOI.AbstractVectorFunction, <:MOI.AbstractVectorSet}) = 1:optimizer.cone.nrows[constroffset(optimizer, ci)]

output_index(t::MOI.VectorAffineTerm) = t.output_index
variable_index_value(t::MOI.ScalarAffineTerm) = t.variable_index.value
variable_index_value(t::MOI.VectorAffineTerm) = variable_index_value(t.scalar_term)
coefficient(t::MOI.ScalarAffineTerm) = t.coefficient
coefficient(t::MOI.VectorAffineTerm) = coefficient(t.scalar_term)

function MOIU.load_constraint(optimizer::Optimizer, ci::CI, f::MOI.VectorAffineFunction, s::MOI.AbstractVectorSet)
    A = sparse(output_index.(f.terms), variable_index_value.(f.terms), coefficient.(f.terms))
    # sparse combines duplicates with + but does not remove zeros created so we call dropzeros!
    dropzeros!(A)
    I, J, V = findnz(A)

    offset = constroffset(optimizer, ci)
    rows = constrrows(s)
    optimizer.cone.nrows[offset] = length(rows)
    i = offset .+ rows

    b = f.constants
    optimizer.data.b[i] = b
    append!(optimizer.data.I, offset .+ I)
    append!(optimizer.data.J, J)
    append!(optimizer.data.V, -V)
    return
end

function MOIU.allocate_variables(optimizer::Optimizer, nvars::Integer)
    optimizer.cone = ConeData()
    return VI.(1:nvars)
end

function MOIU.allocate_constrained_variables(optimizer::Optimizer, set::MOI.Nonnegatives)
    offset = length(optimizer.varmap)
    ci = MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives}(offset + 1)
    return [MOI.VariableIndex(i) for i in offset .+ (1:MOI.dimension(set))], ci
end

function MOIU.allocate_constrained_variables(optimizer::Optimizer, set::MOI.Zeros)
    offset = length(optimizer.varmap)
    ci = MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64}, MOI.Zeros}(offset + 1)
    return [MOI.VariableIndex(i) for i in offset .+ (1:MOI.dimension(set))], ci
end

function MOIU.load_constrained_variables(
    optimizer::Optimizer, vis::Vector{MOI.VariableIndex},
    ci::MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64}},
    set::SS)
end

# allocate block
function MOIU.allocate(::Optimizer, ::MOI.VariablePrimalStart,
                       ::MOI.VariableIndex, ::Union{Nothing, Float64})
end
function MOIU.allocate(::Optimizer, ::MOI.ConstraintPrimalStart,
                       ::MOI.ConstraintIndex,
                       ::Union{Nothing, AbstractVector{Float64}})
end
function MOIU.allocate(::Optimizer, ::MOI.ConstraintDualStart,
                       ::MOI.ConstraintIndex,
                       ::Union{Nothing, AbstractVector{Float64}})
end
function MOIU.allocate(optimizer::Optimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    optimizer.maxsense = sense == MOI.MAX_SENSE
end
function MOIU.allocate(::Optimizer, ::MOI.ObjectiveFunction,
                       ::MOI.Union{MOI.SingleVariable,
                                   MOI.ScalarAffineFunction{Float64}})
end

function MOIU.load_variables(optimizer::Optimizer, nvars::Integer)
    cone = optimizer.cone
    m = cone.f + cone.l
    I = Int[]
    J = Int[]
    V = Float64[]
    b = zeros(m)
    c = zeros(nvars)
    optimizer.data = ModelData(m, nvars, I, J, V, b, 0., c)
    return
end

function MOIU.load(::Optimizer, ::MOI.ObjectiveSense, ::MOI.OptimizationSense)
end
function MOIU.load(optimizer::Optimizer, ::MOI.ObjectiveFunction,
                   f::MOI.SingleVariable)
    MOIU.load(optimizer,
              MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
              MOI.ScalarAffineFunction{Float64}(f))
end
function MOIU.load(optimizer::Optimizer, ::MOI.ObjectiveFunction,
                   f::MOI.ScalarAffineFunction)
    c0 = Vector(sparsevec(variable_index_value.(f.terms), coefficient.(f.terms),
                          optimizer.data.n))
    optimizer.data.obj_constant = f.constant
    optimizer.data.c = optimizer.maxsense ? -c0 : c0
    return nothing
end

# MOI getters
MOI.get(::Optimizer, ::MOI.SolverName) = "LPIP"
MOI.get(optimizer::Optimizer, ::MOI.SolveTime) = optimizer.sol.solve_time
MOI.get(optimizer::Optimizer, ::MOI.PrimalStatus) = optimizer.sol.status in (1, 3) ? MOI.FEASIBLE_POINT : MOI.INFEASIBLE_POINT
MOI.get(optimizer::Optimizer, ::MOI.DualStatus) = optimizer.sol.status in (1, 2) ? MOI.FEASIBLE_POINT : MOI.INFEASIBLE_POINT

function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus)
    s = optimizer.sol.status
    @assert 0 <= s <= 5
    if s == 0
        return MOI.OPTIMIZE_NOT_CALLED
    elseif s == 1
        return MOI.OPTIMAL
    elseif s == 2
        return MOI.INFEASIBLE
    elseif s == 3
        return MOI.DUAL_INFEASIBLE
    elseif s == 4
        return MOI.TIME_LIMIT
    elseif s == 5
        return MOI.ITERATION_LIMIT
    end
end
function MOI.get(optimizer::Optimizer, ::MOI.ResultCount)
    if MOI.get(optimizer, MOI.TerminationStatus()) in (MOI.INFEASIBLE, MOI.DUAL_INFEASIBLE)
        return 0
    else
        return 1
    end
end

function MOI.get(optimizer::Optimizer, ::MOI.VariablePrimal, vi::VI)
    return optimizer.sol.primal[vi.value]
end
function MOI.get(optimizer::Optimizer, a::MOI.VariablePrimal, vi::Vector{VI})
    return MOI.get.(optimizer, a, vi)
end
function MOI.get(optimizer::Optimizer, ::MOI.ObjectiveValue)
    return optimizer.maxsense ? -optimizer.sol.obj_val : optimizer.sol.obj_val
end
function MOI.get(optimizer::Optimizer, ::MOI.ConstraintPrimal,
                 ci::CI{<:MOI.AbstractFunction, S}) where S <: MOI.AbstractSet
    offset = constroffset(optimizer, ci)
    rows = constrrows(optimizer, ci)
    primal = optimizer.sol.slack[offset .+ rows]
    return primal
end
function MOI.get(optimizer::Optimizer, ::MOI.ConstraintDual,
                 ci::CI{<:MOI.AbstractFunction, S}) where S <: MOI.AbstractSet
    offset = constroffset(optimizer, ci)
    rows = constrrows(optimizer, ci)
    return -optimizer.sol.dual[offset .+ rows]
end

# MOI setters

# MOI optimize
function MOI.optimize!(optimizer::Optimizer)
    cone = optimizer.cone
    A = sparse(optimizer.data.I, optimizer.data.J, optimizer.data.V)
    b = optimizer.data.b
    obj_constant = optimizer.data.obj_constant
    c = optimizer.data.c
    # optimizer.data = nothing # Allows GC to free optimizer.data before A is loaded to LPIP

    lpip_pb = LPIP.LPIPLinearProblem{Float64}(A, b, c, obj_constant, cone.l)

    t0 = time()
    sol = LPIP.interior_points(lpip_pb, optimizer.params)
    solve_time = time() - t0

    # Query solutions
    status = sol.status
    primal = sol.variables.x
    slack = sol.variables.x[end - cone.l:end]
    dual = sol.variables.p
    obj_val = sol.obj_val
    dual_obj_val = sol.dual_obj_val
    iter = sol.iter
    optimizer.sol = MOISolution(status, primal, slack, dual, obj_val, dual_obj_val, solve_time, iter)
    return
end


text:
using SafeTestsets
@safetestset "MyExample.jl" begin include("my_f_tests.jl") end

max_stars_repo_name: DomagojGalic/MyExample

max_issues_repo_name: DomagojGalic/MyExample

max_forks_repo_name: DomagojGalic/MyExample


text:
#!/bin/julia
using NetCDF
using PyPlot
using Dierckx
using ..constants:mp,me
using ..gs2ls: GS2specs, GS2geometry

"""
plotSummary(file::String)

Displays a set of time-dependent plots of relevant plasma quantities. Useful for determining the overall health of the run and ideal time samples.
"""
function plotTranspSummary(file::String)

   t = ncread(file,"TIME")
   Ti0 = ncread(file,"TI0")
   Te0 = ncread(file,"TE0")
   q95 = ncread(file,"Q95")
   ne0 = ncread(file,"NE")[1,:]'

   #TODO: Add some kind of global diagnostic of fast ions

end

"""
plotSummary(file::String,time::Float64)

Displays a set of radius-dependent plots of relevant plasma quantities at the given time.
"""
function plotTranspSummary(file::String,time::Float64)

   # TODO: Make radial plots
   kappa = ncread(file,"ELONG")[1,:]
   tri = ncread(file,"TRIANG")[1,:]

end

function getValFrom2Darray(rho::Float64,time::Float64,rgrid::Vector,tgrid::Vector,y::Matrix)
   yfunc = Spline2D(rgrid,tgrid,y,bc="nearest")
   return evaluate(yfunc,rho,time)
end

function getGradient(rho::Float64,time::Float64,rgrid::Vector,tgrid::Vector,y::Matrix)
   dr = 0.01
   yfunc = Spline2D(rgrid,tgrid,y,bc="nearest")
   y2 = evaluate(yfunc,rho+dr,time)
   y1 = evaluate(yfunc,rho-dr,time)
   return (y2-y1)/(2.0*dr)
end

function getArrayAtTime(t::Float64,tgrid::Vector,ygrid::Vector)
   nx = length(ygrid[:,1])
   x = zeros(Float64,nx)
   for i in 1:nx
      yfunc = Spline1D(tgrid,y[i,:]',bc="nearest")
      x[i]  = evaluate(yfunc,t)
   end
   return x
end

function getLocalData(file::String,rho::Float64,t0::Float64,Zi::Float64,Ai::Float64;Zz::Float64=-999.9,Az::Float64=-999.9)

   t = ncread(file,"TIME")
   x = getArrayAtTime(t0,t,ncread(file,"X"))
   xb = getArrayAtTime(t0,t,ncread(file,"XB"))

   beta_e = getValFrom2Darray(rho,t0,x,t,ncread(file,"BTE"))

   if (Zz == -999.9) && (Az == -999.9)
      include_impurity = false
   else
      include_impurity = true
   end

# Params and species:
# To return: beta_e, [Ti, Te, Tz, ni, ne, nz, + gradients]
# Corresponding: BTE, NI, NE, XDENS, TI, TE, TX

# Geometry:
# qinp, shat, [r_geo, shift], rmaj, akappa,akappri,tri,tripri,
# Q, SHAT, [RMJMP,gradient], RAXIS, [ELONG, TRIANG + gradients]

# Other relevant quantities:
# BDENS: beam ion density. How does this affect quasineutrality

   # TODO: Define nuii, nuee, nuzz

   specs=GS2species[]
   push!(specs,GS2species(dens=ni/ne,temp=1.0,z=Zi,mass=1.0,vnewk=nuii*a/vti,fprim=gradni*a,trpim=gradTi*a))
   if include_impurity
      push!(specs,GS2species(dens=nz/ne,temp=Tz/Ti,z=Zimp,mass=Aimp/Ai,vnewk=nuzz*a/vti,fprim=gradnz*a,trpim=gradTz*a))
   end
   push!(specs,GS2species(dens=1.0,temp=Te/Ti,z=-1.0,mass=me/(Ai*mp),vnewk=nuee*a/vti,fprim=gradne*a,trpim=gradTe*a))

   geom = GS2geometry(equilibrium_option="eik",rhoc=rho,irho=1,qinp=q,shat=shat,
                       shift=Rprime,akappa=kappa,akappri=kappri,tri=tri,tripri=tripri)

   return beta_e, geom, specs
end

max_stars_repo_name: gjwilkie/PlasmaTools

max_issues_repo_name: gjwilkie/PlasmaTools

max_forks_repo_name: gjwilkie/PlasmaTools


text:
using Printf
using Statistics
using LinearAlgebra
using DelimitedFiles
using IterativeSolvers
using SparseArrays
using TickTock

# rearrange the edge sequence in the increasing order of layer id and starting node id
function edge_seq_rearrangement(edge_seq::Array{Int64,2})
    layer_list = setdiff(edge_seq[:,4],[])
    L = maximum(layer_list)
    if length(setdiff([i for i = 1:L],layer_list)) != 0
        println("please label layers in the order 1, 2, ..., L")
    end
    node_list = setdiff(edge_seq[:,1],[])
    N = maximum(edge_seq[:,1])
    if length(setdiff([i for i = 1:N],node_list)) != 0
        println("please label nodes in the order 1, 2, ..., N")
    end

    edge_seq_rearranged = zeros(Int64,size(edge_seq))
    id = 1
    for ell in layer_list
        edge_list_ell = findall(isequal(ell),edge_seq[:,4])
        M = sparse(edge_seq[edge_list_ell,1],edge_seq[edge_list_ell,2],edge_seq[edge_list_ell,3])
        inds = findall(!iszero,M)
        a = getindex.(inds, 1)
        b = getindex.(inds, 2)
        edge_seq_rearranged[id:id+length(a)-1,1] = b
        edge_seq_rearranged[id:id+length(a)-1,2] = a
        edge_seq_rearranged[id:id+length(a)-1,3] = M[inds]
        edge_seq_rearranged[id:id+length(a)-1,4] = ell*ones(Int64,length(a))
        id = id+length(a)
    end

    return edge_seq_rearranged
end

# solve a set of equations to obtain beta_i, beta_ij, and gamma_ij, for any number of layers
# this function is used for rhoA > (rhoA)^o in layer one under a given strategy configuration
# input---
# multilayer network in the form of an edge sequanece: edge_seq [node id, node id, edge weight, layer id], where edge weights here are set to be any positive integer
# initial strategy configuration: xi_seq [node id, strategy, layer id]
# output---
# the first column of "solution" gives beta_ij
# the other columns of "solution" gives gamma_ij, where the ll_th column is the strategy assortment between layer 1 and layter ll
function beta_gamma(edge_seq::Array{Int64,2}, xi_seq::Array{Int64,2})
    N = maximum(edge_seq[:,1])
    L = maximum(edge_seq[:,4])
    M = spzeros(Float64,N*L,N)
    for ell = 1:L
        edge_list_ell = findall(isequal(ell),edge_seq[:,4])
        M_ell = sparse(edge_seq[edge_list_ell,1],edge_seq[edge_list_ell,2],edge_seq[edge_list_ell,3],N,N)
        M[(ell-1)*N+1:ell*N,:] = M_ell
    end
    pi = spzeros(Float64,L*N)
    for ell = 1:L
        pi_2 = sum(M[(ell-1)*N+1:ell*N,:], dims = 2)
        pi_2 = pi_2[:,1]
        pi_transient = spzeros(Float64, N)
        presence = findall(!iszero,pi_2)
        pi_transient[presence] = pi_2[presence]/sum(pi_2)
        pi[(ell-1)*N+1:ell*N] = pi_transient

        M_transient = M[(ell-1)*N+1:ell*N,:]
        M_transient[presence,:] = M_transient[presence,:]./pi_2[presence]
        M[(ell-1)*N+1:ell*N,:] = M_transient
    end

    xi = spzeros(Float64,L*N)
    for ell = 1:L
        xi_transient = spzeros(Int64,N)
        presence = findall(isequal(ell),xi_seq[:,3])
        xi_transient[xi_seq[presence,1]] = xi_seq[presence,2]
        xi[(ell-1)*N+1:ell*N] = xi_transient
    end

    # obtain beta_i
    presence_list = findall(!iszero,pi[1:N])
    N1 = length(presence_list)
    absence_list = findall(iszero,pi[1:N])
    MatA_i = M[1:N,:]-sparse(Matrix(1.0I, N, N))
    pi_list = findall(!isequal(presence_list[1]),[i for i = 1:N])
    MatA_i[:,pi_list] = MatA_i[:,pi_list] - MatA_i[:,presence_list[1]]*transpose(pi[pi_list]/pi[presence_list[1]])
    MatA_i_reduced = MatA_i[pi_list,pi_list]

    xi_hat = sum(pi[1:N].*xi[1:N])
    MatB_i = (xi_hat*ones(Float64,N)-xi[1:N])*N1
    MatB_i_reduced = MatB_i[pi_list]

    beta_i_solution_reduced = idrs(MatA_i_reduced, MatB_i_reduced)
    beta_i_solution = zeros(Float64, N)
    beta_i_solution[pi_list] = beta_i_solution_reduced
    beta_i_solution[presence_list[1]] = -sum(pi[1:N].*beta_i_solution)/pi[presence_list[1]]

    # obtain beta_ij
    Mat1 = M[1:N,:]
    inds = findall(!iszero,Mat1)
    a = getindex.(inds, 1)
    b = getindex.(inds, 2)
    vec1 = [i for i = 1:N]
    vec1 = transpose(vec1)
    X1 = N*(a.-1)*ones(Int64,1,N)+ones(Int64,length(a))*vec1
    Y1 = N*(b.-1)*ones(Int64,1,N)+ones(Int64,length(b))*vec1
    W1 = Mat1[inds]*ones(Int64,1,N)

    vec2 = [(i-1)*N for i = 1:N]
    vec2 = transpose(vec2)
    X2 = a*ones(Int64,1,N)+ones(Int64,length(a))*vec2
    Y2 = b*ones(Int64,1,N)+ones(Int64,length(b))*vec2
    W2 = Mat1[inds]*ones(Int64,1,N)

    X11 = reshape(X1,:,1)
    Y11 = reshape(Y1,:,1)
    W11 = reshape(W1,:,1)
    X22 = reshape(X2,:,1)
    Y22 = reshape(Y2,:,1)
    W22 = reshape(W2,:,1)
    MatA_ij = sparse(X11[:,1],Y11[:,1],W11[:,1],N^2,N^2)/2+sparse(X22[:,1],Y22[:,1],W22[:,1],N^2,N^2)/2
    MatA_ij = MatA_ij - sparse(Matrix(1.0I, N*N, N*N))
    beta_obtained = [(i-1)*N+i for i = 1:N]
    beta_missed = setdiff([i for i = 1:N^2], beta_obtained)
    MatA_ij_reduced = MatA_ij[beta_missed, beta_missed]

    MatB_ij = -sum(MatA_ij[:,beta_obtained].*transpose(beta_i_solution), dims = 2)
    Mat = xi_hat*ones(Float64,N,N)*N1/2-xi[1:N]*transpose(xi[1:N])*N1/2
    Mat[absence_list,:] = spzeros(Float64,length(absence_list),N)
    Mat[:,absence_list] = spzeros(Float64,N,length(absence_list))
    MatB_ij = MatB_ij+reshape(transpose(Mat),:)
    MatB_ij_reduced = MatB_ij[beta_missed]

    beta_ij_solution_reduced = idrs(MatA_ij_reduced, MatB_ij_reduced)
    beta_ij_solution = zeros(Float64, N*N)
    beta_ij_solution[beta_missed] = beta_ij_solution_reduced
    beta_ij_solution[beta_obtained] = beta_i_solution

    solution = spzeros(Float64, N*N, L)
    solution[:,1] = beta_ij_solution
    # obtain gamma_ij
    if L > 1
        for ell = 2:L
            xi_ell = xi[(ell-1)*N+1:ell*N]
            pi_ell = pi[(ell-1)*N+1:ell*N]
            presence_ell = findall(!iszero, pi_ell)
            N_ell = length(presence_ell)
            GammaA_ij = sparse(X11[:,1],Y11[:,1],W11[:,1],N^2,N^2)*(N_ell-1)/(N1+N_ell-1)

            Mat_ell = M[(ell-1)*N+1:ell*N,:]
            inds_ell = findall(!iszero,Mat_ell)
            a_ell = getindex.(inds_ell, 1)
            b_ell = getindex.(inds_ell, 2)
            vec_ell = [(i-1)*N for i = 1:N]
            vec_ell = transpose(vec_ell)
            X_ell = a_ell*ones(Int64,1,N)+ones(Int64,length(a_ell))*vec_ell
            Y_ell = b_ell*ones(Int64,1,N)+ones(Int64,length(b_ell))*vec_ell
            W_ell = Mat_ell[inds_ell]*ones(Int64,1,N)
            X_ellell = reshape(X_ell,:,1)
            Y_ellell = reshape(Y_ell,:,1)
            W_ellell = reshape(W_ell,:,1)
            GammaA_ij = GammaA_ij+sparse(X_ellell[:,1],Y_ellell[:,1],W_ellell[:,1],N^2,N^2)*(N1-1)/(N1+N_ell-1)

            vec_gamma1 = reshape(ones(Int64, N)*vec1,:)
            vec_gamma2 = reshape(transpose(vec1)*ones(Int64,1,N),:)
            X_gamma1 = N*(a.-1)*ones(Int64,1,N*N)+ones(Int64,length(a))*transpose(vec_gamma1)
            Y_gamma1 = N*(b.-1)*ones(Int64,1,N*N)+ones(Int64,length(b))*transpose(vec_gamma2)
            W_gamma1 = Mat1[inds]*ones(Int64,1,N*N)

            vec_gamma3 = [(i-1)*N for i in vec_gamma1]
            vec_gamma4 = [(i-1)*N for i in vec_gamma2]
            X_gamma2 = a_ell*ones(Int64,1,N*N)+ones(Int64,length(a_ell))*transpose(vec_gamma3)
            Y_gamma2 = b_ell*ones(Int64,1,N*N)+ones(Int64,length(b_ell))*transpose(vec_gamma4)
            W_gamma2 = Mat_ell[inds_ell]*ones(Int64,1,N*N)

            X_gamma_11 = reshape(X_gamma1,:,1)
            Y_gamma_11 = reshape(Y_gamma1,:,1)
            W_gamma_11 = reshape(W_gamma1,:,1)
            X_gamma_22 = reshape(X_gamma2,:,1)
            Y_gamma_22 = reshape(Y_gamma2,:,1)
            W_gamma_22 = reshape(W_gamma2,:,1)
            GammaA_ij = GammaA_ij + sparse(X_gamma_11[:,1],Y_gamma_11[:,1],W_gamma_11[:,1],N^2,N^2).*sparse(X_gamma_22[:,1],Y_gamma_22[:,1],W_gamma_22[:,1],N^2,N^2)/(N1+N_ell-1)

            GammaA_ij = GammaA_ij - sparse(Matrix(1.0I, N*N, N*N))
            pi_transient = pi[1:N].*pi_ell
            presence_list = findall(!iszero,pi_transient)
            absence_list = findall(iszero,pi_transient)
            GammaA_delete = (presence_list.-1)*N+presence_list
            GammaA_ij[:,GammaA_delete] = GammaA_ij[:,GammaA_delete] - GammaA_ij[:,GammaA_delete[1]]*transpose(pi[presence_list])/pi[presence_list[1]]
            pi_list = findall(!isequal(GammaA_delete[1]),[i for i = 1:N*N])
            GammaA_ij_reduced = GammaA_ij[pi_list,pi_list]

            xi_ell_hat = sum(pi_ell.*xi_ell)
            GammaB_ij = (xi_hat*xi_ell_hat*ones(Float64,N,N) - xi[1:N]*transpose(xi_ell))*N1*N_ell/(N1+N_ell-1)
            presence1 = findall(iszero, pi[1:N])
            presence1 = presence1[:,1]
            presence_ell = findall(iszero, pi_ell)
            presence_ell = presence_ell[:,1]
            GammaB_ij[presence1,:] = spzeros(Float64, length(presence1), N)
            GammaB_ij[:,presence_ell] = spzeros(Float64, N, length(presence_ell))
            GammaB_ij  = reshape(transpose(GammaB_ij), :)
            GammaB_ij_reduced = GammaB_ij[pi_list]

            gamma_ij_solution_reduced = idrs(GammaA_ij_reduced, GammaB_ij_reduced)
            gamma_ij_solution = zeros(Float64, N*N)
            gamma_ij_solution[pi_list] = gamma_ij_solution_reduced
            gamma_ij_solution[GammaA_delete[1]] = -sum(pi[presence_list].*gamma_ij_solution[GammaA_delete])/pi[presence_list[1]]
            solution[:,ell] = gamma_ij_solution
        end
    end

    return solution
end

# provide the value of coefficient of Eq. (2), namely theta and phi
# input---
# multilayer network in the form of an edge sequanece: edge_seq [node id, node id, edge weight, layer id]
# initial strategy configuration: xi_seq [node id, strategy, layer id]
# for an edge between node i and j in layer ll with weight wij_ll, two terms should be added into the edge sequence, i.e. [i j wij_ll ll] and [j i wij_ll ll]
# the code here requires wij_ll to be a positive integer
# output---
# strategy assortment:  theta_phi[:,1] = [theta0, theta1, theta2, theta3], strategy correlation within layer 1
#                       theta_phi[:,ll] = [phi00, phi01, phi20, phi21], strategy correlation between layer 1 and layer ll
# if there are only two layers, in Eq. (2) in the main text, [theta0, theta1, theta2, theta3] = theta_phi[:,1] and [phi00, phi01, phi20, phi21] = theta_phi[:,2]
function bc_multilayer_DB(edge_seq::Array{Int64,2}, xi_seq::Array{Int64,2})
    edge_seq = edge_seq_rearrangement(edge_seq)
    solution = beta_gamma(edge_seq, xi_seq)
    N = maximum(edge_seq[:,1])
    L = maximum(edge_seq[:,4])
    M = spzeros(Float64,N*L,N)
    for ell = 1:L
        edge_list_ell = findall(isequal(ell),edge_seq[:,4])
        M_ell = sparse(edge_seq[edge_list_ell,1],edge_seq[edge_list_ell,2],edge_seq[edge_list_ell,3],N,N)
        M[(ell-1)*N+1:ell*N,:] = M_ell
    end
    pi = spzeros(Float64,L*N)
    for ell = 1:L
        pi_2 = sum(M[(ell-1)*N+1:ell*N,:], dims = 2)
        pi_2 = pi_2[:,1]
        pi_transient = spzeros(Float64, N)
        presence = findall(!iszero,pi_2)
        pi_transient[presence] = pi_2[presence]/sum(pi_2)
        pi[(ell-1)*N+1:ell*N] = pi_transient

        M_transient = M[(ell-1)*N+1:ell*N,:]
        M_transient[presence,:] = M_transient[presence,:]./pi_2[presence]
        M[(ell-1)*N+1:ell*N,:] = M_transient
    end

    M1 = M[1:N,:]
    beta = transpose(reshape(solution[:,1], N, :))
    theta1 = sum(pi[1:N].*M1.*beta)
    M2 = M1*M1
    theta2 = sum(pi[1:N].*M2.*beta)
    M3 = M2*M1
    theta3 = sum(pi[1:N].*M3.*beta)

    theta_phi = zeros(Float64, 4, L)
    theta_phi[:,1] = [0.0,theta1,theta2,theta3]

    if L > 1
        for ell = 2:L
            gamma = transpose(reshape(solution[:,ell], N, :))
            M_ell = M[(ell-1)*N+1:ell*N,:]
            phi01 = sum(pi[1:N].*M_ell.*gamma)
            phi20 = sum(pi[1:N].*M2.*gamma)
            M_ell21 = M2*M_ell
            phi21 = sum(pi[1:N].*M_ell21.*gamma)
            theta_phi[:,ell] = [0,phi01,phi20,phi21]
        end
    end

    return theta_phi
end

# an example about the using of above functions
# see Fig. 5A in the main text, where nodes are labelled and initial strategy configuration is presented in SFig. 2
# the edge sequence is given in the form [node id, node id, edge weight, layer id]
edge_seq = [1 6 1 1;
            2 6 1 1;
            3 6 1 1;
            4 6 1 1;
            5 6 1 1;
            6 1 1 1;
            6 2 1 1;
            6 3 1 1;
            6 4 1 1;
            6 5 1 1;
            1 2 1 2;
            1 3 1 2;
            1 4 1 2;
            1 5 1 2;
            1 6 1 2;
            2 1 1 2;
            3 1 1 2;
            4 1 1 2;
            5 1 1 2;
            6 1 1 2]
# the strategy configuration is given in the form [node id, strategy, layer id]
xi_seq = [1 1 1;
          2 0 1;
          3 0 1;
          4 0 1;
          5 0 1;
          6 0 1;
          1 0 2;
          2 0 2;
          3 0 2;
          4 0 2;
          5 0 2;
          6 1 2]
# see Eq.(2) in the main text
# theta_phi[:,1] = [theta0, theta1, theta2, theta3], strategy correlation within layer 1
# theta_phi[:,2] = [phi00, phi01, phi20, phi21], strategy correlation between layer 1 and layer 2
theta_phi = bc_multilayer_DB(edge_seq, xi_seq)


text:
using Revise
using Dice
using Dice: num_flips, num_nodes, to_dice_ir

code = @dice begin
    function discrete(p::Vector{Float64})
        mb = length(p)
        v = Vector(undef, mb)
        sum = 1
        for i=1:mb
            v[i] = p[i]/sum
            sum = sum - p[i]
        end

        # println(v)
        ans = DistInt(dicecontext(), mb-1)
        for i=mb-1:-1:1
            ans = if flip(v[i]) DistInt(dicecontext(), i-1) else ans end
        end
        return ans
    end

    function sendChar(key, observation)
        gen = discrete([0.5, 0.25, 0.125, 0.125])

        enc = key + gen
        d = prob_equals(observation, enc[1])
        return d
    end

    key = discrete([0.25, 0.25, 0.25, 0.25])
    d = true
    for i=1:4
        d&=sendChar(key, 2)
    end
    key, d
end

# BDD analysis
bdd = compile(code)
infer(code, :bdd)

# IR analysis
# to_dice_ir(code)
# has_dice_binary() && rundice(code)
# has_dice_binary() && infer(code, :ocaml)


text:
# Basic API for iterated maps

abstract AbstractIteratedMap

type IteratedMap{S,D} <: AbstractIteratedMap
    func::Function
end

# Iterated maps are iterators

# Perhaps one way to do this is to prime the iterator
abstract AbstractIteratedMapIterator

type IteratedMapIterator{S,D} <: AbstractIteratedMapIterator
    start
    func::Function
end

# When you pass a starting value to an IteratedMap,
# you get an IteratedMapIterator ready to go

function Base.call{S,D}(f::IteratedMap{S,D},x)
    @assert length(x) == D
    IteratedMapIterator{S,D}(x,f.func)
end

function Base.call(f::IteratedMapIterator)
    f.func(f.start)
end

# start needs something to start with
Base.start(f::IteratedMapIterator) = (f.start,1)

# This one is easy
Base.next(f::IteratedMapIterator,state) = (state[1],(f.func(state[1]),state[2]+1))

# done implies that there is a stopping point for the map
Base.done(f::IteratedMapIterator,state) = false

Base.eltype{S,D}(::Type{IteratedMapIterator{S,D}}) = D>1 ? Vector{Float64} : Float64

max_stars_repo_name: wkearn/IteratedMaps

max_issues_repo_name: wkearn/IteratedMaps

max_forks_repo_name: wkearn/IteratedMaps


text:
export IHTVariable, IHTVariables, use_A2_as_minor_allele, standardize_snpmatrix
export use_A2_as_minor_allele, project_k!

"""
Object to contain intermediate variables and temporary arrays. Used for cleaner code in L0_reg
"""
struct IHTVariable
    b    :: Vector{Float64}     # the statistical model, most will be 0
    b0   :: Vector{Float64}     # previous estimated model in the mm step
    xb   :: Vector{Float64}     # vector that holds x*b
    xb0  :: Vector{Float64}     # previous xb in the mm step
    xk   :: Matrix{Float64}     # the n by k subset of the design matrix x corresponding to non-0 elements of b
    gk   :: Vector{Float64}     # Numerator of step size μ.   gk = df[idx] is a temporary array of length `k` that arises as part of the gradient calculations. I avoid doing full gradient calculations since most of `b` is zero.
    xgk  :: Vector{Float64}     # Demonimator of step size μ. x * gk also part of the gradient calculation
    idx  :: BitArray{1}         # BitArray indices of nonzeroes in b for A_mul_B
    idx0 :: BitArray{1}         # previous iterate of idx
    r    :: Vector{Float64}     # n-vector of residuals
    df   :: Vector{Float64}     # the gradient: df = -x' * (y - xb)
end

function IHTVariables(
    x :: SnpData,
    y :: Vector{Float64},
    k :: Int64
)
    n, p = x.people, x.snps + 1 #adding 1 for p because we need an intercept

    #check if k is sensible
    if k > p;  throw(ArgumentError("k cannot exceed the number of SNPs")); end
    if k <= 0; throw(ArgumentError("k must be positive integer")); end

    b    = zeros(Float64, p)
    df   = zeros(Float64, p)
    xb   = zeros(Float64, n)
    r    = zeros(Float64, n)
    b0   = zeros(Float64, p)
    xb0  = zeros(Float64, n)
    xk   = zeros(Float64, n, k)
    xgk  = zeros(Float64, n)
    gk   = zeros(Float64, k)
    idx  = falses(p)
    idx0 = falses(p)
    return IHTVariable(b, b0, xb, xb0, xk, gk, xgk, idx, idx0, r, df)
end

immutable IHTResult
    time :: Float64
    loss :: Float64
    iter :: Int64
    beta :: Vector{Float64}
end

function IHTResults(
    time :: Float64,
    loss :: Float64,
    iter :: Int64,
    beta :: Vector{Float64})
    return IHTResult(time, loss, iter, beta)
end

# """
#     project_k!(x, k)

# This function projects a vector `x` onto the set S_k = { y in R^p : || y ||_0 <= k }.
# It does so by first finding the pivot `a` of the `k` largest components of `x` in magnitude.
# `project_k!` then thresholds `x` by `abs(a)`, sending small components to 0.

# Arguments:

# - `b` is the vector to project.
# - `k` is the number of components of `b` to preserve.
# """
# function project_k!(
#     x :: Vector{Float64},
#     k :: Int;
# )
#     a = select(x, k, by = abs, rev = true)
#     for i in eachindex(x)
#         if abs(x[i]) < abs(a)
#             x[i] = 0.0
#         end
#     end
#     return nothing
# end

"""
Returns ω, a constant we need to bound the step size μ to guarantee convergence.
"""

function compute_ω!(
    v         :: IHTVariable,
    snpmatrix :: Matrix{Float64},
)
    #update v.xb
    A_mul_B!(v.xb, snpmatrix, v.b)

    #calculate ω efficiently (old b0 and xb0 have been copied before calling iht!)
    return sqeuclidean(v.b, v.b0) / sqeuclidean(v.xb, v.xb0)
end

"""
This function is needed for testing purposes only.

Converts a SnpArray to a matrix of float64 using A2 as the minor allele. We want this function
because SnpArrays.jl uses the less frequent allele in each SNP as the minor allele, while PLINK.jl
always uses A2 as the minor allele, and it's nice if we could cross-compare the results.
"""
function use_A2_as_minor_allele(snpmatrix :: SnpArray)
    n, p = size(snpmatrix)
    matrix = zeros(n, p)
    for i in 1:p
        for j in 1:n
            if snpmatrix[j, i] == (0, 0); matrix[j, i] = 0.0; end
            if snpmatrix[j, i] == (0, 1); matrix[j, i] = 1.0; end
            if snpmatrix[j, i] == (1, 1); matrix[j, i] = 2.0; end
            if snpmatrix[j, i] == (1, 0); matrix[j, i] = missing; end
        end
    end
    return matrix
end

"""
A function for determining whether or not to backtrack. If all conditions are satisfied,
then we DONT do line search, which means _iht_backtrack need to return TRUE.
"""
function _iht_backtrack(
    v :: IHTVariable,
    ω :: Float64,
    μ :: Float64
)
    μ < 0.99*ω
    #&& sum(v.idx) != 0 &&
    #sum(xor.(v.idx,v.idx0)) != 0
end

"""
This function computes the gradient step v.b = P_k(β + μ∇f(β)), and updates v.idx.
Recall calling axpy! implies v.b = v.b + μ*v.df, but v.df stores an extra negative sign.
"""
function _iht_gradstep(
    v  :: IHTVariable,
    μ  :: Float64,
    k  :: Int;
)
    BLAS.axpy!(μ, v.df, v.b) # take the gradient step: v.b = b + μ∇f(b) (which is a plus since df stores X(-1*(Y-Xb)))
    project_k!(v.b, k)       # P_k( β - μ∇f(β) ): preserve top k components of b
    _iht_indices(v, k)       # Update idx. (find indices of new beta that are nonzero)

    # If the k'th largest component is not unique, warn the user.
    sum(v.idx) <= k || warn("More than k components of b is non-zero! Need: VERY DANGEROUS DARK SIDE HACK!")
end

"""
this function updates finds the non-zero index of b, and set v.idx = 1 for those indices.
"""
function _iht_indices(
    v :: IHTVariable,
    k :: Int
)
    # set v.idx[i] = 1 if v.b[i] != 0 (i.e. find components of beta that are non-zero)
    v.idx .= v.b .!= 0

    # if idx is the 0 vector, v.idx[i] = 1 if i is one of the k largest components
    # of the gradient (stored in v.df), and set other components of idx to 0.
    if sum(v.idx) == 0
        a = select(v.df, k, by=abs, rev=true)
        v.idx[abs.(v.df) .>= abs(a)-2*eps()] .= true
        v.gk .= zeros(sum(v.idx))
    end

    return nothing
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMCoupling.jl/blob/master/LICENSE

using Base.Test

@testset "test_plain_strain_kinematic_coupling.jl" begin
    include("test_plain_strain_kinematic_coupling.jl")
end

@testset "test_dcoupling_2d.jl" begin
    include("test_dcoupling_2d.jl")
end

@testset "test_dcoupling_3d.jl" begin
    include("test_dcoupling_3d.jl")
end


text:
module Papers

export Model, Mapper, Resource, Service

include("Model.jl")
using .Model

include("Mapper.jl")
using .Mapper

include("Service.jl")
using .Service

include("Resource.jl")
using .Resource

end # module


text:
function readin!(cfgdata::cfgStruct,datafileName)
    fn = basename( datafileName )
    extensionName1  = uppercase( splitext( fn)[2]  )

    rawdatafileObj= open( datafileName,"r")

    if  extensionName1  in ( ".RAW3"  ,".RAW4"  ,".RAW1",".RAW2")
            seek( rawdatafileObj,512)
            cfgdata.meterperSeg = read( rawdatafileObj, UInt16 )/5000
            cfgdata.chn2SegBegIdx = read( rawdatafileObj, UInt16 )
            cfgdata.reflectorFactor = read( rawdatafileObj, UInt16 )/10000
            cfgdata.attenutionFactor = read( rawdatafileObj, UInt16 )/10000
            cfgdata.scanrate = read( rawdatafileObj, UInt16 )
            cfgdata.calSamples = read( rawdatafileObj, UInt16 )
            cfgdata.SegNumber = read( rawdatafileObj, UInt16 )
    end
    close( rawdatafileObj )
end

function show( dataObj::RNG )
    println( "RNG Data :-----------------------------------" )
    println( " AlrmProcStartID:$(dataObj.AlrmProcStartID)")
    println( " AlrmProcEndtID1:$(dataObj.AlrmProcEndtID1)")
    println( " startID1:$(dataObj.startID1)")
    println( " endID1:$(dataObj.endID1)")
    println( " chnid:$(dataObj.chnid)")
    println( " ProtectSpaceWidth:$(dataObj.ProtectSpaceWidth)")
    println("----------------------")

end
###################
#############
function seg2Meters( CfgInfoObj::cfgStruct,segID::Int  )
    meter = CfgInfoObj.meterperSeg
    chn2Pos = CfgInfoObj.chn2SegBegIdx
    ## write the title s: 通道号下的光缆 时间、各段的米标 。
    # print( "时间" ,"," )
    # for segID = CumSumDataobj.AlrmProcStartID : CumSumDataobj.AlrmProcEndtID
    if segID  >= chn2Pos;
        return ( round( UInt,Float64(segID-chn2Pos)* meter) ,2 )
    else
        return (round( UInt,Float64( segID  )* meter ),1)
    end
end

function getRNG!( x::RNG,cfgObj::cfgStruct, startMeters::Int ,endmeters::Int ,ProtectMeters::Int,chnid::Int )
    ProtectSegs =  round(Int, Float64( ProtectMeters ) / cfgObj.meterperSeg)+1

    AlrmProcStartID =  round(Int, Float64( startMeters ) / cfgObj.meterperSeg)
    AlrmProcEndtID1 =   round(Int, Float64( endmeters ) / cfgObj.meterperSeg)

    if chnid == 1 ;
        AlrmProcStartID < 1 ? AlrmProcStartID = 1 : 0  ;
        AlrmProcEndtID1 >=  cfgObj.chn2SegBegIdx  ? AlrmProcEndtID1 = cfgObj.chn2SegBegIdx -1 : 0  ;

        startID1 = AlrmProcStartID - ProtectSegs
        endID1 = AlrmProcEndtID1 + ProtectSegs

        startID1 < 1 ? startID1 = 1 :  0 ;
        endID1 >=  cfgObj.chn2SegBegIdx  ? endID1 = cfgObj.chn2SegBegIdx : 0 ;
    else
        AlrmProcEndtID1 +=  cfgObj.chn2SegBegIdx
        AlrmProcStartID += cfgObj.chn2SegBegIdx

        AlrmProcEndtID1 >  cfgObj.SegNumber  ? AlrmProcEndtID1 = cfgObj.SegNumber :  0 ;
        AlrmProcStartID >  cfgObj.SegNumber  ? AlrmProcStartID = cfgObj.SegNumber :  0 ;

        startID1 = AlrmProcStartID - ProtectSegs
        endID1 = AlrmProcEndtID1 + ProtectSegs

        startID1 < cfgObj.chn2SegBegIdx ? startID1 = cfgObj.chn2SegBegIdx : 0 ;
        endID1 >  cfgObj.SegNumber  ? endID1 = cfgObj.SegNumber :  0 ;
        # println( AlrmProcEndtID1 ,"=======$AlrmProcStartID===$endID1 ===$startID1====")
    end

    x.ProtectSpaceWidth = ProtectSegs

    x.AlrmProcStartID = AlrmProcStartID
    x.AlrmProcEndtID1 = AlrmProcEndtID1
    x.startID1 = startID1
    x.endID1 = endID1
    x.chnid = chnid >= 2 ? chnid =2 : chnid =1;

end

function setValue!( x::cfgStruct  , y::CfgDataStruct)
    x.meterperSeg     = y.SpatialResolution
    x.chn2SegBegIdx   = y.LenCh1 +1
    x.reflectorFactor = y.Refractivity
    x.attenutionFactor=y.ATT
    x.scanrate        =y.ScanRate
    x.calSamples      = y.EnergyCnt
    x.SegNumber       = y.LenCh1 +y.LenCh2
end
function show( dataObj::CfgDataStruct   )
    println( "Config Data :-----------------------------------" )
    println( " SpatialResolution:$(dataObj.SpatialResolution)")
    println( " Refractivity:$(dataObj.Refractivity)")
    println( " OffsetM:$(dataObj.OffsetM)")
    println( " ATT:$(dataObj.ATT)")
    println( " ScanRate:$(dataObj.ScanRate)")
    println( " EnergyCnt:$(dataObj.EnergyCnt)")
    println( " LcLimitHigh:$(dataObj.LcLimitHigh)")
    println( " LcLimitLow:$(dataObj.LcLimitLow)")
    println( " StartCh1:$(dataObj.StartCh1)")
    println( " LenCh1:$(dataObj.LenCh1)")
    println( " StartCh2:$(dataObj.StartCh2)")
    println( " LenCh2:$(dataObj.LenCh2)")

    # println("----------------------")
end
function show(cfgdata::cfgStruct)
    println(" meterPerSeg = $(cfgdata.meterperSeg) ")
    println(" chn2SegBegIdx = $(cfgdata.chn2SegBegIdx )")
    println(" reflectorFactor = $(cfgdata.reflectorFactor) ")
    println(" attenutionFactor = $(cfgdata.attenutionFactor) ")
    println(" calSamples = $(cfgdata.calSamples) ")
    println(" scanrate = $(cfgdata.scanrate) ")
    println(" SegNumber = $(cfgdata.SegNumber) ")
end
function readCfgFileData(filename)
    cfgdataObj = CfgDataStruct()
    cfgfileObj= open( filename,"r")

            num  =  read( cfgfileObj, UInt16 )
            cfgdataObj.SpatialResolution =  num*1.0 / 5000.0;			# result.SpatialResolution = (double)array[0] / 5000.0;

            num  =  read( cfgfileObj, UInt16 ) ;
            cfgdataObj.Refractivity =  num*1.0 / 10000.0;	# result.Refractivity = (double)array[1] / 10000.0;

            num  =  read( cfgfileObj, UInt16 ) ;
            cfgdataObj.OffsetM =  num*1.0 ;            		# result.OffsetM = (double)array[2];

            num  =  read( cfgfileObj, UInt16 ) ;
            cfgdataObj.ATT =  num*1.0 / 10000.0; 			# result.ATT = (double)array[3] / 10000.0;

            num  =  read( cfgfileObj, UInt16 ) ;
            cfgdataObj.ScanRate =  num*1.0  ; 	          # result.ScanRate = (double)array[4];

            num  =  read( cfgfileObj, UInt16 ) ;
            cfgdataObj.EnergyCnt =  num   ;# result.EnergyCnt = (int)array[5];

            num  =  read( cfgfileObj, UInt8 ) ;
            cfgdataObj.LcLimitHigh =  num  *1.0 ;			# result.LcLimitHigh = (double)(array[6] >> 8);

			num  =  read( cfgfileObj, UInt8 ) ;
            cfgdataObj.LcLimitLow =  num   *1.0;	# result.LcLimitLow = (double)((array[6] & 255) - 256);

			num  =  read( cfgfileObj, UInt16 ) ;
            cfgdataObj.StartCh1 =  num    ;			# result.StartCh1 = (int)array[7];

			num  =  read( cfgfileObj, UInt16 ) ;
            cfgdataObj.LenCh1 =  num    ;				# result.LenCh1 = (int)array[8];

			num  =  read( cfgfileObj, UInt16 ) ;
            cfgdataObj.StartCh2 =  num			# result.StartCh2 = (int)array[9];

            num  =  read( cfgfileObj, UInt16 ) ;
            cfgdataObj.LenCh2 =  num          # result.LenCh2 = (int)array[10];

    close( cfgfileObj)
    return cfgdataObj

end
function read!(x::cfgStruct,filename)
    cfgdataObj  = readCfgFileData(filename)
    setValue!(x,cfgdataObj )

end

function TestRead(filename ="/media/wang/SSDP-128G/code/20190212T131313.conf")
    filename ="/media/wang/SSDP-128G/code/20190212T131313.conf"
    cfgdataObj  = readCfgFileData(filename)
    show(cfgdataObj )
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
using sem, Arrow, ModelingToolkit, LinearAlgebra, SparseArrays, DataFrames, Optim, LineSearches

cd("test")

## Observed Data
dat = DataFrame(Arrow.Table("comparisons/data_growth.arrow"))
par_ml = DataFrame(Arrow.Table("comparisons/par_growth_ml.arrow"))
par_ls = DataFrame(Arrow.Table("comparisons/par_growth_ls.arrow"))

# observed
semobserved = SemObsCommon(data = Matrix{Float64}(dat))

#diff
diff_fin = SemFiniteDiff(
    BFGS(),
    Optim.Options())

## Model definition
@ModelingToolkit.variables x[1:22]

#x = rand(31)
S = zeros(Num,12,12)
S[diagind(S)] .= x[1:12]
S[7,8] = x[13]; S[8,7] = x[13]
S[10,12] = x[15]; S[12,10] = x[15]
S[11,12] = x[16]; S[12,11] = x[16]
S[10,11] = x[14]; S[11,10] = x[14]

F = zeros(9, 12)
F[diagind(F)] .= 1.0

A =[0  0  0  0  0  0  0  0  0     1     0     0
    0  0  0  0  0  0  0  0  0     x[17] 0     0
    0  0  0  0  0  0  0  0  0     x[18] 0     0
    0  0  0  0  0  0  0  0  0     0     1     0
    0  0  0  0  0  0  0  0  0     0     x[19] 0
    0  0  0  0  0  0  0  0  0     0     x[20] 0
    0  0  0  0  0  0  0  0  0     0     0     1
    0  0  0  0  0  0  0  0  0     0     0     x[21]
    0  0  0  0  0  0  0  0  0     0     0     x[22]
    0  0  0  0  0  0  0  0  0     0     0     0
    0  0  0  0  0  0  0  0  0     0     0     0
    0  0  0  0  0  0  0  0  0     0     0     0]

S = sparse(S)

#F
F = sparse(F)

#A
A = sparse(A)

par_order = [collect(11:22); collect(10); collect(23:25); 2;3;5;6;8;9]

start_val_ml = Vector{Float64}(par_ml.start[par_order])
start_val_ls = Vector{Float64}(par_ls.start[par_order])
start_val_snlls = Vector{Float64}(par_ls.start[par_order][17:22])

# loss
loss_ml = Loss([SemML(semobserved, [0.0], similar(start_val_ml))])
loss_ls = Loss([sem.SemWLS(semobserved)])
loss_snlls = Loss([sem.SemSWLS(semobserved)])

#start_val_ml = [fill(1.0, 14); fill(0, 6); fill(1, 8); fill(0, 3)]
#start_val_ls = [fill(1.0, 14); fill(0, 6); fill(1, 8); fill(0, 3)]
#start_val_snlls = [fill(1.0, 8); fill(0, 3)]

# imply
imply_ml = ImplySymbolic(A, S, F, x, start_val_ml)
imply_ls = sem.ImplySymbolicWLS(A, S, F, x, start_val_ls)
imply_snlls = sem.ImplySymbolicSWLS(A, S, F, x[17:22], start_val_snlls)

# model
model_ml = Sem(semobserved, imply_ml, loss_ml, diff_fin)
model_ls = Sem(semobserved, imply_ls, loss_ls, diff_fin)
model_snlls = Sem(semobserved, imply_snlls, loss_snlls, diff_fin)

#checks
imply_ml(start_val_ml, model_ml)
imply_ls(start_val_ls, model_ls)

ind = CartesianIndices(imply_ml.imp_cov)
ind = filter(x -> (x[1] >= x[2]), ind)
s = imply_ml.imp_cov[ind]

imply_ls.imp_cov ≈ s

imply_ml(x, model_ml)
imp_cov = F*inv(I-A)*S*transpose(inv(I-A))*transpose(F)
imp_cov ≈ imply_ml.imp_cov

imply_snlls(start_val_snlls, model_snlls)
t = imply_snlls.G*start_val_ls[[1, 2, 3, 4, 15, 5, 16, 17, 6, 18, 7, 19, 8, 9, 20, 10, 11, 12, 13, 14]]
t ≈ imply_ls.imp_cov

V = I
s = cov(Matrix(dat))[tril(trues(size(cov(Matrix(dat)))))]
G = imply_snlls.G

transpose(s)*V*s - transpose(s)*V*G*inv(transpose(G)*V*G)*(transpose(G)*V*s)

# fit
solution_ml = sem_fit(model_ml)
solution_ls = sem_fit(model_ls)
solution_snlls = sem_fit(model_snlls)

all(#
        abs.(solution_ml.minimizer .- par_ml.est[par_order]
            ) .< 0.05*abs.(par_ml.est[par_order]))

all(#
            abs.(solution_ls.minimizer .- par_ls.est[par_order]
                ) .< 0.05*abs.(par_ls.est[par_order]))

all(#
            abs.(solution_snlls.minimizer .- par_ls.est[par_order][17:22]
                ) .< 0.05*abs.(par_ls.est[par_order][17:22]))

using BenchmarkTools

@benchmark solution_ml = sem_fit(model_ml)
@benchmark solution_ls = sem_fit(model_ls)
@benchmark solution_snlls = sem_fit(model_snlls)

@benchmark model_ml(start_val_ml)
@benchmark model_ls(start_val_ls)
@benchmark model_snlls(start_val_snlls)

a = rand(100,100)
a = a*a'
b = rand(100,1)
b'*a*b

function inv_dot(a, b)
    u = cholesky(a)
    z = u.L\b
    w = u.U\z
    sol = b'*w
    return sol
end

function inv_dot(a, b)
    u = cholesky(a)
    z = u.L\b
    w = u.U\z
    sol = b'*w
    return sol
end

function inv_dot_naive(a, b)
    sol = b'*inv(a)*b
    return sol
end

function inv_dot_dot(a, b)
    sol = dot(b', inv(a), b)
    return sol
end

function inv_solve(a, b)
    right = a\b
    sol = b'*right
    return sol
end

function inv_chol(a, b)
    c = cholesky(a)
    right = c\b
    sol = b'*right
    return sol
end

inv_dot(a, b) ≈ inv_dot_naive(a, b)
inv_dot_dot(a, b) ≈ inv_dot_naive(a, b)[1]
inv_solve(a, b) ≈ inv_dot_naive(a, b)
inv_chol(a, b) ≈ inv_dot_naive(a, b)

@benchmark inv_dot($a, $b)
@benchmark inv_dot_naive($a, $b)
@benchmark inv_dot_dot($a, $b)
@benchmark inv_solve($a, $b)
@benchmark inv_chol($a, $b)


text:
# compute R²
function rsq(Y, X)
    β = X\Y
    errors = Y - X*β
    R² = 1.0 - sum(errors.^2.0)/(size(X,1)-size(X,2))
end

max_stars_repo_name: manuelmontesinos/Econometrics

max_issues_repo_name: manuelmontesinos/Econometrics

max_forks_repo_name: manuelmontesinos/Econometrics


text:
module JPEGs

include("color_space_transformation.jl")
include("downsampling.jl")
include("block_splitting.jl")
include("dct.jl")
include("quantization.jl")
include("entropycoding.jl")
include("segments.jl")

function encode_image(img::Array{ColorTypes.RGB{FixedPointNumbers.Normed{UInt8,8}},2})
    i2 = JPEGs.rgb_to_ycbcr(img)
    y,cb,cr = JPEGs.split_channels(i2)
    #dcb = JPEGs.downsample_4_2_0(cb)
    #dcr = JPEGs.downsample_4_2_0(cr)
    y = JPEGs.split_into_blocks(y)
    cb = JPEGs.split_into_blocks(cb)
    cr = JPEGs.split_into_blocks(cr)
    y = JPEGs.dct.(y)
    cb = JPEGs.dct.(cb)
    cr = JPEGs.dct.(cr)

    y = JPEGs.quantize.(y)
    cb = JPEGs.quantize.(cb)
    cr = JPEGs.quantize.(cr)
    y = JPEGs.zigzag.(y)
    cb = JPEGs.zigzag.(cb)
    cr = JPEGs.zigzag.(cr)
    encoded_y = JPEGs.encode_zigzagged_block.(y)
    encoded_cb = JPEGs.encode_zigzagged_block.(cb)
    encoded_cr = JPEGs.encode_zigzagged_block.(cr)
    interlaced = [vcat(i...) for i in zip(encoded_y, encoded_cb, encoded_cr)]
    flattened = reshape(interlaced, *(size(interlaced)...))
    encoded_blocks = vcat(flattened...)
    serialized = JPEGs.serialize_jpeg_scan(encoded_blocks)
    return serialized
end

function decode_image(serialized_jpeg::Vector{UInt8})#::Array{ColorTypes.RGB{FixedPointNumbers.Normed{UInt8,8}},2}
    allbits = JPEGs.deserialize_into_jpeg_scan(serialized_jpeg)
    decoded_blocks = []
    while length(decoded_blocks) < (64*64*3)-1
        if length(decoded_blocks) == 12287
            println(allbits)
        end
        push!(decoded_blocks, decode_zigzagged_block(allbits))
    end
    y = []
    cb = []
    cr = []
    for b in 1:3:length(decoded_blocks)
        push!(y, JPEGs.unzigzag(decoded_blocks[b]))
        push!(cb, JPEGs.unzigzag(decoded_blocks[b+1]))
        push!(cr, JPEGs.unzigzag(decoded_blocks[b+2]))
    end
    y = JPEGs.unquantize.(y)
    cb = JPEGs.unquantize.(cb)
    cr = JPEGs.unquantize.(cr)
    y = JPEGs.invdct.(y)
    cb = JPEGs.invdct.(cb)
    cr = JPEGs.invdct.(cr)
    y = reshape(y, (64,64))
    cb = reshape(cb, (64,64))
    cr = reshape(cr, (64,64))
    y = JPEGs.join_blocks(y)
    cb = JPEGs.join_blocks(cb)
    cr = JPEGs.join_blocks(cr)
    i2 = JPEGs.join_channels((y,cb,cr))
    img = JPEGs.ycbcr_to_rgb(i2)
    return img
end

end


text:
using Quiz

PROB_ID = qProb([
    qT("A rigid tank contains " )
    qV(:m1, 1.5:0.5:4.5, u"kg")
    qT(" of Oxygen gas (M = ")
    qC(:MO2, 31.999, u"kg/kmol")
    qT(") at ")
    qV(:T1, 300.0:25.0:450.0, u"K")
    qT(" and ")
    qV(:P1, 100.0:50.0:500.0, u"kPa")
    qT(". Assuming ideal gas behavior for the Oxygen gas, determine:")
    qA(3, "The specific volume of the gas is _.", :v1, u"m^3/kg")
    qA(1, "The volume of the tank is _.", :V1, u"m^3")
])

# Sampling
SPL = sample(PROB_ID, 20; mode=:max | :all)

# Solving
SOL = solve(SPL)

# Backend output
OUT = apply(PROB_ID, MOO, SOL)


text:
"""
    PyObjectArray(undef, dims...)
    PyObjectArray(array)

An array of `Py`s which supports the Python buffer protocol.

Internally, the objects are stored as an array of pointers.
"""
mutable struct PyObjectArray{N} <: AbstractArray{Py,N}
    ptrs::Array{C.PyPtr,N}
    function PyObjectArray{N}(::UndefInitializer, dims::NTuple{N,Integer}) where {N}
        x = new{N}(fill(C.PyNULL, dims))
        finalizer(pyobjectarray_finalizer, x)
    end
end
const PyObjectVector = PyObjectArray{1}
const PyObjectMatrix = PyObjectArray{2}
export PyObjectVector, PyObjectMatrix, PyObjectArray

PyObjectArray{N}(undef::UndefInitializer, dims::Vararg{Integer,N}) where {N} = PyObjectArray(undef, dims)
PyObjectArray(undef::UndefInitializer, dims::NTuple{N,Integer}) where {N} = PyObjectArray{N}(undef, dims)
PyObjectArray(undef::UndefInitializer, dims::Vararg{Integer,N}) where {N} = PyObjectArray{N}(undef, dims)
PyObjectArray{N}(x::AbstractArray{T,N}) where {T,N} = copyto!(PyObjectArray{N}(undef, size(x)), x)
PyObjectArray(x::AbstractArray{T,N}) where {T,N} = PyObjectArray{N}(x)

pyobjectarray_finalizer(x::PyObjectArray) = begin
    if C.CTX.is_initialized
        C.with_gil(false) do
            for ptr in x.ptrs
                C.Py_DecRef(ptr)
            end
        end
    end
end

Base.IndexStyle(x) = Base.IndexStyle(x.ptrs)

Base.length(x::PyObjectArray) = length(x.ptrs)

Base.size(x::PyObjectArray) = size(x.ptrs)

@propagate_inbounds function Base.isassigned(x::PyObjectArray, i::Vararg{Integer})
    @boundscheck checkbounds(Bool, x, i...) || return false
    return @inbounds x.ptrs[i...] != C_NULL
end

@propagate_inbounds function Base.getindex(x::PyObjectArray, i::Integer...)
    @boundscheck checkbounds(x, i...)
    @inbounds ptr = x.ptrs[i...]
    ptr == C_NULL && throw(UndefRefError())
    return pynew(incref(ptr))
end

@propagate_inbounds function Base.setindex!(x::PyObjectArray, v, i::Integer...)
    @boundscheck checkbounds(x, i...)
    v_ = Py(v)
    @inbounds decref(x.ptrs[i...])
    @inbounds x.ptrs[i...] = getptr(v_)
    pystolen!(v_)
    return x
end

@propagate_inbounds function Base.deleteat!(x::PyObjectVector, i::Integer)
    @boundscheck checkbounds(x, i)
    @inbounds decref(x.ptrs[i])
    deleteat!(x.ptrs, i)
    return x
end

function pyjlarray_array_interface(x::PyObjectArray)
    return pydict(
        shape = size(x),
        typestr = "|O",
        data = (UInt(pointer(x.ptrs)), false),
        strides = strides(x.ptrs) .* sizeof(C.PyPtr),
        version = 3,
    )
end

# C._pyjlarray_get_buffer(o, buf, flags, x::PyObjectArray) = C.pyjl_get_buffer_impl(
#     o,
#     buf,
#     flags,
#     pointer(x.ptrs),
#     sizeof(CPyPtr),
#     length(x),
#     ndims(x),
#     "O",
#     size(x),
#     strides(x.ptrs) .* sizeof(CPyPtr),
#     true,
# )


text:
module TestEnsemblePTSampler

using Distributed

using Test: @test, @testset

using Ensemble

using Statistics

function test_brewer()
    nd = 5
    nw = 128
    nt = 16

    function loglike(x)
        v = 0.1
        u = 0.01
        m = 0.031

        logl1 = -nd*log(v) - nd/2*log(2*pi) - 0.5*sum(x.*x/(v*v))
        logl2 = -nd*log(u) - nd/2*log(2*pi) - 0.5*sum((x.-m).*(x.-m)/(u*u))

        return Stats.logsumexp(logl1, log(99.0) + logl2)
    end

    function logprior(x)
        if any(x .> 0.5) || any(x .< -0.5)
            return -Inf
        else
            return 0.0
        end
    end

    pts = rand(nd,nw,nt) .- 0.5

    betas = collect(range(1, stop=0, length=nt))

    chain, lnlikes, lnpriors, betas = EnsemblePTSampler.run_mcmc(pts, loglike, logprior, betas, 8192, thin=4)

    @testset "Brewer phase transition distribution" begin
        for i in 1:nd
            @test isapprox(mean(chain[i,:,1,:]), 0.99*0.031, atol=0.1*sqrt(0.99*0.01^2 + 0.01*0.1^2))
        end

        lnZ, dlnZ = EnsemblePTSampler.lnZ(lnlikes, betas)
        @test isapprox(lnZ, log(100.0), atol=3*dlnZ)
    end

    chain, lnlikes, lnpriors, betas
end

function testall()
    @testset "EnsemblePTSampler test suite" begin
        test_brewer()
    end
end

end


text:
## Custom serialization

# The default, write a type as itself
writeas(T::Type) = T

# wconvert and rconvert do type conversion before reading and writing,
# respectively. These fall back to convert.
wconvert(T, x) = convert(T, x)
rconvert(T, x) = convert(T, x)

# Select an ODR, incorporating custom serialization only if the types do not
# match
CustomSerialization(::Type{WrittenAs}, ::Type{WrittenAs}, odr) where {WrittenAs} = odr
CustomSerialization(::Type{WrittenAs}, ::Type{ReadAs}, odr) where {WrittenAs,ReadAs} =
    CustomSerialization{WrittenAs,odr}

odr_sizeof(::Type{CustomSerialization{T,ODR}}) where {T,ODR} = odr_sizeof(ODR)

# Usually we want to convert the object and then write it.
@inline h5convert!(out::Pointers, ::Type{CustomSerialization{T,ODR}}, f::JLDFile,
                   x, wsession::JLDWriteSession) where {T,ODR} =
    h5convert!(out, ODR, f, wconvert(T, x)::T, wsession)

# When writing as a reference, we don't want to convert the object first. That
# should happen automatically after write_dataset is called so that the written
# object gets the right written_type attribute.
@inline h5convert!(out::Pointers, odr::Type{CustomSerialization{T,RelOffset}},
                   f::JLDFile, x, wsession::JLDWriteSession) where {T} =
    h5convert!(out, RelOffset, f, x, wsession)

# When writing as a reference to something that's being custom-serialized as an
# array, we have to convert the object first.
@inline h5convert!(out::Pointers, odr::Type{CustomSerialization{T,RelOffset}},
            f::JLDFile, x, wsession::JLDWriteSession) where {T<:Array} =
    h5convert!(out, RelOffset, f, wconvert(T, x)::T, wsession)

h5convert_uninitialized!(out::Pointers, odr::Type{CustomSerialization{T,ODR}}) where {T,ODR} =
    h5convert_uninitialized!(out, ODR)

jlconvert_canbeuninitialized(::ReadRepresentation{T,CustomSerialization{S,ODR}}) where {T,S,ODR} =
    jlconvert_canbeuninitialized(ODR)
jlconvert_isinitialized(::ReadRepresentation{T,CustomSerialization{S,ODR}}, ptr::Ptr) where {T,S,ODR} =
    jlconvert_isinitialized(ReadRepresentation{S,ODR}(), ptr)
jlconvert(::ReadRepresentation{T,CustomSerialization{S,ODR}},
          f::JLDFile, ptr::Ptr, header_offset::RelOffset) where {T,S,ODR} =
    rconvert(T, jlconvert(ReadRepresentation{S,ODR}(), f, ptr, header_offset))::T


text:
module TempleBigEyeball

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Temple Big Eyeball" => Ahorn.EntityPlacement(
        Maple.TempleBigEyeball
    )
)

bodySprite = "danger/templeeye/body00.png"
pupilSprite = "danger/templeeye/pupil.png"

function Ahorn.selection(entity::Maple.TempleBigEyeball)
    x, y = Ahorn.position(entity)

    return Ahorn.coverRectangles([
        Ahorn.getSpriteRectangle(bodySprite, x, y),
        Ahorn.getSpriteRectangle(pupilSprite, x, y)
    ])
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.TempleBigEyeball, room::Maple.Room)
    Ahorn.drawSprite(ctx, bodySprite, 0, 0)
    Ahorn.drawSprite(ctx, pupilSprite, 0, 0)
end

end


text:
"""
SPUNG: State Research Program for Utilization of Natural Gas

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

The idea is to use a "shape model" that provides a corresponding states parameters
and a "reference model" that implements a helmholtz energy function.

this SPUNG by default uses the propane reference equation of state (`PropaneRef`) as the reference model
and SRK for the shape model.
"""
struct SPUNG{S<:EoSModel,E<:EoSModel} <: EoSModel
    shape_model::S
    shape_ref::S
    model_ref::E
end

function SPUNG(components::Vector{String},refmodel=PropaneRef(),shapemodel=SRK(components),shaperef = SRK(component_names(refmodel)))
    model = SPUNG(shapemodel,shaperef,refmodel)
    return model
end

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

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

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

mw(model::SPUNG) = mw(model.shape_model)
molecular_weight(model::SPUNG,z=SA[1.0]) = molecular_weight(model.shape_model,z)

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

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

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

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

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

#=
ideally we could perform SPUNG only providing x0, but i cant find the error here

=#
#overloading sat_pure for SPUNG directly seems to be the way
function sat_pure(model::SPUNG,T::Real,v0=[zero(T)/zero(T),zero(T)/zero(T)])
    lb_v0 = lb_volume(model.model_ref)
    f,h = shape_factors(model,lb_v0,T,SA[1.0]) #h normaly should be independent of temperature
    T0 = T/f
    if isnan(v0[1]) && isnan(v0[2])
        v0 = x0_sat_pure(model.model_ref,T0)
    else
        vl = exp10(v0[1])
        vv = exp10(v0[2])
        v0 = [log10(vl*h),log10(vv*h)]
    end
    psat0,vl0,vv0 = sat_pure(model.model_ref,T0,v0)
    p = pressure(model,vl0*h,T)
    return (p,vl0*h,vv0*h)
end

#============
uncomment when sat_pure_p is ready
====================
function sat_pure_p(model::SPUNG,p::Real)
    return naive_sat_pure_p(model,p)
end
=#

# function split_model(model::SPUNG)
#     #only the shape model is splittable
#     shape_model = split_model(model.shape_model)
#     len = length(shape_model)
#     shape_ref = fill(model.shape_ref,len)
#     model_ref = fill(model.model_ref,len)
#     return SPUNG.(shape_model,shape_ref,model_ref)
# end

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

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

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

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

=#

export SPUNG


text:
"""
Function creates writeable object and creates & writes -∞ Array to it. Array Dimmensions are not stored.
# Arguments
- `dimsize::NTuple{N,:<Int}`: Dimmensions of an array to write
- `dir_name::String = "/tmp/io.mmap": Path to write a file.
"""
function create_mmap_array(dimsize::NTuple{N,<:Integer}, dir_name::String = "/tmp/io.mmap") where N

    s = Mmap.open(dir_name, "w+")
    Mmap.write(s, repeat([-Inf::Float64], dimsize... ) )
    close(s)

end

"""
Function opens with r+ permission a readable file with an array and returns IO Array (Mmap)
# Arguments
- `dimsize::NTuple{N,:<Int}`: Ntuple of dimensions of Array
- `dir_name::String = "/tmp/io.mmap": Path to a file.
"""
function read2write_mmap_array(dimsize::NTuple{N,<:Integer}, dir_path::String = "/tmp/io.mmap") where N

    s = Mmap.open(dir_path, "r+")
    A = Mmap.mmap(s, Array{Float64,length(dimsize)}, dimsize)
    close(s)
    return(A)

end

"""
Function opens with r permission a readable file with an array and returns IO Array (Mmap)
# Arguments
- `dimsize::NTuple{N,:<Int}`: Ntuple of dimensions of Array
- `dir_name::String = "/tmp/io.mmap": Path to a file.
"""
function read_mmap_array(dimsize::NTuple{N,<:Integer}, dir_path::String = "/tmp/io.mmap") where N

    s = Mmap.open(dir_path, "r")
    A = Mmap.mmap(s, Array{Float64,length(dimsize)}, dimsize)
    close(s)
    return(A)

end

"""
Function that returns Utility for given prices and coordinates. Requires precompilation of Utility Cube. MMap (hard drive) version. Requires precompilation of utils cube on disk. Accesses values directly. Requires few GB of space.
# Arguments
- `ind_vec::NTuple{M,<:Integer}`: coordinates of Utility.
- `mdl::Model`: mutable struct Model."""

function cube_Mmap( ind_vec::Vector{<:Integer}, mdl::Model, additional_args=nothing)

    grid = mdl.GridDef
    index, status = ndim2lin(ind_vec, mdl.GridDef)
    return(mdl.Utils.qb[ status ][ index ])

end

"""
Function building cube. GPU version. Requires precompilation of utils cube in VRAM. Accesses values directly.Requires as many GPUs as there are statuses.
# Arguments
- `mdl::ModelInit`: struct ModelInit.
# Examples```julia-repl\njulia> execute_cube(my_model, "CPU")\n```
"""
function cube_Mmap_qb(mdl::Model)

    n_dim = tuple( length.(mdl.GridDef.Grids)[1:(end-1)]... )
    cube = Vector{Array{Float64}}()

    for status = mdl.GridDef.Grids[end]

        path_ = joinpath(pwd(), "/tmp/utils_" * string(status) * ".io")
        create_mmap_array(n_dim, path_)
        append!( cube, [ read2write_mmap_array( n_dim, path_ ) ] )

    end

    Threads.@threads for idx = CartesianIndices(n_dim)

        for status = mdl.GridDef.Grids[end]

            ind_vec = CartesianIndex( idx, status )
            cube[status][idx] =  mdl.HAutility( mdl.HAconstraint( Tuple(ind_vec) , mdl )... )

        end

    end

    cube_out = Vector{Array{Float64}}()

    for status = mdl.GridDef.Grids[end]

        Mmap.sync!(cube[status])
        path_ = "/tmp/utils_" * string(status) * ".io"

        append!( cube_out, [read_mmap_array(n_dim, path_)])

    end

    return(cube_out)

end


text:
export Abstract_Solver_Result
export converged, iteration_count, objective_value, solution

# Result abstraction for unconstrained problem ================
#
@doc raw"""
```julia
abstract type Abstract_Solver_Result end
```

This is the base type returned by the [`solve`](@ref) method. It
contains the information related to the founded solution.

# Interface
- [`converged`](@ref)
- [`iteration_count`](@ref)
- [`objective_value`](@ref)
- [`solution`](@ref)

# Implementations
- [`LevenbergMarquardt_Result`](@ref)
- [`LevenbergMarquardt_BC_Result`](@ref)

"""
abstract type Abstract_Solver_Result end

"""
    converged(::Abstract_Solver_Result)

Return `true` if the solver converged

See: [`Abstract_Solver_Result`](@ref)

"""
converged(::Abstract_Solver_Result) = @assert(false,"To implement")

"""
    iteration_count(::Abstract_Solver_Result)

Return the number of consumed iteration

See: [`Abstract_Solver_Result`](@ref)

"""
iteration_count(::Abstract_Solver_Result) = @assert(false,"To implement")

"""
    objective_value(::Abstract_Solver_Result)

Returns objective value at the point [`solution`](@ref).

See: [`Abstract_Solver_Result`](@ref)

"""
objective_value(r::Abstract_Solver_Result) =  @assert(false,"To implement")

"""
    solution(::Abstract_Solver_Result)

Returns the founded solution

See: [`Abstract_Solver_Result`](@ref)

"""
solution(r::Abstract_Solver_Result) =  @assert(false,"To implement")


text:
Nfile = 8

myFiles = 1:Nfile

Aout = EmptyAssoc()
Ain = EmptyAssoc()

for i = myFiles
    tic()
    fname = "data/"* string(i)

    fidRow = open(fname * "r.txt","r")
    fidCol = open(fname * "c.txt","r")
    fidVal = open(fname * "v.txt","r")

    rowStr = readline(fidRow)
    colStr = readline(fidCol)
    valStr = readline(fidVal)

    close(fidRow)
    close(fidCol)
    close(fidVal)

    A = Assoc(rowStr,colStr,1,(sum))
    Aout += sum(A,2)
    Ain  += sum(A,1)

#TODO update Time
    pTime = toq()
    println("Sum Time: ", pTime)
    println("Edges/sec: ", string(NumStr(rowStr)/pTime))
end

#= #PyPlot
figure()
loglog(full(sparse(Adj(Aout),1,1)),"o")
xlabel("out degree")

figure()
loglog(full(sparse(Adj(Ain),1,1)),"o")
xlabel("in degree")

=#


text:
include("../common.jl")
build_julia(v"1.5.1")


text:
const tolerance = 1e-16

const QSpins = Union{AbstractVector{Int64}, AbstractVector{NTuple{2, Int64}}, CartesianRange}

function build_quantum_links(L::Lattice, quantum_spins::AbstractVector{NTuple{2, Int64}})
    q_links = Vector{Tuple{Int64,Int64,Float64}}()
    num_of_spins = length(quantum_spins)
    Dh = 2^num_of_spins
    for i in 1:num_of_spins-1
        for j in i+1:num_of_spins
            val = get_value(L, quantum_spins[i], quantum_spins[j])*0.25
            if val!=0.
                push!(q_links, (i, j, val))
            end
        end
    end
    return q_links
end

#function to calculate number of nnz matrix elements
#arising from $S_i^\sigma S_j^\sigma$ for $\sigma = x,y0$
#if J_x == J_y
function calculate_nnz_xx_yy(Dh::Int64, q_links::Vector{Tuple{Int64,Int64,Float64}})
    nnz_xx_yy = zero(Dh)
    for vec in 1:Dh
        for link in q_links
            nnz_xx_yy += div(1-getSign(vec, link[1], link[2]),2)
        end
    end
    return nnz_xx_yy
end

function build_Hamiltonian(L::Lattice, quantum_spins::AbstractVector{NTuple{2, Int64}}, T::Type)
    q_links = build_quantum_links(L, quantum_spins)
    num_of_links = length(q_links)
    num_of_spins = length(quantum_spins)
    Dh = 2^num_of_spins
    #calculate number of nonzero elements to store in the matrix
    assert((L.Js |> collect |> norm) + (L.hs |> collect |> norm) > 0)
    num_of_nnzs = zero(Dh)
    if ([L.Js[3], L.hs[3]] |> norm) > 0
        num_of_nnzs += Dh
    end
    if L.Js[1] == L.Js[2] == 0
        #do nothing
    elseif L.Js[1] == L.Js[2]
        num_of_nnzs += calculate_nnz_xx_yy(Dh, q_links)
    elseif L.Js[1] == -L.Js[2]
        num_of_nnzs += Dh*num_of_links - calculate_nnz_xx_yy(Dh, q_links)
    else
        num_of_nnzs += Dh*num_of_links
    end
    num_of_nnzs += Dh*num_of_spins*((L.hs[1]!=0) + (L.hs[2]!=0))
    #define arrays to store data
    data = SharedVector{T}(num_of_nnzs)
    rowind = SharedVector{Int64}(num_of_nnzs)
    colptr = SharedVector{Int64}(Dh + 1)
    #main cycle
    index = 1
    for vec in 1:Dh
        colptr[vec] = index
        data_temp = 0.
        for pos in 1:num_of_spins
            val, vec2 = pauliZ(vec, pos)
            data_temp -= 0.5*L.hs[3]*val
        end
        for link in q_links
            sign = getSign(vec, link[1], link[2])
            data_temp -= link[3]*L.Js[3]*sign
        end
        if abs(data_temp)>tolerance
            rowind[index] = vec
            data[index] = data_temp
            index+=1
        end
        for link in q_links
            vec2 = invert(vec, link[1], link[2])
            sign = getSign(vec, link[1], link[2])
            rowind_temp, data_temp = vec2, (L.Js[2]*sign-L.Js[1])*link[3]
            if abs(data_temp)>tolerance
                rowind[index] = rowind_temp
                data[index] = data_temp
                index += 1
            end
        end
        if L.hs[1]!=0.
            for pos in 1:num_of_spins
                val, rowind[index] = pauliX(vec, pos)
                data[index] = -0.5*val*L.hs[1]
                index += 1
            end
        end
        if L.hs[2]!=0.
            for pos in 1:num_of_spins
                val, rowind[index] = pauliY(vec, pos)
                data[index] = vec, -0.5*val*L.hs[1]
                index += 1
            end
        end
        colind = @view rowind[colptr[vec]:index-1]
        coldata = @view data[colptr[vec]:index-1]
        perm = sortperm(colind)
        permute!(colind, perm)
        permute!(coldata, perm)
    end
    colptr[end] = index
    #if(index<num_of_nnzs+1)
    #    resize!(rowind, index-1)
    #    resize!(data, index-1)
    #end
    return SharedSparseMatrixCSC(Dh, Dh, colptr, rowind, data)
end
build_Hamiltonian(L::Lattice{D,F,T,true}, quantum_spins::AbstractVector{NTuple{2, Int64}}) where {D,F,T} = build_Hamiltonian(L, quantum_spins, Complex{Float64})
build_Hamiltonian(L::Lattice{D,F,T,false}, quantum_spins::AbstractVector{NTuple{2, Int64}}) where {D,F,T} = build_Hamiltonian(L, quantum_spins, Float64)
build_Hamiltonian(A::Approx) where Approx<:AbstractQuantumApproximation = build_Hamiltonian(A.L, get_q_spins(A))

function build_Spin_Operator(Dh::Int64, spin::Int64, sigma::Int64)
    assert(1<=sigma<=3)
    num_of_nnzs = Dh
    #define arrays to store data
    data = SharedVector{Complex{Float64}}(num_of_nnzs)
    rowind = SharedVector{Int64}(num_of_nnzs)
    colptr = SharedVector{Int64}(Dh + 1)
    #main cycle
    index = 1
    for vec in 1:Dh
        colptr[vec] = index
        data[index], rowind[index] = pauli[sigma](vec, spin)
        index += 1
    end
    colptr[end] = index
    return SharedSparseMatrixCSC(Dh, Dh, colptr, rowind, data)
end
function build_Spin_Operator(Dh::Int64, spins::Vector{Int64}, sigma::Int64)
    if sigma==3
        num_of_nnzs = Dh
        #define arrays to store data
        data = SharedVector{Complex{Float64}}(num_of_nnzs)
        rowind = SharedVector{Int64}(num_of_nnzs)
        colptr = SharedVector{Int64}(Dh + 1)
        #main cycle
        index = 1
        for vec in 1:Dh
            colptr[vec] = index
            data_temp = 0.
            for pos in spins
                val, vec2 = pauliZ(vec,pos)
                data_temp += val
            end
            rowind[index] = vec
            data[index] = data_temp
            index+=1
        end
        colptr[end] = index
        return SharedSparseMatrixCSC(Dh, Dh, colptr, rowind, data)
    else
        num_of_nnzs = Dh*length(spins)
        #define arrays to store data
        data = SharedVector{Complex{Float64}}(num_of_nnzs)
        rowind = SharedVector{Int64}(num_of_nnzs)
        colptr = SharedVector{Int64}(Dh + 1)
        #main cycle
        index = 1
        for vec in 1:Dh
            colptr[vec] = index
            for pos in spins
                data[index], rowind[index] = pauli[sigma](vec, pos)
                index += 1
            end
            colind = @view rowind[colptr[vec]:index-1]
            coldata = @view data[colptr[vec]:index-1]
            perm = sortperm(colind)
            permute!(colind, perm)
            permute!(coldata, perm)
        end
        colptr[end] = index
        return SharedSparseMatrixCSC(Dh, Dh, colptr, rowind, data)
    end
end

function build_ClassicalInteractions(L::Lattice{D,typeof(nearest_neighbours),T,C}, cl_spins::Vector{Tuple{Int64,Int64}}) where {D,T,C}
    cl_num = length(cl_spins)
    IM = zeros((cl_num, cl_num))
    for i in 1:cl_num-1 for j in i+1:cl_num
        IM[i,j] = get_value(L, cl_spins[i], cl_spins[j])
        IM[j,i] = IM[i,j]
    end end
    return shared_sparse(IM)::SharedSparseMatrixCSC{Float64, Int64}
end
function build_ClassicalInteractions(L::Lattice{D,MDI{D2},T,C}, cl_spins::Vector{Tuple{Int64,Int64}}) where {D,D2,T,C}
    cl_num = length(cl_spins)
    IM = SharedMatrix{Float64}((cl_num, cl_num))
    fill!(IM,0.0)
    for i in 1:cl_num-1 for j in i+1:cl_num
        IM[i,j] = get_value(L, cl_spins[i], cl_spins[j])
        IM[j,i] = IM[i,j]
    end end
    print(issymmetric(IM),"\n")
    return Symmetric(IM)::Symmetric{Float64, SharedMatrix{Float64}}
end
function build_ClassicalInteractions(L::Lattice{D,MDI{D2},T,C}, cl_spins::Vector{Tuple{Int64,Int64}}, factor::Float64) where {D,D2,T,C}
    cl_num = length(cl_spins)
    IM = SharedMatrix{Float64}((cl_num, cl_num))
    fill!(IM,0.0)
    for i in 1:cl_num-1 for j in i+1:cl_num
        IM[i,j] = get_value(L, cl_spins[i], cl_spins[j])*factor
        IM[j,i] = IM[i,j]
    end end
    print(issymmetric(IM),"\n")
    return Symmetric(IM)::Symmetric{Float64, SharedMatrix{Float64}}
end
build_ClassicalInteractions(A::Approx) where Approx<:Union{PureClassicalApprox, HybridApprox} = build_ClassicalInteractions(A.L, get_cl_spins(A))
build_ClassicalInteractions(A::PureClassicalApprox, factor::Float64) = build_ClassicalInteractions(A.L, get_cl_spins(A), factor)

@inline build_InterClusterInteractions(A::ClusteredApprox) = _build_InterClusterInteractions(A)
function build_InterClusterInteractions(A::ClusteredApprox{Lattice{D,typeof(nearest_neighbours),T,C}}) where {D,T,C}
    edge_spins, IM = _build_InterClusterInteractions(A)
    return edge_spins::Vector{Int64}, shared_sparse(IM)::SharedSparseMatrixCSC{Float64,Int64}
end

function _build_InterClusterInteractions(A::ClusteredApprox)
    edge_spins = Set{Int64}()
    cluster_spins = get_first_cluster(A)
    links = Vector{Tuple{Int64,Int64,Int64,Float64}}()
    cluster_pos1 = ind2sub(A.outer_cluster_dims, 1)
    cl_pos1 = collect(cluster_pos1)
    for sp1 in eachindex(cluster_spins)
        index1 = (cluster2lin(A, cluster_spins[sp1][1],cluster_pos1), cluster_spins[sp1][2])
        for cluster in 2:A.cluster_num
            cluster_pos2 = ind2sub(A.outer_cluster_dims, cluster)
            for sp2 in eachindex(cluster_spins)
                index2 = (cluster2lin(A, cluster_spins[sp2][1], cluster_pos2), cluster_spins[sp2][2])
                val = get_value(A, index1, index2)
                if val!=0.0
                    push!(edge_spins, sp1)
                    push!(links, (sp1, sp2, cluster, val))
                end
            end
        end
    end
    io = open("test.txt", "w")
    #print(io, "$edge_spins\n$links\n")
    eedge_spins = Int64[edge_spins...]
    sort!(eedge_spins)
    len = length(eedge_spins)*A.cluster_num
    IM = SharedMatrix{Float64}(len, len)
    fill!(IM, 0.0)
    for cluster in 1:A.cluster_num
        cl_pos = ind2sub(A.outer_cluster_dims, cluster)
        delta = collect(cl_pos) - cl_pos1
        #print(io,"$delta\n\n")
        for link in links
            p1 = sub2ind((length(eedge_spins), A.cluster_num), findfirst(eedge_spins, link[1]), cluster)
            cl_p2 = (collect(ind2sub(A.outer_cluster_dims, link[3])) .+ delta .- 1).%collect(A.outer_cluster_dims) .+ 1
            #print(io,"$cl_p2\n")
            #print(io,"$link\n")
            p2 = sub2ind((length(eedge_spins), A.cluster_num), findfirst(eedge_spins, link[2]), sub2ind(A.outer_cluster_dims, cl_p2...))
            #print(io,p1,"\t",p2,"\n")
            IM[p1,p2] = link[4]
        end
        #print(io,"\n\n\n")
    end
    close(io)
    IM .= .-IM.*sqrt(2^(length(cluster_spins))+1).*0.25
    return eedge_spins::Vector{Int64}, IM::SharedMatrix{Float64}
end

function _build_QuantumClassicalInteractions(A::HybridApprox)
    edge_spins = Set{Float64}()
    q_spins = get_q_spins(A)
    cl_spins = get_cl_spins(A)
    links = Vector{Tuple{Int64,Int64,Float64}}()
    for q in eachindex(q_spins)
        for cl in eachindex(cl_spins)
            val = get_value(A, q_spins[q], cl_spins[cl])*0.5
            if val!= 0.0
                push!(edge_spins, q)
                push!(links, (q,cl,val))
            end
        end
    end
    edge_spins = Int64[edge_spins...]
    sort!(edge_spins)
    q2cl = SharedMatrix{Float64}(length(cl_spins), length(edge_spins))
    fill!(q2cl, 0.0)
    for link in links
        q2cl[link[2], findfirst(edge_spins,link[1])] += link[3]
    end
    cl2q = SharedMatrix{Float64}(length(edge_spins), length(cl_spins))
    transpose!(cl2q,q2cl)
    scale!(cl2q,-1.0)
    q2cl .*= sqrt(2^length(q_spins)+1)
    return edge_spins, q2cl, cl2q
end
@inline build_QuantumClassicalInteractions(A::HybridApprox{Lattice{D,MDI{D2},T,C,B}}) where {D,D2,T,C,B} = _build_QuantumClassicalInteractions(A)

function build_QuantumClassicalInteractions(A::HybridApprox{Lattice{D,typeof(nearest_neighbours),T,C,B}}) where {D,T,C,B}
    edge_spins = Set{Float64}()
    q_spins = get_q_spins(A)
    cl_spins = get_cl_spins(A)
    links = Vector{Tuple{Int64,Int64,Float64}}()
    for q in eachindex(q_spins)
        for cl in eachindex(cl_spins)
            val = get_value(A, q_spins[q], cl_spins[cl])*0.5
            if val!= 0.0
                push!(edge_spins, q)
                push!(links, (q,cl,val))
            end
        end
    end
    edge_spins = Int64[edge_spins...]
    sort!(edge_spins)
    q2cl = Matrix{Float64}(length(cl_spins), length(edge_spins))
    fill!(q2cl, 0.0)
    for link in links
        q2cl[link[2], findfirst(edge_spins,link[1])] += link[3]
    end
    cl2q = -transpose(q2cl)
    q2cl .*= sqrt(2^length(q_spins)+1)
    return edge_spins, shared_sparse(q2cl)::SharedSparseMatrixCSC{Float64,Int64}, shared_sparse(cl2q)::SharedSparseMatrixCSC{Float64,Int64}
end

#####################################
#What is written here works only for long-range coupling.

function build_InterLattice_Couplings(A::CompositeApproximation{Approx1,Approx2}, edge_spins_len::Int64) where {Approx1<:HybridApprox, Approx2<:PureClassicalApprox}
    q1tocl2, cl2toq1 = build_InterLattice_q2cl(A.A1, A.A2, edge_spins_len, A.gamma)
    cl1tocl2, cl2tocl1 = build_InterLattice_cl2cl(A.A1, A.A2, A.gamma)
    return q1tocl2, cl2toq1, cl1tocl2, cl2tocl1
end

function build_InterLattice_q2cl(A1::Approx1, A2::Approx2, edge_spins_len::Int64, gamma::Float64) where {Approx1<:HybridApprox, Approx2<:PureClassicalApprox}
    q_spins = get_q_spins(A1)
    assert(length(q_spins)==edge_spins_len)
    cl_spins = get_cl_spins(A2)
    L1 = A1.L
    L2 = A2.L
    q1tocl2 = SharedMatrix{Float64}(length(cl_spins), edge_spins_len)
    fill!(q1tocl2, 0.0)
    for q in eachindex(q_spins)
        for cl in eachindex(cl_spins)
            q1tocl2[cl, q] = get_value(L1, L2, q_spins[q], cl_spins[cl])*0.5*gamma
        end
    end
    cl2toq1 = SharedMatrix{Float64}(edge_spins_len, length(cl_spins))
    transpose!(cl2toq1,q1tocl2)
    scale!(cl2toq1,-1.0)
    q1tocl2 .*= sqrt(2^edge_spins_len+1)
    return q1tocl2, cl2toq1
end

function build_InterLattice_cl2cl(A1::Approx1, A2::Approx2, gamma::Float64) where {Approx1<:HybridApprox, Approx2<:PureClassicalApprox}
    cl_spins1 = get_cl_spins(A1)
    cl_spins2 = get_cl_spins(A2)
    L1 = A1.L
    L2 = A2.L
    cl1tocl2 = SharedMatrix{Float64}(length(cl_spins2), length(cl_spins1))
    fill!(cl1tocl2, 0.0)
    for cl2 in eachindex(cl_spins2)
        for cl1 in eachindex(cl_spins1)
            cl1tocl2[cl2, cl1] = get_value(L1, L2, cl_spins1[cl1], cl_spins2[cl2])*gamma
        end
    end
    cl2tocl1 = SharedMatrix{Float64}(length(cl_spins1), length(cl_spins2))
    transpose!(cl2tocl1, cl1tocl2)
    return cl1tocl2, cl2tocl1
end