text:
module CompositeStructures

using ..Prerequisites: rawtype
using ..TypeTraits: efficientoperations

export CompositeTuple,CompositeNTuple,CompositeVector,CompositeDict,NamedContainer

"""
    CompositeTuple{T<:Tuple}

A composite tuple is a tuple that is implemented by including an ordinary `Tuple` as one of its attributes with the name `:contents`.
"""
abstract type CompositeTuple{T<:Tuple} end
Base.length(::CompositeTuple{T}) where T=fieldcount(T)
Base.length(::Type{<:CompositeTuple{T}}) where T=fieldcount(T)
Base.eltype(::CompositeTuple{T}) where T=eltype(T)
Base.eltype(::Type{<:CompositeTuple{T}}) where T=eltype(T)
@generated function Base.hash(ct::CompositeTuple,h::UInt)
    ts=Expr(:tuple,[:(getfield(ct,$i)) for i=1:fieldcount(ct)]...)
    return :(hash($ts,h))
end
Base.:(==)(ct1::CompositeTuple,ct2::CompositeTuple) = ==(efficientoperations,ct1,ct2)
Base.isequal(ct1::CompositeTuple,ct2::CompositeTuple)=isequal(efficientoperations,ct1,ct2)
Base.getindex(ct::CompositeTuple,i::Union{<:Integer,CartesianIndex})=getfield(ct,:contents)[i]
@generated function Base.getindex(ct::CompositeTuple,inds)
    exprs=[name==:contents ? :(getfield(ct,:contents)[inds]) : :(getfield(ct,$i)) for (i,name) in enumerate(ct|>fieldnames)]
    return :(rawtype(typeof(ct))($(exprs...)))
end
Base.lastindex(ct::CompositeTuple)=lastindex(getfield(ct,:contents))
Base.iterate(ct::CompositeTuple)=iterate(getfield(ct,:contents))
Base.iterate(ct::CompositeTuple,state)=iterate(getfield(ct,:contents),state)
Base.iterate(rv::Iterators.Reverse{<:CompositeTuple},state=length(rv.itr))=state<1 ? nothing : (rv.itr[state],state-1)
Base.keys(ct::CompositeTuple)=keys(getfield(ct,:contents))
Base.values(ct::CompositeTuple)=values(getfield(ct,:contents))
Base.pairs(ct::CompositeTuple)=pairs(getfield(ct,:contents))
@generated function Base.reverse(ct::CompositeTuple)
    exprs=[name==:contents ? :(reverse(getfield(ct,:contents))) : :(getfield(ct,$i)) for (i,name) in enumerate(ct|>fieldnames)]
    return :(rawtype(typeof(ct))($(exprs...)))
end
Base.convert(::Type{Tuple},ct::CompositeTuple)=getfield(ct,:contents)

"""
    CompositeNTuple{N,T}

A composite ntuple is a ntuple that is implemented by including an ordinary `NTuple` as one of its attributes with the name `:contents`.

Alias for `CompositeTuple{NTuple{N,T}}`.
"""
const CompositeNTuple{N,T}=CompositeTuple{NTuple{N,T}}

"""
    CompositeVector{T}

A composite vector is a vector that is implemented by including an ordinary `Vector` as one of its attributes with the name `:contents`.
"""
abstract type CompositeVector{T} <:AbstractVector{T} end
Base.size(cv::CompositeVector)=size(getfield(cv,:contents))
Base.size(cv::CompositeVector,i)=size(getfield(cv,:contents),i)
Base.length(cv::CompositeVector)=length(getfield(cv,:contents))
Base.:(==)(cv1::CompositeVector,cv2::CompositeVector) = ==(efficientoperations,cv1,cv2)
Base.isequal(cv1::CompositeVector,cv2::CompositeVector)=isequal(efficientoperations,cv1,cv2)
Base.getindex(cv::CompositeVector,i::Union{<:Integer,CartesianIndex})=getfield(cv,:contents)[i]
@generated function Base.getindex(cv::CompositeVector,inds)
    exprs=[name==:contents ? :(getfield(cv,:contents)[inds]) : :(getfield(cv,$i)) for (i,name) in enumerate(cv|>fieldnames)]
    return :(rawtype(typeof(cv))($(exprs...)))
end
Base.lastindex(cv::CompositeVector)=lastindex(getfield(cv,:contents))
Base.setindex!(cv::CompositeVector,value,inds)=(getfield(cv,:contents)[inds]=value)
Base.push!(cv::CompositeVector,values...)=(push!(getfield(cv,:contents),values...);cv)
Base.pushfirst!(cv::CompositeVector,values...)=(pushfirst!(getfield(cv,:contents),values...);cv)
Base.insert!(cv::CompositeVector,index::Integer,value)=(insert!(getfield(cv,:contents),index,value);cv)
Base.append!(cv::CompositeVector,values)=(append!(getfield(cv,:contents),values);cv)
Base.prepend!(cv::CompositeVector,values)=(prepend!(getfield(cv,:contents),values);cv)
Base.splice!(cv::CompositeVector,index::Integer,replacement=Base._default_splice)=splice!(getfield(cv,:contents),index,replacement)
@generated function Base.splice!(cv::CompositeVector,range::UnitRange{<:Integer},replacement=Base._default_splice)
    exprs=[name==:contents ? :(splice!(getfield(cv,:contents),range,replacement)) : :(getfield(cv,$i)) for (i,name) in enumerate(cv|>fieldnames)]
    return :(rawtype(typeof(cv))($(exprs...)))
end
Base.deleteat!(cv::CompositeVector,indices)=(deleteat!(getfield(cv,:contents),indices);cv)
Base.pop!(cv::CompositeVector)=pop!(getfield(cv,:contents))
Base.popfirst!(cv::CompositeVector)=popfirst!(getfield(cv,:contents))
Base.empty!(cv::CompositeVector)=(empty!(getfield(cv,:contents));cv)
@generated function Base.empty(cv::CompositeVector)
    exprs=[name==:contents ? :(empty(getfield(cv,:contents))) : :(getfield(cv,$i)) for (i,name) in enumerate(cv|>fieldnames)]
    return :(rawtype(typeof(cv))($(exprs...)))
end
@generated function Base.reverse(cv::CompositeVector)
    exprs=[name==:contents ? :(reverse(getfield(cv,:contents))) : :(getfield(cv,$i)) for (i,name) in enumerate(cv|>fieldnames)]
    return :(rawtype(typeof(cv))($(exprs...)))
end
@generated function Base.similar(cv::CompositeVector,dtype::Type{T}=eltype(cv),dims::NTuple{N,Int}=size(cv)) where {T,N}
    exprs=[name==:contents ? :(similar(getfield(cv,:contents),dtype,dims)) : :(getfield(cv,$i)) for (i,name) in enumerate(cv|>fieldnames)]
    return :(rawtype(typeof(cv))($(exprs...)))
end
Base.iterate(cv::CompositeVector,state=1)=iterate(getfield(cv,:contents),state)
Base.keys(cv::CompositeVector)=keys(getfield(cv,:contents))
Base.values(cv::CompositeVector)=values(getfield(cv,:contents))
Base.pairs(cv::CompositeVector)=pairs(getfield(cv,:contents))
Base.convert(::Type{Vector},cv::CompositeVector)=getfield(cv,:contents)

"""
    CompositeDict{K,V}

A composite dict is a dict that is implemented by including an ordinary `Dict` as one of its attributes with the name `:contents`.
"""
abstract type CompositeDict{K,V} <: AbstractDict{K,V} end
Base.isempty(cd::CompositeDict)=isempty(getfield(cd,:contents))
Base.length(cd::CompositeDict)=length(getfield(cd,:contents))
Base.haskey(cd::CompositeDict,key)=haskey(getfield(cd,:contents),key)
Base.in(p::Pair,cd::CompositeDict,valcmp=(==))=in(p,getfield(cd,:contents),valcmp)
Base.:(==)(cd1::CompositeDict,cd2::CompositeDict) = ==(efficientoperations,cd1,cd2)
Base.isequal(cd1::CompositeDict,cd2::CompositeDict)=isequal(efficientoperations,cd1,cd2)
Base.get(cd::CompositeDict,key,default)=get(getfield(cd,:contents),key,default)
Base.get(f::Base.Callable,cd::CompositeDict,key)=get(f,getfield(cd,:contents),key)
Base.getkey(cd::CompositeDict,key,default)=getkey(getfield(cd,:contents),key,default)
Base.getindex(cd::CompositeDict{K,V},index::K) where {K,V}=getfield(cd,:contents)[index]
Base.push!(cd::CompositeDict,ps::Pair...)=(push!(getfield(cd,:contents),ps...);cd)
Base.get!(cd::CompositeDict,key,default)=get!(getfield(cd,:contents),key,default)
Base.get!(default,cd::CompositeDict,key)=get!(default,getfield(cd,:contents),key)
Base.setindex!(cd::CompositeDict{K,V},value::V,index::K) where {K,V}=(getfield(cd,:contents)[index]=value)
Base.pop!(cd::CompositeDict)=pop!(getfield(cd,:contents))
Base.pop!(cd::CompositeDict,key)=pop!(getfield(cd,:contents),key)
Base.pop!(cd::CompositeDict,key,default)=pop!(getfield(cd,:contents),key,default)
Base.delete!(cd::CompositeDict,key)=(delete!(getfield(cd,:contents),key);cd)
Base.empty!(cd::CompositeDict)=(empty!(getfield(cd,:contents));cd)
Base.merge(cd::CD,others::CD...) where CD<:CompositeDict=merge!(empty(cd),cd,others...)
Base.merge(combine::Function,cd::CD,others::CD...) where CD<:CompositeDict=merge!(combine,empty(cd),cd,others...)
@generated function Base.empty(cd::CompositeDict)
    exprs=[name==:contents ? :(Dict{cd|>keytype,cd|>valtype}()) : :(getfield(cd,$i)) for (i,name) in enumerate(cd|>fieldnames)]
    return :(rawtype(typeof(cd))($(exprs...)))
end
Base.iterate(cd::CompositeDict)=iterate(getfield(cd,:contents))
Base.iterate(cd::CompositeDict,state)=iterate(getfield(cd,:contents),state)
Base.keys(cd::CompositeDict)=keys(getfield(cd,:contents))
Base.values(cd::CompositeDict)=values(getfield(cd,:contents))
Base.pairs(cd::CompositeDict)=pairs(getfield(cd,:contents))
Base.convert(::Type{Dict},cd::CompositeDict)=getfield(cd,:contents)

"""
    NamedContainer{T,Names}=NamedTuple{Names,<:Tuple{Vararg{T}}}

NamedContainer is just a wrapper of Julia NamedTuple, but not a composite type.
"""
const NamedContainer{T,Names}=NamedTuple{Names,<:Tuple{Vararg{T}}}

"""
    NamedContainer{Names}(contents) where Names -> NamedTuple{Names,typeof(contents)}

Construct a named container.
"""
@generated function NamedContainer{Names}(contents::Tuple) where Names
    @assert length(Names)==fieldcount(contents) "NamedContainer error: dismatched length between names and contents."
    fieldcount(contents)==0 && return NamedTuple()
    return Expr(:tuple,[:($name=contents[$i]) for (i,name) in enumerate(Names)]...)
end

end #module


text:
function get_raw_sections()
    return ["CASE IDENTIFICATION", "BUS", "LOAD", "FIXED SHUNT",
            "GENERATOR", "BRANCH", "TRANSFORMER", "AREA INTERCHANGE",
            "TWO-TERMINAL DC", "VOLTAGE SOURCE CONVERTER",
            "IMPEDANCE CORRECTION", "MULTI-TERMINAL DC",
            "MULTI-SECTION LINE", "ZONE", "INTER-AREA TRANSFER", "OWNER",
            "FACTS CONTROL DEVICE", "SWITCHED SHUNT", "GNE DEVICE",
            "INDUCTION MACHINE"]
end

function get_raw_section_info()
    case_identification = [["SBASE", 2, Float64]]
    bus = [["I", 1, Int], ["AREA", 5, Int], ["VM", 8, Float64],
           ["VA", 9, Float64], ["NVHI", 10, Float64], ["NVLO", 11, Float64],
           ["EVHI", 12, Float64], ["EVLO", 13, Float64], ["TYPE", 4, Int]]
    load = [["I", 1, Int], ["ID", 2, String], ["STATUS", 3, Int],
            ["PL", 6, Float64], ["QL", 7, Float64]]
    fixed_shunt = [["I", 1, Int], ["ID", 2, String], ["STATUS", 3, Int],
                   ["GL", 4, Float64], ["BL", 5, Float64]]
    generator = [["I", 1, Int], ["ID", 2, String], ["PG", 3, Float64],
                 ["QG", 4, Float64], ["QT", 5, Float64], ["QB", 6, Float64],
                 ["STAT", 15, Int], ["PT", 17, Float64], ["PB", 18, Float64]]
    branch = [["I", 1, Int], ["J", 2, Int], ["CKT", 3, String],
              ["R", 4, Float64], ["X", 5, Float64], ["B", 6, Float64],
              ["RATEA", 7, Float64], ["RATEC", 9, Float64], ["ST", 14, Int]]
    transformer = [["I", 1, Int], ["J", 2, Int], ["CKT", 4, String],
                   ["MAG1", 8, Float64], ["MAG2", 9, Float64],
                   ["STAT", 12, Int], ["R12", 22, Float64],
                   ["X12", 23, Float64], ["WINDV1", 25, Float64],
                   ["ANG1", 27, Float64], ["RATA1", 28, Float64],
                   ["RATC1", 30, Float64], ["WINDV2", 42, Float64]]
    switched_shunt = [["I", 1, Int], ["STAT", 4, Int], ["BINIT", 10, Float64],
                      ["N1", 11, Float64], ["B1", 12, Float64],
                      ["N2", 13, Float64], ["B2", 14, Float64],
                      ["N3", 15, Float64], ["B3", 16, Float64],
                      ["N4", 17, Float64], ["B4", 18, Float64],
                      ["N5", 19, Float64], ["B5", 20, Float64],
                      ["N6", 21, Float64], ["B6", 22, Float64],
                      ["N7", 23, Float64], ["B7", 24, Float64],
                      ["N8", 25, Float64], ["B8", 26, Float64]]

    section_info = Dict{String,Array}(
        "CASE IDENTIFICATION" => case_identification,
        "BUS" => bus,
        "LOAD" => load,
        "FIXED SHUNT" => fixed_shunt,
        "GENERATOR" => generator,
        "BRANCH" => branch,
        "TRANSFORMER" => transformer,
        "AREA INTERCHANGE" => [],
        "TWO-TERMINAL DC" => [],
        "VOLTAGE SOURCE CONVERTER" => [],
        "IMPEDANCE CORRECTION" => [],
        "MULTI-TERMINAL DC" => [],
        "MULTI-SECTION LINE" => [],
        "ZONE" => [],
        "INTER-AREA TRANSFER" => [],
        "OWNER" => [],
        "FACTS CONTROL DEVICE" => [],
        "SWITCHED SHUNT" => switched_shunt,
        "GNE DEVICE" => [],
        "INDUCTION MACHINE" => [])

    return section_info
end

function idx_bus()
    BUS_B = 1
    BUS_AREA = 2
    BUS_VM = 3
    BUS_VA = 4
    BUS_NVHI = 5
    BUS_NVLO = 6
    BUS_EVHI = 7
    BUS_EVLO = 8
    BUS_TYPE = 9

    return BUS_B, BUS_AREA, BUS_VM, BUS_VA, BUS_NVHI, BUS_NVLO, BUS_EVHI, BUS_EVLO, BUS_TYPE
end

function idx_gen()
    GEN_BUS = 1
    GEN_ID = 2
    GEN_PG = 3
    GEN_QG = 4
    GEN_QT = 5
    GEN_QB = 6
    GEN_STAT = 7
    GEN_PT = 8
    GEN_PB = 9

    return GEN_BUS, GEN_ID, GEN_PG, GEN_QG, GEN_QT, GEN_QB, GEN_STAT, GEN_PT, GEN_PB
end

function idx_load()
    LOAD_BUS = 1
    LOAD_ID = 2
    LOAD_STAT = 3
    LOAD_PL = 4
    LOAD_QL = 5
    return LOAD_BUS, LOAD_ID, LOAD_STAT, LOAD_PL, LOAD_QL
end

function idx_branch()
    BR_FR = 1
    BR_TO = 2
    BR_ID = 3
    BR_R = 4
    BR_X = 5
    BR_B = 6
    BR_RATEA = 7
    BR_RATEC = 8
    BR_STAT = 9

    return BR_FR, BR_TO, BR_ID, BR_R, BR_X, BR_B, BR_RATEA, BR_RATEC, BR_STAT
end

function idx_fshunt()
    FSH_BUS = 1
    FSH_ID = 2
    FSH_STAT = 3
    FSH_G = 4
    FSH_B = 5

    return FSH_BUS, FSH_ID, FSH_STAT, FSH_G, FSH_B
end

function idx_transformer()
    TR_FR = 1
    TR_TO = 2
    TR_ID = 3
    TR_MAG1 = 4
    TR_MAG2 = 5
    TR_STAT = 6
    TR_R = 7
    TR_X = 8
    TR_WINDV1 = 9
    TR_ANG = 10
    TR_RATEA = 11
    TR_RATEC = 12
    TR_WINDV2 = 13

    return TR_FR, TR_TO, TR_ID, TR_MAG1, TR_MAG2, TR_STAT, TR_R, TR_X,
           TR_WINDV1, TR_ANG, TR_RATEA, TR_RATEC, TR_WINDV2
end

function idx_sshunt()
    SSH_BUS = 1
    SSH_STAT = 2
    SSH_BINIT = 3
    SSH_N1 = 4
    SSH_B1 = 5

    return SSH_BUS, SSH_STAT, SSH_BINIT, SSH_N1, SSH_B1
end

function raw_xfmr(lines::Array{String}, pos::Int, info::Array;
                  delim=",")
    start = last = pos
    while strip(lines[last])[1] ∉ [ '0', 'Q' ]
        last += 4
    end

    ndata = Int((last - start)/4)
    data = Array{Any}(undef, ndata, length(info))

    for i=1:ndata
        # Each transformer record consists of 4 lines.
        # We read them all and make a single line.
        offset = 4*(i-1)
        line = replace(lines[start + offset], "'" => "")
        for j=1:3
            line *= ","*replace(lines[start + offset + j], "'" => "")
        end

        fields = split(line, delim)
        for j=1:length(info)
            if info[j][3] == String
                data[i,j] = strip(fields[info[j][2]])
            else
                data[i,j] = parse(info[j][3], fields[info[j][2]])
            end
        end
    end

    pos = last + 1
    return pos, data
end

function parse_raw(filename::AbstractString; delim=",")
    f = open(filename, "r")
    lines = readlines(f)
    close(f)

    sections = get_raw_sections()
    section_info = get_raw_section_info()
    pos = 1

    rawdata = Dict{String,Array}()
    for s in sections
        info = section_info[s]

        if isempty(info)
            pos = parse_skipsection(lines, pos) + 1
            continue
        end

        if s == "CASE IDENTIFICATION"
            line = lines[pos]
            data = [parse(Float64, split(line, delim)[info[1][2]])]
            pos += 3
        elseif s == "TRANSFORMER"
            pos, data = raw_xfmr(lines, pos, info; delim=delim)
        else
            pos, data = parse_data(lines, pos, info; delim=delim)
        end

        rawdata[s] = data
    end

    return rawdata
end

function raw_to_exapf(rawdata::Dict{String, Array})

    BUS_B, BUS_AREA, BUS_VM, BUS_VA, BUS_NVHI, BUS_NVLO, BUS_EVHI,
    BUS_EVLO, BUS_TYPE = idx_bus()
    GEN_BUS, GEN_ID, GEN_PG, GEN_QG, GEN_QT, GEN_QB, GEN_STAT,
    GEN_PT, GEN_PB = idx_gen()
    LOAD_BUS, LOAD_ID, LOAD_STATUS, LOAD_PL, LOAD_QL = idx_load()
    BR_FR, BR_TO, BR_CKT, BR_R, BR_X, BR_B, BR_RATEA, BR_RATEC,
        BR_STAT = idx_branch()
    TR_FR, TR_TO, TR_CKT, TR_MAG1, TR_MAG2, TR_STAT, TR_R, TR_X, TR_WINDV1,
        TR_ANG, TR_RATEA, TR_RATEC, TR_WINDV2 = idx_transformer()
    FSH_BUS, FSH_ID, FSH_STAT, FSH_G, FSH_B = idx_fshunt()

    data = Dict{String,Array}()
    bus = rawdata["BUS"]
    gen = rawdata["GENERATOR"]
    load = rawdata["LOAD"]
    branch = rawdata["BRANCH"]
    trans = rawdata["TRANSFORMER"]
    fsh = rawdata["FIXED SHUNT"]
    SBASE = rawdata["CASE IDENTIFICATION"][1]

    nbus = size(bus, 1)
    ngen = size(gen, 1)
    nload = size(load, 1)
    nbranch = size(branch, 1)
    ntrans = size(trans, 1)
    nfsh = size(fsh, 1)

    # create new data structure
    data = Dict{String,Array}()

    # perunit base
    data["baseMVA"] = [SBASE]

    # buses
    bus_array = zeros(nbus, 17)
    for i in 1:nbus
        bus_array[i, 1] = bus[i, BUS_B]
        bus_array[i, 2] = bus[i, BUS_TYPE]
        bus_array[i, 3] = 0.0
        bus_array[i, 4] = 0.0
        bus_array[i, 5] = 0.0
        bus_array[i, 6] = 0.0
        bus_array[i, 7] = bus[i, BUS_AREA]
        bus_array[i, 8] = bus[i, BUS_VM]
        bus_array[i, 9] = bus[i, BUS_VA]
    end

    # fixed shunt
    for i in 1:nfsh
        if fsh[i, FSH_STAT] == 0
            break
        end
        bus_idx = fsh[i, FSH_BUS]
        bus_array[bus_idx, 5] += fsh[i, FSH_G]
        bus_array[bus_idx, 6] += fsh[i, FSH_B]
    end

    # add load
    for i in 1:nload
        if load[i, LOAD_STATUS] == 0
            break
        end
        bus_idx = load[i, LOAD_BUS]
        Pd = load[i, LOAD_PL]
        Qd = load[i, LOAD_QL]
        bus_array[bus_idx, 3] += Pd
        bus_array[bus_idx, 4] += Qd
    end
    data["bus"] = bus_array

    # gens
    gen_array = zeros(ngen, 25)
    for i in 1:ngen
        gen_array[i, 1] = gen[i, GEN_BUS]
        gen_array[i, 2] = gen[i, GEN_PG]
        gen_array[i, 3] = gen[i, GEN_QG]
        gen_array[i, 8] = gen[i, GEN_STAT]
    end
    data["gen"] = gen_array

    # branches
    branch_array = zeros(nbranch + ntrans, 21)
    for i in 1:nbranch
        branch_array[i, 1] = branch[i, BR_FR]
        branch_array[i, 2] = branch[i, BR_TO]
        branch_array[i, 3] = branch[i, BR_R]
        branch_array[i, 4] = branch[i, BR_X]
        branch_array[i, 5] = branch[i, BR_B]
        branch_array[i, 9] = 0.0
        branch_array[i, 10] = 0.0
        branch_array[i, 11] = branch[i, BR_STAT]
    end

    for i in 1:ntrans
        branch_array[nbranch + i, 1] = trans[i, TR_FR]
        branch_array[nbranch + i, 2] = trans[i, TR_TO]
        branch_array[nbranch + i, 3] = trans[i, TR_R]
        branch_array[nbranch + i, 4] = trans[i, TR_X]
        branch_array[nbranch + i, 5] = 0.0
        tap = trans[i, TR_WINDV1] / trans[i, TR_WINDV2]
        branch_array[nbranch + i, 9] = tap
        branch_array[nbranch + i, 10] = trans[i, TR_ANG]
        branch_array[nbranch + i, 11] = trans[i, TR_STAT]
    end
    data["branch"] = branch_array

    return data
end

ext: jl


text:
export Runtime

export assert_runtime_alive

const runtime_ref_count = Int[0]

type Runtime
    is_alive::Bool

    function Runtime()
        test_status(hsa_init())
        runtime_ref_count[1] += 1
        val = new(true)
        finalizer(val, finalize_runtime)
        return val
    end
end

function finalize_runtime(r::Runtime)
    if r.is_alive
        err = hsa_shut_down()
        r.is_alive = false
        runtime_ref_count[1] -= 1

        test_status(err)
    end
end

function assert_runtime_alive()
    @assert runtime_ref_count[1] > 0
end


text:
function mark_cross_with_barriers!(r)
    for side in 0:3
        counters = move_through_barrier!(r, side)
        putmarker!(r)
        return_to_start!(r, counters, side)
    end
    putmarker!(r)
end

function move_through_barrier!(r, side)
    counters = []
    while (length(counters)==0 || counters[length(counters)]!=0)
        count1 = 0
        count2 = 0
        count3 = 0
        while !isborder(r, HorizonSide(side))
            move!(r, HorizonSide(side))
            putmarker!(r)
            count1+=1
        end
        while (isborder(r, HorizonSide(side)) && !isborder(r, HorizonSide((side+1)%4)))
            move!(r, HorizonSide((side+1)%4))
            count2 += 1
        end
        if !isborder(r, HorizonSide(side))
            move!(r, HorizonSide(side))
            count3 += 1
            while isborder(r, HorizonSide((side+3)%4))
                move!(r, HorizonSide(side))
                count3 += 1
            end
            for i in 1:count2
                move!(r, HorizonSide((side+3)%4))
            end
        else
            for i in 1:count2
                move!(r, HorizonSide((side+3)%4))
            end
            count2 = 0
        end
        putmarker!(r)
        push!(counters,count1)
        push!(counters,count2)
        push!(counters,count3)
        push!(counters,count2)
    end
    return counters
end

function move_forward!(r, side, count)
    for i in 0:count-1
        move!(r, HorizonSide(side))
    end
end

function return_to_start!(r, counters, side)
    counters = reverse(counters)
    for i in 1:length(counters)
        rside = side
        if i%2 == 1
            rside = (rside+i)%4
        else
            rside = (rside+2)%4
        end
        move_forward!(r, rside, counters[i])
    end
end

#mark_cross_with_barriers!(r)


text:
function calculate_residual(flux_vector,STM,W)

    # compute the species error -
    species_error_vector = STM*flux_vector

    # compute the residual -
    residual = transpose(species_error_vector)*W*species_error_vector

    # return -
    return residual
end

function calculate_gradient!(G,flux_vector,STM,W)

    # compute the residual -
    residual = STM*flux_vector

    # compute the grad -
    tmp_vector = transpose(residual)*(transpose(W)+W)*STM
    index = 1
    for value in tmp_vector
        G[index] = value
        index = index + 1
    end
end

function sample_flux_space(solution_bounds_array::Array{Float64,2}, data_dictionary::Dict{String,Any}, number_of_samples::Int64)

    # Declare a progress meter for user feedback -
    p = Progress(number_of_samples,color=:yellow)

    # get some stuff from the data dictionary -
    STM = data_dictionary["stoichiometric_matrix"]

    # Compute the W matrix -
    (number_of_species,number_of_reactions) = size(STM)
    W = 10000*Matrix{Float64}(I,number_of_species,number_of_species)

    # setup the call to the optimizer --
    objective_function(x) = calculate_residual(x,STM,W)
    gradient_function!(G,x) = calculate_gradient!(G,x,STM,W)

    # get the bounds array -
    lower_bound_array = solution_bounds_array[:,1]
    upper_bound_array = solution_bounds_array[:,2]

    # initialize the ensemble -
    number_of_fluxes = length(lower_bound_array)
    flux_ensemble = zeros(number_of_fluxes,1)

    sobol_sequence = SobolSeq(number_of_fluxes)
    good_sample_count = 1
    q_vector = Float64[]
    flux = Float64[]
    for sample in sobol_sequence

        # generate a q_vector -
        q_vector = lower_bound_array.*sample + (1 .- sample).*upper_bound_array

        # make a call to the optim package -
        result = optimize(objective_function,gradient_function!,lower_bound_array, upper_bound_array, q_vector, Fminbox(LBFGS()))

        # get the flux -
        flux = Optim.minimizer(result)

        # check the error?
        converged_flag = Optim.converged(result)
        local_error = STM*flux

        if (converged_flag == true && norm(local_error)<1e-8)

            # user message -
            msg = "Completed $(good_sample_count) of $(number_of_samples) trials ..."

            # update the progress bar -
            ProgressMeter.next!(p; showvalues = [(:status,msg)])

            # package -
            flux_ensemble = [flux_ensemble flux]

            # update the sample count -
            good_sample_count = good_sample_count + 1
        end

        # should we stop?
        if (good_sample_count>number_of_samples)
            @info "Completed ...\r"
            break
        end
    end

    # return -
    return flux_ensemble[:,2:end]
end


text:
type FileMonitor
    handle::Ptr{Void}
    cb::Callback
    open::Bool
    notify::Condition
    function FileMonitor(cb, file)
        handle = c_malloc(_sizeof_uv_fs_event)
        err = ccall(:jl_fs_event_init,Int32, (Ptr{Void}, Ptr{Void}, Ptr{Uint8}, Int32), eventloop(),handle,file,0)
        if err < 0
            c_free(handle)
            throw(UVError("FileMonitor",err))
        end
        this = new(handle,cb,false,Condition())
        associate_julia_struct(handle,this)
        finalizer(this,uvfinalize)
        this
    end
    FileMonitor(file) = FileMonitor(false,file)
end

function close(t::FileMonitor)
    if t.handle != C_NULL
        ccall(:jl_close_uv,Void,(Ptr{Void},),t.handle)
    end
end

immutable FileEvent
    renamed::Bool
    changed::Bool
    timedout::Bool
end
# libuv file watching event flags
const UV_RENAME = 1
const UV_CHANGE = 2
const FE_TIMEDOUT = 4
FileEvent() = FileEvent(false,false,false)
FileEvent(flags::Integer) = FileEvent((flags & UV_RENAME) != 0,
                                      (flags & UV_CHANGE) != 0,
                                      (flags & FE_TIMEDOUT) != 0)

immutable FDEvent
    readable::Bool
    writable::Bool
    timedout::Bool
end
# libuv file descriptor event flags
const UV_READABLE = 1
const UV_WRITABLE = 2
const FD_TIMEDOUT = 4

isreadable(f::FDEvent) = f.readable
iswritable(f::FDEvent) = f.writable
FDEvent() = FDEvent(false,false,false)
FDEvent(flags::Integer) = FDEvent((flags & UV_READABLE) != 0,
                                  (flags & UV_WRITABLE) != 0,
                                  (flags & FD_TIMEDOUT) != 0)
fdtimeout() = FDEvent(false,false,true)

#Wrapper for an OS file descriptor (for Windows)
@windows_only immutable WindowsRawSocket
    handle::Ptr{Void}   # On Windows file descriptors are HANDLE's and 64-bit on 64-bit Windows...
end

@windows_only export WindowsRawSocket

abstract UVPollingWatcher

type PollingFileWatcher <: UVPollingWatcher
    handle::Ptr{Void}
    file::String
    open::Bool
    notify::Condition
    cb::Callback
    function PollingFileWatcher(cb, file)
        handle = c_malloc(_sizeof_uv_fs_poll)
        err = ccall(:uv_fs_poll_init,Int32,(Ptr{Void},Ptr{Void}),eventloop(),handle)
        if err < 0
            c_free(handle)
            throw(UVError("PollingFileWatcher",err))
        end
        this = new(handle, file, false, Condition(), cb)
        associate_julia_struct(handle,this)
        finalizer(this,uvfinalize)
        this
    end
    PollingFileWatcher(file) =  PollingFileWatcher(false,file)
end

@unix_only typealias FDW_FD RawFD
@windows_only typealias FDW_FD WindowsRawSocket

@unix_only _get_osfhandle(fd::RawFD) = fd
@windows_only _get_osfhandle(fd::RawFD) = WindowsRawSocket(ccall(:_get_osfhandle,Ptr{Void},(Cint,),fd.fd))
@windows_only _get_osfhandle(fd::WindowsRawSocket) = fd

type FDWatcher <: UVPollingWatcher
    handle::Ptr{Void}
    fd::FDW_FD
    open::Bool
    notify::Condition
    cb::Callback
    events::FDEvent
    FDWatcher(handle::Ptr,fd,open::Bool,notify::Condition,cb::Callback,events::FDEvent) =
        new(handle,_get_osfhandle(fd),open,notify,cb,events)
end
function FDWatcher(fd::RawFD)
    handle = c_malloc(_sizeof_uv_poll)
    @unix_only if ccall(:jl_uv_unix_fd_is_watched,Int32,(Int32,Ptr{Void},Ptr{Void}),fd.fd,handle,eventloop()) == 1
        c_free(handle)
        error("file descriptor $(fd.fd) is already being watched by another watcher")
    end
    err = ccall(:uv_poll_init,Int32,(Ptr{Void},Ptr{Void},Int32),eventloop(),handle,fd.fd)
    if err < 0
        c_free(handle)
        throw(UVError("FDWatcher",err))
    end
    this = FDWatcher(handle,fd,false,Condition(),false,FDEvent())
    associate_julia_struct(handle,this)
    finalizer(this,uvfinalize)
    this
end
@windows_only function FDWatcher(fd::WindowsRawSocket)
    handle = c_malloc(_sizeof_uv_poll)
    err = ccall(:uv_poll_init_socket,Int32,(Ptr{Void},   Ptr{Void}, Ptr{Void}),
                                            eventloop(), handle,    fd.handle)
    if err < 0
        c_free(handle)
        throw(UVError("FDWatcher",err))
    end
    this = FDWatcher(handle,fd,false,Condition(),false,FDEvent())
    associate_julia_struct(handle,this)
    finalizer(this,uvfinalize)
    this
end

function fdw_wait_cb(fdw::FDWatcher, events::FDEvent, status)
    if status == -1
        notify_error(fdw.notify,UVError("FDWatcher",status))
    else
        notify(fdw.notify,events)
    end
end

function _wait(fdw::FDWatcher,readable,writable)
    if !readable && !writable
        error("must watch for at least one event")
    end
    events = FDEvent(readable | fdw.events.readable,
                     writable | fdw.events.writable,
                     fdw.events.timedout)
    if !fdw.open || (events != fdw.events)
        # (re)initialize fdw
        start_watching(fdw_wait_cb,fdw,events)
    end
    while true
        events = wait(fdw.notify)
        if isa(events, FDEvent) &&
            ((readable && isreadable(events)) || (writable && iswritable(events)))
            break
        end
    end
    if isempty(fdw.notify.waitq)
        stop_watching(fdw)
    end
    events
end

# On Unix we can only have one watcher per FD, so we need to keep an explicit
# list of them. On Windows, I think it is techincally possible to have more than one
# watcher per FD, but in order to keep compatibility, we do the same on windows as we do
# on unix

let
    global fdwatcher_init, wait
    @unix_only begin
        local fdwatcher_array
        function fdwatcher_init()
            fdwatcher_array = Array(FDWatcher,0)
        end

        function wait(fd::RawFD; readable=false, writable=false)
            old_length = length(fdwatcher_array)
            if fd.fd+1 > old_length
                resize!(fdwatcher_array,fd.fd+1)
            end
            if !isdefined(fdwatcher_array,fd.fd+1)
                fdwatcher_array[fd.fd+1] = FDWatcher(fd)
            end
            _wait(fdwatcher_array[fd.fd+1],readable,writable)
        end
    end
    @windows_only begin
        local fdwatcher_array
        function fdwatcher_init()
            fdwatcher_array = Dict{WindowsRawSocket,FDWatcher}()
        end

        function wait(fd::RawFD; readable=false, writable=false)
            wait(_get_osfhandle(fd); readable=readable, writable=writable)
        end

        function wait(socket::WindowsRawSocket; readable=false, writable=false)
            if !haskey(fdwatcher_array,socket.handle)
                fdwatcher_array[socket] = FDWatcher(socket)
            end
            _wait(fdwatcher_array[socket],readable,writable)
        end
    end
end

function pfw_wait_cb(pfw::PollingFileWatcher, prev, cur, status)
    if status < 0
        notify_error(pfw.notify,UVError("PollingFileWatcher",status))
    else
        notify(pfw.notify,(prev,cur))
    end
end

function wait(pfw::PollingFileWatcher; interval=2.0)
    if !pfw.open
        start_watching(pfw_wait_cb,pfw,interval)
    end
    prev,curr = stream_wait(pfw,pfw.notify)
    if isempty(pfw.notify.waitq)
        stop_watching(pfw)
    end
    (prev,curr)
end

function wait(m::FileMonitor)
    err, filename, events = stream_wait(m,m.notify)
    filename, events
end

close(t::UVPollingWatcher) = ccall(:jl_close_uv,Void,(Ptr{Void},),t.handle)

function start_watching(t::FDWatcher, events::FDEvent)
    associate_julia_struct(t.handle, t)
    @unix_only if ccall(:jl_uv_unix_fd_is_watched,Int32,(Int32,Ptr{Void},Ptr{Void}),t.fd,t.handle,eventloop()) == 1
        error("cannot watch an FD more than once on Unix")
    end
    uv_error("start_watching (FD)",
        ccall(:jl_poll_start, Int32, (Ptr{Void},Int32), t.handle,
              events.readable*UV_READABLE + events.writable*UV_WRITABLE + events.timedout*FD_TIMEDOUT))
end
start_watching(f::Function, t::FDWatcher, events::FDEvent) = (t.cb = f; start_watching(t,events))

function start_watching(t::PollingFileWatcher, interval)
    associate_julia_struct(t.handle, t)
    uv_error("start_watching (File)",
             ccall(:jl_fs_poll_start, Int32, (Ptr{Void},Ptr{Uint8},Uint32),
                   t.handle, t.file, iround(interval*1000)))
end
start_watching(f::Function, t::PollingFileWatcher, interval) = (t.cb = f;start_watching(t,interval))

function stop_watching(t::FDWatcher)
    disassociate_julia_struct(t.handle)
    uv_error("stop_watching (FD)",
        ccall(:uv_poll_stop,Int32,(Ptr{Void},),t.handle))
end

function stop_watching(t::PollingFileWatcher)
    disassociate_julia_struct(t.handle)
    uv_error("stop_watching (File)",
        ccall(:uv_fs_poll_stop,Int32,(Ptr{Void},),t.handle))
end

function _uv_hook_fseventscb(t::FileMonitor,filename::Ptr,events::Int32,status::Int32)
    fname = filename == C_NULL ? "" : bytestring(convert(Ptr{Uint8},filename))
    fe = FileEvent(events)
    if isa(t.cb,Function)
        t.cb(fname, fe, status)
    end
    if status < 0
        notify_error(t.notify,(UVError("FileMonitor",status), fname, fe))
    else
        notify(t.notify,(status, fname, fe))
    end
end

function _uv_hook_pollcb(t::FDWatcher,status::Int32,events::Int32)
    if isa(t.cb,Function)
        t.cb(t, FDEvent(events), status)
    end
end

function _uv_hook_fspollcb(t::PollingFileWatcher,status::Int32,prev::Ptr,cur::Ptr)
    if isa(t.cb,Function)
        t.cb(t, StatStruct(convert(Ptr{Uint8},prev)), StatStruct(convert(Ptr{Uint8},cur)), status)
    end
end

_uv_hook_close(uv::FileMonitor) = (uv.handle = 0; nothing)
_uv_hook_close(uv::UVPollingWatcher) = (uv.handle = 0; nothing)

function poll_fd(s, seconds::Real; readable=false, writable=false)
    wt = Condition()

    @schedule (args = wait(s; readable=readable, writable=writable); notify(wt,(:poll,args)))
    @schedule (sleep(seconds); notify(wt,(:timeout,fdtimeout())))

    _, ret = wait(wt)

    return ret
end

function poll_file(s, interval_seconds::Real, seconds::Real)
    wt = Condition()
    pfw = PollingFileWatcher(s)

    @schedule (wait(pfw;interval=interval_seconds); notify(wt,(:poll)))
    @schedule (sleep(seconds); notify(wt,(:timeout)))

    result = wait(wt)
    if result == :timeout
        stop_watching(pfw)
    end
    result == :poll
end

watch_file(s; poll=false) = watch_file(false, s, poll=poll)
function watch_file(cb, s; poll=false)
    if poll
        pfw = PollingFileWatcher(cb,s)
        start_watching(pfw)
        return pfw
    else
        return FileMonitor(cb,s)
    end
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
using Parameters
import Base: length

abstract type AbstractPhase end
export AbstractPhase

abstract type IdealPhase <: AbstractPhase end
export IdealPhase

struct EmptyPhase <: AbstractPhase end
export EmptyPhase

@with_kw struct IdealGas{Q<:AbstractReaction} <: IdealPhase
    name::String = ""
    species::Array{Species,1}
    reactions::Array{Q,1}
    spcdict::Dict{String,Int64}
    diffusionlimited::Bool = false
end
IdealGas(species,reactions; name="",diffusionlimited=false) = IdealGas(species=species,reactions=reactions,name=name,
diffusionlimited=diffusionlimited,spcdict=Dict([sp.name=>sp.index for sp in species]))
export IdealGas

@with_kw struct IdealDiluteSolution{Q<:AbstractReaction} <: IdealPhase
    name::String = ""
    species::Array{Species,1}
    reactions::Array{Q,1}
    solvent::Solvent
    spcdict::Dict{String,Int64}
    diffusionlimited::Bool = true
end
IdealDiluteSolution(species,reactions,solvent; name="",diffusionlimited=true) = IdealDiluteSolution(species=species,reactions=reactions,
solvent=solvent,name=name,diffusionlimited=diffusionlimited,spcdict=Dict([sp.name=>sp.index for sp in species]))
export IdealDiluteSolution

@with_kw struct HomogeneousCatalyst{Q<:AbstractReaction} <: AbstractPhase
    name::String = ""
    species::Array{Species,1}
    reactions::Array{Q,1}
    spcdict::Dict{String,Int64}
end
export HomogeneousCatalyst

length(p::T) where {T<:AbstractPhase} = 1
export length

iterate(p::T) where {T<:AbstractPhase} = p
export iterate

Broadcast.broadcastable(p::T) where {T<:AbstractPhase} = Ref(p)
export broadcastable


text:
function open_externally(res_dir)
    MyTools.run_command("open $(res_dir)");
end

max_stars_repo_name: josePereiro/MyResultsManager

max_issues_repo_name: josePereiro/MyResultsManager

max_forks_repo_name: josePereiro/MyResultsManager


text:
############################################################
## get methods
## - reorganise result_type here?

function Base.get(
    parser::W, sequence, till,
    after, i, state
) where {W <: WrappedParser}
    get(parser.parser, sequence, till, after, i, state)
end

function Base.get(
    parser::W, sequence, till,
    after, i, state
) where {W <: LookAround}
    after_ = nextind(sequence,i,parser.parser,state)
    get(parser.parser, sequence, till,
        after_, i, state)
end

"""
    Base.get(parser::Bytes{T}, sequence::Vector{UInt8})

!!! note
    The endian type is not configurable.  Might depend on the arch.
    Your hint and help would be much appreciated!
"""
function Base.get(parser::Bytes{T}, sequence::Vector{UInt8}, till,
                  after, i, state) where {T}
    if isbitstype(T)
        reinterpret(T,sequence[i:after-1])[1]
    else
        T(sequence[i:after-1])
    end
end

function Base.get(parser::AbstractToken,
                  sequence, till,
                  after, i, state)
    parser isa CombinedParser && error("define Base.get(parser::$(typeof(parser)), sequence, till, after, i, state)")
    state.state
end

Base.get(
    parser::Union{Always,NegativeLookahead,AtStart,AtEnd,NegativeLookbehind},
    sequence, till, after,
    i, state) =
    parser

"""
    Base.get(parser::AbstractToken{Nothing}, sequence, till, after, i, state)

Default method for parser types returning nothing
"""
Base.get(::AbstractToken{<:Nothing},
         sequence, till, after,
         i, state) =
             nothing

function Base.get(t::PositiveLookbehind,
                  str, till,
                  after, i, state)
    rseq = revert(str)
    i_ = reverse_index(rseq,prevind(rseq,i))
    after_ = nextind(rseq,i_,t.parser,state)
    get(t.parser, rseq, till,
        after_, i_, state)
end

export JoinSubstring
"""
    JoinSubstring(x)
    (!)(x::AbstractToken)

Parser Transformation getting the matched SubString.
"""
@auto_hash_equals struct JoinSubstring{P,S} <: WrappedParser{P,S,SubString}
    parser::P
    JoinSubstring(x) =
        new{typeof(x),state_type(x)}(x)
end
Base.map(f::Type{<:JoinSubstring}, p::AbstractToken) = JoinSubstring(p)
revert(x::JoinSubstring) = JoinSubstring(revert(x.parser))
function print_constructor(io::IO,x::JoinSubstring)
    print_constructor(io,x.parser)
    printstyled(io," |> !", color=:blue)
end

deepmap_parser(f::Function,mem::AbstractDict,x::JoinSubstring,a...;kw...) =
    get!(mem,x) do
        JoinSubstring(
            deepmap_parser(f,mem,x.parser,a...;kw...))
    end

export map_match
map_match(f::Function,p_) =
    map(f, JoinSubstring(parser(p_)))

function Base.get(
    x::Union{AbstractString,JoinSubstring,ConstantParser{<:Any,<:AbstractString}},
    sequence, till,
    after, i, state)
    li = prevind(sequence,after)
    li<i ? "" : @inbounds SubString(sequence,i,li)
end

function Base.get(
    parser::Union{Char,NIndexParser{1,Char},ConstantParser{1,Char}},
    sequence, till,
    after, i, state)
    @inbounds sequence[i]
end

# Base.get(parser::ConstantParser, sequence, till, after, i, state) = parser.parser

function Base.get(parser::JoinSubstring{<:NamedParser}, sequence, till, after, i, state)
    parser.parser.name => get(parser.parser.parser, sequence, till, after, i, state)

end

function Base.get(parser::FlatMap, sequence, till, after, i, state)
    li = nextind(sequence,i,parser.left,tuple_pos(state))
    get(right_parser(state),
        sequence, till,
        after, li,
        right_state(state))
end

function Base.get(parser::Either,
                  sequence, till,
                  after, i, state)
    get(parser.options[either_state_option(state)],
        sequence, till,
        after, i, either_state_state(state))
end

function Base.get(parser::Either{<:Trie},
                  sequence, till,
                  after, i, state::NCodeunitsState)
    get(state.state)
end

function Base.get(parser::Optional,
                  sequence, till,
                  after, i, state)
    if state === None()
        parser.default
    else
        get(parser.parser,sequence, till, after, i, state)
    end
end

function Base.get(parser::Repeat,
                  sequence, till,
                  after, i, state::Vector)
    r = Vector{result_type(parser.parser)}(undef,length(state))
    i_=i
    for (p,s) in enumerate(state)
        after_=nextind(sequence,i_,parser.parser,s)
        @inbounds r[p] = get(parser.parser,sequence, till, after_, i_, s)
        i_=after_
    end
    r
end

function Base.get(parser::Repeat,
                  sequence, till,
                  after, i, state::Int)
    r = Vector{result_type(parser.parser)}(undef,state)
    i_=i
    s=MatchState()
    for p in 1:state
        after_=nextind(sequence,i_,parser.parser,s)
        @inbounds r[p] = get(parser.parser,sequence, till, after_, i_, s)
        i_=after_
    end
    r
end

# Base.get(parser::Sequence, sequence, till, after, i, state::MatchState) =
#     get(parser, sequence, till, after, i, ( MatchState() for i in 1:length(parser.parts)) )
# function Base.get(parser::Sequence, sequence, till::Int, after::Int, i::Int, state)
#     r = Vector{Any}(undef,length(parser.parts))
#     i_::Int = i
#     for (p,s) in enumerate(state)
#         after_=nextind(sequence,i_,parser.parts[p],s)
#         r[p] = get(parser.parts[p],sequence, till, after_, i_, s)
#         i_=after_
#     end
#     1
#     tuple(r...)
# end
@generated function get(parser::Sequence{pts,sts}, sequence, till::Int, after::Int, posi::Int, states) where {pts,sts}
    fpts = fieldtypes(pts)
    spts = Type[ Union{Nothing,state_type(t)} for t in fpts ]
    n = length(fpts)
    subresult = Symbol[ gensym(:r) for i in fpts ]
    part = Symbol[ gensym(:part) for i in fpts ]
    pposi = Symbol[ gensym(:pos) for i in 1:(n+1) ]
    substate = Symbol[ gensym(:s) for i in fpts ]
    init = if states<:MatchState
        [
            quote
            $(substate[i])::MatchState = MatchState()
            @inbounds $(part[i])::$p = parser.parts[$i]
            $(pposi[i])::Int = 0
            end
            for (i,(t,p)) in enumerate(zip(spts,fpts))
        ]
    else
        [
            quote
            $(substate[i])::$t = states[$i]
            @inbounds $(part[i])::$p = parser.parts[$i]
            $(pposi[i])::Int = 0
            end
            for (i,(t,p)) in enumerate(zip(spts,fpts))
        ]
    end
    parseparts = [
        quote
        $(pposi[i+1]) = nextind(sequence,$(pposi[i]),$(part[i]),$(substate[i]))
        $(subresult[i]) = get($(part[i]),sequence, till, $(pposi[i+1]), $(pposi[i]), $(substate[i]))
        end
        for (i,t) in enumerate(fpts)
    ]
    R = quote
        $(pposi[end])::Int = after
        $(init...)
        $(pposi[1]) = posi
        $(parseparts...)
        tuple($(subresult...))
    end
    R
end

## Base.get(x::Palimdrome, str, state) = SubString(str,state...)
## TODO: dispatch Base.get(x::CombinedParser, str, till::Int, after::Int, posi::Int, state) = Base.get(x, str, after, posi, state)
## TODO: dispatch Base.get(x::CombinedParser, str, after::Int, posi::Int, state) = Base.get(x, str, posi, state)
## TODO: dispatch Base.get(x::CombinedParser, str, posi::Int, state) = Base.get(x, str, state)

export Transformation
"""
    Transformation{T}(f::Function, p_) where {T}
    Base.map(f::Function, Tc::Type, p::ParserTypes, a...)
    Base.map(f::Function, p::ParserTypes, a...)

Parser transforming result of a wrapped parser.
`a...` is passed as additional arguments to `f` (at front .
"""
@auto_hash_equals struct Transformation{F,P,S,T} <: WrappedParser{P,S,T}
    transform::F
    parser::P
    Transformation(T::Type, p_) =
        let p = parser(p_)
            new{Type,typeof(p),state_type(p),T}(T, p)
        end
    Transformation{T}(transform, p_) where {T} =
        let p = parser(p_)
            new{typeof(transform),typeof(p),state_type(p),T}(transform, p)
        end
    Transformation{T}(transform, p_::NamedParser) where {T} =
        let p = p_.parser
            tp = new{typeof(transform),typeof(p),state_type(p),T}(transform, p)
            with_name(p_.name,tp)
        end
end
deepmap_parser(f::Function,mem::AbstractDict,x::Transformation,a...;kw...) =
    get!(mem,x) do
        Transformation{result_type(x)}(
            x.transform,
            deepmap_parser(f,mem,x.parser,a...;kw...))
    end

"""
A Transformation{<:Constant} skips evaluation of get(.parser).
"""
struct Constant{T}
    value::T
end
Base.show(io::IO, x::Constant) = show(io,x.value)
function map_constant(transform, p::AbstractToken)
    T=typeof(transform)
    Transformation{T}(Constant(transform), p)
end
"""
    parser(constant::Pair{<:ParserTypes})
    map_constant(constant, parser::AbstractToken)

A parser mapping matches of `parser` `x.first` to `constant` `x.second`.
"""
parser(constant::Pair{<:ParserTypes}) =
    map_constant(constant.second, parser(constant.first))
"""
    Base.get(parser::Transformation{<:Constant}, a...)

does not evaluate `get(parser.transform,...)`.
"""
function Base.get(parser::Transformation{<:Constant}, sequence, till, after, i, state)
    parser.transform.value
end

function _string(io::IO,x::Constant)
    print(io,"Constant(")
    show(io,x.value)
    print(io,")")
end
_string(io::IO,x::Function) = print(io,x)
_string(io::IO,x::Type) = print(io,x)
children(x::Transformation) = children(x.parser)
function print_constructor(io::IO,x::Transformation)
    print_constructor(io,x.parser)
    print(io," |> map(")
    printstyled(io,x.transform, color=:blue)
    print(io,")")
end

function print_constructor(io::IO,x::Transformation{Type, <:JoinSubstring})
    print_constructor(io,x.parser)
    print(io,"!")
end

function print_constructor(io::IO,x::Transformation{<:Constant})
    print_constructor(io,x.parser)
    printstyled(io," => ",string(x.transform), color=:blue)

end

"""
    Base.get(parser::Transformation{<:Function}, a...)
    Base.get(parser::Transformation{<:Type}, a...)

Function call `parser.transform(get(parser.parser,a...))`.
"""
function Base.get(parser::Transformation{<:Function}, sequence, till, after, i, state)
    v = get(parser.parser, sequence, till, after, i, state)
    parser.transform(v)
end
function Base.get(parser::Transformation{<:Type}, sequence, till, after, i, state)
    v = get(parser.parser, sequence, till, after, i, state)
    parser.transform(v)
end

export IndexAt
"""
Struct for fast access to an index of a `Transformation`.

```jldoctest
julia> using CombinedParsers.Regexp

julia> p = re"(?:a|b*)."[1]
🗄 Sequence[1]
├─ |🗄... Either
│  ├─ a
│  └─ b*  |> Repeat
└─ [^\\n] CharNotIn
::Union{Char, Array{Char,1}}
```

See also [`getindex`](@ref), [`Sequence`](@ref).
"""
struct IndexAt{I}
    i::I
end

"""
    Base.get(parser::Transformation{<:IndexAt}, a...)

`getindex(get(parser.parser,a...).parser.transform)`
"""
function Base.get(parser::Transformation{IndexAt{I}}, sequence, till, after, i, state) where {I <: Integer}
    v = get(parser.parser,sequence, till, after, i, state)
    v[parser.transform.i]
end
function Base.get(parser::Transformation{IndexAt{Is}}, sequence, till, after, i, state) where {Is <: Union{Tuple, Vector, UnitRange}}
    tuple(get(parser.parser,sequence, till, after, i, state)[parser.transform.i]...)
end

function print_constructor(io::IO,x::Transformation{<:IndexAt})
    print_constructor(io,x.parser)
    printstyled(io,"[",x.transform.i,"]", color=:blue)
end

function _iterate(parser::Transformation, sequence, till, posi, next_i, state)
    r = _iterate(parser.parser, sequence, till, posi, next_i, state )
end

"""
    infer_result_type(f::Function,Tc::Type,p::ParserTypes,onerror::AbstractString,ts::Type...; throw_empty_union=true)

Used by Parser Transformations to infer result type of a parser.
Throws error if type inference fails, if throw_empty_union=true.
"""
function infer_result_type(f::Function,Tc::Type,p::ParserTypes,onerror::AbstractString,ts::Type...; throw_empty_union=true)
    Ts = Base.return_types(f, tuple(result_type(p),ts...))
    isempty(Ts) && error("transformation type signature mismatch $f$(tuple(result_type(p),ts...))::$Ts<:$Tc")
    ( length(Ts) > 1 || Any <: first(Ts) ) && return Tc ##error(onerror*"  $f$(tuple(result_type(p),ts...))::$Ts<:$Tc")
    T = first(Ts)
    if throw_empty_union && T <: Union{}
        error("transformation type signature mismatch $f$(tuple(result_type(p),ts...))::$Ts<:$Tc")
    elseif T <: Tc
        T
    else
        @warn "type mismatch $f$(tuple(result_type(p),ts...))::$T<:$Tc"
        Tc
    end
end

# export Index
# Index = map_at((v,i) -> i, Always())

import Base: map
"""
    map(f::Function, p::AbstractToken, a...)

Parser matching `p`, transforming parsing results (`x`) with function `f(x,a...)`.

See also: [`get`](@ref), [`Transformation`](@ref)
"""
function Base.map(f::Function, p::AbstractToken, a...;
                  throw_empty_union=true)
    T = infer_result_type(f,Any,p,"call seq(function,type,parts...)",typeof.(a)...;
                          throw_empty_union=throw_empty_union)
    Transformation{T}(isempty(a) ? f : v -> f(v, a...), p)
end

"""
    map(T::Type, p::AbstractToken, a...)

Parser matching `p`, transforming `p`s parsing result with constructor `T(x,a...)`.

See also: [`get`](@ref), [`Transformation`](@ref)
"""
function Base.map(Tc::Type, p::AbstractToken, a...)
    Transformation{Tc}(isempty(a) ? Tc : v -> Tc(a..., v), p)
end

"""
    map(index::IndexAt, p::AbstractToken, a...)
    map(constant, p::AbstractToken, a...)

Parser matching `p`, transforming `p`s parsing results to `getindex(x,index)` or `constant`.

See also: [`get`](@ref), [`Transformation`](@ref)

"""
function Base.map(index::IndexAt{<:Integer}, p::AbstractToken)
    T=result_type(p)
    Transformation{fieldtype(T,index.i)}(index, p)
end
function Base.map(index::IndexAt{<:UnitRange}, p::AbstractToken)
    T=Tuple{fieldtypes(result_type(p))[index.i]...}
    Transformation{T}(index, p)
end
##instance(Tc::Type, p::ParserTypes, a...)
##    instance(Tc::Type, p::ParserTypes)
function instance(Tc::Type, p::ParserTypes, a...)
    Transformation{Tc}((v) -> Tc(a..., v), p)
end
function instance(Tc::Type, p::ParserTypes)
    Transformation(Tc, p)
end
function Base.map(f::Function, Tc::Type, p::AbstractToken, a...)
    T = infer_result_type(f,Tc,p,"call seq(function,type,parts...)",typeof.(a)...)
    Transformation{Tc}(isempty(a) ? f : v -> f(v, a...), p)
end
Base.map(f::typeof(identity), p::AbstractToken) = p

@deprecate map(T::Type, f::Function, p::AbstractToken, a...) map(f,T,p,a...)
@deprecate instance(f::Function,p,a...) map(f,parser(p),a...)


text:
Base.@kwdef mutable struct Sobol <: GSAMethod
    N::Int=1000
    order::Array{Int}=[0,1]
    nboot::Int=0
    conf_int::Float64=0.95
end

function give_rand_p!(p_range,p,p_fixed=nothing,indices=nothing)
    if p_fixed === nothing
        for j in 1:length(p_range)
            p[j] = (p_range[j][2] -p_range[j][1])*rand() + p_range[j][1]
        end
    else
        j = 1
        for i in 1:length(p_range)
            if i in indices
                p[i] = p_fixed[j]
                j += 1
            else
                p[i] = (p_range[i][2] -p_range[i][1])*rand() + p_range[i][1]
            end
        end
    end
end

function calc_mean_var(f,p_range,N)
    p = Array{eltype(p_range[1])}(undef, length(p_range))
    give_rand_p!(p_range,p)
    y1 = f(p)
    y0 = zero(y1)
    v = zero(y1)
    for i in 1:N
        give_rand_p!(p_range,p)
        y1 = f(p)
        @. y0 += y1
        @. v += y1^2
    end
    @. y0 = y0/N
    @. v = v/(N-1) - (N*y0^2)/(N-1)
    y0,v
end

function first_order_var(f,p_range,N,y0,v,p1,p2,p3)
    ys = Array{typeof(y0)}(undef,length(p_range))
    for i in 1:length(p_range)
        y = zero(y0)
        indices = [k for k in 1:length(p_range) if k != i]
        i_arr = [i]
        for j in 1:N
            give_rand_p!(p_range,p2)
            give_rand_p!(p_range,p1,@view(p2[i:i]),i_arr)
            give_rand_p!(p_range,p3,@view(p1[indices]),indices)
            y .+= (f(p2)) .* (f(p1) .- f(p3))
        end
        ys[i] = y/N
    end
    for i in 1:length(ys)
        @. ys[i] = ys[i] / v
    end
    ys
end

function second_order_var(f,p_range,N,y0,v,p1,p2,p3)
    ys = Array{typeof(y0)}(undef,Int((length(p_range)*(length(p_range)-1))/2))
    curr = 1
    for i in 1:length(p_range)
        for j in i+1:length(p_range)
            y = zero(y0)
            i_arr = [l for l in 1:length(p_range) if l != i]
            j_arr = [l for l in 1:length(p_range) if l != j]
            for k in 1:N
                give_rand_p!(p_range,p2)
                give_rand_p!(p_range,p1,@view(p2[i_arr]),i_arr)
                give_rand_p!(p_range,p3,@view(p2[j_arr]),j_arr)
                y .+= (f(p1) .- f(p3)).^2
            end
            ys[curr] = y/(2*N)
            curr += 1
        end
    end
    for i in 1:length(ys)
        ys[i] = @. ys[i] / v
    end
    ys
end

function total_var(f,p_range,N,y0,v,p1,p2,p3)
    ys = Array{typeof(y0)}(undef,length(p_range))
    for i in 1:length(p_range)
        y = zero(y0)
        indices = [k for k in 1:length(p_range) if k != i]
        for j in 1:N
            give_rand_p!(p_range,p2)
            give_rand_p!(p_range,p1,@view(p2[indices]),indices)
            y .+= (f(p2) .- f(p1)).^2
        end
        ys[i] = y/(2*N)
    end
    for i in 1:length(ys)
        @. ys[i] = ys[i] / v
    end
    ys
end

function first_total_var(f,p_range,N,y0,v,p1,p2,p3)
    ys_first = Array{typeof(y0)}(undef,length(p_range))
    ys_tot = Array{typeof(y0)}(undef,length(p_range))
    for i in 1:length(p_range)
        y_first = zero(y0)
        y_tot = zero(y0)
        indices = [k for k in 1:length(p_range) if k != i]
        i_arr = [i]
        for j in 1:N
            give_rand_p!(p_range,p2)
            give_rand_p!(p_range,p1,@view(p2[i:i]),i_arr)
            give_rand_p!(p_range,p3,@view(p1[indices]),indices)
            f_p1 = f(p1)
            f_p3 = f(p3)
            y_first .+=  (f(p2)) .* (f_p1 .- f_p3)
            y_tot .+=  (f_p1 .- f_p3).^2
        end
        ys_first[i] = y_first/N
        ys_tot[i] = y_tot/(2*N)
    end
    for i in 1:length(ys_tot)
        @. ys_first[i] = ys_first[i] / v
        @. ys_tot[i] = ys_tot[i] / v
    end
    [ys_first,ys_tot]
end

mutable struct SobolResult{T1,T2}
    S1::T1
    S1_Conf_Int::T2
    S2::T1
    S2_Conf_Int::T2
    ST::T1
    ST_Conf_Int::T2
end

function calc_ci(f,p_range,N,y0,v,nboot,conf_int,sa_func,p1,p2,p3)
    conf_int_samples = [sa_func(f,p_range,N,y0,v,p1,p2,p3) for i in 1:nboot]
    elems_ = []
    for i in 1:length(conf_int_samples[1])
        elems = []
        for k in 1:length(conf_int_samples[1][1])
            elem = eltype(conf_int_samples[1][1])[]
            for j in 1:length(conf_int_samples)
                push!(elem,conf_int_samples[j][i][k])
            end
            push!(elems,elem)
        end
        push!(elems_,elems)
    end
    z = -quantile(Normal(), (1-conf_int)/2)
    S1_Conf_Int = [[z*std(sample) for sample in el] for el in elems_]
end

function gsa(f,p_range::AbstractVector,method::Sobol)
    @unpack N, order, nboot, conf_int = method
    y0,v = calc_mean_var(f,p_range,N)
    p2 = Array{eltype(p_range[1])}(undef, length(p_range))
    p1 = Array{eltype(p_range[1])}(undef, length(p_range))
    p3 = Array{eltype(p_range[1])}(undef, length(p_range))
    sobol_sens = SobolResult(Array{T where T}[],Array{Array{T where T},1}[],Array{T where T}[],Array{Array{T where T},1}[],Array{T where T}[],Array{Array{T where T},1}[])
    if 0 in order && 1 in order
        first_total = first_total_var(f,p_range,N,y0,v,p1,p2,p3)
        sobol_sens.S1 = first_total[1]
        sobol_sens.ST = first_total[2]
        if nboot > 0
            ci = calc_ci(f,p_range,N,y0,v,nboot,conf_int,first_total_var,p1,p2,p3)
            sobol_sens.S1_Conf_Int = [vec.(first_total[1]) - ci[1], vec.(first_total[1]) + ci[1]]
            sobol_sens.ST_Conf_Int = [vec.(first_total[2]) - ci[2], vec.(first_total[2]) + ci[2]]
        end
    elseif 1 in order
        first_order = first_order_var(f,p_range,N,y0,v,p1,p2,p3)
        sobol_sens.S1 = first_order
        if nboot > 0
            ci = calc_ci(f,p_range,N,y0,v,nboot,conf_int,first_order_var,p1,p2,p3)
            sobol_sens.S1_Conf_Int = [vec.(first_order) - ci, vec.(first_order) + ci]
        end
    elseif 0 in order
        total_indices = total_var(f,p_range,N,y0,v,p1,p2,p3)
        sobol_sens.ST = total_indices
        if nboot > 0
            ci = calc_ci(f,p_range,N,y0,v,nboot,conf_int,total_var)
            sobol_sens.ST_Conf_Int = [vec.(total_indices) - ci, vec.(total_indices) + ci]
        end
    end
    if 2 in order
        second_order = second_order_var(f,p_range,N,y0,v,p1,p2,p3)
        sobol_sens.S2 = second_order
        if nboot > 0
            ci = calc_ci(f,p_range,N,y0,v,nboot,conf_int,second_order_var,p1,p2,p3)
            sobol_sens.S2_Conf_Int = [vec.(second_order) - ci, vec.(second_order) + ci]
        end
    end
    sobol_sens
end

function gsa(prob::DiffEqBase.DEProblem,alg::DiffEqBase.DEAlgorithm,t,p_range::AbstractVector,method::Sobol)
    f = function (p)
        prob1 = remake(prob;p=p)
        Array(solve(prob1,alg;saveat=t))
    end
    @assert length(prob.p) == length(p_range)
    gsa(f,p_range,method)
end


text:
# Primal infeasible test problems
# Problems with m1 feasible equality constraints, one infeasible m2-sized second-order-cone constraint and one feasible psd constraint
# constrain s[m1] = t of [t;x] to be -1, which makes the second-order-cone constraint infeasible

# This test is precision agnostic and will be run with T precision
if !@isdefined(UnitTestFloats)
    UnitTestFloats = [Float64] #if not run in full test setup, just do it for one float type
end

for T in UnitTestFloats

  # COSMO currently doesn't support PSD constraints with BigFloat precision
  if precision(T) >= precision(Float64) && T != BigFloat

    @testset "Primal infeasible QP problems - Testset 3" begin
      # for iii = 1:nn

        # nn = 1
        rng = Random.MersenneTwister(1313)

        # choose size of problem
        n = rand(rng,10:50)
        m1 = rand(rng,2:10) #equality constraint
        m2 = rand(rng,3:10) # second order cone constraint
        r = rand(rng,4:10)
        m3 = r^2 # one sdp cone constraint of size r
        m = m1+m2+m3

        A = sprand(rng,T, m,n,0.8)*T(50)
        xtrue = rand(rng,T, n,1)*T(50)
        s1 = zeros(T, m1,1)
        s2 = rand(rng, T, (m2, 1))
        s3 = vec(generate_pos_def_matrix(rng, r, MT = T))
        s = [s1;s2;s3]
        b = A*xtrue+s

        # enforce the t in ||x|| < t of second order cone constraint to be t < 0 --> infeasible
        A[m1+1,:] = zeros(T, n)
        b[m1+1] = -one(T)
        b = vec(b)

        # create dual feasible problem Px+q+A'y = 0, and y in K*
        P = generate_pos_def_matrix(rng, n, MT = T)
        ytrue_1 = rand(rng, T, m1,1)*T(50)
        ytrue_2 = rand(rng, T, m2-1,1)*T(50)
        ytrue_2 = [norm(ytrue_2,2)+one(T);ytrue_2]
        Ytrue_3 = generate_pos_def_matrix(rng, r, T(0.1), T(5.), MT = T)
        ytrue_3 = vec(Ytrue_3)
        ytrue = [ytrue_1;ytrue_2;ytrue_3]
        q = (-P*xtrue - A'*ytrue)[:]

        A1 = -A[1:m1,:]
        A2 = -A[m1+1:m1+m2,:]
        A3 = -A[m1+m2+1:end,:]
        b1 = b[1:m1,:]
        b2 = b[m1+1:m1+m2]
        b3 = b[m1+m2+1:end]

        cs1 = COSMO.Constraint(A1,b1,COSMO.ZeroSet)
        cs2 = COSMO.Constraint(A2,b2,COSMO.SecondOrderCone)
        cs3 = COSMO.Constraint(A3,b3,COSMO.PsdCone)

        settings = COSMO.Settings{T}(max_iter=10000,eps_abs = 1e-5,eps_rel=1e-5)

        model = COSMO.Model{T}()
        assemble!(model,P,q,[cs1;cs2;cs3], settings = settings)
        res = COSMO.optimize!(model);

        @test res.status == :Primal_infeasible
      end
    # end
  end
end
nothing


text:
# using MatrixCompletion

const TEST_CHAINED_ALM_SMALL_RANDOM = true
const TEST_ONESHOT_ALM_SMALL_RANDOM = false

import LinearAlgebra: norm
import Random

function run_chained_alm_randomized(;m = nothing, n = nothing, k = nothing,
                                    impute_round = 10, sample_rate = nothing, data = nothing, sampled_data = nothing,
                                    prox_params = nothing,
                                    init_X = nothing, init_Y = nothing)
    local truth_matrix
    if isnothing(data)
        truth_matrix = randn(m, k) * randn(k, n)
    else
        truth_matrix = data
        m, n = Base.size(data)
    end
    if isnothing(sample_rate) && !isnothing(sampled_data)
        input_matrix = sampled_data
    else
        sample_model = provide(Sampler{BernoulliModel}(), rate = sample_rate / 100)
        input_matrix = sample_model.draw(truth_matrix)
    end
    manual_type_matrix = Array{Symbol}(undef, m, n)
    manual_type_matrix .= :Gaussian
    imputed, X, Y, tracker = complete(ChainedALM(),
                                      A = input_matrix,
                                      type_assignment = manual_type_matrix,
                                      block_size = Int64(500 * 500 / 10),
                                      rx = MatrixCompletion.LowRankModels.QuadReg(0),
                                      ry = MatrixCompletion.LowRankModels.QuadReg(0),
                                      target_rank = k,
                                      imputation_round = impute_round,
                                      initialX = init_X,
                                      initialY = init_Y,
                                      proximal_params = prox_params)

    return truth_matrix, imputed, X, Y, tracker
end

function run_one_shot_alm_randomized(;m = nothing, n = nothing, k = nothing,
                                     sample_rate = nothing, data = nothing, prox_params = nothing,
                                     sampled_data = nothing,
                                     init_X = nothing, init_Y = nothing)
    local truth_matrix, input_matrix
    if isnothing(data)
        truth_matrix = randn(m, k) * randn(k, n)
    else
        truth_matrix = data
        m, n = Base.size(data)
    end
    if isnothing(sample_rate) && !isnothing(sampled_data)
        input_matrix = sampled_data
    else
        sample_model = provide(Sampler{BernoulliModel}(), rate = sample_rate / 100)
        input_matrix = sample_model.draw(truth_matrix)
    end
    manual_type_matrix = Array{Symbol}(undef, m, n)
    manual_type_matrix .= :Gaussian
    imputed, X, Y, tracker = complete(OneShotALM(),
                                      A = input_matrix,
                                      type_assignment = manual_type_matrix,
                                      target_rank = k,
                                      initialX = init_X,
                                      initialY = init_Y,
                                      proximal_params = prox_params)
    return truth_matrix, imputed, X, Y, tracker
end

function get_diagnostic(A, A_imputed, X, Y, tracker)
    ret = Dict{Symbol, Any}()
    ret[:L2_total_error] = norm(A[tracker[:Missing][:Total]] - A_imputed[tracker[:Missing][:Total]]) ^ 2
    ret[:L2_relative_error] = ret[:L2_total_error] / norm(A[tracker[:Missing][:Total]]) ^ 2
    ret[:MissingEntries] = tracker[:Missing][:Total]
    ret[:Truth] = A
    ret[:Imputed] = Base.convert(Array{Float64, 2}, A_imputed)
    ret[:X] = X
    ret[:Y] = Y
    return ret
end

# @testset "$(format("Algorithm: OneShotALM[Randomized, Small]"))" begin
#     for i in 1:1
#         @test get_diagnostic(run_one_shot_alm_randomized(m = 200, n = 200, k = 10, sample_rate = 80)...)[:L2_relative_error] < 0.05
#     end
# end

# @testset "$(format("Algorithm: ChainedALM[Randomized, Small]"))" begin
#     for i in 1:1
#        @test get_diagnostic(run_chained_alm_randomized(m = 200, n = 200, k = 10, sample_rate = 80)...)[:L2_relative_error] < 0.05
#     end
# end

let
    Random.seed!(65536)
    m = 500
    n = 500
    for k in collect(10:10:500)
        for sample_rate in collect(10:1:99)
            RESULTS_DIR  = GLOBAL_SIMULATION_RESULTS_DIR *
                "random_continuous/small_500x500/" *
                "rank" * string(k) * "/"  *
                "sample" * string(sample_rate) * "/"
            DATA_FILE_NAME_ONESHOT = "oneshot_saved_variables.h5"
            DATA_FILE_PATH_ONESHOT = RESULTS_DIR * DATA_FILE_NAME_ONESHOT
            DATA_FILE_NAME_CHAINED = "chained_saved_variables.h5"
            DATA_FILE_PATH_CHAINED = RESULTS_DIR * DATA_FILE_NAME_CHAINED
            Base.Filesystem.mkpath(RESULTS_DIR)

            truth_matrix = randn(m, k) * randn(k, n)
            sample_rate = 80
            sample_model = provide(Sampler{BernoulliModel}(), rate = sample_rate / 100)
            input_matrix = sample_model.draw(truth_matrix)
            initX = randn(k, m)
            initY = randn(k, n)
            param_oneshot = ProxGradParams(max_iter = 200)
            param_chained = ProxGradParams(max_iter = 200)

            result_oneshot_alm = get_diagnostic(run_one_shot_alm_randomized(data = deepcopy(truth_matrix),
                                                                            sampled_data = deepcopy(input_matrix),
                                                                            k = k,
                                                                            init_X = deepcopy(initX),
                                                                            init_Y = deepcopy(initY),
                                                                            prox_params = param_oneshot)...)

            result_chained_alm = get_diagnostic(run_chained_alm_randomized(data = deepcopy(truth_matrix),
                                                                           sampled_data = deepcopy(input_matrix),
                                                                           impute_round = 5,
                                                                           init_X = deepcopy(initX),
                                                                           init_Y = deepcopy(initY),
                                                                           k = k,
                                                                           prox_params = param_chained)...)
            pickle(DATA_FILE_PATH_ONESHOT, result_oneshot_alm)
            pickle(DATA_FILE_PATH_CHAINED, result_chained_alm)
        end
    end
end


text:
import LinearAlgebra
import FileIO
import GeometryBasics
import Observables
import Makie
import GLMakie
import DataFrames
import CSV
import ColorTypes
import FixedPointNumbers
import LinearAlgebra
import DifferentialEquations

using Porta

startframe = 85
frames = 3600
FPS = 60
resolution = (3840, 2160)
segments = 36
exportmode = ["gif", "frames", "video"]
exportmode = exportmode[2]
modelname = "spin3_alpha1"

# The scene object that contains other visual objects
#Makie.reasonable_resolution() = (800, 800)
scene = Makie.Scene(backgroundcolor = :white,
                    show_axis = false,
                    resolution = resolution)

geocoordinates = Dict()
geocoordinates["Antarctica"] = deg2rad.([135.00; 82.86])
geocoordinates["Japan"] = deg2rad.([138.25; -36.20])
geocoordinates["Italy"] = deg2rad.([12.56; -41.87])
geocoordinates["Germany"] = deg2rad.([10.45; -51.16])
geocoordinates["Canada"] = deg2rad.([-106.34; -56.13])
geocoordinates["United States of America"] = deg2rad.([-95.71; -37.09])
geocoordinates["United Kingdom"] = deg2rad.([-3.43; -55.37])
geocoordinates["France"] = deg2rad.([2.21; -46.22])

"""
    Represents a line segment.

fields: p₁ and p₂.
"""
struct Line
    p₁::ℝ³
    p₂::ℝ³
end

function getdistance(point::ℝ³, line::Line)
    p₀ = point
    p₁, p₂ = line.p₁, line.p₂
    norm(cross(p₂ - p₁, p₁ - p₀)) / norm(p₂ - p₁)
end

function decimate(points::Array{Geographic,1}, ϵ::Float64)
    # Find the point with the maximum distance
    dmax = 0
    index = 1
    number = length(points)
    for i in 2:number-1
        line = Line(ℝ³(Cartesian(points[begin])), ℝ³(Cartesian(points[end])))
        point = ℝ³(Cartesian(points[i]))
        d = getdistance(point, line)
        if d > dmax
            index = i
            dmax = d
        end
    end

    array = []

    # If max distance is greater than epsilon, recursively simplify
    if dmax > ϵ
        # Recursive call
        array1 = decimate(points[begin:index], ϵ)
        array2 = decimate(points[index:end], ϵ)

        # Build the result list
        array = [array1[begin:end-1]; array2]
    else
        array = [points[begin]; points[end]]
    end
    # Return the result
    array
end

function getcolor(name::String, points::Array{Geographic,1}, color::Array{ColorTypes.RGBA{FixedPointNumbers.Normed{UInt8, 8}}, 2}, α::Float64)
    height, width = size(color)
    margin = 25
    if name in keys(geocoordinates)
        ϕ, θ = geocoordinates[name]
    else
        ϕ = map(x -> x.ϕ, points)
        θ = map(x -> x.θ, points)
        number = length(points)
        ϕ = sum(ϕ) / number
        θ = sum(θ) / number
    end
    ϕ, θ = (ϕ + π) / 2π, (θ + π / 2) / π
    x = Int(floor(ϕ * (width - 1))) + 1
    y = Int(floor(θ * (height - 1))) + 1
    colors = Dict()
    for i in -margin:margin
        for j in -margin:margin
            c = color[y + i, x + j]
            if isapprox(c.r, 0) && isapprox(c.g, 0) && isapprox(c.b, 0) continue end
            colors[c] = get(colors, c, 0) + 1
        end
    end
    array = []
    for (key, value) in colors
        push!(array, (value, key))
    end
    lt(x, y) = x[1] < y[1]
    array = sort(array, lt = lt, rev = true)
    c = array[begin][end]
    r, g, b, a = c.r, c.g, c.b, c.alpha
    Makie.RGBAf(r, g, b, α)
end

A(p::ComplexPlane) = begin
    z₀, z₁ = vec(p)
    r = vec(Quaternion(p)) # radial
    f = [-r[2]; r[1]; -r[4]; r[3]] # fiber tangent
    X₀, X₁ = vec(ComplexPlane(Quaternion(f)))
    U1(1 / 2 * (conj(z₀) * X₀ - z₀ * conj(X₀) + conj(z₁) * X₁ - z₁ * conj(X₁)))
end

# Use QGIS to design a geo map
colorref = FileIO.load("data/gdp/map_color_reference.png")
color = FileIO.load("data/gdp/map_transparent.png")
# Made with Natural Earth.
# Free vector and raster map data @ naturalearthdata.com.

# 1. Download cultural data admin 0 from natural earth data
# 2. Install qgis desktop
# 3. Install mmqgis plugin for exporting to CSV

# The path to the dataset
attributespath = "data/gdp/geometry-attributes.csv"
attributes = DataFrames.DataFrame(CSV.File(attributespath))
attributes = DataFrames.sort(attributes, :shapeid, rev = true)
nodespath = "data/gdp/geometry-nodes.csv"
nodes = DataFrames.DataFrame(CSV.File(nodespath))

attributesgroup = DataFrames.groupby(attributes, :NAME)
nodesgroup = DataFrames.groupby(nodes, :shapeid)
number = length(attributesgroup)
ϵ = 1e-3
countries = Dict("shapeid" => [], "name" => [], "gdpmd" => [],
                 "gdpyear" => [], "economy" => [], "partid" => [], "nodes" => [])
selection = collect(keys(geocoordinates))
for i in 1:number
    shapeid = attributesgroup[i][!, :shapeid][1]
    name = attributesgroup[i][!, :NAME][1]
    if !(name in selection) continue end
    gdpmd = attributesgroup[i][!, :GDP_MD][1]
    gdpyear = attributesgroup[i][!, :GDP_YEAR][1]
    economy = attributesgroup[i][!, :ECONOMY][1]
    subdataframe = nodes[nodes.shapeid .== shapeid, :]
    uniquepartid = unique(subdataframe[!, :partid])
    # if economy == "1. Developed region: G7"
    #     push!(selection, name)
    # end
    histogram = Dict()
    for id in uniquepartid
        sub = subdataframe[subdataframe.partid .== id, :]
        ϕ = sub.x ./ 180 .* π
        θ = sub.y ./ 180 .* π
        coordinates = map(x -> Geographic(1, x[1], x[2]), eachrow([ϕ θ]))[begin:end-1]
        coordinates = decimate(coordinates, ϵ)
        histogram[id] = length(coordinates)
    end
    partsnumber = max(values(histogram)...)
    index = findfirst(x -> histogram[x] == partsnumber, uniquepartid)
    partid = uniquepartid[index]
    subdataframe = subdataframe[subdataframe.partid .== partid, :]
    ϕ = subdataframe.x ./ 180 .* π
    θ = subdataframe.y ./ 180 .* π
    coordinates = map(x -> Geographic(1, x[1], x[2]), eachrow([ϕ θ]))[begin:end-1]
    println("Length of points: $name : $(length(coordinates))")
    coordinates = decimate(coordinates, ϵ)
    println("Length of points: $name : $(length(coordinates))")
    push!(countries["shapeid"], shapeid)
    push!(countries["name"], name)
    push!(countries["gdpmd"], gdpmd)
    push!(countries["gdpyear"], gdpyear)
    push!(countries["economy"], economy)
    push!(countries["partid"], partid)
    push!(countries["nodes"], coordinates)
end

for i in 1:length(countries["nodes"])
    println(length(countries["nodes"][i]))
end

countries["gdpmd"] = countries["gdpmd"] ./ (max(countries["gdpmd"]...) - min(countries["gdpmd"]...))

configuration = Biquaternion(ℝ³(0, 0, 0))
solidwhirls = []
ghostwhirls = []
number = length(selection)
for i in 1:number
    name = countries["name"][i]
    points = countries["nodes"][i]
    gdpmd = countries["gdpmd"][i]
    α₁ = name == "Antarctica" ? 0.1 : 0.7
    α₂ = name == "Antarctica" ? 0.05 : 0.3
    solidcolor = getcolor(name, points, colorref, α₁)
    ghostcolor = getcolor(name, points, colorref, α₂)
    points = map(x -> σmap(x), points)
    solidgauge1 = [U1(0) for j in 1:length(points)]
    solidgauge2 = [U1(gdpmd * 2π) for j in 1:length(points)]
    solidgauge1 = convert(Array{U1,1}, solidgauge1)
    solidgauge2 = convert(Array{U1,1}, solidgauge2)
    ghostgauge1 = solidgauge2
    ghostgauge2 = solidgauge2
    solidwhirl = Whirl(scene,
                       points,
                       solidgauge1,
                       solidgauge2,
                       configuration = configuration,
                       segments = segments,
                       color = solidcolor,
                       transparency = true)
    push!(solidwhirls, solidwhirl)
    ghostwhirl = Whirl(scene,
                       points,
                       ghostgauge1,
                       ghostgauge2,
                       configuration = configuration,
                       segments = segments,
                       color = ghostcolor,
                       transparency = true)
    push!(ghostwhirls, ghostwhirl)
end

framesprites = []
framesprite1 = Frame(scene,
                     σmap,
                     color,
                     configuration = configuration,
                     segments = 3segments,
                     transparency = true)
framesprite2 = Frame(scene,
                     σmap,
                     color,
                     configuration = configuration,
                     segments = 3segments,
                     transparency = true)
push!(framesprites, framesprite1)
push!(framesprites, framesprite2)

group2algebra(A::SU2, r::Float64) = begin
    τ₃ = -0.5 .* [im 0; 0 -im]
    (-2 * r) .* (A.a * τ₃ * LinearAlgebra.inv(A.a))
end

exponentiate(M::Array{<:Complex,2}, t::Float64) = begin
    eigenvalues = LinearAlgebra.eigvals(M)
    eigenvectors = LinearAlgebra.eigvecs(M)
    S = eigenvectors
    S⁻¹ = LinearAlgebra.inv(S)
    e = LinearAlgebra.Diagonal(eigenvalues) .* t
    S = convert.(Complex{BigFloat}, S)
    S⁻¹ = convert.(Complex{BigFloat}, S⁻¹)
    e = convert.(Complex{BigFloat}, e)
    G = S * ℯ.^e * S⁻¹
    z₁, z₂ = G[1,1], G[1,2]
    v = LinearAlgebra.normalize([real(z₁); imag(z₁); real(z₂); imag(z₂)])
    z₁, z₂ = Complex{Float64}(v[1] + im * v[2]), Complex{Float64}(v[3] + im * v[4])
    G = [z₁ -conj(z₂); z₂ conj(z₁)]
    SU2(G)
end

L₁ = [Complex(0) Complex(1); Complex(1) Complex(0)]
L₂ = [Complex(0) Complex(-im); Complex(im) Complex(0)]
L₃ = [Complex(1) Complex(0); Complex(0) Complex(-1)]
rotation = 4π
operator = exponentiate(L₃, 0.0)
ϵ = float(π / 2)

"""
    animate(i)

Update the state of observables with the given frame number `i`.
"""
function animate(i::Int)
    step = i / frames
    τ = step * rotation
    ψ = Int(τ ÷ ϵ)
    θ = τ % ϵ

    q₁ = [exponentiate(group2algebra(operator, 1.0), ϵ) for x in 1:ψ]
    q = exponentiate(group2algebra(operator, 1.0), θ)
    for item in q₁
        q = SU2(normalize(Quaternion(item * q)))
    end

    for index in 1:number
        # gdpmd = countries["gdpmd"][index]
        points = countries["nodes"][index]
        p = map(x -> ComplexPlane(normalize(Quaternion(SU2(σmap(x)) * q))), points)
        # α = gdpmd * 2π
        g₁ = map(x -> U1(0), p)
        g₂ = A.(p)
        g₃ = map(x -> U1(2π), p)
        update(solidwhirls[index], p, g₁, g₂)
        update(ghostwhirls[index], p, g₂, g₃)
    end

    σ₁(x) = begin
        p = σmap(x)
        p = ComplexPlane(normalize(Quaternion(SU2(p) * q)))
        α = A(p)
        _g = transformg(p, U1(0), α, 3segments)[end]
        S¹action(p, _g)
    end
    σ₂(x) = begin
        p = σmap(x)
        p = ComplexPlane(normalize(Quaternion(SU2(p) * q)))
        α = A(p)
        S¹action(p, α)
    end
    update(framesprite1, σ₁)
    update(framesprite2, σ₂)
end

function saveframe(filepath, scene, resolution)
    FileIO.save(filepath,
                scene;
                resolution = resolution,
                pt_per_unit = 400.0,
                px_per_unit = 400.0)
end

n = ℝ³(0, 0, 1)
v = ℝ³(-1, -1, 0) * (pi / 2)
#v = rotate(v, Quaternion(π / 8, ℝ³(0, 0, 1)))
# v = rotate(v, Quaternion(-π / 8, ℝ³(0, 0, 1)))
# update eye position
# scene.camera.eyeposition.val
upvector = GeometryBasics.Vec3f0(vec(n)...)
eyeposition = GeometryBasics.Vec3f0(vec(v)...)
lookat = GeometryBasics.Vec3f0(0, 0, 0)
#camera = Makie.Camera3D(scene)
# Makie.update_cam!(scene, camera, eyeposition, lookat, upvector)
Makie.update_cam!(scene, eyeposition, lookat, upvector)
scene.center = false # prevent scene from recentering on display
if exportmode ∈ ["gif", "video"]
    outputextension = exportmode == "gif" ? "gif" : "mkv"
    Makie.record(scene, "gallery/$modelname.$outputextension",
                            framerate = FPS) do io
        for i in startframe:frames
            animate(i) # animate the scene
            Makie.recordframe!(io) # record a new frame
            step = (i - 1) / frames
            println("Completed step $(100step).\n")
        end
    end
elseif exportmode == "frames"
    directory = joinpath("gallery", modelname)
    !isdir(directory) && mkdir(directory)
    for i in startframe:frames
        start = time()
        animate(i) # animate the scene
        elapsed = time() - start
        sleep(0.1)
        println("Generating frame $i took $elapsed (s).")

        start = time()

        paddingnumber = length(digits(frames))
        imageid = lpad(i, paddingnumber, "0")
        imagename = "$(modelname)_$(resolution[2])p_$(FPS)fps_$imageid.jpeg"
        filepath = joinpath(directory, imagename)
        saveframe(filepath, scene, resolution)
        elapsed = time() - start
        sleep(0.1)
        println("Saving file $filepath took $elapsed (s).")

        step = (i - 1) / frames
        println("Completed step $(100step).\n")
        sleep(1)

        stitch() = begin
            part = i ÷ (frames / 3)
            exportdir = joinpath(directory, "export")
            !isdir(exportdir) && mkdir(exportdir)
            WxH = "$(resolution[1])x$(resolution[2])"
            #WxH = "1920x1080"
            commonpart = "$(modelname)_$(resolution[2])p_$(FPS)fps"
            inputname = "$(commonpart)_%0$(paddingnumber)d.jpeg"
            inputpath = joinpath(directory, inputname)
            outputname = "$(commonpart)_$(part).mp4"
            outputpath = joinpath(exportdir, outputname)
            command1 = `ffmpeg -y -f image2 -framerate $FPS -i $inputpath -s $WxH -pix_fmt yuvj420p $outputpath`
            run(command1)
        end

        if isapprox(i % (frames / 3), 0)
            sleep(1)
            stitch()
            sleep(1)
        end
        # break
    end
end

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
# Add Julia kernel to Jupyter Notebook and/or JupyterLab IDEs.
using Pkg
Pkg.add("IJulia")

# In case the kernel is not installed, another method is to use the following code.
using IJulia
installkernel("Julia")

max_stars_repo_name: a-mhamdi/Eg-Softs

max_issues_repo_name: a-mhamdi/Eg-Softs

max_forks_repo_name: a-mhamdi/Eg-Softs


text:
module Hydrogen

using LinearAlgebra
using SparseArrays

using AtomicLevels
using AngularMomentumAlgebra
import AngularMomentumAlgebra: powneg1
using AngularMomentumAlgebra.Dipoles
using HypergeometricFunctions
using Polynomials: Polynomial
using SpecialPolynomials: Laguerre
using SpecialFunctions: gamma

include("sparse_builder.jl")
include("factorial_ratio.jl")

include("constants.jl")
include("energies.jl")
include("dipoles.jl")
include("orbitals.jl")

end