text:
include("DADA\\errormodels.jl")
ρₚₒᵢₛ = StatsFuns.poispdf
get_pval(nⱼλᵢⱼ,aⱼ) = poisccdf(nⱼλᵢⱼ,aⱼ-1)

function get_norm(nⱼ,TAIL_APPROX_CUTOFF =1e-7)
    norm = (1.0 - exp(-nⱼ))
    if norm < TAIL_APPROX_CUTOFF
      norm = nⱼ - 0.5*nⱼ*nⱼ
    end
    return norm
end

get_normed_pval(nⱼ,λᵢⱼ,aⱼ) = get_pval(nⱼ*λᵢⱼ,aⱼ)/get_norm(nⱼ*λᵢⱼ) #get_norm(nⱼ*λᵢⱼ)?

function compare!(seqs,quals,counts,ind,centres,partition,pvals,n,reads,compI,compV, err,mers,hdists,E_minmax; band_size = 16,kdist_cut = 0.42,gapless = false)
    affinegap = AffineGapScoreModel(match=5,
                                mismatch=-4,
                                 gap_open=-0,
                                 gap_extend=-8)
    base_vals = Base.ImmutableDict(DNA_A=>1,DNA_C=>2,DNA_G=>3,DNA_T=>4)
    l = keys(base_vals)
    mer_dist = 0.0
    λᵢⱼ = 0.0
    nⱼ = sum(counts[partition .== ind])
    N = sum(counts)
    #Threads.@threads for i in 1:n
     for i in 1:n
            hdist = 0.0
                        if !in(i,centres)
                            mer_dist = kmer_distance(mers[:,ind],mers[:,i],length(seqs[ind]),length(seqs[i]))
                            if mer_dist < kdist_cut
                                if gapless & (kord_distance(seqs[ind], seqs[i]) == mer_dist)
                                    aln = Tuple([(seqs[ind][j],seqs[i][j]) for j in 1:length(seqs[i])])
                                else
                                    if band_size >0
                                        pa = pairalign(OverlapAlignment(), seqs[ind], seqs[i],affinegap,true,band_size,band_size)
                                    else
                                        pa = pairalign(OverlapAlignment(), seqs[ind], seqs[i],affinegap)
                                    end
                                aln = (collect(pa.aln))
                                end
                                qs = quals[i]
                                q = 0
                                λᵢⱼ = 1.0
                                for j in 1:length(aln)

                                    if in(aln[j][2],l)
                                        q+=1
                                            j_ind = base_vals[aln[j][2]]
                                            i_ind = in(aln[j][1],l) ? base_vals[aln[j][1]] : base_vals[aln[j][2]]
                                            λᵢⱼ*=  err[(i_ind-1)*4+j_ind,qs[q]+1]
                                        hdist += 1
                                    end
                                end
                            else
                                λᵢⱼ = 0.0
                                hdist = Inf
                            end
                            if λᵢⱼ * N >= E_minmax[i]
                                push!(compI[i],length(centres))
                                push!(compV[i],λᵢⱼ)
                            end
                            if λᵢⱼ * counts[ind] > E_minmax[i]
                                E_minmax[i] = λᵢⱼ * counts[ind]
                            end
                            #λ[i] = λᵢⱼ
                            hdists[i] = hdist
                        end
                    end
end

boolfunction(pvals,n,ωₐ) = any(pvals .* n .< ωₐ)

function minp(pvals,counts,centres,partition,n = length(pvals))
    minval = 1.0
    min_ind = 1
        @inbounds for i in 1:n
        if !in(i,centres)
            if pvals[i] < minval
                minval = pvals[i]
                min_ind = i
            elseif (pvals[i] == minval) & (counts[i]> counts[min_ind])
                minval = pvals[i]
                min_ind = i
            end
        end
    end
    return minval, min_ind
end

function shuffle(centres,partition,compI,compV,reads,pvals,counts,n,shuff,max_shuffle)
    change = false
    @inbounds for i in (1:n)
        if !in(i,centres)
            _,k = findmax(compV[i] .* reads[compI[i]])
            j = compI[i][k]
            prev_j = findfirst(centres .== [partition[i]])
            if j !== prev_j
                change = true
                reads[j] += counts[i]
                reads[prev_j] -= counts[i]
                partition[i] = centres[j]
            end
            if shuff == max_shuffle
                pvals[i] = get_normed_pval(reads[j],compV[i][k],counts[i])
            end
        end
    end
    return change
end

function clustering(seqs, counts, quals,mers,err,ωₐ;band_size =16,log_p = false,log_λ=false, kdist_cut = 0.42,
    gapless = false, max_shuffle = 10)

    n = length(counts)
    pvals = zeros(n)
    centres = []
    E_minmax = zeros(n)
    hdists = zeros(n)
    compI = [Vector{Int}(undef,0) for i in 1:n]
    compV = [Vector{Float64}(undef,0) for i in 1:n]
    ind = findmax(counts)[2]
    centres = [ind]
    pvals[ind] = 1.0
    partition = fill(ind,n)
    compare!(seqs,quals,counts,ind,centres,partition,pvals,n,[],compI,compV, err,mers,hdists,E_minmax,band_size = band_size,kdist_cut = Inf,gapless = gapless)
    base_vals = Base.ImmutableDict(DNA_A=>1,DNA_C=>2,DNA_G=>3,DNA_T=>4)
    l = keys(base_vals)
    reads=[sum(counts)]

    for i in 1:n
        if !in(i,centres)
            partition[i] = centres[1]
            pvals[i] = get_normed_pval(reads[1],compV[i][1],counts[i])
        end
    end
    while true

        val, ind = minp(pvals,counts,centres,partition)
        if  boolfunction(val,n,ωₐ) & !in(ind, centres)
            partition[ind] = ind
            pvals[ind] = 1.0
            partition[centres] .= centres
            push!(centres,ind)
            hdists = hcat(hdists,zeros(n))
            compare!(seqs,quals,counts,ind,centres,partition,pvals,n,reads,compI,compV,
            err,mers,view(hdists,:,length(centres)),E_minmax,
            band_size = band_size,kdist_cut = kdist_cut,gapless = gapless)
        else
            break
        end
        reads = [sum(counts[partition .== i]) for i in centres]'
        change = true
        shuff = 1
        while shuff < max_shuffle
            shuff +=1
            if change ==false
               shuff = max_shuffle
            end
            change = false
            change = shuffle(centres,partition,compI,compV,reads,pvals,counts,n,shuff,max_shuffle)
        end
    end
    reads = [sum(counts[partition .== i]) for i in centres]'
    trans_mat = get_trans(centres,seqs,quals,counts,partition,n)
    cluster_map = Vector{Int}(undef,n)
    pc = copy(partition)
    center_raws = []
    for i in 1:length(centres)
        boolmask = partition .== centres[i]
        raws = seqs[boolmask]
        abunds = counts[boolmask]
        inds = collect(1:n)[boolmask]
        push!(center_raws,raws[findmax(abunds)[2]])
        centres[i] = inds[findmax(abunds)[2]]
        partition[boolmask] .= centres[i]
    end
    partition_map = Dict([centres[i] => i for i in 1:length(centres)]...)
    for i in 1:n
        cluster_map[i] = partition_map[partition[i]]
    end
    #trans_mat = get_trans(centres,seqs,quals,counts,partition,n)

    return (df =DataFrame([center_raws,vec(reads)],[:sequence,:abundance]),trans =  trans_mat,pvals = pvals,cluster_map =cluster_map,
    new_partition = partition,old_partition = pc)
end

function get_trans(centres,seqs,quals,counts,partition,n)
    affinegap = AffineGapScoreModel(match=5,
                                mismatch=-4,
                                 gap_open=-0,
                                 gap_extend=-8)
    base_vals = Base.ImmutableDict(DNA_A=>1,DNA_C=>2,DNA_G=>3,DNA_T=>4)
    l = keys(base_vals)
    trans_mat = zeros(Int,16,41)
    for i in 1:n
        if 1==1#!in(i,centres)
            pa = pairalign(OverlapAlignment(), seqs[partition[i]], seqs[i],affinegap)
            #qs2 = quals[partition[i]]#Int.(round.(derep.quality[i]))
            qs = quals[i]#Int.(round.(derep.quality[i]))
            q = 0
            aln = (collect(pa.aln))
            for j in 1:length(aln)# need to make this start from start of seq[i] in alignment e.g. firt !== DNA_GAP
                    if in(aln[j][2],l)
                        q+=1
                        if in(aln[j][1],l)
                            i_ind = base_vals[aln[j][1]]
                            j_ind =  base_vals[aln[j][2]]
                            #i_ind = in(aln[j][1],l) ? base_vals[aln[j][1]] : base_vals[aln[j][2]]
                            #Q = !in(i,centres) ? qs[q]+1 : qs[q]
                            #trans_mat[(i_ind-1)*4+j_ind, Q] +=counts[i]
                            trans_mat[(i_ind-1)*4+j_ind, qs[q]+1] +=counts[i]
                            #trans_mat[(i_ind-1)*4+j_ind,Int(round(mean([qs[q],qs2[q]])))+1] +=counts[i]
                        end
                    end
            end
        end
    end
    return trans_mat
end
function learnErrors(dereps ::Vector{DataFrame}, max_tries = 10;nbases = 1e8, band_size = 16,ωₐ = 1e-40,kdist_cut = 0.42)

    st = time()
    base_count = 0
    dereps_for_errs = []
    mers = []
    for (j,derep) in enumerate(dereps)
        derep.sequence =LongDNASeq.(derep.sequence)
        push!(dereps_for_errs, derep)
        push!(mers,kmer_count(derep))
        base_count += sum(derep.count .* length.(derep.sequence))
        if base_count > nbases
            break
        end
    end
    println("$(base_count) bases in $(length(dereps_for_errs)) samples")
    #counts = derep.count
    err = ones(16,41)
    #err = fill(.25,16,41)
    trans_mat = zeros(16,41)
    errs = [err]
    trans_mats =Vector{Array}(undef,length(dereps_for_errs))
    for i in 2:12
        #Threads.@threads for j in 1:length(dereps_for_errs)
            for j in 1:length(dereps_for_errs)
            derep = dereps_for_errs[j]
            trans_mats[j]= clustering(derep.sequence,derep.count,derep.quality,mers[j],err,ωₐ,
             band_size = band_size,log_p = true, kdist_cut = kdist_cut).trans
        end
        trans_mat = sum(trans_mats)
        err = loess_errors(trans_mat)
        push!(errs,err)
        println("Time elapsed: $(time() - st) seconds")
       println(UnicodePlots.heatmap(errs[i] .- errs[i-1], title = "Change from previous error estimate"))
       if any([errs[i][.!isnan.(errs[i])] ≈ errs[j][.!isnan.(errs[j])] for j in 1:i-1])
            println("Successfully converged")
            break
       end
    end
    return errs,trans_mat, trans_mats
end
function learnErrors(dereps ::Vector{String};nbases = 1e8, band_size = 16,ωₐ = 1e-40,kdist_cut = 0.42)
    dereps = dereplicate.(dereps)
    return learnErrors(dereps,nbases = nbases , band_size = band_size,ωₐ = ωₐ,kdist_cut = kdist_cut)
end

function learnErrors(derep ::DataFrame;nbases = 1e8, band_size = 16,ωₐ = 1e-40,kdist_cut = 0.42)
    return learnErrors([derep],nbases = nbases , band_size = band_size,ωₐ = ωₐ,kdist_cut = kdist_cut)
end

function dada(derep,err; band_size = 16,ωₐ = 1e-40,kdist_cut = 0.42,log_λ = false, log_p = true, gapless = false)
    println(derep.sample)
    counts = derep.count
    derep.df.sequence =LongDNASeq.(derep.sequence)
    seqs = LongDNASeq.(derep.sequence)
    qs = derep.quality
    partition = Vector{Int}(undef,length(counts))
    mers = kmer_count(derep.df)
    seq_abund = []
    pvals = []
    partition=[]
    λ= []
    hdists = []
    reads = []
    trans_mat = []
        clust = clustering(seqs,counts,qs,mers,err,ωₐ,
        band_size = band_size,log_λ = log_λ, log_p = log_p, kdist_cut = kdist_cut, gapless = gapless)

    return clust
end

function mergePairs(dadaF,derepF,dadaR,derepR)

    rf = dadaF.cluster_map[derepF.read_map]
    rr = dadaR.cluster_map[derepR.read_map]
    pairsArr = hcat(rf,rr)
    unqs = unique(pairsArr, 1)

    Fseqs = dadaF.df.sequence[unqs[:,1]]
    Rseqs = dadaR.df.sequence[unqs[:,2]]
    affinegap = AffineGapScoreModel(match=5,
    mismatch=-4,
     gap_open=-64,
     gap_extend=-64)
     alignments = pairalign.(OverlapAlignment(), Fseqs,Rseqs,affinegap)
    counts = countmap(eachrow(unqs))
    abundance = [counts[unq] for unq in eachrow(unqs)]
    df = Dataframe(hcat(pairsArr,aln.(alignments),abundance),[:rev,:fwd,:aln,:abundance])
    return df
end

function pairalign(::OverlapAlignment, a::S1, b::S2, score::AffineGapScoreModel{T}, banded::Bool,lower_offset,upper_offset;
    score_only::Bool=false) where {S1,S2,T}
    m = length(a)
    n = length(b)
    if banded
        if m > n
            L = m - n + lower_offset
            U = upper_offset
        else
            L = lower_offset
            U = n - m + upper_offset
        end
        bnw = BioAlignments.BandedNeedlemanWunsch{T}(m, n, L, U)
        score = BioAlignments.run!(bnw, a, b, score.submat,
        T(0), T(0), score.gap_open, score.gap_extend, T(0), T(0),
        T(0), T(0), score.gap_open, score.gap_extend, T(0), T(0),
        )
        if score_only
            return BioAlignments.PairwiseAlignmentResult{S1,S2}(score, true)
        else
            a′ = BioAlignments.traceback(bnw, a, b, (m, n))
            return BioAlignments.PairwiseAlignmentResult(score, true, a′, b)
        end
    end
end


text:
function solve(filename)
    #load txt file
    values = readlines(filename*".txt")
    count0 = 0
    count1 = 0
    γ = ""
    ϵ = ""
    for i in 1:length(values[1])
        for j in 1:length(values)
            if values[j][i] == '0'
                count0 += 1
            else
                count1 += 1
            end
        end
        if count0 > count1
            γ = γ * "0"
            ϵ = ϵ * "1"
        else
            γ = γ * "1"
            ϵ = ϵ * "0"
        end
        count0 = 0
        count1 = 0
    end
    println(γ)
    println(parse(Int, γ, base = 2))
    γ = parse(Int, γ, base = 2)
    println(ϵ)
    println(parse(Int, ϵ, base = 2))
    ϵ = parse(Int, ϵ, base = 2)
    println("Result: "*string(γ * ϵ))
end

solve("smalltest")
solve("input")

function solve2(filename)
    #load txt file
    values = readlines(filename*".txt")
    count0 = 0
    count1 = 0
    count0_co2 = 0
    count1_co2 = 0
    currMostCommon = '0'
    currLeastCommon = '0'
    values_CO2 = copy(values)
    for i in 1:length(values[1])
        #println("i: "*string(i))
        #println("size of values: "*string(length(values)))
        #println("size of CO2vals: "*string(length(values_CO2)))
        if length(values) > 1 || length(values_CO2) > 1
            #compute most common bit
            for j in 1:length(values)
                if values[j][i] == '0'
                    count0 += 1
                else
                    count1 += 1
                end
            end
            for j in 1:length(values_CO2)
                if values_CO2[j][i] == '0'
                    count0_co2 += 1
                else
                    count1_co2 += 1
                end
            end

            #filter out all non-fitting entries
            if length(values) > 1
                if count0 > count1
                    currMostCommon = '0'
                else
                    currMostCommon = '1'
                end
                values = copy(filter!(x -> x[i] == currMostCommon, values))
                count0 = 0
                count1 = 0
                #println(values)
            end
            if length(values_CO2) > 1
                if count0_co2 > count1_co2
                    currLeastCommon = '1'
                else
                    currLeastCommon = '0'
                end
                values_CO2 = copy(filter!(x -> x[i] == currLeastCommon, values_CO2))
                count0_co2 = 0
                count1_co2 = 0
                #println(values_CO2)
            end
        else
            break
        end
    end
    println("Oxygen: "*string(parse(Int, values[1], base = 2)))
    println(values[1])
    println("CO2: "*string(parse(Int, values_CO2[1], base = 2)))
    println(values_CO2[1])
    println("Product: "*string(parse(Int, values[1], base = 2) * parse(Int, values_CO2[1], base = 2)))
end

solve2("smalltest")
solve2("input")
#solve2("addtest")

ext: jl


text:
const VALUE_TYPE_UNKNOWN = 0
const VALUE_TYPE_NULL = 1
const VALUE_TYPE_STRING = 2
const VALUE_TYPE_INTEGER = 3
const VALUE_TYPE_BOOLEAN = 4
const VALUE_TYPE_DOUBLE = 5
const VALUE_TYPE_ARRAY = 6
const VALUE_TYPE_EDGE = 7
const VALUE_TYPE_NODE = 8
const VALUE_TYPE_PATH = 9
const VALUE_TYPE_MAP = 10
const VALUE_TYPE_POINT = 11

struct VALUE_TYPE{x}
end

VALUE_TYPE(x) = VALUE_TYPE{x}()

function parsevalue(g::Graph, ::VALUE_TYPE{VALUE_TYPE_UNKNOWN}, raw_value::String) @assert false "Unknown value type" end
function parsevalue(g::Graph, ::VALUE_TYPE{VALUE_TYPE_NULL}, raw_value::Nothing) return nothing end
function parsevalue(g::Graph, ::VALUE_TYPE{VALUE_TYPE_STRING}, raw_value::String) return raw_value end
function parsevalue(g::Graph, ::VALUE_TYPE{VALUE_TYPE_INTEGER}, raw_value::Int) return raw_value end
function parsevalue(g::Graph, ::VALUE_TYPE{VALUE_TYPE_BOOLEAN}, raw_value::String) return parse(Bool, raw_value) end
function parsevalue(g::Graph, ::VALUE_TYPE{VALUE_TYPE_DOUBLE}, raw_value::String) return parse(Float64, raw_value) end
function parsevalue(g::Graph, ::VALUE_TYPE{VALUE_TYPE_ARRAY}, raw_entry::Vector{T} where T) return parsearray(g, raw_entry) end
function parsevalue(g::Graph, ::VALUE_TYPE{VALUE_TYPE_EDGE}, raw_entry::Vector{T} where T) return parseedge(g, raw_entry) end
function parsevalue(g::Graph, ::VALUE_TYPE{VALUE_TYPE_NODE}, raw_entry::Vector{T} where T) return parsenode(g, raw_entry) end
function parsevalue(g::Graph, ::VALUE_TYPE{VALUE_TYPE_PATH}, raw_entry::Vector{T} where T) return parsepath(g, raw_entry) end
function parsevalue(g::Graph, ::VALUE_TYPE{VALUE_TYPE_MAP}, raw_entry::Vector{T} where T) return parsemap(g, raw_entry) end
function parsevalue(g::Graph, ::VALUE_TYPE{VALUE_TYPE_POINT}, raw_entry::Vector{T} where T) return parsepoint(g, raw_entry) end

function parsepoint(g::Graph, raw_entry::Vector{T} where T)
    point_value = Dict()
    latitude, longitude = raw_entry
    # The Point data type is a set of latitude/longitude coordinates, stored within RedisGraph as a pair of 32-bit floats. (c) https://oss.redislabs.com/redisgraph/datatypes/
    point_value["latitude"] = parse(Float32, latitude)
    point_value["longitude"] = parse(Float32, longitude)
    return point_value
end

function parsemap(g::Graph, raw_entry::Vector{T} where T)
    map = Dict()
    reshaped_entry = reshape(raw_entry, 2, :)
    for entry in eachcol(reshaped_entry)
        map_key, map_value = entry
        entry_type, entry_value = map_value
        map[map_key] = parsevalue(g, VALUE_TYPE(entry_type), entry_value)
    end
    return map
end

function parsepath(g::Graph, raw_entry::Vector{T} where T)
    nodes_entry_type, nodes_entry = raw_entry[1]
    nodes = parsevalue(g, VALUE_TYPE(nodes_entry_type), nodes_entry)
    edges_entry_type, edges_entry = raw_entry[2]
    edges = parsevalue(g, VALUE_TYPE(edges_entry_type), edges_entry)
    return Path(nodes, edges)
end

function parsearray(g::Graph, raw_entry::Vector{T} where T)
    array = []
    for entry in raw_entry
        entry_type, entry_value = entry
        push!(array, parsevalue(g, VALUE_TYPE(entry_type), entry_value))
    end
    return array
end

function parseprops(g::Graph, raw_props::Vector{T} where T)
    props = Dict()
    for raw_prop in raw_props
        prop_name = getprop(g, raw_prop[1])
        prop_value = parsevalue(g, VALUE_TYPE(raw_prop[2]), raw_prop[3])
        props[prop_name] = prop_value
    end
    return props
end

function parsenode(g::Graph, raw_entry::Vector{T} where T)
    node_id = raw_entry[1]
    label = NaN
    if length(raw_entry[2]) != 0
        label = getlabel(g, raw_entry[2][1])
    end
    properties = parseprops(g, raw_entry[3])
    return Node(node_id, "", label, properties)
end

function parseedge(g::Graph, raw_entry::Vector{T} where T)
    edge_id = raw_entry[1]
    relation = getrelation(g, raw_entry[2])
    src_node_id = Node(raw_entry[3])
    dst_node_id = Node(raw_entry[4])
    properties = parseprops(g, raw_entry[5])
    return Edge(edge_id, relation, src_node_id, dst_node_id, properties)
end

function parsestatistics(raw_stats::Vector{String})
    statistics = Dict()
    for stat in raw_stats
        stat_name, stat_value = split(stat, ": ")
        statistics[stat_name] = parse(Float64, split(stat_value, " ")[1])
    end
    return statistics
end

function parseresults(g::Graph, raw_results::Vector{Vector{T} where T})
    if length(raw_results) == 0
        return nothing, nothing
    end

    header = ResultHeaderEntry.(raw_results[1])
    results = Vector{Any}()
    for row in raw_results[2]
        for (header_entry, raw_result_entry) in zip(header, row)
            entry = parsevalue(g, VALUE_TYPE(raw_result_entry[1]), raw_result_entry[2])
            push!(results, entry)
        end
    end
    return header, results
end

struct ResultHeaderEntry
    r_entry_type::Int
    r_entry_name::String
    ResultHeaderEntry(raw_header_entry) = new(raw_header_entry[1], raw_header_entry[2])
end

struct QueryResult
    statistics::Dict{String, Float64}
    header::Union{Vector{ResultHeaderEntry}, Nothing}
    results::Union{Vector{Any}, Nothing}
    function QueryResult(g::Graph, raw_query_result::Vector{Vector{String}})
        stats = parsestatistics(raw_query_result[length(raw_query_result)])
        new(stats, nothing, nothing)
    end
    function QueryResult(g::Graph, raw_query_result::Vector{Vector{T} where T})
        stats = parsestatistics(raw_query_result[length(raw_query_result)])
        header, results = parseresults(g, raw_query_result[1:length(raw_query_result)-1])
        new(stats, header, results)
    end
end

function isempty(result::QueryResult)
    length(result.header) == 0
end


text:
import Redis: flatten, flatten_command, convert_response

@testset "Flatten" begin
    @test flatten("simple") == "simple"
    @test flatten(1) == "1"
    @test flatten(2.5) == "2.5"
    @test flatten([1, "2", 3.5]) == ["1", "2", "3.5"]

    s = Set([1, 5, "10.9"])
    @test Set(flatten(s)) == Set(["1", "5", "10.9"])

    d = Dict{Any, Any}(1 => 2, 3 => 4)
    @test Set(flatten(d)) == Set(["1", "2", "3", "4"])
end

@testset "Commands" begin
    result = flatten_command(1, 2, ["4", "5", 6.7], 8)
    @test result == ["1", "2", "4", "5", "6.7", "8"]
end

@testset "Convert" begin
    @test convert_response(Dict{AbstractString, AbstractString}, ["1","2","3","4"]) == Dict("1" => "2", "3" => "4")
    @test convert_response(Dict{AbstractString, AbstractString}, []) == Dict()
    @test convert_response(Float64, "12.3") ≈ 12.3
    @test convert_response(Float64, 10) ≈ 10.
    @test convert_response(Bool, "OK")
    @test !convert_response(Bool, "f")
    @test convert_response(Bool, 1)
    @test !convert_response(Bool, 0)
end


text:
function get_indentation(nspaces::Int)
    res = ""
    for i = 1:nspaces
        res *= " "
    end
    return res
end

function is_leafnode(obj)
    if isa(obj, Array{Any, 1})
        n = length(obj)
        non_dict = true
        i = 1
        while non_dict
            if !isa(obj[i], ASCIIString) && !(typeof(obj[i]) <: Number)
                non_dict = false
                res = false
            elseif i == n
                res = true
                break
            end
            i += 1
        end
    else
        res = false
    end
    return res
end


text:
using DataFrames, XLSX

# Definimos la estructura de BasicInfo
struct BasicInfo
    match_id::String
    basic_info::DataFrame
    team_id::String
    team_info::DataFrame
end

# Funcion de lectura del fichero BasicInfo
function leer_BasicInfo(path::String)
    # Leemos fichero
    xf = XLSX.readxlsx(path)
    sps_names = XLSX.sheetnames(xf)

    # Leemos primera hoja
    info_match = xf[sps_names[1]]
    # Convertimos a DataFrame
    df_info_match = DataFrame(info_match[:], :auto)
    # Eliminamos primera fila
    delete!(df_info_match, 1)
    # Reasignamos los nombres a las columnas
    df_info_match = rename(df_info_match, ["Basic Info", "Data"])
    # Obtenemos id del match
    match_id = df_info_match[2,2]

    # Leemos segunda hoja
    info_players = xf[sps_names[2]]
    # Convertimos a DataFrame
    df_info_players = DataFrame(info_players[:], :auto)
    # Guardamos los nombres de la primera fila
    column_names = Vector{String}(df_info_players[1, :])
    # Renombramos los nombres de las columnas
    df_info_players = rename(df_info_players, column_names)
    # Eliminamos primera fila
    delete!(df_info_players, 1)
    # Obtenemos id del equipo
    team_id = df_info_match[5, 2]

    # Devolvemos la informacion pertinente
    BasicInfo(match_id, df_info_match, team_id, df_info_players)

end

max_stars_repo_name: vntr-CC/Forward-Football

max_issues_repo_name: vntr-CC/Forward-Football

max_forks_repo_name: vntr-CC/Forward-Football


text:
# TODO place this doc
# """
# Type assertions are supported in the left hand of rules
# to match and access literal values both when using classic
# rewriting and EGraph based rewriting.
# To use a type assertion pattern, add `::T` after
# a pattern variable in the `left_hand` of a rule.
# """

# TODO HASH CACHING!

using Parameters

import Base.==

abstract type Rule end
# Must override
==(a::Rule, b::Rule) = false

abstract type SymbolicRule <: Rule end

"""
Rules defined as `left_hand => right_hand` are
called *symbolic rewrite* rules. Application of a *rewrite* Rule
is a replacement of the `left_hand` pattern with
the `right_hand` substitution, with the correct instantiation
of pattern variables. Function call symbols are not treated as pattern
variables, all other identifiers are treated as pattern variables.
Literals such as `5, :e, "hello"` are not treated as pattern
variables.

```julia
Rule(:(a * b => b * a))
```
"""
@auto_hash_equals struct RewriteRule <: SymbolicRule
    left::Pattern
    right::Pattern
    patvars::Vector{Symbol}
    function RewriteRule(l,r)
        pvars = patvars(l) ∪ patvars(r)
        # sort!(pvars)
        setindex!(l, pvars)
        setindex!(r, pvars)
        new(l,r,pvars)
    end
end

# =============================================================================

abstract type BidirRule <: SymbolicRule end

"""
This type of *anti*-rules is used for checking contradictions in the EGraph
backend. If two terms, corresponding to the left and right hand side of an
*anti-rule* are found in an [`EGraph`], saturation is halted immediately.
"""
@auto_hash_equals struct UnequalRule <: BidirRule
    left::Pattern
    right::Pattern
    patvars::Vector{Symbol}
    function UnequalRule(l,r)
        pvars = patvars(l) ∪ patvars(r)
        # sort!(pvars)
        setindex!(l, pvars)
        setindex!(r, pvars)
        new(l,r,pvars)
    end
end

"""
```julia
Rule(:(a * b == b * a))
```
"""
@auto_hash_equals struct EqualityRule <: BidirRule
    left::Pattern
    right::Pattern
    patvars::Vector{Symbol}
    function EqualityRule(l,r)
        pvars = patvars(l) ∪ patvars(r)
        # sort!(pvars)
        setindex!(l, pvars)
        setindex!(r, pvars)
        new(l,r,pvars)
    end
end

"""
Rules defined as `left_hand |> right_hand` are
called `dynamic` rules. Dynamic rules behave like anonymous functions.
Instead of a symbolic substitution, the right hand of
a dynamic `|>` rule is evaluated during rewriting:
matched values are bound to pattern variables as in a
regular function call. This allows for dynamic computation
of right hand sides.

Dynamic rule
```julia
Rule(:(a::Number * b::Number |> a*b))
```
"""
@auto_hash_equals struct DynamicRule <: Rule
    left::Pattern
    right::Any
    patvars::Vector{Symbol} # useful set of pattern variables
    function DynamicRule(l, r)
        pvars = patvars(l)
        # sort!(pvars)
        setindex!(l, pvars)
        new(l, r, pvars)
    end
end


text:
#
# Defining mutually-recursive types.
#

macro rectypes(blk)
    if Meta.isexpr(blk, :block)
        names = Symbol[]
        for def in blk.args
            if Meta.isexpr(def, :macrocall) && length(def.args) > 1 && def.args[1] == GlobalRef(Core, Symbol("@doc"))
                def = def.args[end]
            end
            if Meta.isexpr(def, :struct, 3) && def.args[2] isa Symbol
                push!(names, def.args[2])
            end
        end
        name2deps = Dict{Symbol,Vector{Symbol}}()
        for def in blk.args
            if Meta.isexpr(def, :macrocall) && length(def.args) > 1 && def.args[1] == GlobalRef(Core, Symbol("@doc"))
                def = def.args[end]
            end
            if Meta.isexpr(def, :struct, 3) && def.args[2] isa Symbol
                name = def.args[2]
                deps = name2deps[name] = Symbol[]
                collect_deps!(name, name2deps, names, def)
                rewrite_deps!(name2deps, def)
            end
        end
        for name in reverse(names)
            repl = Symbol("$(name)_")
            deps = name2deps[name]
            push!(blk.args, Expr(:const, Expr(:(=), name, !isempty(deps) ? Expr(:curly, repl, deps...) : repl)))
            push!(blk.args, :(Base.show_datatype(io::IO, ::Type{$name}) = print(io, $(QuoteNode(name)))))
        end
    end
    return esc(blk)
end

collect_deps!(name, name2deps, names, ex) =
    if ex isa Expr
        for arg in ex.args
            collect_deps!(name, name2deps, names, arg)
        end
    elseif ex isa Symbol && ex in names && !(ex in name2deps[name])
        if ex in keys(name2deps)
            for dep in name2deps[ex]
                collect_deps!(name, name2deps, names, dep)
            end
        elseif ex != name
            push!(name2deps[name], ex)
        end
    end

rewrite_deps!(name2deps, ex) =
    if ex isa Expr
        for (i, arg) in enumerate(ex.args)
            if arg isa Symbol && arg in keys(name2deps)
                repl = Symbol("$(arg)_")
                deps = name2deps[arg]
                ex.args[i] = !isempty(deps) ? Expr(:curly, repl, deps...) : repl
                rewrite_deps!(name2deps, ex.args[i])
            else
                fn = Meta.isexpr(ex, (:function, :(=)), 2) &&
                     i == 1 &&
                     Meta.isexpr(arg, :call) &&
                     length(arg.args) >= 1 &&
                     !Meta.isexpr(arg.args[1], :curly)
                rewrite_deps!(name2deps, arg)
                if fn && Meta.isexpr(arg.args[1], :curly)
                    ex.args[i] = Expr(:where, arg, arg.args[1].args[2:end]...)
                end
            end
        end
    end


text:
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = test_page_model()
    include("../src/components/ClimateTemperature.jl")

    climatetemperature = addclimatetemperature(m, use_seaice)

    climatetemperature[:y_year_0] = 2015.
    climatetemperature[:y_year] = Mimi.dim_keys(m.md, :time)
    climatetemperature[:fant_anthroforcing] = readpagedata(m, "test/validationdata/$valdir/fgt.csv")

    p = load_parameters(m)
    set_leftover_params!(m, p)

    ##running Model
    run(m)

    rtl = m[:ClimateTemperature, :rtl_realizedtemperature]
    rtl_compare = readpagedata(m, "test/validationdata/$valdir/rtl_realizedtemperature.csv")

    @test rtl ≈ rtl_compare rtol = 1e-3

    rto = m[:ClimateTemperature, :rto_g_oceantemperature]
    rto_compare = readpagedata(m, "test/validationdata/$valdir/rto_g_oceantemperature.csv")

    @test rto ≈ rto_compare rtol = 1e-3
end


text:
struct ProxALEvaluator <: AbstractNLPEvaluator
    problem::ProxALProblem
    modelinfo::ModelInfo
    algparams::AlgParams
    opfdata::OPFData
    rawdata::RawData
    space::AbstractBackend
    comm::Union{MPI.Comm,Nothing}
end

"""
    ProxALEvaluator(
        case_file::String,
        load_file::String,
        modelinfo::ModelInfo,
        algparams::AlgParams,
        space::AbstractBackend=JuMPBackend(),
        comm::MPI.Comm = MPI.COMM_WORLD
    )

Instantiate multi-period ACOPF
specified in `case_file` with loads in `load_file` with model parameters
`modelinfo`, algorithm parameters `algparams`, modeling backend `space`, and
a MPI communicator `comm`.

"""
function ProxALEvaluator(
    case_file::String,
    load_file::String,
    modelinfo::ModelInfo,
    algparams::AlgParams,
    space::AbstractBackend=JuMPBackend(),
    opt_sol::Dict = Dict(),
    lyapunov_sol::Dict = Dict(),
    comm::Union{MPI.Comm,Nothing} = MPI.COMM_WORLD
)
    rawdata = RawData(case_file, load_file)
    opfdata = opf_loaddata(
        rawdata;
        time_horizon_start = 1,
        time_horizon_end = modelinfo.num_time_periods,
        load_scale = modelinfo.load_scale,
        ramp_scale = modelinfo.ramp_scale,
        corr_scale = modelinfo.corr_scale
    )
    if modelinfo.time_link_constr_type != :penalty
        @warn("ProxAL is guaranteed to converge only when time_link_constr_type = :penalty\n"*
              "         Forcing time_link_constr_type = :penalty\n")
        modelinfo.time_link_constr_type = :penalty
    end
    if modelinfo.num_ctgs > 1 && algparams.decompCtgs
        if modelinfo.ctgs_link_constr_type ∉ [:frequency_penalty, :preventive_penalty, :corrective_penalty]
            str = "ProxAL is guaranteed to converge only when "*
                  "ctgs_link_constr_type ∈ [:frequency_penalty, :preventive_penalty, :corrective_penalty]\n"
            if modelinfo.ctgs_link_constr_type == :preventive_equality
                @warn(str * "         Forcing ctgs_link_constr_type = :preventive_penalty\n")
                modelinfo.ctgs_link_constr_type = :preventive_penalty
            elseif modelinfo.ctgs_link_constr_type ∈ [:corrective_equality, :corrective_inequality]
                @warn(str * "         Forcing ctgs_link_constr_type = :corrective_penalty\n")
                modelinfo.ctgs_link_constr_type = :corrective_penalty
            else
                @warn(str * "         Forcing ctgs_link_constr_type = :frequency_penalty\n")
                modelinfo.ctgs_link_constr_type = :frequency_penalty
            end
        end
    end

    # ctgs_arr = deepcopy(rawdata.ctgs_arr)
    problem = ProxALProblem(opfdata, rawdata, modelinfo, algparams, space, comm, opt_sol, lyapunov_sol)
    return ProxALEvaluator(problem, modelinfo, algparams, opfdata, rawdata, space, comm)
end

"""
    optimize!(nlp::ProxALEvaluator)

Solve problem using the `nlp` evaluator
of the decomposition algorithm.

"""
function optimize!(
    nlp::ProxALEvaluator;
    print_timings::Bool=false,
    θ_t_initial::Union{Real,Nothing}=nothing,
    θ_c_initial::Union{Real,Nothing}=nothing,
    ρ_t_initial::Union{Real,Nothing}=nothing,
    ρ_c_initial::Union{Real,Nothing}=nothing,
    τ_initial::Union{Real,Nothing}=nothing
)
    algparams = nlp.algparams
    modelinfo = nlp.modelinfo
    runinfo   = nlp.problem
    opfdata   = nlp.opfdata
    comm      = nlp.comm

    has_ctgs(modelinfo, algparams) = (algparams.decompCtgs && modelinfo.num_ctgs > 0)
    τ_default(modelinfo, algparams) = has_ctgs(modelinfo, algparams) ? 2.0*max(algparams.ρ_t, algparams.ρ_c) : 2.0*algparams.ρ_t
    maxθ(modelinfo, algparams) = has_ctgs(modelinfo, algparams) ? max(algparams.θ_t, algparams.θ_c) : algparams.θ_t

    algparams.θ_t = algparams.θ_c = (1/algparams.tol^2)
    algparams.ρ_t = algparams.ρ_c = modelinfo.obj_scale
    algparams.τ = τ_default(modelinfo, algparams)
    !isnothing(θ_t_initial) && (algparams.θ_t = θ_t_initial)
    !isnothing(θ_c_initial) && (algparams.θ_c = θ_c_initial)
    !isnothing(ρ_t_initial) && (algparams.ρ_t = ρ_t_initial)
    !isnothing(ρ_c_initial) && (algparams.ρ_c = ρ_c_initial)
    !isnothing(τ_initial) && (algparams.τ = τ_initial)
    if (!isnothing(ρ_t_initial) || !isnothing(ρ_c_initial)) && isnothing(τ_initial)
        algparams.τ = τ_default(modelinfo, algparams)
    end
    runinfo.initial_solve &&
        (algparams_copy = deepcopy(algparams))
    opfBlockData = runinfo.opfBlockData
    nlp_opt_sol = runinfo.nlp_opt_sol
    nlp_soltime = runinfo.nlp_soltime
    ngen = length(opfdata.generators)
    nbus = length(opfdata.buses)
    T = modelinfo.num_time_periods
    K = modelinfo.num_ctgs + 1 # base case counted separately
    x = runinfo.x
    λ = runinfo.λ
    # number of contingency per blocks
    k_per_block = (algparams.decompCtgs) ? 1 : K

    function transfer!(blk, opt_sol, solution)
        # Pg
        fr = 1 ; to = ngen * k_per_block
        @views opt_sol[fr:to, blk] .= solution.pg[:]
        # Qg
        fr = to + 1 ; to += ngen * k_per_block
        @views opt_sol[fr:to, blk] .= solution.qg[:]
        # vm
        fr = to +1 ; to = fr + nbus * k_per_block - 1
        @views opt_sol[fr:to, blk] .= solution.vm[:]
        # va
        fr = to + 1 ; to = fr + nbus * k_per_block - 1
        @views opt_sol[fr:to, blk] .= solution.va[:]
        # wt
        fr = to +1  ; to = fr + k_per_block -1
        @views opt_sol[fr:to, blk] .= solution.ωt[:]
        # St
        fr = to +1  ; to = fr + ngen - 1
        @views opt_sol[fr:to, blk] .= solution.st[:]
        # zt
        fr = to +1  ; to = fr + ngen -1
        if haskey(solution, :zt)
            @views opt_sol[fr:to, blk] .= solution.zt[:]
        end
        # sk
        fr = to +1  ; to = fr + ngen * k_per_block - 1
        if haskey(solution, :sk)
            @views opt_sol[fr:to, blk] .= solution.sk[:]
        end
        # zk
        fr = to +1  ; to = fr + ngen * k_per_block - 1
        if haskey(solution, :zk)
            @views opt_sol[fr:to, blk] .= solution.zk[:]
        end
    end
    #------------------------------------------------------------------------------------
    function blocknlp_copy(blk, x_ref, λ_ref, alg_ref)
        model = opfBlockData.blkModel[blk]
        # Update objective
        set_objective!(model, alg_ref, x_ref, λ_ref)
        x0 = @view opfBlockData.colValue[:, blk]
        solution = optimize!(model, x0, alg_ref)
        transfer!(blk, nlp_opt_sol, solution)
    end
    #------------------------------------------------------------------------------------
    function blocknlp_recreate(blk, x_ref, λ_ref, alg_ref)
        model = opfBlockData.blkModel[blk]
        init!(model, alg_ref)
        # Update objective
        set_objective!(model, alg_ref, x_ref, λ_ref)
        x0 = @view opfBlockData.colValue[:, blk]
        solution = optimize!(model, x0, alg_ref)
        transfer!(blk, nlp_opt_sol, solution)
    end
    #------------------------------------------------------------------------------------
    function primal_update()
        runinfo.wall_time_elapsed_actual += @elapsed begin
            # Primal update except penalty vars
            nlp_soltime .= 0.0
            nlp_soltime_local = 0.0
            for blk in runinfo.par_order
                if ismywork(blk, comm)
                    nlp_opt_sol[:,blk] .= 0.0
                    # nlp_soltime[blk] = @elapsed blocknlp_copy(blk; x_ref = x, λ_ref = λ, alg_ref = algparams)
                    nlp_soltime[blk] = @elapsed blocknlp_recreate(blk, x, λ, algparams)
                    nlp_soltime_local += nlp_soltime[blk]
                end
            end
            if comm_rank(comm) == 0
              print_timings && println("Solve subproblems(): $nlp_soltime_local")
            end

            print_timings && comm_barrier(comm)
            elapsed_t = @elapsed begin
                requests = comm_neighbors!(nlp_opt_sol, opfBlockData, runinfo, comm)
                # Every worker sends his contribution
                comm_sum!(nlp_soltime, comm)
                comm_wait!(requests)
                print_timings && comm_barrier(comm)
            end
            if comm_rank(comm) == 0
                print_timings && println("Comm primals: $elapsed_t")
            end

            # Update primal values
            elapsed_t = @elapsed begin
                for blk in runinfo.par_order
                    block = opfBlockData.blkIndex[blk]
                    k = block[1]
                    t = block[2]
                    if ismywork(blk, comm)
                        # Updating my own primal values
                        opfBlockData.colValue[:,blk] .= nlp_opt_sol[:,blk]
                        update_primal_nlpvars(x, opfBlockData, blk, modelinfo, algparams)
                        for blkn in runinfo.par_order
                            blockn = opfBlockData.blkIndex[blkn]
                            kn = blockn[1]
                            tn = blockn[2]
                            # Updating the received neighboring primal values
                            if (tn == t-1 || tn == t+1) && !ismywork(blkn, comm)
                                opfBlockData.colValue[:,blkn] .= nlp_opt_sol[:,blkn]
                                update_primal_nlpvars(x, opfBlockData, blkn, modelinfo, algparams)
                            end
                        end
                    end
                end
                print_timings && comm_barrier(comm)
            end

            if comm_rank(comm) == 0
              print_timings && println("update_primal_nlpvars(): $elapsed_t")
            end

            # Primal update of penalty vars
            elapsed_t = @elapsed begin
                for blk in runinfo.par_order
                    if ismywork(blk, comm)
                        update_primal_penalty(x, opfdata, opfBlockData, blk, x, λ, modelinfo, algparams)
                    end
                end
                print_timings && comm_barrier(comm)
            end
            if comm_rank(comm) == 0
              print_timings && println("update_primal_penalty(): $elapsed_t")
            end
        end

        print_timings && comm_barrier(comm)
        elapsed_t = @elapsed begin
            requests = comm_neighbors!(x.Zt, opfBlockData, runinfo, comm)
            comm_wait!(requests)
            print_timings && comm_barrier(comm)
        end
        if comm_rank(comm) == 0
            print_timings && println("Comm penalty: $elapsed_t")
        end
        runinfo.wall_time_elapsed_ideal += isempty(runinfo.par_order) ? 0.0 : maximum([nlp_soltime[blk] for blk in runinfo.par_order])
        runinfo.wall_time_elapsed_ideal += elapsed_t
    end
    #------------------------------------------------------------------------------------
    function dual_update()
        elapsed_t = @elapsed begin
            maxviol_t = 0.0; maxviol_c = 0.0
            for blk in runinfo.par_order
                block = opfBlockData.blkIndex[blk]
                k = block[1]
                t = block[2]
                if ismywork(blk, comm)
                    lmaxviol_t, lmaxviol_c = update_dual_vars(λ, opfdata, opfBlockData, blk, x, modelinfo, algparams)
                    maxviol_t = max(maxviol_t, lmaxviol_t)
                    maxviol_c = max(maxviol_c, lmaxviol_c)
                end
            end
            print_timings && comm_barrier(comm)
        end
        if comm_rank(comm) == 0
          print_timings && println("update_dual_vars(): $elapsed_t")
        end
        elapsed_t = @elapsed begin
            requests = comm_neighbors!(λ.ramping, opfBlockData, runinfo, comm)
            maxviol_t = comm_max(maxviol_t, comm)
            maxviol_c = comm_max(maxviol_c, comm)
            push!(runinfo.maxviol_t, maxviol_t)
            push!(runinfo.maxviol_c, maxviol_c)
            comm_wait!(requests)
            print_timings && comm_barrier(comm)
        end
        if comm_rank(comm) == 0
            print_timings && println("Comm duals: $elapsed_t")
        end
        runinfo.wall_time_elapsed_actual += elapsed_t
        runinfo.wall_time_elapsed_ideal += elapsed_t
    end
    #------------------------------------------------------------------------------------
    function proximal_parameter_update()
        elapsed_t = @elapsed begin
            if algparams.updateτ && runinfo.iter > 1
                delta = (runinfo.lyapunov[end-1] - runinfo.lyapunov[end])/abs(runinfo.lyapunov[end])
                if delta < -1e-4 && algparams.τ < 320.0*maxθ(modelinfo, algparams)
                    algparams.τ *= 2.0
                end
            end
        end
        runinfo.wall_time_elapsed_actual += elapsed_t
        runinfo.wall_time_elapsed_ideal += elapsed_t
    end
    #------------------------------------------------------------------------------------
    function penalty_parameter_update()
        elapsed_t = @elapsed begin
            if max(runinfo.maxviol_t[end], runinfo.maxviol_d[end]) <= algparams.tol
                if runinfo.maxviol_t_actual[end] > algparams.tol && algparams.θ_t < 1e+12
                    algparams.θ_t *= 10.0
                    if algparams.verbose > 0 && comm_rank(comm) == 0
                        (algparams.θ_t >= 1e+12) &&
                            (@warn "penalty parameter too large. problem will likely not converge" maxlog=1)
                    end
                end
            end
            if (algparams.decompCtgs && modelinfo.num_ctgs > 0 &&
                modelinfo.ctgs_link_constr_type ∈ [:preventive_penalty, :corrective_penalty] &&
                max(runinfo.maxviol_c[end], runinfo.maxviol_d[end]) <= algparams.tol
            )
                if runinfo.maxviol_c_actual[end] > algparams.tol && algparams.θ_c < 1e+12
                    algparams.θ_c *= 10.0
                    if algparams.verbose > 0 && comm_rank(comm) == 0
                        (algparams.θ_c >= 1e+12) &&
                            (@warn "penalty parameter too large. problem will likely not converge" maxlog=1)
                    end
                end
            end
        end
        runinfo.wall_time_elapsed_actual += elapsed_t
        runinfo.wall_time_elapsed_ideal += elapsed_t
    end
    #------------------------------------------------------------------------------------
    function auglag_parameter_update()
        elapsed_t = @elapsed begin
            # ρ_t update
            if algparams.updateρ_t
                if runinfo.maxviol_t[end] > 10.0*runinfo.maxviol_d[end] && algparams.ρ_t < 32.0*algparams.θ_t
                    algparams.ρ_t = min(2.0*algparams.ρ_t, 32.0*algparams.θ_t)
                    algparams.τ = τ_default(modelinfo, algparams)
                elseif runinfo.maxviol_d[end] > 10.0*runinfo.maxviol_t[end]
                    algparams.ρ_t *= 0.5
                    algparams.τ = τ_default(modelinfo, algparams)
                end
            end

            # ρ_c update
            if algparams.updateρ_c && algparams.decompCtgs && modelinfo.num_ctgs > 0
                @info "tau-update not finalized when decompCtgs = true" maxlog=1
                if runinfo.maxviol_c[end] > 10.0*runinfo.maxviol_d[end] && algparams.ρ_c < 32.0*algparams.θ_c
                    algparams.ρ_c = min(2.0*algparams.ρ_c, 32.0*algparams.θ_c)
                    algparams.τ = 2.0*max(algparams.ρ_t, algparams.ρ_c)
                elseif runinfo.maxviol_d[end] > 10*runinfo.maxviol_c[end]
                    algparams.ρ_c *= 0.5
                    algparams.τ = 2.0*max(algparams.ρ_t, algparams.ρ_c)
                end
            end
        end
        runinfo.wall_time_elapsed_actual += elapsed_t
        runinfo.wall_time_elapsed_ideal += elapsed_t
    end

    # Initialization of Pg, Qg, Vm, Va via a OPF solve
    algparams.init_opf && opf_initialization!(nlp)

    function iteration()

        if runinfo.initial_solve
            if runinfo.iter == 1
                algparams.ρ_t = 0.0
                algparams.ρ_c = 0.0
                algparams.τ = 0.0
            elseif runinfo.iter == 2
                algparams = deepcopy(algparams_copy)
            end
        end

        # use this to compute the KKT error at the end of the loop
        runinfo.xprev = deepcopy(x)
        runinfo.λprev = deepcopy(λ)

        # Primal update
        primal_update()

        # Dual update
        dual_update()

        # Update counters and write output
        runinfo_update(runinfo, opfdata, opfBlockData, modelinfo, algparams, comm)

        # Prox update
        proximal_parameter_update()

        # Penalty update
        penalty_parameter_update()

        # Auglag parameter update
        auglag_parameter_update()
    end

    for runinfo.iter=1:algparams.iterlim
        iteration()

        # Check convergence
        if max(
            runinfo.maxviol_t[end],
            runinfo.maxviol_c[end],
            runinfo.maxviol_t_actual[end],
            runinfo.maxviol_c_actual[end],
            runinfo.maxviol_d[end]
        ) <= algparams.tol
            break
        end
    end
    return runinfo
end

function opf_initialization!(nlp::ProxALEvaluator)
    runinfo   = nlp.problem
    modelinfo = nlp.modelinfo
    opfdata   = nlp.opfdata
    rawdata   = nlp.rawdata

    modelinfo_single = deepcopy(modelinfo)
    modelinfo_single.num_time_periods = 1
    primal = ProxAL.OPFPrimalSolution(opfdata, modelinfo_single)
    dual = ProxAL.OPFDualSolution(opfdata, modelinfo_single)
    algparams = AlgParams()
    algparams.mode = :coldstart
    algparams.optimizer = optimizer_with_attributes(Ipopt.Optimizer,
            "print_level" => Int64(algparams.verbose > 0)*5)
    blockmodel = ProxAL.JuMPBlockBackend(1, opfdata, rawdata, algparams, modelinfo_single, 1, 1, 0)
    ProxAL.init!(blockmodel, algparams)
    ProxAL.set_objective!(blockmodel, algparams, primal, dual)
    n = JuMP.num_variables(blockmodel.model)
    x0 = zeros(n)
    solution = ProxAL.optimize!(blockmodel, x0, algparams)
    nbus = length(opfdata.buses)
    ngen = length(opfdata.generators)
    K = modelinfo.num_ctgs
    T = modelinfo.num_time_periods
    for i in 1:T*(K+1)
        for i=1:ngen
            runinfo.x.Pg[i,:,:] .= solution.pg[i]
            runinfo.x.Qg[i,:,:] .= solution.qg[i]
        end
        for i=1:nbus
            runinfo.x.Vm[i,:,:] .= solution.vm[i]
            runinfo.x.Va[i,:,:] .= solution.va[i]
        end
    end
end


text:
module CP_VB

include("../Misc.jl")
using .Misc

using Distributions, SpecialFunctions

export standard_VB, online_VB

function standard_VB(X::Array{ℜ,3}, R::Ƶ ; μ::ℜ=3.5, γ::ℜ=0.1, EPOCHS::Ƶ=1, ϵ::ℜ=1e-16) where {ℜ<:Real, Ƶ<:Int}
    I::Ƶ, J::Ƶ, K::Ƶ = size(X)

    a::ℜ, b::ℜ = I*J*K*γ, γ/μ
    α_R, α_IR, α_JR, α_KR = fill(a/R,R), fill(a/(I*R),I,R), fill(a/(J*R),J,R), fill(a/(K*R),K,R)

    log_λ = log(rand(Gamma(a,1.0/b)))
    log_θ_R = rand(Dirichlet(α_R .+ 1.0/R))
    log_θ_IR = reshape([log(θ_ir) for r=1:R for θ_ir=rand(Dirichlet(α_IR[:,r] .+ 1.0/I))],I,R)
    log_θ_JR = reshape([log(θ_jr) for r=1:R for θ_jr=rand(Dirichlet(α_JR[:,r] .+ 1.0/J))],J,R)
    log_θ_KR = reshape([log(θ_kr) for r=1:R for θ_kr=rand(Dirichlet(α_KR[:,r] .+ 1.0/K))],K,R)

    s, S_R, S_IR, S_JR, S_KR, S⁺ = zeros(ℜ,R), zeros(ℜ,R), zeros(ℜ,I,R), zeros(ℜ,J,R), zeros(ℜ,K,R), sum(X)
    q, log_ρ, log_ρ_ijk::ℜ = Array{ℜ}(undef,R), Array{ℜ}(undef,R), 0.0

    ELBO::Array{ℜ} = zeros(ℜ,EPOCHS)
    X_full = similar(X)

    for eph=1:EPOCHS
        S_R .= 0.0
        S_IR .= 0.0
        S_JR .= 0.0
        S_KR .= 0.0
        S⁺ = 0.0

        for k=1:K, j=1:J, i=1:I #order of traversal is important
            log_ρ .= log_λ .+ log_θ_R .+ log_θ_IR[i,:] .+ log_θ_JR[j,:] .+ log_θ_KR[k,:]
            log_ρ_ijk = logsumexp(log_ρ)
            q .= exp.(log_ρ .- log_ρ_ijk)

            X_full[i,j,k] = isnan(X[i,j,k]) ? exp(log_ρ_ijk) : X[i,j,k]
            s .= X_full[i,j,k] .* q

            S_R .+= s
            S_IR[i,:] .+= s
            S_JR[j,:] .+= s
            S_KR[k,:] .+= s
            S⁺ += X_full[i,j,k]

            ELBO[eph] += isnan(X[i,j,k]) ? X_full[i,j,k] : X[i,j,k]*log_ρ_ijk - lgamma(X[i,j,k] + 1.0)
        end

        ELBO[eph] += a*log(b) - (S⁺+a)*log(b+1.0) - sum(lgamma.(α_IR)) + sum(lgamma.(α_IR .+ S_IR))
        ELBO[eph] += sum(lgamma.(α_JR .+ S_JR)) - sum(lgamma.(α_JR)) + sum(lgamma.(α_KR .+ S_KR)) - sum(lgamma.(α_KR))
        ELBO[eph] += 2.0*sum(lgamma.(α_R)) - 2.0*sum(lgamma.(α_R .+ S_R))
        ELBO[eph] -= S⁺ * log_λ + sum(S_R .* log_θ_R) + sum(S_IR .* log_θ_IR) + sum(S_JR .* log_θ_JR) + sum(S_KR .* log_θ_KR)

        log_λ = digamma(S⁺+a) - log(b+1.0)
        log_θ_R .= digamma.(S_R.+α_R) .- digamma(S⁺+a)
        log_θ_IR .= digamma.(S_IR.+α_IR) .- digamma.(S_R.+α_R)'
        log_θ_JR .= digamma.(S_JR.+α_JR) .- digamma.(S_R.+α_R)'
        log_θ_KR .= digamma.(S_KR.+α_KR) .- digamma.(S_R.+α_R)'

    end
    return ELBO, X_full, log_λ, (log_θ_R, log_θ_IR, log_θ_JR, log_θ_KR)
end

function online_VB(X::Array{ℜ,3}, R::Ƶ, Ω::AbstractArray{ℜ,1}; μ::ℜ=3.5, γ::ℜ=0.1, ϵ::ℜ=1e-16) where {ℜ<:Real, Ƶ<:Int}
    I::Ƶ, J::Ƶ, K::Ƶ = size(X)

    a::ℜ, b::ℜ = I*J*K*γ, γ/μ
    α_R, α_IR, α_JR, α_KR = fill(a/R,R), fill(a/(I*R),I,R), fill(a/(J*R),J,R), fill(a/(K*R),K,R)

    log_λ = log(rand(Gamma(a,1.0/b)))
    log_θ_R = rand(Dirichlet(α_R .+ 1.0/R))
    log_θ_IR = reshape([log(θ_ir) for r=1:R for θ_ir=rand(Dirichlet(α_IR[:,r] .+ 1.0/I))],I,R)
    log_θ_JR = reshape([log(θ_jr) for r=1:R for θ_jr=rand(Dirichlet(α_JR[:,r] .+ 1.0/J))],J,R)
    log_θ_KR = reshape([log(θ_kr) for r=1:R for θ_kr=rand(Dirichlet(α_KR[:,r] .+ 1.0/K))],K,R)

    s, S_R, S_IR, S_JR, S_KR, S⁺ = zeros(ℜ,R), zeros(ℜ,R), zeros(ℜ,I,R), zeros(ℜ,J,R), zeros(ℜ,K,R), sum(X)
    q, log_ρ, log_ρ_ijk::ℜ = Array{ℜ}(undef,R), Array{ℜ}(undef,R), 0.0

    ELBO::Array{ℜ} = zeros(ℜ,length(Ω))
    X_full = similar(X)

    for (eph,ω) in enumerate(Ω)
        S_R .= 0.0
        S_IR .= 0.0
        S_JR .= 0.0
        S_KR .= 0.0
        S⁺ = 0.0

        for k=1:K, j=1:J, i=1:I #order of traversal is important
            log_ρ .= log_λ .+ log_θ_R .+ log_θ_IR[i,:] .+ log_θ_JR[j,:] .+ log_θ_KR[k,:]
            log_ρ_ijk = logsumexp(log_ρ)
            q .= exp.(log_ρ .- log_ρ_ijk)

            X_full[i,j,k] = isnan(X[i,j,k]) ? ω*exp(log_ρ_ijk) : rand(Binomial(Ƶ(X[i,j,k]),ω))
            s .= X_full[i,j,k] .* q

            S_R .+= s
            S_IR[i,:] .+= s
            S_JR[j,:] .+= s
            S_KR[k,:] .+= s
            S⁺ += X_full[i,j,k]

            ELBO[eph] += isnan(X[i,j,k]) ? X_full[i,j,k] : X[i,j,k]*log_ρ_ijk - lgamma(X[i,j,k] + 1.0)

        end

        ELBO[eph] += a*log(b) - (S⁺+a)*log(b+ω) - sum(lgamma.(α_IR)) + sum(lgamma.(α_IR .+ S_IR))
        ELBO[eph] += sum(lgamma.(α_JR .+ S_JR)) - sum(lgamma.(α_JR)) + sum(lgamma.(α_KR .+ S_KR)) - sum(lgamma.(α_KR))
        ELBO[eph] += 2.0*sum(lgamma.(α_R)) - 2.0*sum(lgamma.(α_R .+ S_R))
        ELBO[eph] -= S⁺ * log_λ + sum(S_R .* log_θ_R) + sum(S_IR .* log_θ_IR) + sum(S_JR .* log_θ_JR) + sum(S_KR .* log_θ_KR)

        log_λ = digamma(S⁺ + a) - log(b + ω) ## should be checked
        log_θ_R .= digamma.(S_R.+α_R) .- digamma(S⁺+a)
        log_θ_IR .= digamma.(S_IR.+α_IR) .- digamma.(S_R.+α_R)'
        log_θ_JR .= digamma.(S_JR.+α_JR) .- digamma.(S_R.+α_R)'
        log_θ_KR .= digamma.(S_KR.+α_KR) .- digamma.(S_R.+α_R)'

    end
    return ELBO, X_full, log_λ, (log_θ_R, log_θ_IR, log_θ_JR, log_θ_KR)
end

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
using SequentialMonteCarlo
using RNGPool
import SMCExamples.SMCSampler.defaultSMCSampler
using StaticArrays

include("test.jl")

setRNGs(0)

model, ltarget = defaultSMCSampler()

## just run the algorithm a few times

testSMC(model, 1024*1024, 2, false)
testSMC(model, 1024*1024, 2, false, 0.5)
testSMCParallel(model, 1024*1024, 2, false)
testSMCParallel(model, 1024*1024, 2, false, 0.5)


text:
module DS

using UUIDs
import Base.append!
import Base.iterate
import Base.push!

"""
Node element in LList-type structure. Maintains forward and reverse links, data, and UUID.
"""
mutable struct Node{T}
    data::Union{Nothing, T}
    next::Union{Nothing, Node{T}}
    prev::Union{Nothing, Node{T}}
    uuid::UUID
    Node{T}() where {T} = new(nothing, nothing, nothing, UUIDs.uuid4())
    Node{T}(data) where {T} = new(data, nothing, nothing, UUIDs.uuid4())
    Node{T}(data, prev) where {T} = new(data, nothing, prev, UUIDs.uuid4())
    Node{T}(data, next, prev) where {T} = new(data, next, prev, UUIDs.uuid4())
end

"""
DLList structure for representing GateCalls. Maintains start, end and position pointer.
"""
mutable struct DLList{T}
    len::Int
    head::Union{Nothing, Node{T}}
    tail::Union{Nothing, Node{T}}
    current::Union{Nothing, Node{T}}
    DLList{T}() where {T} = new(0, nothing, nothing, nothing)
end

function push!(cl::DLList{T}, data::T) where T
    if cl.len == 0
        n = Node{T}(data)
        cl.head = n
        cl.tail = n
        cl.current = n
        cl.current = cl.head
    else
        cl.tail.next = Node{T}(data, cl.tail)
        cl.tail.next.prev = cl.tail
        cl.tail = cl.tail.next
    end
    cl.len += 1
end

function clear!(l::DLList{T}) where T
    l.head = nothing
    l.tail = nothing
    l.current = nothing
    l.len = 0
end

function append!(l1::DLList{T}, l2::DLList{T}) where T
    if l2.current === nothing
        return
    end
    l2c = deepcopy(l2)
    l1.tail.next = l2c.head
    l2c.head.prev = l1.tail
    l1.len += l2c.len
    clear!(l2)
    return l1
end

function copy(l::DLList{T}) where T
    lc = DLList{T}()
    for val in l
        push!(lc, val)
    end
    return lc
end

function iterate(list::DLList{T}) where T
    if list.len == 0
        return nothing
    else
        list.current = list.head
        return (list.current.data, list.current.next)
    end
end

function iterate(list::DLList{T}, n::Union{Node{T}, Nothing}) where T
    if list.current.next === nothing
        return nothing
    else
        res = (list.current.next.data, list.current.next)
        list.current = list.current.next
        return res
    end
end

function delete!(list::DLList{T}, index::Int) where T
    if index < list.len
        throw(BoundsError(list, index))
    else
        #for i in
    end
end

function delete!(list::DLList{T}, uuid::UUID) where T
    if index < list.len
        throw(BoundsError(list, index))
    else
        for node in list
            if node.uuid === uuid
                #remove this node
                break;
            end
        end
    end
end

end


text:
# ------------------------------------------------------------------------------------------
#
#
# This notebook accompanies the [workshop paper](http://jiahao.github.io/parallel-prefix/)
# by Jiahao Chen and Alan Edelman entitled "Parallel Prefix Polymorphism Permits
# Parallelization, Presentation & Proof" and will appear in the proceedings of the [First
# Workshop for High Performance Technical Computing in Dynamic
# Languages](http://jiahao.github.io/hptcdl-sc14/), held in conjunction with [SC14: The
# International Conference on High Performance Computing, Networking, Storage and
# Analysis](http://sc14.supercomputing.org/).
#
#
# ------------------------------------------------------------------------------------------

# using Pkg; Pkg.add(["Compose", "Gadfly", "Interact"])
using Compose, Gadfly, Interact

# ------------------------------------------------------------------------------------------
# # `reduce()`
#
# Reduction applies a binary operator to a vector repeatedly to return a scalar. Thus +
# becomes sum, and * becomes prod.
#
# It is considered a basic parallel computing primitive.
#
#
# ------------------------------------------------------------------------------------------

reduce(+, 1:8), sum(1:8)  # triangular numbers

reduce(*, 1:8), prod(1:8) # factorials

boring(a,b)=a
@show reduce(boring, 1:8)
boring2(a,b)=b
@show reduce(boring2, 1:8)

# ------------------------------------------------------------------------------------------
# You can also use reduce to compute Fibonacci numbers using their recurrences.
#
# $$\begin{pmatrix} f_2 \\f_1 \end{pmatrix} = \begin{pmatrix} f_1 + f_0 \\ f_1 \end{pmatrix}
# = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} f_1 \\ f_0 \end{pmatrix} $$
#
# $$\begin{pmatrix} f_{n+1} \\ f_n \end{pmatrix} = \dots
# = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^n \begin{pmatrix} f_1 \\ f_0 \end{pmatrix}
# $$
#
# From this, you can show that
#
# $$\begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^n  = \begin{pmatrix} f_{n+1} & f_n \\ f_n &
# f_{n-1} \end{pmatrix} $$
#
# (this applies reduce to the same argument over and over again -- there are of course other
# ways)
# ------------------------------------------------------------------------------------------

M=[1 1; 1 0]
reduce(*,fill(M,3))
prod(fill(M,3))

@manipulate for n=1:100
    prod(fill(big.(M),n))
end

fib(j)=reduce(*, fill(M,j))
fib.([4,7])

# ------------------------------------------------------------------------------------------
# You can solve recurrences of any complexity using `reduce`. For example, `reduce` can
# compute a Hadamard matrix from its definition in terms of its submatrices:
#
# $$H_2 = \begin{pmatrix} H_1 & H_1 \\ H_1 & -H_1 \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1
# & -1 \end{pmatrix} \otimes H_1$$
#
# and so on.
#
# (Note: this is just using reduce to compute a matrix power.
# One can think of alternative ways for sure.)
# ------------------------------------------------------------------------------------------

# [A B]
# If A is m x n
# If B is p x q
# then kron(A,B) is mp x nq and has all the elements of A times all of the elements of B

A=[1 2;3 4]
B=[10 100; 1 -10]
⊗(A,B)=kron(A,B)

M=[ 1 1;1 -1]
H=⊗(⊗(⊗(M,M),M),M)
H*H'

Hadamard(n)=reduce(⊗, fill(M,n))
H=Hadamard(3)

H'H

H*H' #This is a legitimate Hadamard matrix

using StatsBase

rolldice() = rand(1:6)+rand(1:6)
rolls(n)=fit(Histogram,[rolldice() for i=1:n ],2:12,closed=:left)

rolls(1000)

[rolls(100) for i=1:10]

reduce(merge,[rolls(100) for i=1:10])

#Reduction of matrix multiplications
M=[rand(-3:3,2,2) for i=1:4]

display(M[4]*M[3]*M[2]*M[1])
display(reduce((A,B)->B*A, M)) #backward multiply
display(reduce(*, M))          #forward multiply

# ------------------------------------------------------------------------------------------
# In the following example we apply `reduce()` to  function composition:
# ------------------------------------------------------------------------------------------

(sin ∘ cos)(1), sin(cos(1))

((-) ∘ sin)(1), -sin(1)

( (!) ∘ )

# h=reduce( (f,g) -> (x->f(g(x))), [sin cos tan])
h=reduce(∘, [sin cos tan])
[h(π) sin(cos(tan(π)))]

cumsum(1:8)  # It is useful to know that cumsum is a linear operator
# You can use power method! Below is the underlying matrix
A=tril(ones(Int,8,8))

[A*(1:8),cumsum(1:8)]

# ------------------------------------------------------------------------------------------
# # `prefix`
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Having discussed `reduce`, we are now ready for the idea behind prefix sum.
# Prefix or scan is long considered an important parallel
# primitive as well.
#
# Suppose you wanted to compute the partial sums of a vector, i.e. given
# `y[1:n]`, we want to overwrite the vector `y` with the vector of partial sums
#
# ```julia
# new_y[1] = y[1]
# new_y[2] = y[1] + y[2]
# new_y[3] = y[1] + y[2] + y[3]
# ...
# ```
#
# At first blush, it seems impossible to parallelize this, since
#
# ```julia
# new_y[1] = y[1]
# new_y[2] = new_y[1] + y[2]
# new_y[3] = new_y[2] + y[3]
# ...
# ```
#
# which appears to be an intrinsically serial process.
# ------------------------------------------------------------------------------------------

function prefix_serial!(y, ⊕)
    for i=2:length(y)
        y[i] = y[i-1] ⊕ y[i]
        end
    y
end

prefix_serial!([1:8;],*)

prefix_serial!([rand(1:5,2,2) for i=1:4],*)

# ------------------------------------------------------------------------------------------
# However, it turns out that because addition (`+`) is associative, we can regroup the
# _order_ of how these sums are carried out. (This of course extends to other associative
# operations such as multiplication.) Another ordering of 8 associative operations is
# provided by `prefix8!`:
# ------------------------------------------------------------------------------------------

# Eight only for pedagogy
function prefix8!(y, ⋅)
    length(y)==8 || error("length 8 only")
    for i in [2,4,6,8]; y[i] = y[i-1] ⋅ y[i]; end
    for i in [  4,  8]; y[i] = y[i-2] ⋅ y[i]; end
    for i in [      8]; y[i] = y[i-4] ⋅ y[i]; end
    for i in [    6  ]; y[i] = y[i-2] ⋅ y[i]; end
    for i in [ 3,5,7 ]; y[i] = y[i-1] ⋅ y[i]; end
    y
end

# Generalization to any n
function prefix!(y, ⋅)
    l=length(y)
    k=ceil(Int, log2(l))
    @inbounds for j=1:k, i=2^j:2^j:min(l, 2^k)              #"reduce"
        y[i] = y[i-2^(j-1)] ⋅ y[i]
    end
    @inbounds for j=(k-1):-1:1, i=3*2^(j-1):2^j:min(l, 2^k) #"broadcast"
        y[i] = y[i-2^(j-1)] ⋅ y[i]
    end
    y
end

prefix!([1:12;],*)

# ------------------------------------------------------------------------------------------
# # Polymorphism for visualization
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We can visualize the operations with a little bit of trickery. In Julia, arrays are simply
# types that expose the array protocol. In particular, they need to implement  methods for
# the generic functions `length`, `getindex` and `setindex!`. The last two are used in
# indexing operations, since statements
#
#     y = A[1]
#     A[3] = y
#
# get desugared to
#
#     y = getindex(A, 1)
#     setindex!(A, y, 3)
#
# respectively.
#
# We can trace through the iterable by introduce a dummy array type, `DummyArray`, which
# stores no useful information but records every access to `getindex` and `setindex!`.
#
# Specifically:
#
# - `length(A::DummyArray)` returns an integer that is stored internally in the `A.length`
# field
# - `getindex(A::DummyArray, i)` records read access to the index `i` in the `A.read` field
# and always retuns `nothing`.
# - `setindex!(A::DummyArray, x, i)` records write access to the index `i`. The `A.history`
# field is appended with a new tuple consisting of the current `A.read` field and the index
# `i`.
#
# The way `DummyArray` works, it assumes an association between a single `setindex!` call
# and and all the preceding `getindex` calls since the previous `setindex!` call, which is
# sufficient for the purposes of tracing through prefix calls.
# ------------------------------------------------------------------------------------------

import Base: getindex, setindex!, length

isdefined(:DummyArray) || type DummyArray
    length :: Int
    read :: Vector
    history :: Vector{Any}
    DummyArray(length, read=[], history=[])=new(length, read, history)
end

length(A::DummyArray)=A.length

function getindex(A::DummyArray, i)
    push!(A.read, i)
    nothing
end

function setindex!(A::DummyArray, x, i)
    push!(A.history, (A.read, [i]))
    A.read = []
end

M = DummyArray(4)

M[7] = M[3] + M[2]
M

import Base.+
+(a::Void, b::Void)=a

@manipulate for i=1:100
    M = DummyArray(4)
    M[i] = M[i-2] + M[2*i]
    M
end

# ------------------------------------------------------------------------------------------
# We also need a dummy associative operator to pass to the prefix call.
# ------------------------------------------------------------------------------------------

A=prefix8!(DummyArray(8),+)

A.history

# ------------------------------------------------------------------------------------------
# Now let's visualize this! Each entry in `A.history` is rendered by a gate object:
# ------------------------------------------------------------------------------------------

mutable struct gate
    ins :: Vector
    outs:: Vector
end

import Gadfly.render

function render(G::gate, x₁, y₁, y₀; rᵢ=0.1, rₒ=0.25)
    ipoints = [(i, y₀+rᵢ) for i in G.ins]
    opoints = [(i, y₀+0.5) for i in G.outs]
    igates  = [Compose.circle(i..., rᵢ) for i in ipoints]
    ogates  = [Compose.circle(i..., rₒ) for i in opoints]
    lines = [line([i, j]) for i in ipoints, j in opoints]
    compose(context(units=UnitBox(0.5,0,x₁,y₁+1)),
    compose(context(), stroke(colorant"black"), fill(colorant"white"),
            igates..., ogates...),
    compose(context(), linewidth(0.3mm), stroke(colorant"black"),
            lines...))
end

A=gate([1,2],[2])
render(A,2,0,0)

# ------------------------------------------------------------------------------------------
# Now we render the whole algorithm. We have to scan through the trace twice; the first time
# merely calculates the maximum depth that needs to be drawn and the second time actually
# generates the objects.
# ------------------------------------------------------------------------------------------

function render(A::DummyArray)
    #Scan to find maximum depth
    olast = depth = 0
    for y in A.history
        (any(y[1] .≤ olast)) && (depth += 1)
        olast = maximum(y[2])
    end
    maxdepth = depth

    olast = depth = 0
    C = []
    for y in A.history
        (any(y[1] .≤ olast)) && (depth += 1)
        push!(C, render(gate(y...), A.length, maxdepth, depth))
        olast = maximum(y[2])
    end

    push!(C, compose(context(units=UnitBox(0.5,0,A.length,1)),
      [line([(i,0), (i,1)]) for i=1:A.length]...,
    linewidth(0.1mm), stroke(colorant"grey")))
    compose(context(), C...)
end

render(prefix!(DummyArray(8), +))

# ------------------------------------------------------------------------------------------
# Now we can see that `prefix!` rearranges the operations to form two spanning trees:
# ------------------------------------------------------------------------------------------

render(prefix!(DummyArray(120),+))

function prefix!(y, +)
    l=length(y)
    k=ceil(Int, log2(l))
    @inbounds for j=1:k, i=2^j:2^j:min(l, 2^k)              #"reduce"
        y[i] = y[i-2^(j-1)] + y[i]
    end
    @inbounds for j=(k-1):-1:1, i=3*2^(j-1):2^j:min(l, 2^k) #"broadcast"
        y[i] = y[i-2^(j-1)] + y[i]
    end
    y
end
@show prefix!(DummyArray(8),+)
render(prefix!(DummyArray(9),+))

#Look at the compose tree
#render(prefix8!(AccessArray(8),+)) |> introspect

# ------------------------------------------------------------------------------------------
# as contrasted with the serial code:
# ------------------------------------------------------------------------------------------

render(prefix_serial!(DummyArray(8),+))

@manipulate for npp=1:180
    render(prefix!(DummyArray(npp),+))
end

@manipulate for npp=1:180
    render(prefix_serial!(DummyArray(npp),+))
end

draw(PDF("serial.pdf", 3.5inch, 3.5inch), render(prefix_serial!(AccessArray(8),+)))
draw(PDF("tree.pdf", 3.5inch, 3.5inch), render(prefix8!(AccessArray(8),+)))

# ------------------------------------------------------------------------------------------
# # Polymorphism for parallel operations
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Now let's run `prefix` in parallel on every core on the CPU. Use `addprocs` to attach
# additional worker processes.
# ------------------------------------------------------------------------------------------

import Base: fetch, length, +, *
fetch(t::Vector) = map(fetch, t) #Vectorize fetch

#Define elementary operations on remote data
length(r1::Future)=length(fetch(r1))
+(r1::Future,r2::Future)=@spawnat r2.where fetch(r1)+fetch(r2)
*(r1::Future,r2::Future)=@spawnat r2.where fetch(r1)*fetch(r2)

# ------------------------------------------------------------------------------------------
# For 8 processes, the serial version requires 7 operations. The parallel version uses 11
# operations, but they are grouped into 5 simultaneous chunks of operations. Hence we expect
# that the parallel version runs in 5/7 the time needed for the naïve serial code.
# ------------------------------------------------------------------------------------------

#Before timing, force Julia to JIT compile the functions
[f([randn(3,3) for i=1:2], *) for f in (prefix_serial!, prefix!)]

# addprocs(max(0, Sys.CPU_CORES-nprocs())) #Add all possible processors
using JuliaRunClient
# initializeCluster(2)
initializeCluster(16) # more workers can be added with a larger subscription

timings=Dict{Int,Tuple{Float64,Float64}}()
for np in 2:nprocs()
    gc_enable(false) #Disable garbage collector for timing purposes
    n=1024#4096#3072

    #Create RemoteRefs on np processes
    r=[@spawnat p randn(n,n) for p in procs()[1:np]]
    s=fetch(r) #These are in-memory matrices

    #Repeat timings a couple of times and take smallest
    t_ser = minimum([(tic(); prefix_serial!(s, *); toc()) for i=1:2])
    t_par = minimum([(tic(); @sync prefix!(r, *) ; toc()) for i=1:2])
    timings[np] = (t_ser, t_par)
    gc_enable(true)
end
timings

#timings = #Data from julia.mit.edu
#Dict(
#     1 => ( 15.585,  15.721),
#     2 => ( 16.325,  15.893),
#     3 => ( 32.079,  31.678),
#     4 => ( 46.037,  33.257),
#     5 => ( 61.947,  51.412),
#     6 => ( 77.247,  50.769),
#     7 => ( 92.588,  67.796),
#     8 => (108.057,  65.185),
#     9 => (123.556,  68.197),
#    10 => (138.408,  66.658),
#    11 => (153.392,  82.580),
#    12 => (168.874,  83.373),
#    13 => (184.227,  82.875),
#    14 => (199.741,  83.098),
#    15 => (215.881, 101.271),
#    16 => (230.531, 104.860),
#    17 => (246.576, 103.229),
#    18 => (261.412, 102.497),
#    19 => (276.926, 103.367),
#    20 => (294.404, 104.162),
#    21 => (308.995, 104.944),
#    22 => (323.287, 104.838),
#    23 => (340.173, 121.495),
#    24 => (353.717, 119.729),
#    25 => (369.113, 120.281),
#    26 => (384.638, 118.710),
#    27 => (402.484, 119.237),
#    28 => (417.036, 119.980),
#    29 => (431.273, 120.463),
#    30 => (446.288, 120.560),
#    31 => (461.756, 137.209),
#    32 => (476.653, 135.873),
#    33 => (493.989, 138.150),
#    34 => (508.089, 137.998),
#    35 => (523.703, 139.059),
#    36 => (539.305, 135.201),
#    37 => (554.201, 138.242),
#    38 => (569.276, 137.470),
#    39 => (585.268, 137.261),
#    40 => (599.921, 137.871),
#    41 => (615.718, 137.887),
#    42 => (632.463, 137.579),
#    43 => (647.355, 138.181),
#    44 => (664.344, 139.917),
#    45 => (678.899, 138.382),
#    46 => (695.815, 138.354),
#    47 => (717.068, 154.826),
#    48 => (730.902, 155.603),
#    49 => (745.608, 154.549),
#    50 => (756.012, 154.244),
#    51 => (769.611, 154.000),
#    52 => (793.350, 152.893),
#    53 => (801.568, 155.276),
#    54 => (818.174, 153.436),
#    55 => (832.994, 155.923),
#    56 => (849.526, 155.721),
#    57 => (863.337, 155.289),
#    58 => (881.577, 156.188),
#    59 => (893.735, 157.889),
#    60 => (911.597, 155.557),
#    61 => (925.117, 157.048),
#    62 => (940.897, 153.338),
#    63 => (954.825, 173.942),
#    64 => (990.757, 173.236),
#    65 => (987.457, 172.005),
#    66 => (1009.688, 173.869),
#    67 => (1022.357, 173.303),
#    68 => (1043.536, 175.044),
#    69 => (1052.156, 173.361),
#    70 => (1067.174, 172.230),
#    71 => (1079.543, 173.491),
#    72 => (1101.516, 174.311),
#    73 => (1119.993, 174.057),
#    74 => (1136.092, 175.804),
#    75 => (1154.758, 175.075),
#    76 => (1170.732, 173.453),
#    77 => (1184.598, 171.936),
#    78 => (1190.922, 172.809),
#    79 => (1204.645, 172.878),
#)

speedup(p::Integer) = ((p-1)/(floor(Int, log2(p)) + 1 + floor(Int, log2(p/3))))
#for (np, (t_ser, t_par)) in sort!([x for x in timings])
#    @printf("np: %3d Serial: %.3f sec  Parallel: %.3f sec  speedup: %.3fx (theory=%.3fx)\n",
#        np, t_ser, t_par, t_ser/t_par, speedup(np))
#end

xend = 1+maximum(keys(timings))

theplot=plot(layer(x=2:xend, y=map(speedup, 2:xend), Geom.line(), Theme(default_color=colorant"darkkhaki")),
layer(x=[np+1 for (np, (t_ser, t_par)) in timings],
      y=[t_ser/t_par for (np, (t_ser, t_par)) in timings],
    Geom.point, Theme(default_color=colorant"firebrick",
        #=default_point_size=1.25pt, =#highlight_width=0pt)),
    Guide.xlabel("Number of processors, p"), Guide.ylabel("Speedup, r"))

labeledplot = compose(render(theplot),
compose(context(), stroke(colorant"purple"), #Powers of 2
text(0.21, 0.73, "2"), text(0.24, 0.72, "4"), text(0.27, 0.7, "8"), text(0.33, 0.63, "16"),
text(0.47, 0.53, "32"), text(0.76, 0.35, "64")),
compose(context(), stroke(colorant"green"), #3x powers of 2
text(0.25, 0.62, "6"), text(0.28, 0.57, "12"), text(0.39, 0.49, "24"), text(0.63, 0.32, "48"),
)
)
draw(SVG(3.5inch, 3.2inch), labeledplot)

draw(SVG("scaling.svg", 3.5inch, 3.2inch), labeledplot)
#draw(PDF("scaling.pdf", 3.5inch, 3inch), theplot)

# ------------------------------------------------------------------------------------------
# The exact same serial code now runs in parallel thanks to multiple dispatch!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Formal verification
#
# Julia allows us to implement very easily the interval of summations monoid technique for
# verifying the correctness of prefix sums
# ([doi:10.1145/2535838.2535882](http://dx.doi.org/10.1145/2535838.2535882)).
# ------------------------------------------------------------------------------------------

#Definition 4.3 of Chong et al., 2014
abstract type ⊤ end
abstract type Id end
const 𝕊 = Union{UnitRange, Type{⊤}, Type{Id}}

⊕(I::UnitRange, J::UnitRange) = I.stop+1==J.start ? (I.start:J.stop) : ⊤
⊕(::Type{Id}, ::Type{Id}) = Id
⊕(I::𝕊, ::Type{Id}) = I
⊕(::Type{Id}, I::𝕊) = I
⊕(I::𝕊, J::𝕊) = ⊤

for kernel in (prefix_serial!, prefix!), n=1:1000
    @assert kernel([k:k for k=1:n], ⊕) == [1:k for k=1:n]
end

# ------------------------------------------------------------------------------------------
# Wrong kernels produce type errors!
# ------------------------------------------------------------------------------------------

function prefix_wrong!(y, .+)
    l = length(y)
    offset = 1
    z = copy(y)
    while offset < l
        for i=offset+1:l
            z[i] = y[i-offset] .+ z[i]
        end
        offset += 2
        y = copy(z)
    end
    y
end

prefix_wrong!(Any[k:k for k=1:8], ⊕) #Should produce conversion error

@manipulate for npp=1:10
    render(prefix_wrong!(AccessArray(npp),+))
end

# ------------------------------------------------------------------------------------------
# # Experimental code which doesn't work
#
# In this section I play with other parallel prefix algorithms. All the ones here suffer
# from read antidependencies and they break the logic I used to draw the gates.
# ------------------------------------------------------------------------------------------

#Kogge-Stone variant
#Works but doesn't render properly
function prefix_ks!(y, .+)
    l = length(y)
    offset = 1
    z = copy(y)
    while offset < l
        for i=offset+1:l
            z[i] = y[i-offset] .+ z[i]
        end
        offset *= 2
        y = copy(z)
    end
    y
end

@assert prefix_ks!(ones(8), +) == [1:8;]

Base.copy(A::AccessArray) = A

render(prefix_ks!(AccessArray(8),+))

#Hillis-Steele = Kogge-Stone
function prefix_hs!(y, +)
    l = length(y)
    lk = floor(Int, log2(l))
    for d=1:lk
        for k=2d:l
            y[k] = y[k-2^(d-1)] + y[k]
        end
    end
    y
end
prefix_hs!(ones(4), +)

render(prefix_hs!(AccessArray(4), +))

#Sklansky
function prefix_s!(y, +)
    n = length(y)
    for d=1:log2(n)
        for k=1:n÷2+1
            block = 2 * (k - (mod(k, 2^d)))
            me = block + mod(k, 2^d) + 2^d
            spine = block + 2^d - 1
            if spine ≤ n && me ≤ n
                y[me] = y[me] + y[spine]
            end
        end
    end
    y
end
prefix_s!(ones(4), +)

function prefix(v,⊕)
    @show v
    k=length(v)
    k≤1 && return v
    v[2:2:k] = @show v[1:2:k-1]⊕v[2:2:k]  # Line 1
    v[2:2:k] = @show prefix(v[2:2:k],⊕)   # Line 2
    v[3:2:k] = @show v[2:2:k-1]⊕v[3:2:k]  # Line 3
    v
end

prefix(ones(8),+)

function ⊕{T<:Range}(v1::Vector{T}, v2::Vector{T})
    for v in (v1, v2), i in 1:length(v)-1
        @assert v[i].stop + 1 == v[i+1].start
    end
    if length(v1)>0 && length(v2)>0
        @assert v1[end].stop + 1 == v2[1].start
        v1[1].start:v2[end].stop
    elseif length(v1)>0 #&& length(v2)==0
        v1[1].start:v1[end].stop
    elseif length(v2)>0 #&& length(v1)==0
        v2[1].start:v2[end].stop
    else #Both are empty
        v1
    end
end

for kernel in (prefix,), n=1:100#(prefix_serial!, prefix!, prefix_ks!), n=1:100
    @assert kernel([k:k for k=1:n], ⊕) == [1:k for k=1:n]
end

max_stars_repo_path: .nbexports/introductory-tutorials/intro-to-julia/Reduce and Parallel Prefix - JuliaBox.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
# --------------------------------------------------------------------------
# ACE1.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

module Transforms

import Base:   ==
import JuLIP:  cutoff, AtomicNumber
import JuLIP.FIO: read_dict, write_dict
import ACE1: _allfieldsequal

abstract type DistanceTransform end

export PolyTransform, IdTransform, MorseTransform, AgnesiTransform

# fall-backs from species dependent to independent transforms
transform(t::DistanceTransform, r::Number, z::AtomicNumber, z0::AtomicNumber) = transform(t, r)
transform_d(t::DistanceTransform, r::Number, z::AtomicNumber, z0::AtomicNumber) = transform_d(t, r)
inv_transform(t::DistanceTransform, x::Number, z::AtomicNumber, z0::AtomicNumber) = inv_transform(t, x)

poly_trans(p, r0, r) = @fastmath(((1+r0)/(1+r))^p)

poly_trans_d(p, r0, r) = @fastmath((-p/(1+r0)) * ((1+r0)/(1+r))^(p+1))

poly_trans_inv(p, r0, x) = ( (1+r0)/(x^(1/p)) - 1 )

# TODO: generalise the distance transform to allow
#       ((a + r0) / (a + r) )^p

@doc raw"""
Implements the distance transform
```math
   x(r) = \Big(\frac{1 + r_0}{1 + r}\Big)^p
```

Constructor:
```
PolyTransform(p, r0)
```
"""
struct PolyTransform{TP, T} <: DistanceTransform
   p::TP
   r0::T
end

PolyTransform(; p = 2, r0 = 2.5) = PolyTransform(p, r0)

write_dict(T::PolyTransform) =
   Dict("__id__" => "ACE1_PolyTransform", "p" => T.p, "r0" => T.r0)

PolyTransform(D::Dict) = PolyTransform(D["p"], D["r0"])

read_dict(::Val{:ACE1_PolyTransform}, D::Dict) = PolyTransform(D)

transform(t::PolyTransform, r::Number) = poly_trans(t.p, t.r0, r)

transform_d(t::PolyTransform, r::Number) = poly_trans_d(t.p, t.r0, r)

inv_transform(t::PolyTransform, x::Number) = poly_trans_inv(t.p, t.r0, x)

(t::PolyTransform)(x) = transform(t, x)

"""
`IdTransform`: Implements the distance transform `z -> z`;
Primarily used for the z-coordinate for the EnvPairPots

Constructor: `IdTransform()`
"""
struct IdTransform <: DistanceTransform
end

write_dict(T::IdTransform) =  Dict("__id__" => "ACE1_IdTransform")
IdTransform(D::Dict) = IdTransform()
read_dict(::Val{:ACE1_IdTransform}, D::Dict) = IdTransform(D)
transform(t::IdTransform, z::Number) = z
transform_d(t::IdTransform, r::Number) = one(r)
inv_transform(t::IdTransform, x::Number) = x

@doc raw"""
Implements the distance transform
```math
   x(r) = \exp( - \lambda (r/r_0) )
```

Constructor:
```
MorseTransform(lambda, r0)
```
"""
struct MorseTransform{T} <: DistanceTransform
   lambda::T
   r0::T
end

write_dict(T::MorseTransform) =
   Dict("__id__" => "ACE1_MorseTransform", "lambda" => T.p, "r0" => T.r0)
MorseTransform(D::Dict) = MorseTransform(D["lambda"], D["r0"])
read_dict(::Val{:ACE1_MorseTransform}, D::Dict) = MorseTransform(D)
transform(t::MorseTransform, r::Number) = exp(- t.lambda * (r/t.r0))
transform_d(t::MorseTransform, r::Number) = (-t.lambda/t.r0) * exp(- t.lambda * (r/t.r0))
inv_transform(t::MorseTransform, x::Number) = - t.r0/t.lambda * log(x)
(t::MorseTransform)(x) = transform(t, x)

@doc raw"""
Implements the distance transform
```math
   x(r) = \frac{1}{1 + a (r/r_0)^p}
```
with default $a = (p-1)/(p+1)$. That default is chosen such that
$|x'(r)|$ is maximised at $r = r_0$. Default for $p$ is $p = 2$. Any value
$p > 1$ is permitted.

Constructor:
```
AgnesiTransform(r0, [p, [, a]])
```
"""
struct AgnesiTransform{T, TP} <: DistanceTransform
   r0::T
   p::TP
   a::T
   c::T
end

AgnesiTransform(r0, p=2, a=(p-1)/(p+1)) = (@assert p > 1;
                             AgnesiTransform(r0, p, a, - a * p / r0))

write_dict(T::AgnesiTransform) =
Dict("__id__" => "ACE1_AgnesiTransform", "r0" => T.r0, "p" => T.p, "a" => T.a)
AgnesiTransform(D::Dict) = AgnesiTransform(D["r0"], D["p"], D["a"])
read_dict(::Val{:ACE1_AgnesiTransform}, D::Dict) = AgnesiTransform(D)
transform(t::AgnesiTransform, r::Number) = @fastmath 1 / (1 + t.a * (r/t.r0)^t.p)
transform_d(t::AgnesiTransform, r::Number) = (s1 = (r/t.r0); s2 = s1^(t.p-1);
                                         @fastmath t.c * s2 / (1+t.a * s2*s1)^2)
inv_transform(t::AgnesiTransform, x::Number) = t.r0 * ( (1/x-1)/t.a )^(1/t.p)
(t::AgnesiTransform)(x) = transform(t, x)

# --------- AnalyticTransform

import JuLIP
import JuLIP.Potentials: ScalarFun

"""
`AnalyticTransform`: implements a distance transform that can be specified
by an analytic expression.

Constructor:
```julia
AnalyticTransform(forwardmap, inversmap)
```
For `forwardmap` and `inversemap` must both be of type `AnalyticFunction`.
(cf. `JuLIP.@analytic`).

Example:
```julia
using ACE1: AnalyticTransform
using JuLIP: @analytic
trans = AnalyticTransform( :(r -> exp( - 2 * r )),
                           :(x -> -0.5 * log(x)) )
```
"""
struct AnalyticTransform{T} <: DistanceTransform
   f::ScalarFun{T}
   df::ScalarFun{T}
   finv::ScalarFun{T}
   str_f::String
   str_finv::String
end

==(T1::AnalyticTransform, T2::AnalyticTransform) = T1.str_f == T2.str_f

Base.show(io::IO, trans::AnalyticTransform) =
      print(io, "AnalyticTransform($(trans.str_f))")

function AnalyticTransform(str_f::String, str_finv::String; T=Float64)
   ex_f = Meta.parse(str_f)
   ex_df = JuLIP.Potentials.fdiff( ex_f, 1 )
   ex_finv = Meta.parse(str_finv)
   f = ScalarFun{T}(Meta.eval(ex_f))
   df = ScalarFun{T}(Meta.eval(ex_df))
   finv = ScalarFun{T}(Meta.eval(ex_finv))
   return AnalyticTransform(f, df, finv, str_f, str_finv)
end

write_dict(T::AnalyticTransform) =
         Dict("__id__" => "ACE1_AnalyticTransform",
              "f" => T.str_f, "finv" => T.str_finv)
AnalyticTransform(D::Dict) = AnalyticTransform(D["f"], D["finv"])
read_dict(::Val{:ACE1_AnalyticTransform}, D::Dict) = AnalyticTransform(D)

transform(t::AnalyticTransform, r::Number) = t.f(r)
transform_d(t::AnalyticTransform, r::Number) = t.df(r)
inv_transform(t::AnalyticTransform, x::Number) = t.finv(x)
(t::AnalyticTransform)(r) = transform(t, r)

import ForwardDiff
transform(t::AnalyticTransform, r::ForwardDiff.Dual) = t.f.obj.x(r)

# --------- Utility function - affine transform

"""
`AffineT` : wraps another transform and then applies an affine transformation.

Constructor:
```julia
AffineT(transform, x1, x2, y1, y2)
```
then `x = t(r)` and then `x -> y` with `xi -> yi`.
"""
struct AffineT{T, TT} <: DistanceTransform
   t::TT   # the inner transform  : t(r) = x
   x1::T   # intervals to be transformed x1 -> y1 etc...
   x2::T
   y1::T
   y2::T
end

transform(t::AffineT, r) = t.y1 + (transform(t.t, r) - t.x1) * (t.y2-t.y1)/(t.x2-t.x1)
transform_d(t::AffineT, r) = ((t.y2-t.y1)/(t.x2-t.x1)) * transform_d(t.t, r)
inv_transform(t::AffineT, y) = inv_transform(t.t, t.x1 + (y - t.y1) * (t.x2-t.x1)/(t.y2-t.y1))

==(t1::AffineT, t2::AffineT) = _allfieldsequal(t1, t2)

write_dict(T::AffineT) =
      Dict("__id__" => "ACE1_AffineT",
           "t" => write_dict(T.t),
           "xy" => [T.x1, T.x2, T.y1, T.y2] )

read_dict(::Val{:ACE1_AffineT}, D::Dict) = AffineT(D)

AffineT(D::Dict) = AffineT(read_dict(D["t"]),
                         D["xy"]... )

# --------- Multi-transform: species-dependent transform

import JuLIP: chemical_symbol
import JuLIP.Potentials: ZList, SZList, i2z, z2i
using StaticArrays: SMatrix
struct MultiTransform{NZ, TT} <: DistanceTransform
   zlist::SZList{NZ}
   transforms::SMatrix{NZ, NZ, TT}
end

cutoff_extrema(T::MultiTransform) =
   minimum( inv_transform(t, -1.0) for t in T.transforms ),
   maximum( inv_transform(t,  1.0) for t in T.transforms )

# FIO

==(T1::MultiTransform, T2::MultiTransform) =  _allfieldsequal(T1, T2)

write_dict(T::MultiTransform) =
      Dict("__id__" => "ACE1_MultiTransform",
           "zlist" => write_dict(T.zlist),
           "transforms" => write_dict.(T.transforms[:]))

read_dict(::Val{:ACE1_MultiTransform}, D::Dict) = MultiTransform(D)

function MultiTransform(D::Dict)
   zlist = read_dict(D["zlist"])
   NZ = length(zlist)
   transforms = SMatrix{NZ, NZ}( read_dict.(D["transforms"])... )
   return MultiTransform(zlist, transforms)
end

#  Constructor

function multitransform(D::Dict; rin=nothing, rcut=nothing, cutoffs = nothing)

   if rin != nothing && rcut != nothing && cutoffs == nothing
      cutoffs = Dict( [key => (rin, rcut) for key in keys(D)]... )
   end

   species = Symbol[]
   for key in keys(D)
      append!(species, [key...])
   end
   species = unique(species)
   zlist = ZList(species, static=true)
   NZ = length(zlist)
   transforms = Matrix{Any}(undef, NZ, NZ)
   for i = 1:NZ, j = 1:NZ
      Si = chemical_symbol(i2z(zlist, i))
      Sj = chemical_symbol(i2z(zlist, j))
      if haskey(D, (Si, Sj))
         key = (Si, Sj)
      else
         key = (Sj, Si)
      end
      # get the transform from the dict
      t = D[key]
      # apply another affine transform so all transforms have the same range
      if cutoffs != nothing
         rin, rcut = cutoffs[key]
         x1 = transform(t, rin)
         x2 = transform(t, rcut)
         transforms[i, j] = AffineT(t, x1, x2, -1.0, 1.0)
      else
         transforms[i, j] = t
      end
   end
   transforms = identity.(transforms)  # infer the type
   return MultiTransform(zlist, SMatrix{NZ, NZ}(transforms...))
end

transform(t::MultiTransform, r::Number, z::AtomicNumber, z0::AtomicNumber) =
      transform(t.transforms[z2i(t.zlist, z), z2i(t.zlist, z0)], r)

transform_d(t::MultiTransform, r::Number, z::AtomicNumber, z0::AtomicNumber) =
      transform_d(t.transforms[z2i(t.zlist, z), z2i(t.zlist, z0)], r)

inv_transform(t::MultiTransform, x::Number, z::AtomicNumber, z0::AtomicNumber) =
      inv_transform(t.transforms[z2i(t.zlist, z), z2i(t.zlist, z0)], r)

# # NOTE: This is a bit of a hack, I'm checking whether
# #       the transforms have been transformed to the domain [-1, 1]
# #       then I'm checking whether r is either rin or rcut and only
# #       then do I return the value
# function transform(t::MultiTransform{NZ, TT}, r::Number) where {NZ, TT}
#    @assert (TT <: AffineT) "transform(::MultiTransfrom, r) is only defined if rin, rcut are specified during construction"
#    x = sum( transform(_t, r)  for _t in t.transforms ) / NZ^2
#    @assert (abs(abs(x) - 1) <= 1e-7) "transform(::MultiTransfrom, r) is only defined for r = rin, rcut"
#    return x
# end

end


text:
module Pkg3

include("Types.jl")
include("Operations.jl")
include("Loading.jl")

import .Types: VersionSpec, @vs_str
import .Loading: Loader, LoadInstalled
import .Operations:
    add, find_config, find_manifest, user_depot, depots, registries

export @vs_str

end # module


text:
"""Validates the contents of a Storage."""
function validate(storage::Storage)::Bool
    is_valid = true

    @debug "Storage validation" storage.name is_valid
    return is_valid
end


text:
using FractionalDiffEq
using Plots, LaTeXStrings
using MittagLeffler

#Analytical solution
target = []

#MittagLeffler.jl doesn't support array argument
for i in 0:0.01:20
    push!(target, i^1.8*mittleff(1.8,2.8,-i^1.8))
end

s="\$D^{1.8}y(x)=1-y(x),\\ y(0)=0\$"

#Numerical solution
fun(x, y) = 1-y
prob = FDEProblem(fun, 1.8, 0, 20, 0.01)
result=solve(prob, PECE())
tspan=collect(0:0.01:20)

gr()

plot(tspan, result, title=s, linewidth=5, label="Numerical", legend=:bottomright)
plot!(tspan, target, lw=3, ls=:dash, label="Analytical")


text:
using PGFPlotsX
using CSV
using FinEtools

results = [
(1.00000e-02, Any[2.52679e-07, 2.52672e-07, 2.52662e-07, 2.52651e-07, 2.52671e-07, 2.52651e-07, 1.90903e-07, 1.34989e-07])
(2.00000e-03, Any[2.52633e-07, 2.52675e-07, 2.52699e-07, 2.52659e-07, 2.41211e-07, 2.15746e-07, 1.52555e-07, 1.07873e-07])
(1.00000e-03, Any[2.24037e-07, 2.24064e-07, 2.24032e-07, 2.24037e-07, 2.06195e-07, 1.84427e-07, 1.30409e-07, 9.22134e-08])
(1.00000e-04, Any[1.93031e-07, 1.93028e-07, 1.93035e-07, 1.93027e-07, 1.84096e-07, 1.64660e-07, 1.16432e-07, 8.23301e-08])
(1.00000e-05, Any[1.92667e-07, 1.92665e-07, 1.92671e-07, 1.92664e-07, 1.83822e-07, 1.64415e-07, 1.16259e-07, 8.22077e-08])
]

drilling_stiffness_scales = [0.01, 0.1, 1.0, 1.5, 2.0, 2.5, 5.0, 10.0]

objects = []
marks = ["x", "o", "square", "+", "diamond"]

for i in eachindex(results)
    xdata, ydata = drilling_stiffness_scales, results[i][2] ./ maximum(results[i][2]);
    @pgf p = PGFPlotsX.Plot(
    {
    color = "black",
    mark = "$(marks[i])",
    line_width  = 1.0,
    },
    Coordinates([v for v in  zip(xdata, ydata)])
    )
    push!(objects, p)
    push!(objects, LegendEntry("t/R=$(round(results[i][1]/0.1, digits=4))"))
end

@pgf ax = Axis(
{
xlabel = "Drilling Stiffness Multiplier [ND]",
ylabel = "Normalized time step [ND]",
# ymin = -1.1e-5,
# ymax = 1.1e-5,
# xmin = 0.0,
# xmax = 1.0,
xmode = "log",
ymode = "linear",
yminorgrids = "true",
grid = "both",
legend_style = {
at = Coordinate(0.35, 0.05),
anchor = "south",
# legend_columns = -1
},
},
objects...
)

display(ax)
pgfsave("clamped_cylinder_vibration-time_step.pdf", ax)


text:
lines = readlines("d03.in")
const height, width = size(lines, 1), length(lines[1])

function trees(right, down)
    x, y, tree = 0, 0, 0
    while y < height
        if lines[y+1][x+1] == '#'
            tree += 1
        end
        y += down
        x = (x + right) % width
    end
    return tree
end

P1 = trees(3, 1)
P2 = P1 * trees(1, 1) * trees(5, 1) * trees(7, 1) * trees(1, 2)

println("P1: ", P1)
println("P2: ", P2)

max_stars_repo_name: eduellery/adventofcode

max_issues_repo_name: eduellery/adventofcode

max_forks_repo_name: eduellery/adventofcode


text:
module Sbatch
export Cluster, CaltechCentral
export sbatch, modules

abstract type Cluster end
struct CaltechCentral <: Cluster end

function sbatch end
function modules end

function modules(::CaltechCentral)
    ["julia/1.5.4 ", "cuda/10.2 ", "openmpi/4.0.4_cuda-10.2 "]
end

function sbatch(
    cluster::Cluster,
    script::String,
    project_dir = "",
    wait = true;
    exclusive = true,
    time = Time(0, 30, 0),
    nodes = 1,
    ntasks = 1,
    cpus_per_task = 1,
    gpu = 0,
    kwargs...,
)
    project = isempty(project_dir) ? "" : "--project=$(project_dir)"
    args = String[]
    exclusive && push!(args, "--exclusive")
    append!(args, [
        "--export=ALL",
        "--time=$(string(time))",
        "--nodes=$(nodes)",
        "--ntasks=$(ntasks)",
        "--cpus-per-task=$(cpus_per_task)",
    ])
    if gpu > 0
        push!(args, "--gres=gpu:$(gpu)")
        extra_arg = ""
    else
        extra_arg = "--disable-gpu"
    end
    append!(args, ["--$kw=$val" for (kw, val) in kwargs])

    ss, io = mktemp()
    println(io, "set -euo pipefail")
    println(io, "set -x")
    println(io, "export TICTOC_PRINT_RESULTS=1")
    println(io, "export OPENBLAS_NUM_THREADS=1")
    println(io, "module load $(modules(cluster)...)")
    println(io, "time mpiexec julia --color=no $(project) $(script) $(extra_arg)")
    close(io)

    sbatchcmd = `echo sbatch`
    append!(sbatchcmd.exec, args)
    append!(sbatchcmd.exec, `$ss`)

    # XXX: debug
    for wd in sbatchcmd.exec
        print(wd, " ")
    end
    println()

    _, jobid = rsplit(chomp(String(read(sbatchcmd))), limit = 2)
    return String(jobid)
end

end


text:
using Makie
using Bridge
using Bridge: increment
using StaticArrays
using Trajectories
using Colors
import Trajectories: Trajectory
using Random
Trajectory(X::SamplePath) = Trajectory(X.tt, X.yy)
Wnr = Wiener()
P = P̂
R = 𝕏(1.5,6.0)
rebirth(α, R) = x -> (rand() > α  ? x : (2rand(typeof(x)) .- 1).*R)

k = 250
n = 250
pts = [rebirth(1.0, R)(x0[1]) for i in 1:n]
xraw = [copy(pts) for i in 1:k]
x = [Node(xraw[i]) for i in 1:k]
col = [Node(RGBA{Float32}(0.0, 0.0, 0.0, 0.0)) for i in 1:k]
c = 1
ms = 0.08
i = 1;

limits = FRect(-R[1], -R[2], 2R[1], 2R[2])
p = Scene(resolution=(800,800), limits=limits, backgroundcolor = RGB{Float32}(0.04, 0.11, 0.22))
for i in randperm(k)
    scatter!(p, x[i], color = col[i], markersize = ms, #=show_axis = true,limits=limits,=#  glowwidth = 0.005, glowcolor = :white)
end
#update_cam!(p, limits)
axis = p[Axis]
axis[:grid, :linewidth] =  (1, 1)
axis[:grid, :linecolor] = (RGBA{Float32}(0.5, 0.7, 1.0, 0.5),RGBA{Float32}(0.5, 0.7, 1.0, 0.5))
axis[:names][:textsize] = (0.0,0.0)
axis[:ticks, :textcolor] = (RGBA{Float32}(0.5, 0.7, 1.0, 0.8),RGBA{Float32}(0.5, 0.7, 1.0, 0.8))
r = range(-R[1], 2R[1], length=100)

display(p)

function update!(t, x, y, dt, P, W, jump)
    for i in eachindex(x)
        y[i] = jump(x[i])
        y[i] = y[i] + BSI.b(t, y[i], P)*dt + BSI.σ(t, y[i], P)*rand(increment(dt, W))
    end
    y
end
dtf = dt
update!(x, y, dt, P, W, jump = indentity) = update!(NaN, x, y, dt, P, W, jump)
sleep(1)
N = 4000
jump = rebirth(0.0001, R)
for i in 1:N
    global c
    cnew = mod1(c+1, k)
    update!(x[c][], x[cnew][], dtf, P, Wnr, jump)
    c = cnew
    x[c][] = x[cnew][]

    for i in 0:k-1
        f = (i + 5*(i!=0))/(k+5)
        col[mod1(c-i, k)][] = RGBA{Float32}(0.5, 0.7, 1.0, (1-sqrt(f))/3)
    end
    sleep(1e-10)
end