text:
function compInit2WayGains(bins::Vector{Int}, fixed::Vector{Int}, H::Hypergraph, incidence::Incidence)
    e = H.e
    n = H.n
    netDegs = [zeros(Int, 2) for i in 1:e]
    gains = zeros(Int, n)
    initEdgeDegs!(H, bins, netDegs)
    initNodeGains!(H, fixed, bins, gains, netDegs, incidence)

    return (netDegs, gains)
end

function initEdgeDegs!(H::Hypergraph, bins::Vector{Int}, netDegs::Vector)
    partCost = 0
    edgeDegree = 0
    n = H.n
    e = H.e
    hedges = H.hedges
    eptr = H.eptr

    for i in 1:e
        u_loc = eptr[i]
        v_loc = eptr[i+1]
        for j in u_loc:v_loc-1
            v = hedges[j]
            netDegs[i][bins[v]+1] += 1
        end
    end
end

function initNodeGains!(H::Hypergraph, fixed::Vector{Int}, bins::Vector{Int}, gains::Vector{Int}, netDegs::Vector, incidence::Incidence)
    n_gain = 0
    eGain = 0
    pId = 0
    oId = 0
    fromDeg = 0
    toDeg = 0
    n = H.n
    i_hedges = incidence.hedges
    i_eptr = incidence.eptr
    h_hedges = H.hedges
    h_eptr = H.eptr

    if H.weighted == true
        w_ = H.w_
    else
        w_ = ones(Int, H.e)
    end

    for i in 1:n
        if fixed[i] == 0
            pId = bins[i]
            oId = pId == 0 ? 1 : 0
            u_loc = i_eptr[i]
            v_loc = i_eptr[i+1]

            for j in u_loc:v_loc-1
                h = i_hedges[j]
                h_size = h_eptr[h+1]-h_eptr[h]
                eWt = w_[h]
                #=if h_size > _FM_max_fanout
                    continue
                end=#
                fromDeg = netDegs[h][pId+1]
                toDeg = netDegs[h][oId+1]
                if toDeg == 0 && fromDeg > 0
                    eGain = -eWt
                elseif fromDeg == 1 && toDeg >= 0
                    eGain = eWt
                else
                    eGain = 0
                end
                n_gain += eGain
            end
            gains[i] = n_gain
            n_gain = 0
        end
    end
end

function updAfterMove!(gains::Vector{Int}, fixed::Vector{Int}, moved::Vector{Int}, marked::Vector{Int},
                       netDegs::Vector{Vector{Int}}, nIdx::Int, bins::Vector{Int}, H::Hypergraph, incidence::Incidence, order::Int, pass::Int)
    n = H.n
    e = H.e
    i_hedges = incidence.hedges
    i_eptr = incidence.eptr
    h_hedges = H.hedges
    h_eptr = H.eptr
    w_ = H.w_

    if H.weighted == false
        w_ = ones(Int, H.e)
    end

    moved[nIdx] = pass
    toPId = bins[nIdx]
    fromPId = toPId == 0 ? 1 : 0
    gains[nIdx] = -gains[nIdx]
    nbrs = Vector{Int}()

    for i in i_eptr[nIdx]:i_eptr[nIdx+1]-1
        h = i_hedges[i]
        eWt = w_[h]
        u_loc = h_eptr[h]
        v_loc = h_eptr[h+1]
        h_size = v_loc-u_loc
        netDegs[h][toPId+1] += 1
        netDegs[h][fromPId+1] -= 1

        for j in u_loc:v_loc-1
            v = h_hedges[j]
            if v == nIdx || fixed[v] > 0
                continue
            end
            nodePId = bins[v]
            fromDeg = netDegs[h][fromPId+1]
            toDeg = netDegs[h][toPId+1]

            if nodePId == fromPId
                if fromDeg == 1
                    gains[v] += eWt
                end
                if toDeg == 1
                    gains[v] += eWt
                end
            else
                if fromDeg == 0
                    gains[v] -= eWt
                end
                if toDeg == 2
                    gains[v] -= eWt
                end
            end

            if moved[v] != pass && marked[v] != order
                push!(nbrs, v)
                marked[v] = order
            end
        end
    end

    return nbrs
end

function updateGainHeap!(perm::Vector{Int}, bins::Vector{Int}, gainHeaps::Vector, gainHeapLocators::Vector, numNodes::Vector{Int}, gains::Vector{Int})
    nId = Int(0)
    pId = Int(0)
    heaploc = Int(0)

    for pIdx in 1:length(perm)
        nId = perm[pIdx]
        pId = bins[nId]
        push!(gainHeaps[pId+1], [nId, gains[nId]])
        numNodes[pId+1] += 1
        heaploc = numNodes[pId+1]
        gainHeapLocators[nId] = [pId+1, numNodes[pId+1]]
        heapifyUp!(heaploc, gainHeapLocators, gainHeaps, gains[nId], pId+1)
    end
end

function rollBack!(minCostOrder::Int, fixed::Vector{Int}, binsArea::Array{Int}, gains::Vector{Int}, swaps::Vector{Int}, bins::Vector{Int}, area::Array{Int},
                   moved::Vector{Int}, marked::Vector{Int}, netDegs::Vector, H::Hypergraph, incidence::Incidence, pass::Int, order::Int, cut_size::Int)
    node = 0
    fromId = 0
    toId = 0
    nbrs = Vector{Int}()

    @inbounds for oIdx in minCostOrder+1:length(swaps)
        node = swaps[oIdx]
        fromId = bins[node]
        toId = fromId == 0 ? 1 : 0
        binsArea[:, fromId+1] -= area[node, :]
        binsArea[:, toId+1] +=  area[node, :]
        bins[node] = toId
        cut_size += gains[node]
        nbrs = updAfterMove!(gains, fixed, moved, marked, netDegs, node, bins, H, incidence, order, pass)
    end
end

function updateNbrs!(nbrs::Vector{Int}, gains::Vector{Int}, gainHeapLocators::Vector, gainHeaps::Vector, numNodes::Vector{Int}, bins::Vector{Int})
    for nbrIdx in 1:length(nbrs)
        nbr = nbrs[nbrIdx]
        newGain = gains[nbr]
        locT = gainHeapLocators[nbr]
        oldGain = gainHeaps[locT[1]][locT[2]][2]
        heapify!(nbr, numNodes, newGain, oldGain, bins, gainHeapLocators, gainHeaps)
    end
end

function heapify!(node::Int, numNodes::Vector{Int}, newGain::Int, oldGain::Int, bins::Vector{Int}, gainHeapLocators::Vector, gainHeaps::Vector)
    partId = bins[node]+1
    heapLoc = gainHeapLocators[node][2]
    gainHeaps[partId][heapLoc][2] = newGain

    if oldGain > newGain
        heapifyDown!(heapLoc, numNodes[partId], gainHeapLocators, gainHeaps, newGain, partId)
    elseif oldGain < newGain
        heapifyUp!(heapLoc, gainHeapLocators, gainHeaps, newGain, partId)
    end
end

function heapifyUp!(heapLoc::Int, gainHeapLocators::Vector, gainHeaps::Vector, newGain::Int, partId::Int)
    parent = 0
    gainHeap = Vector()
    h_node = 0
    child = zeros(Int, 2)

    while heapLoc > 1
        parent = heapLoc >> 1
        if newGain > gainHeaps[partId][parent][2]
            child = gainHeaps[partId][heapLoc]
            gainHeaps[partId][heapLoc] = gainHeaps[partId][parent]
            h_node = gainHeaps[partId][parent][1]
            gainHeapLocators[h_node] = [partId, heapLoc]
            gainHeaps[partId][parent] = child
            gainHeapLocators[child[1]] = [partId, parent]
            heapLoc = parent
        else
            break
        end
    end
end

function heapifyDown!(heaploc::Int, numNodes::Int, gainHeapLocators::Vector, gainHeaps::Vector, newGain::Int, partId::Int)
    j = heaploc << 1
    gain_l = Vector()
    gain_r = Vector()
    parent = zeros(Int ,2)
    child = zeros(Int, 2)

    while j <= numNodes
        gain_l = gainHeaps[partId][j]

        if j+1 <= numNodes
            gain_r = gainHeaps[partId][j+1]
            if gain_l[2] < gain_r[2]
                j += 1
            end
        else
            if newGain > gain_l[2]
                break
            end
        end

        if newGain > gain_l[2] && newGain > gain_r[2]
            break
        end

        parent = gainHeaps[partId][heaploc]
        child = gainHeaps[partId][j]
        gainHeaps[partId][heaploc] = child
        gainHeapLocators[child[1]] = [partId, heaploc]
        gainHeaps[partId][j] = parent
        gainHeapLocators[parent[1]] = [partId, j]
        heaploc = j
        j = heaploc << 1
    end
end

function bestNodeToMove(gainHeaps::Vector, gainHeapLocators::Vector, binsArea::Array{Int}, area::Array{Int}, numNodes::Vector{Int}, target_area::Array{Int})
    ptr = 1
    best_node = 0
    best_side = 0

    if numNodes[1] == 0 && numNodes[2] == 0
        return (-1, -1)

    elseif numNodes[1] == 0 && numNodes[2] > 0
        node = gainHeaps[2][ptr][1]
        done = binsArea[:, 1] .+ area[node] <= target_area[:, 1]

        while done == false
            if ptr < numNodes[2]
                ptr += 1
                done = binsArea[:, 1] .+ area[node] <= target_area[:, 1]
            else
                return (-1, -1)
            end
        end

        best_node = gainHeaps[2][ptr][1]
        best_side = 1

    elseif numNodes[2] == 0 && numNodes[1] > 0
        node = gainHeaps[1][ptr][1]
        done = binsArea[:, 2] .+ area[node] <= target_area[:, 2]

        while done == false
            if ptr < numNodes[1]
                ptr += 1
                done = binsArea[:, 2] .+ area[node] <= target_area[:, 2]
            else
                return (-1, -1)
            end
        end

        best_node = gainHeaps[1][ptr][1]
        best_side = 0

    else
        gain_0 = gainHeaps[1][1][2]
        gain_1 = gainHeaps[2][1][2]

        if gain_0 == gain_1
            node_0 = gainHeaps[1][ptr][1]
            node_1 = gainHeaps[2][ptr][1]
            balance_0 = binsArea[:, 2] .+ area[node_0] <= target_area[:, 2]
            balance_1 = binsArea[:, 1] .+ area[node_1] <= target_area[:, 1]

            if balance_0 == true && balance_1 == true
                best_node = node_0 < node_1 ? node_0 : node_1
                best_side = noed_0 < noed_1 ? 0 : 1
            elseif balance_0 == true && balance_1 == false
                best_node = node_0
                best_side = 0
            elseif balance_0 == false && balance_1 == true
                best_node = node_1
                best_side = 1
            else
                return (-1, -1) #need to rework this
            end

        else
            side = gainHeaps[1][1][2] > gainHeaps[2][1][2] ? 0 : 1
            other_side = side == 0 ? 1 : 0
            node = gainHeaps[side+1][ptr][1]
            balance = binsArea[:, other_side+1] .+ area[node] <= target_area[:, other_side+1]

            if balance == false
                if ptr == numNodes[side+1]
                    heap_block = gainHeaps[other_side+1]
                    node = heap_block[ptr][1]
                    balance = binsArea[:, other_side+1] .+ area[node] <= target_area[:, other_side+1]

                    while balance == false
                        if ptr < numNodes[other_side+1]
                            ptr += 1
                            node = heap_block[ptr][1]
                            balance = binsArea[:, other_side+1] .+ area[node] <= target_area[:, other_side+1]
                        else
                            return (-1, -1)
                        end
                    end

                    best_side = other_side
                    best_node = node
                else
                    best_side = other_side
                    heap_block = gainHeaps[other_side+1]
                    best_node = heap_block[ptr][1]
                end
            else
                best_side = side
                heap_block = gainHeaps[side+1]
                best_node = heap_block[ptr][1]
            end
        end
    end

    gain_block = gainHeaps[best_side+1][numNodes[best_side+1]]
    gainHeaps[best_side+1][ptr] = gain_block
    gainHeapLocators[gain_block[1]] = [best_side+1, ptr]
    heapifyDown!(ptr, numNodes[best_side+1]-1, gainHeapLocators, gainHeaps, gain_block[2], best_side+1)
    numNodes[best_side+1] -= 1

    return (best_node, best_side)
end

function pop_heap!(side::Int, gainHeaps::Vector, gainHeapLocators::Vector, numNodes::Vector{Int})
    success_flag = false

    if numNodes[side] > 1
        gain_block = gainHeaps[side][numNodes[side]]
        gainHeaps[side][1] = gain_block
        gainHeapLocators[gain_block[1]] = [side, 1]
        heapifyDown!(1, numNodes[side]-1, gainHeapLocators, gainHeaps, gain_block[2], side)
        numNodes[side] -= 1
        success_flag = true
    elseif numNodes[side] == 1
        numNodes[side] = 0
        success_flag = true
    else
        success_flag = false
    end

    return success_flag
end

function improveAreaCost(bins::Vector{Int}, binsArea::Array{Int}, fixed::Vector{Int}, area::Array{Int}, H::Hypergraph,
                        B::Incidence, cut_size::Int, target_cut::Int, max_cut::Int, target_area::Array{Int}, capacity::Array{Int}, cut_cost::Float64)
    init_cost = deepcopy(cut_cost)
    n = H.n
    side = 0
    base_cut = 1e06
    base_area = 1e08
    done = false
    order = 1
    pass = 1
    keep_going = 0
    moved = zeros(Int, n)
    marked = zeros(Int, n)
    (netDegs, gains) = compInit2WayGains(bins, fixed, H, B)
    perm = sortperm(gains, rev=true)
    gainHeaps = [Vector(), Vector()]
    gainHeapLocators = [zeros(Int, 2) for i in 1:n]
    numNodes = [0, 0]
    updateGainHeap!(perm, bins, gainHeaps, gainHeapLocators, numNodes, gains)

    if binsArea[:, 1] <= target_area[:, 1] && binsArea[:, 2] <= target_area[:, 2]
        side =  binsArea[:, 1] > binsArea[:, 2] ? 0 : 1
    elseif binsArea[:, 1] > target_area[:, 1]
        side = 0
    elseif binsArea[:, 2] > target_area[:, 2]
        side = 1
    end

    while done == false
        block = gainHeaps[side+1][1]
        nId = block[1]
        pop_flag = pop_heap!(side+1, gainHeaps, gainHeapLocators, numNodes)

        if pop_flag == false
            done = true
            continue
        end

        moved[nId] = pass

        if fixed[nId] > 0
            continue
        end

        best_node = nId
        fromP = side
        toP = side == 0 ? 1 : 0
        cut_size -= gains[best_node]
        excess_cut = cut_size - target_cut
        exp_cut = excess_cut/(max_cut-target_cut)
        cut_cost = base_cut^exp_cut
        excess_area_0 = binsArea[1, 1] - target_area[1, 1]
        excess_area_1 = binsArea[1, 2] - target_area[1, 2]
        exp_area_0 = excess_area_0/(capacity[1, 1]-target_area[1, 1])
        exp_area_1 = excess_area_1/(capacity[1, 2]-target_area[1, 2])
        area_cost = base_area^exp_area_0 + base_area^exp_area_1
        total_cost = area_cost + cut_cost

        if binsArea[1, fromP] > target_area[1, fromP] || total_cost < init_cost
            keep_going = 0
            best_node = nId
            bins[best_node] = toP
            binsArea[:, fromP+1] -= area[best_node, :]
            binsArea[:, toP+1] += area[best_node, :]
            nbrs = updAfterMove!(gains, fixed, moved, marked, netDegs, nId, bins, H, B, order, pass)
            updateNbrs!(nbrs, gains, gainHeapLocators, gainHeaps, numNodes, bins)
            order += 1
            init_cost = total_cost
        else
            keep_going += 1

            if keep_going == 2
                cut_size += gains[best_node]
                done = true
                continue
            else
                best_node = nId
                bins[best_node] = toP
                binsArea[:, fromP+1] -= area[best_node, :]
                binsArea[:, toP+1] += area[best_node, :]
                nbrs = updAfterMove!(gains, fixed, moved, marked, netDegs, best_node, bins, H, B, order, pass)
                updateNbrs!(nbrs, gains, gainHeapLocators, gainHeaps, numNodes, bins)
                order += 1
                init_cost = total_cost
            end
        end
    end

    return (cut_size, bins, binsArea)
end

function improveCutCost(bins::Vector{Int}, binsArea::Array{Int}, fixed::Vector{Int}, area::Array{Int}, H::Hypergraph,
                        B::Incidence, cut_size::Int, target_cut::Int, max_cut::Int, target_area::Array{Int}, capacity::Array{Int}, cut_cost::Float64, early_exit::Bool)
    init_cost = deepcopy(cut_size)
    prev_cost = init_cost
    min_cost = init_cost
    max_passes = 10
    n = H.n
    side = 0
    base_cut = 1e06
    base_area = 1e08
    done = false
    order = 1
    order_loc = 0
    pass = 1
    keep_going = 0
    keep_going_threshold = 10
    saturation = 0
    saturation_threshold = 3
    chain_count = 0
    moved = zeros(Int, n)
    marked = zeros(Int, n)
    (netDegs, gains) = compInit2WayGains(bins, fixed, H, B)
    perm = sortperm(gains, rev=true)
    gainHeaps = [Vector(), Vector()]
    gainHeapLocators = [zeros(Int, 2) for i in 1:n]
    numNodes = [0, 0]
    updateGainHeap!(perm, bins, gainHeaps, gainHeapLocators, numNodes, gains)

    if n <= 20_000
        chain_count = n
    else
        chain_count = 10_000
    end

    if early_exit == true
        chain_count = 200
        keep_going_threshold = 5
        saturation_threshold = 2
    end

    while done == false
        swaps = Vector{Int}()
        perm = sortperm(gains, rev=true)
        gainHeaps = [Vector(), Vector()]
        gainHeapLocators = [zeros(Int, 2) for i in 1:n]
        numNodes = [0, 0]
        updateGainHeap!(perm, bins, gainHeaps, gainHeapLocators, numNodes, gains)

        for i in 1:chain_count
            if keep_going == keep_going_threshold
                break
            end

            (best_node, best_side) = bestNodeToMove(gainHeaps, gainHeapLocators, binsArea, area, numNodes, target_area)

            if best_node == -1
                break
            end

            if fixed[best_node] > 0
                continue
            end

            push!(swaps, best_node)
            from = best_side
            to = from == 0 ? 1 : 0
            binsArea[:, from+1] -= area[best_node, :]
            binsArea[:, to+1] += area[best_node, :]
            bins[best_node] = to
            cut_size -= gains[best_node]
            nbrs = updAfterMove!(gains, fixed, moved, marked, netDegs, best_node, bins, H, B, order, pass)
            updateNbrs!(nbrs, gains, gainHeapLocators, gainHeaps, numNodes, bins)
            order += 1

            #=
            excess_cut = cut_size - target_cut
            exp_cut = excess_cut/(max_cut-target_cut)
            cut_cost = base_cut^exp_cut
            excess_area_0 = binsArea[1, 1] - target_area[1, 1]
            excess_area_1 = binsArea[1, 2] - target_area[1, 2]
            exp_area_0 = excess_area_0/(capacity[1, 1]-target_area[1, 1])
            exp_area_1 = excess_area_1/(capacity[1, 2]-target_area[1, 2])
            area_cost = base_area^exp_area_0 + base_area^exp_area_1
            total_cost = area_cost + cut_cost
            =#
            total_cost = cut_size

            if total_cost < min_cost
                min_cost = total_cost
                order_loc = order
                keep_going = 0
            else
                if total_cost > prev_cost
                    keep_going += 1
                else
                    keep_going = 0
                end
            end

            prev_cost = total_cost
        end

        if min_cost >= init_cost
            saturation += 1
        else
            saturation = 0
        end

        if length(swaps) > 0
            init_cost = min_cost
            rollBack!(order_loc, fixed, binsArea, gains, swaps, bins, area, moved, marked, netDegs, H, B, pass, order, cut_size)
        end

        pass += 1

        if pass <= max_passes
            done = true
        end

        if saturation <= saturation_threshold
            done = true
        end
    end

    return (cut_size, bins, binsArea)
end


text:
using Plots
theme(:default)
pyplot()

using DelimitedFiles

push!(LOAD_PATH,"src")
using DalitzPlotAnalysis
using amplitudes_compass
using PWAHelper
using SDMHelper
using FittingPWALikelihood

################################
mass_bin_name = "1540_1560" # "2320_2340"
wavelist = get_wavelist(joinpath(pwd(),"src/wavelist_formated.txt");
    path_to_thresholds=joinpath(pwd(),"src/thresholds_formated.txt"),
    M3pi=Meta.parse(mass_bin_name[1:4])/1000)
Nwaves = size(wavelist,1)
wavenames = get_wavenames(wavelist)
wavebasis = get_wavebasis(wavelist)
################################

function readnlines(f,n)
    lines = ""
    for i=1:n
        lines=lines*"\n"*readline(f)
    end
    return lines
end

readmatrix(f, nlines) = readdlm(IOBuffer(readnlines(f,nlines)))

my_old = open("data/functions_1540_1560_t1_rd.txt") do io
    readmatrix(io, 10)
end

srijan = open("data/srijan_calcs/functions_1540_1560_t1_rd.txt") do io
    readmatrix(io, 10)
end

sort(collect(zip(1:170,(my_old .- srijan)[1,:])), by=x->x[2])

my_dta = open("data/variables_1540_1560_t1_rd.txt") do io
    readmatrix(io, 10)
end

current = [wf(my_dta[1,2:end]..., my_dta[1,1]) for wf in wavebasis]

max(abs.((real.(current) - my_old[1,1:85]))...)

old_comt = [COMPASS_wave(i,my_dta[1,1:end]...) for i in 1:88]

max(abs.((real.(old_comt)[1:4] - my_old[1,1:4]))...)
max(abs.((real.(old_comt)[1:4] - srijan[1,1:4]))...)
max(abs.((real.(  my_old)[1,1:4] - srijan[1,1:4]))...)

my_dta[1:4]

srijan[1,1:4]

my_old[1,1:4]

real.(old_comt)[1:4]

COMPASS_wave(4,my_dta[1,1:end]...)#, my_dta[1,1]

wavenames[1:4]

my_dta[1,1]
my_dta[1,2]
λ(my_dta[1,1], mπ2, my_dta[1,2])

ext: jl


text:
using Weave,Trapz
weave(joinpath(dirname(pathof(Trapz)), "../benchmarks", "benchmark.jmd"), out_path=joinpath(dirname(pathof(Trapz)), "../benchmarks"), doctype="github")


text:
module SafeREPL

Base.Experimental.@optlevel 0

export swapliterals!

using SwapLiterals: SwapLiterals, makedict, floats_use_rationalize!,
                    literalswapper, default_literalswapper, defaultswaps

using REPL

function __init__()
    # condensed equivalent version of `swapliterals!()`, for faster loading
    push!(get_transforms(), default_literalswapper)
end

const LAST_SWAPPER = Ref{Function}(default_literalswapper)

function get_transforms()
    if isdefined(Base, :active_repl_backend) &&
        isdefined(Base.active_repl_backend, :ast_transforms)
        Base.active_repl_backend.ast_transforms::Vector{Any}
    else
        REPL.repl_ast_transforms::Vector{Any}
    end
end

"""
    SafeREPL.swapliterals!(Float64, Int, Int128, BigInt=nothing)

Specify transformations for literals:
argument `Float64` corresponds to literals of type `Float64`, etcetera.

A transformation can be
* a `Symbol`, to refer to a function, e.g. `:big`;
* `nothing` to not transform literals of this type;
* a `String` specifying the name of a string macro, e.g. `"@big_str"`,
  which will be applied to the input. Available only for
  `Int128` and `BigInt`, and experimentally for `Float64`.
"""
function swapliterals!(Float64,
                       Int,
                       Int128,
                       BigInt=nothing)
    @nospecialize
    swapliterals!(; Float64, Int, Int128, BigInt)
end

function swapliterals!(@nospecialize(swaps::AbstractDict))
    swapliterals!(false) # remove previous settings

    LAST_SWAPPER[] = swaps === defaultswaps ? default_literalswapper :
                                              literalswapper(swaps)

    push!(get_transforms(), LAST_SWAPPER[])
    nothing
end

# called only when !isempty(swaps)
swapliterals!(swaps::Pair...) = swapliterals!(makedict(swaps))

# non-public API when !isempty(kwswaps)
function swapliterals!(; kwswaps...)
    swapliterals!(
        isempty(kwswaps) ?
            defaultswaps :
            makedict(Any[getfield(Base, first(sw)) => last(sw) for sw in kwswaps]))
end

function swapliterals!(activate::Bool)
    transforms = get_transforms()
    # first always de-activate
    filter!(f -> parentmodule(f) != SwapLiterals, transforms)
    if activate
        push!(transforms, LAST_SWAPPER[])
    end
    nothing
end

isactive() = any(==(SwapLiterals) ∘ parentmodule, get_transforms())

end # module


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("Parallel_Computing")

# # A brief introduction to Tasks
#
# You're working on a computer that's doing _lots_ of things. It's managing
# inputs, outputs, delegating control of the CPU between Julia and _all_ of
# the other applications you have running. This wasn't always the case — does
# anyone remember the days before you could just switch between applications?
#
# It's not really doing all these things at once, but for the most part it
# gives the _appearance_ of parallelism. We think about our computers as doing
# _lots_ of things simultaneously — but it's not really simultaneous. It's just
# switching between tasks so fast that it feels simultaneous.
#
# This kind of task switching is perfect for situations like an operating system
# where you're just waiting for user input most of the time. The OS multitasking
# you're familiar with is called "preemptive" multitasking — the operating system
# sits at the top and can arbitrarily control who gets to run when. Julia's task
# system uses cooperative multitasking (also known as coroutines or green threads).

#-

# Tasks work best when they're waiting for some _external_ condition to complete
# their work. Let's say we had a directory "results" and wanted to process any
# new files that appeared there:

using FileWatching
isdir("results") || mkdir("results")
watch_folder("results", #= time out in seconds =# 5)

# Julia happily will sit there and wait for something to happen... but it's
# blocking anything else from happening while it's doing so! This is the perfect
# case for a Task. We can say we want a given expression to run asynchronously
# in a Task with the `@async` macro

t = @async watch_folder("results") # no timeout means it will wait forever!

#-

run(`touch results/0.txt`)

#-

file, info = fetch(t)
file # |> process

# We can even bundle this up into a repeating task:

isdone = false
function process_folder(dir)
    !isdir("processed-results") && mkdir("processed-results")
    while !isdone
        file, info = watch_folder(dir)
        path = joinpath(dir, file)
        if isfile(path)
            print("processing $path...")
            run(`cp $path processed-results/$file`) # Or actually do real work...
        end
    end
end

t = @async process_folder("results")

#-

run(`touch results/1.txt`)
sleep(.1)
readdir("processed-results")

#-

run(`touch results/2.txt`)
sleep(.1)
readdir("processed-results")

#-

isdone = true
run(`touch results/3.txt`)
sleep(.1)
readdir("processed-results")

#-

run(`touch results/4.txt`)
sleep(.1)
readdir("processed-results")

#-

rm("results", recursive=true)
rm("processed-results", recursive=true)

# ## Quiz:
#
# How long will this take?

@time for i in 1:10
    sleep(1)
end

# What about this?

@time for i in 1:10
    @async sleep(1)
end

# And finally, this?

@time @sync for i in 1:10
    @async sleep(1)
end

# Now what if I had something that actually did work?

function work(N)
    series = 1.0
    for i in 1:N
        series += (isodd(i) ? -1 : 1) / (i*2+1)
    end
    return 4*series
end
work(1)
@time work(100_000_000)

#-

@time @sync for i in 1:10
    @async work(100_000_000)
end

# # So what's happening here?
#
# `sleep` is nicely cooperating with our tasks

methods(sleep)

# # Fetching values from tasks

#-

# You can even fetch values from tasks

t = @async (sleep(5); rand())

#-

wait(t)

#-

fetch(t)

# # Key takeaways
#
# There is a lot more to tasks, but they form the foundation for reasoning about
# actually _doing_ computation in parallel (and not just hoping that things will
# cooperate for us to emulate parallelism by task switching).
#
# * `@async` creates and starts running a task
# * `@sync` waits for them to all complete
# * We can reason about something that runs asynchronously and may return a value
#   at some point in the future with `fetch`. Or we can just `wait` for it.

max_stars_repo_path: Courses/Parallel_Computing/060 Tasks.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_path: Courses/Parallel_Computing/060 Tasks.jl

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_path: Courses/Parallel_Computing/060 Tasks.jl

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
abstract type Executor end
struct SequentialEx{K} <: Executor
    kwargs::K
end
struct ThreadedEx{K} <: Executor
    kwargs::K
end
struct DistributedEx{K} <: Executor
    kwargs::K
end
struct PreferParallel{K} <: Executor
    kwargs::K
end

"""
    SequentialEx(; kwargs...)
    ThreadedEx(; kwargs...)
    DistributedEx(; kwargs...)

Sequential, threaded, and distributed executor.  An executor specifies
execution strategy and its parameters.

See [`foldxl`](@ref), [`foldxt`](@ref) and [`foldxd`](@ref) for usable
keyword arguments.
"""
(SequentialEx, ThreadedEx, DistributedEx)

"""
    PreferParallel(; simd, basesize)

A "placeholder" executor that indicates preference to parallel execution.

This lets the input data collection decide preferred execution strategy
(e.g., `CUDAEx` for `CuArray` when FoldsCUDA.jl is available), **assuming
that the reducing function is associative**. The default executor is
[`ThreadedEx`](@ref). As an optional feature, some input data collections
support (e.g., `AbstractChannel`) automatically demoting the execution
strategy to [`SequentialEx`](@ref).  An error is thrown if the automatic
detection fails,
"""
PreferParallel

maybe_set_simd(exc::Executor, simd) =
    isdefined(exc.kwargs, :simd) ? exc : @set exc.kwargs = merge(exc.kwargs, (simd = simd,))

(::Type{Ex})(; kwargs...) where {Ex<:Executor} = Ex(values(kwargs))

transduce(xf, rf::RF, init, coll, exc::SequentialEx) where {RF} =
    transduce(xf, rf, init, coll; exc.kwargs...)
transduce(xf, rf, init, coll, exc::ThreadedEx) =
    transduce_assoc(xf, rf, init, coll; exc.kwargs...)
transduce(xf, rf, init, coll, exc::DistributedEx) =
    dtransduce(xf, rf, init, coll; exc.kwargs...)

fold(rf::RF, coll, exec::Executor; init = DefaultInit) where {RF} =
    unreduced(transduce(IdentityTransducer(), Completing(rf), init, coll, exec))
fold(rf::RF, coll; init = DefaultInit, kwargs...) where {RF} =
    fold(rf, coll, PreferParallel(; kwargs...); init = init)

function transduce(xf, rf, init, coll, exc::PreferParallel)
    xf0, coll0 = extract_transducer(coll)
    return transduce(xf ∘ xf0, rf, init, coll0, executor_for(coll0, exc))
end

fold(rf; kw...) = itr -> fold(rf, itr; kw...)

executor_from(::Type{Ex}, exc::PreferParallel) where {Ex<:Executor} = Ex(exc.kwargs)
executor_from(::Type{PreferParallel}, exc::PreferParallel) = ThreadedEx(exc.kwargs)
executor_from(::Type{SequentialEx}, exc::PreferParallel) =
    SequentialEx(; _seq_ex_kwargs(; exc.kwargs...)...)
_seq_ex_kwargs(; simd = Val(false), _...) = (simd = simd,)

executor_for(xs::Any, exc) = executor_from(executor_type(xs), exc)
# * `ThreadedEx` part sould be determined by the property of transducers (and,
#   ideally, reducing function).
# * We can extend the API of `executor_type` that returns a list of comptaible
#   executors rather than only one executor.

# `executor_type(::CustomContainer)` is the main (still hidden) API.
# TODO: should `executor_type` be renamed to `parallel_executor_type` if the
# default type is `PreferParallel`?
@inline executor_type(::Any) = PreferParallel

# Some "famous" containers known to be sequential-only:
@inline executor_type(::AbstractChannel) = SequentialEx
# @inline executor_type(::Distributed.RemoteChannel) = SequentialEx
@inline executor_type(::Iterators.Stateful) = SequentialEx

@inline executor_type(xs::Iterators.Pairs) =
    executor_promote_type(executor_type(keys(xs)), executor_type(values(xs)))

@inline executor_type(xs::Iterators.ProductIterator) =
    mapfoldlargs(executor_type, executor_promote_type, PreferParallel, xs.iterators...)

@inline executor_type(xs::_Zip) =
    mapfoldlargs(executor_type, executor_promote_type, PreferParallel, arguments(xs)...)

@inline executor_promote_type(::Type{T}, ::Type{T}) where {T} = T
@inline function executor_promote_type(::Type{T}, ::Type{S}) where {T,S}
    result =
        executor_promote_result(executor_promote_rule(T, S), executor_promote_rule(S, T))
    result === nothing && throw_promote_failure(T, S)
    return result
end

# `nothing` means promotion failure
@inline executor_promote_result(::Nothing, ::Nothing) = nothing
@inline executor_promote_result(::Type{T}, ::Type{S}) where {T,S} = nothing
@inline executor_promote_result(::Type{T}, ::Type{T}) where {T} = T
@inline executor_promote_result(::Nothing, ::Type{T}) where {T} = T
@inline executor_promote_result(::Type{T}, ::Nothing) where {T} = T

# Only the owner of the first argument should define the rule. Furthermore, a
# custom executor `CustomEx` should define
# `executor_promote_rule(::Type{SequentialEx}, ::Type{OtherEx}) = OtherEx` only
# if it knows how to handle all `xs` such that `executor_type(xs) = OtherEx`.
@inline executor_promote_rule(::Type{PreferParallel}, ::Type{T}) where {T} = T
@inline executor_promote_rule(::Type, ::Type) = nothing

@noinline throw_promote_failure(T, S) = error("failed to promote executors:\n$T\n$S")
# TODO: show information on which container "suggested" which executor?

function Base.show(io::IO, exc::Executor)
    T = typeof(exc.kwargs)
    if (
        fieldtype(typeof(exc), 1) !== typeof(exc.kwargs) ||
        any(i -> typeof(exc.kwargs[i]) != fieldtype(T, i), 1:nfields(exc.kwargs))
    )
        return invoke(show, Tuple{IO,Any}, io, exc)
    end
    Ex = getfield(parentmodule(typeof(exc)), nameof(typeof(exc)))
    print(io, Ex)
    if isempty(exc.kwargs)
        print(io, "()")
    elseif length(exc.kwargs) == 1
        print(io, '(', fieldname(T, 1), " = ", exc.kwargs[1], ')')
    else
        print(io, exc.kwargs)
    end
    return
end


text:
using RobustDMD, Random, Printf

# Generate DMD Synthetic Data
#------------------------------------------------------------------------------
# dimensions
m = 500;    # temporal dimension
n = 1000;    # spatial dimension
k = 3;      # number of modes
T = Float64;
sigman = 1e-3;
# generate data
X, t, at, Bt = simDMD(m, n, k, T; seed=123456);
Xdat = X + sigman*Random.randn!(copy(X))
# specify the loss
function lossFunc(r)
	return sum(abs2, r)
end
#
function lossGrad(r)
	conj!(r)
end
#
# Create object
#------------------------------------------------------------------------------

#
# Initial Guess
#------------------------------------------------------------------------------
# create initial guess using exact dmd
a0, B0 = dmdexactestimate(m, n, k, Xdat, t);

#
# Create Object
#------------------------------------------------------------------------------
params = DMDParams(k, Xdat, t, lossFunc, lossGrad, inner_directl2=true);
 copyto!(params.a, a0);
# copyto!(params.B, B0);
#Random.randn!(params.ar);
#
# Apply Solver
#------------------------------------------------------------------------------
res = solveDMD_withOptimBFGS(params);

@printf "%7.4e relative residual after BFGS\n" sqrt(dmdobjective(params)/lossFunc(X))

params2 = DMDParams(k,Xdat,t,lossFunc,lossGrad);
copyto!(params2.a,at)
@printf "%7.4e relative residual with true evals\n" sqrt(dmdobjective(params2)/lossFunc(X))
copyto!(params2.a,a0)
@printf "%7.4e relative residual with initial guess\n" sqrt(dmdobjective(params2)/lossFunc(X))

@printf "%7.4e abs err evals after BFGS\n" besterrperm(params.a,at)
@printf "%7.4e abs err evals initial guess\n" besterrperm(a0,at)


text:
function GasInflow(pr,J,GasPVT;n=10)

    #Calculate PseudoPressure
    GasPVT.Ps=map((x,y,z)->2*x/(y*z),GasPVT.P,GasPVT.Mug,GasPVT.Z)

    #Integral Matrix
    Integral = DataFrame()

    #Trasnsfer Pressure to Integrate Matrix
    Integral.P= GasPVT.P[2:end]

    #Pressure Differences
    Integral.DiffP=diff(GasPVT.P)

    #PseudoPressure Average
    Integral.AvgPs=0.0
    for i=1:size(GasPVT,1)-1
        Integral.AvgPs[i]=(GasPVT.Ps[i+1]+GasPVT.Ps[i])/2
    end
    Integral.Cum=accumulate(+, Integral.DiffP.*Integral.AvgPs)

    #Interpolation Integral
    mp=LinearInterpolation(Integral.P.*1, Integral.Cum.*1)

    #IPR Estimation

    pres=zeros(n,size(pr,1))
    Qg=zeros(n,size(pr,1)*size(J,1))
    leg=Array{String,1}(undef,size(pr,1)*size(J,1))
    count=0
    for p=1:size(pr,1)

        pres[:,p]=range(14.7,stop=pr[p]-50,length=n)

        for j=1:size(J,1)
    count=count+1
    leg[count]="Res. Pressure=$(pr[p]) psi; J=$(J[j]) Mscf/d*psi*cp"
      Qg[:,count]=map(pwf->J[j]*(mp(pr[p])-mp(pwf)),pres[:,p])

        end
    end

    return pres, Qg, leg

end #End Function


text:
using ProgressMeter
using MATLAB
using Printf
using BSON: @save, @load
using PyCall
pygui(:qt)
using PyPlot
pygui(true)
include("escape_dimer_one_hit.jl")

H=.25

t_end = 1e3
width = 2.25; height = 4.5
n_iter_P = 1001;n_iter_Q =1001;
N = n_iter_P * n_iter_Q;
ArrQ = range(-width, stop = width, length = n_iter_P)
ArrP = range(0, stop = height, length = n_iter_Q)
mesh = [(Q, P) for Q in ArrQ, P in ArrP]
mesh_list = reshape(mesh, 1, :)
t_end_list = t_end * ones(N);
energy_list= H* ones(N);

h = replace(@sprintf("%.13f",H), "." => "_")
nQ = @sprintf("_%d",n_iter_Q)
nP = @sprintf("_%d",n_iter_P)
h_BSON="num_hhits_escape_data"*h*nQ*nP*".bson"
location = "/mnt/bdd38f66-9ece-451a-b915-952523c139d2/Escape/"
println(h_BSON)

h_title = @sprintf("h= %.13f",H)

num_until_exit = @showprogress map(escape_exit_num_dimer, mesh_list, t_end_list, energy_list)

one_hit_exit=[]
for i = 1:N
    Q, P = mesh_list[i]
    if num_until_exit[i] == 1
        push!(one_hit_exit,[Q,P])
    end
end
t_end=1e3
for k=1:3
old_list=copy(one_hit_exit)
energy=.25
radius=.01

@showprogress for (index,point) in enumerate(old_list)
    add_four!(point,radius, one_hit_exit, t_end, H)
end

Q,P=unzipper(one_hit_exit)
Q_one=vcat(Q,-Q); P_one=vcat(P,-P);

plot(Q_one,P_one, ".",markersize=1,c=:blue);
println(length(one_hit_exit))
SAVE_DATA=Dict("one_hit_exit"=>one_hit_exit)
@save h_BSON SAVE_DATA

end


text:
using DashBootstrapComponents

accordion = html_div([
    dbc_accordion(
        [
            dbc_accordionitem(
                "This is the content of the first section. It has a " *
                "default ID of item-0.",
                title = "Item 1: item-0",
            ),
            dbc_accordionitem(
                "This is the content of the second section. It has a " *
                "default ID of item-1.",
                title = "Item 2: item-1",
            ),
            dbc_accordionitem(
                "This is the content of the third section. It has a " *
                "default ID of item-2.",
                title = "Item 3: item-2",
            ),
        ],
        id = "accordion-always-open",
        always_open = true,
    ),
    html_div(id = "accordion-contents-open-ids", className = "mt-3"),
])

callback!(
    app,
    Output("accordion-contents-open-ids", "children"),
    Input("accordion-always-open", "active_item"),
) do item
    return "Item(s) selected: $item"
end;


text:
module CartesianGeometryTests

using Test
using Gridap
using Gridap.DiscreteModels: dict_to_model
using JSON

grid = CartesianGrid(domain=(0.0,1.0,-1.0,2.0),partition=(3,4))
test_grid(grid,20,12)

@test celldim(grid) == 2
@test pointdim(grid) == 2
@test ncells(grid) == 3*4
@test npoints(grid) == (3+1)*(4+1)

x = points(grid)

@test isa(x,CellValue{Point{2,Float64}})

@test length(x) == 20
@test x[13] == Point(0.0,1.25)
@test x[4] == Point(1.0,-1.0)
@test x[1] == Point(0.0,-1.0)
@test x[end] == Point(1.0,2.0)

t = cells(grid)

@test isa(t,CellArray{Int,1})

@test length(t) == 12

@test t[1] == [1,2,5,6]
@test t[11] == [14,15,18,19]

c = celltypes(grid)
@test isa(c,ConstantCellValue{NTuple{2,Int}})
@test c.value == (HEX_AXIS,HEX_AXIS)
@test length(c) == 12

o = cellorders(grid)
@test isa(o,ConstantCellValue{Int})
@test o.value == 1
@test length(o) == 12

trian = Triangulation(grid)

xe = CellPoints(trian)
@test isa(xe,CellPoints{2,Float64})

cb = CellBasis(trian)
@test isa(cb,CellBasis{2,Float64})

model = CartesianDiscreteModel(
  domain=(0.0,1.0,-1.0,2.0,0.0,1.0),
  partition=(3,4,2))

test_discrete_model(model,3)

grid3 = Grid(model,3)
np = (3+1)*(4+1)*(2+1)
nc = 3*4*2
test_grid(grid3,np,nc)

@test isa(grid3,CartesianGrid{3})

model = CartesianDiscreteModel(
  partition=(3,4,2),
  map = (x -> 5*x))

grid4 = Grid(model)
x = points(grid4)
D = 2
@test isa(x,CellValue{Point{3,Float64}})

gridgraph = FullGridGraph(model)

@test isa(gridgraph, FullGridGraph)

grid2 = Grid(model,2)
grid1 = Grid(model,1)
grid0 = Grid(model,0)

labels = FaceLabels(model)

@test isa(labels,FaceLabels)
@test length(labels_on_dim(labels,3)) == ncells(grid3)
@test length(labels_on_dim(labels,2)) == ncells(grid2)
@test length(labels_on_dim(labels,1)) == ncells(grid1)
@test length(labels_on_dim(labels,0)) == ncells(grid0)

@test ntags(labels) == 28
@test name_from_tag(labels,28) == "boundary"

cgrid = CartesianGrid(domain=(0.0,1.0,-1.0,2.0),partition=(3,4))

grid = FlexibleUnstructuredGrid(cgrid)
test_grid(grid,20,12)

grid = UnstructuredGrid(cgrid)
test_grid(grid,20,12)

# Serialization

model = CartesianDiscreteModel(partition=(2,2,2))
s = json(model)

dict = JSON.parse(s)
model = dict_to_model(dict)
test_discrete_model(model,3)

d = mktempdir()
filename = joinpath(d,"model.json")

open(filename,"w") do f
  JSON.print(f,model)
end

model = DiscreteModelFromFile(filename)
test_discrete_model(model,3)

end # module


text:
#!/usr/bin/env julia

function main2()
  lines_read = 0
  res = zeros(Int, 5)
  for l in eachline("./input")
    function calc_index(right_offset)
      1 + ((lines_read * right_offset) % length(l))
    end

    trees = [
             l[calc_index(1)] == '#',
             l[calc_index(3)] == '#',
             l[calc_index(5)] == '#',
             l[calc_index(7)] == '#',
             lines_read % 2 == 0 ? l[calc_index(1)] == '#' : 0
            ]
    res += trees
    lines_read += 1
  end
  println(reduce(*, res))
end

main2()

# part one
function main()
  lines_read = 0
  res = 0
  for l in eachline("./input")
    if (l[1 + ((lines_read * 3) % length(l))] == '#')
      res += 1
    end
    lines_read += 1
  end

  println(res)
end


text:
#=
Simple example of use of montecarlo.jl
From the system prompt, execute
mpirun -np X julia pi.jl
setting X to the number of ranks you'd
like to use, subject to  X-1 being an even divisor
of 1e6, e.g., set X=5.
=#
using Pkg
Pkg.activate(".")
using MPI, Statistics, LinearAlgebra
include("montecarlo.jl")
# this is the function that gets Monte Carlo'd
function pi_wrapper()
    pihat = [4.0*float((norm(rand(2,1)) .< 1.))]
end
# this function reports intermediate results during MC runs
function pi_monitor(sofar, results)
    # examine results every 2.5*10^5 draws
    if mod(sofar,2.5e5)==0.
        m = mean(results[1:sofar,:],dims=1)
        println("reps so far: ", sofar)
        println("pihat: ", m)
        println()
        #if sofar == size(results,1)
        #    writedlm("mcresults.out", results)
        #end
    end
end
# do the monte carlo: 10^6 reps of single draws
function main()
    if !MPI.Initialized()
        MPI.Init()
    end
    comm = MPI.COMM_WORLD
    reps = Int(1e7)  # desired number of MC reps
    nreturns = 1
    pooled = Int(1e5)
    montecarlo(pi_wrapper, pi_monitor, comm, reps, nreturns, pooled)
end
main()

max_issues_repo_name: mcreel/JuliaMPIMonteCarlo

max_forks_repo_name: mcreel/JuliaMPIMonteCarlo


text:
import Pkg
Pkg.activate("examples")
using Revise, GpABC, OrdinaryDiffEq, Distances, Distributions, Plots
import RecipesBase
RecipesBase.debug()
import Logging

# seed!(4)

#
# ABC settings
#

n_particles = 1000
threshold = 2.0
threshold_schedule = [3.0, 2.0, 1.0, 0.5, 0.2]
progress_every = 1000

#
# Emulation settings
#
n_design_points = 100

#
# True parameters
#
true_params =  [2.0, 2.0, 15.0, 1.0, 1.0, 1.0, 100.0, 1.0, 1.0, 1.0]
priors = [Uniform(0.2, 5.), Uniform(0.2, 5.), Uniform(10., 20.),
            Uniform(0.2, 2.), Uniform(0.2, 2.), Uniform(0.2, 2.),
            Uniform(75., 125.),
            Uniform(0.2, 2.), Uniform(0.2, 2.), Uniform(0.2, 2.)]
param_indices = [1,2,3]
n_var_params = length(param_indices)

#
# ODE solver settings
#
Tspan = (0.0, 10.0)
x0 = [0.0, 0.0, 0.0]
solver = RK4()
saveat = 0.1

function ODE_3GeneReg(dx, x, pars, t)
  dx[1] = pars[1]/(1+pars[7]*x[3]) - pars[4]*x[1]
  dx[2] = pars[2]*pars[8]*x[1]/(1+pars[8]*x[1]) - pars[5]*x[2]
  dx[3] = pars[3]*pars[9]*x[1]*pars[10]*x[2]/((1+pars[9]*x[1])*(1+pars[10]*x[2])) - pars[6]*x[3]
end

#
# Returns the solution to the toy model as solved by OrdinaryDiffEq
#
function GeneReg(params::AbstractArray{Float64,1},
    Tspan::Tuple{Float64,Float64}, x0::AbstractArray{Float64,1},
    solver::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm, saveat::Float64)

  prob = ODEProblem(ODE_3GeneReg, x0 ,Tspan, params)
  Obs = solve(prob, solver, saveat=saveat)

  return Array{Float64, 2}(Obs)
end

reference_data = GeneReg(true_params, Tspan, x0, solver, saveat)
# reference_data += randn(size(reference_data)) * 0.1

function simulator_function(var_params::AbstractArray{T, 1}) where {T<:Real}
    params = copy(true_params)
    params[param_indices] .= var_params
    GeneReg(params, Tspan, x0, solver, saveat)
end

logger=Logging.SimpleLogger(stdout,Logging.Debug)
# old_logger = Logging.global_logger()
# Logging.global_logger(logger)
Logging.with_logger(logger) do
    sim_result = SimulatedABCRejection(reference_data, simulator_function, priors[param_indices], threshold, n_particles)
    plot(sim_result)
end


text:
using GraphLaplacians
using LinearAlgebra
using Graphs
using SimpleWeightedGraphs
using Test

tests = [
    "adjmat",
    "sparse",
    "graphs",
]

@testset "GraphLaplacians.jl" begin
    for t in tests
        include("$(t).jl")
    end
end


text:
function allowed_changed_files(::NewPackage, pkg::String)
    first_letter = uppercase(pkg[1:1])
    result = String["Registry.toml",
                    "$(first_letter)/$(pkg)/Compat.toml",
                    "$(first_letter)/$(pkg)/Deps.toml",
                    "$(first_letter)/$(pkg)/Package.toml",
                    "$(first_letter)/$(pkg)/Versions.toml"]
    return result
end

function allowed_changed_files(::NewVersion, pkg::String)
    first_letter = uppercase(pkg[1:1])
    result = String["$(first_letter)/$(pkg)/Compat.toml",
                    "$(first_letter)/$(pkg)/Deps.toml",
                    "$(first_letter)/$(pkg)/Versions.toml"]
    return result
end

const guideline_pr_only_changes_allowed_files =
    Guideline("Only modifies the files that it's allowed to modify",
              data -> pr_only_changes_allowed_files(data.api,
                                                    data.registration_type,
                                                    data.registry,
                                                    data.pr,
                                                    data.pkg;
                                                    auth = data.auth))

function pr_only_changes_allowed_files(api::GitHub.GitHubAPI,
                                       t::Union{NewPackage, NewVersion},
                                       registry::GitHub.Repo,
                                       pr::GitHub.PullRequest,
                                       pkg::String;
                                       auth::GitHub.Authorization)
    _allowed_changed_files = allowed_changed_files(t, pkg)
    _num_allowed_changed_files = length(_allowed_changed_files)
    this_pr_num_changed_files = num_changed_files(pr)
    if this_pr_num_changed_files > _num_allowed_changed_files
        g0 = false
        m0 = "This PR is allowed to modify at most $(_num_allowed_changed_files) files, but it actually modified $(this_pr_num_changed_files) files."
        return g0, m0
    else
        this_pr_changed_files = get_changed_filenames(api, registry, pr; auth = auth)
        if length(this_pr_changed_files) != this_pr_num_changed_files
            g0 = false
            m0 = "Something weird happened when I tried to get the list of changed files"
            return g0, m0
        else
            if issubset(this_pr_changed_files, _allowed_changed_files)
                g0 = true
                m0 = ""
                return g0, m0
            else
                g0 = false
                m0 = string("This pull request modified at least one file ",
                            "that it is not allowed to modify. It is only ",
                            "allowed to modify the following files ",
                            "(or a subset thereof): ",
                            "$(join(_allowed_changed_files, ", "))")
                return g0, m0
            end
        end
    end
end


text:
lines = readlines("input.txt")
numbers = parse.(Int, lines)

println(sum(numbers[2:end] .> numbers[1:end-1]))


text:
export reset_outputs

@defstruct MemoryOutputLayer Layer (
  name :: AbstractString = "memory-output",
  (bottoms :: Vector{Symbol} = [], length(bottoms) > 0)
)
@characterize_layer(MemoryOutputLayer,
  is_sink => true
)

type MemoryOutputLayerState <: LayerState
  layer   :: MemoryOutputLayer
  outputs :: Vector{Vector{Array}}
end

function reset_outputs(state::MemoryOutputLayerState)
  for i = 1:length(state.outputs)
    state.outputs[i] = Array[]
  end
end

function setup(backend::Backend, layer::MemoryOutputLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  outputs = Array(Vector{Array}, length(inputs))
  for i = 1:length(inputs)
    outputs[i] = Array[]
  end

  return MemoryOutputLayerState(layer, outputs)
end

function forward(backend::Backend, state::MemoryOutputLayerState, inputs::Vector{Blob})
  for i = 1:length(inputs)
    push!(state.outputs[i], to_array(inputs[i]))
  end
end

function backward(backend::Backend, state::MemoryOutputLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
end

function shutdown(backend::Backend, state::MemoryOutputLayerState)
end


text:
using Random

function make_jacobian(ws)
    jacobian_orig = zeros(ws.endogenous_nbr,
                          ws.backward_nbr
                          + ws.current_nbr
                          + ws.forward_nbr
                          + ws.exogenous_nbr)
    for i = 1:1000
        Random.seed!(i)
        jacobian = randn(ws.endogenous_nbr,
                          ws.backward_nbr
                          + ws.current_nbr
                          + ws.forward_nbr
                         + ws.exogenous_nbr)
        copy!(jacobian_orig, jacobian)
        LinearRationalExpectations.remove_static!(jacobian, ws)
        LinearRationalExpectations.get_de!(ws, jacobian)
        F = schur(ws.e, ws.d)
        if count(abs.(F.α ./ F.β) .> 1.0) == ws.forward_nbr
            break
        end
    end
    return jacobian_orig
end

algo = "GS"

endogenous_nbr = 10
exogenous_nbr = 3
exogenous_deterministic_nbr = 0
forward_indices = [2, 3, 5, 7, 9]
current_indices = collect(2:10)
backward_indices = [1, 4, 6, 7, 9]
both_indices = [7, 9]
static_indices = [8, 10]

ws = LinearRationalExpectationsWs(algo,
                                  endogenous_nbr,
                                  exogenous_nbr,
                                  exogenous_deterministic_nbr,
                                  forward_indices,
                                  current_indices,
                                  backward_indices,
                                  both_indices,
                                  static_indices)
@test ws.static_indices_j == [12, 14]

results = LinearRationalExpectationsResults(endogenous_nbr,
                                            exogenous_nbr,
                                            ws.backward_nbr)
jacobian = make_jacobian(ws)
jacobian_orig = copy(jacobian)
LinearRationalExpectations.remove_static!(jacobian, ws)
FQ = LinearAlgebra.qr(jacobian_orig[:, [12, 14]])
@test jacobian ≈ transpose(FQ.Q)*jacobian_orig

LinearRationalExpectations.get_de!(ws, jacobian)
@test ws.d[1:8, 1] == zeros(8)
@test ws.d[1:8, 2:10]  == jacobian[3:end, [8, 10, 11, 13, 15, 16, 17, 18, 19]]
@test ws.e[1:8, 1:8]  == -jacobian[3:end, [1, 2, 3, 4, 5, 6, 7, 9]]

targetD = zeros(2, 10)
targetD[1, 4] = 1.0
targetD[2, 5] = 1.0
targetE = zeros(2, 10)
targetE[1, 9] = 1.0
targetE[2, 10] = 1.0
@test ws.d[9:10,:] == targetD
@test ws.e[9:10,:] == targetE

d_orig = copy(ws.d)
e_orig = copy(ws.e)
options = LinearRationalExpectationsOptions()
LinearRationalExpectations.PolynomialMatrixEquations.gs_solver!(ws.solver_ws, ws.d, ws.e, ws.backward_nbr, options.generalized_schur.criterium)
@test d_orig * vcat(I(ws.backward_nbr), ws.solver_ws.g2[:, 1:ws.backward_nbr])*ws.solver_ws.g1 ≈ e_orig * vcat(I(ws.backward_nbr), ws.solver_ws.g2[:, 1:ws.backward_nbr])

results.gs1 .= ws.solver_ws.g1
for i = 1:ws.backward_nbr
    for j = 1:ws.backward_nbr
        x = ws.solver_ws.g1[j,i]
        results.g1[ws.backward_indices[j],i] = x
    end
    for j = 1:(ws.forward_nbr - ws.both_nbr)
        results.g1[ws.purely_forward_indices[j], i] =
            ws.solver_ws.g2[ws.icolsE[ws.backward_nbr + j] - ws.backward_nbr, i]
    end
end

LinearRationalExpectations.add_static!(results, jacobian, ws)
@test ws.b11[:, 2:8] ≈ jacobian[1:ws.static_nbr, [6, 7, 8, 9, 10, 11, 13]]
@test results.g1[ws.forward_indices, 1:ws.backward_nbr] ≈ ws.solver_ws.g2
@test results.g1[ws.backward_indices, 1:ws.backward_nbr] ≈ ws.solver_ws.g1
b10 = jacobian[1:ws.static_nbr, [12, 14]]
@test ws.b10 ≈ b10
target = -b10\(jacobian[1:ws.static_nbr, 15:19]
               *results.g1_1[ws.forward_indices,:]
               *results.gs1
               + jacobian[1:ws.static_nbr, [6, 7, 8, 9, 10, 11, 13]]
               *results.g1_1[[2, 3, 4, 5, 6, 7, 9], :]
               + jacobian[1:ws.static_nbr, 1:ws.backward_nbr])
@test results.g1[ws.static_indices, 1:ws.backward_nbr] ≈ target

A = randn(10, 5)
B = randn(10, 9)
G = randn(endogenous_nbr, length(backward_indices))
LinearRationalExpectations.make_lu_AGplusB!(ws.AGplusB, A, G, B, ws)
AGplusB = zeros(endogenous_nbr, endogenous_nbr)
AGplusB[:, 2:10] = B
AGplusB[:, [1, 4, 6, 7, 9]] += A*G[[2, 3, 5, 7, 9], :]
@test AGplusB ≈ ws.AGplusB
LU = reshape(ws.AGplusB_linsolve_ws.lu, endogenous_nbr, endogenous_nbr)
U = triu(LU)
L = tril(LU, -1) + diagm(ones(10))

function swaprow!(x, i, j)
    for k in axes(x, 2)
        x[i, k], x[j, k] = x[j, k], x[i, k]
    end
end
target = copy(AGplusB)
for i in axes(target, 2)
    swaprow!(target, i, ws.AGplusB_linsolve_ws.ipiv[i])
end
@test target ≈ L*U

LinearRationalExpectations.solve_for_derivatives_with_respect_to_shocks!(results, jacobian, ws)
g1_2 = -AGplusB\jacobian[:, 20:22]
@test results.g1_2 ≈ g1_2

algo = "GS"
ws = LinearRationalExpectationsWs(algo,
                                  endogenous_nbr,
                                  exogenous_nbr,
                                  exogenous_deterministic_nbr,
                                  forward_indices,
                                  current_indices,
                                  backward_indices,
                                  both_indices,
                                  static_indices)
jacobian = make_jacobian(ws)
jacobian_orig = copy(jacobian)
first_order_solver!(results, algo, jacobian, options, ws)
@test d_orig * vcat(I(ws.backward_nbr), ws.solver_ws.g2[:, 1:ws.backward_nbr])*ws.solver_ws.g1 ≈ e_orig * vcat(I(ws.backward_nbr), ws.solver_ws.g2[:, 1:ws.backward_nbr])
A2 = jacobian[:, 15:19]
A1 = jacobian[:, 6:14]
A0 = jacobian[:, 1:5]
@test A2[3:10,:] ≈ d_orig[1:8, 6:10]
@test A1[3:10, [3, 5, 6, 8]] ≈ d_orig[1:8, 2:5]
@test A1[3:10, [1, 2, 4]] ≈ -e_orig[1:8, 6:8]
@test A0[3:10, :] ≈ -e_orig[1:8, 1:5]
@test d_orig[1:8, 1] ≈ zeros(8, 1)
@test e_orig[1:8, 9:10] ≈ zeros(8, 2)
@test results.g1[ws.forward_indices, 1:ws.backward_nbr] ≈ ws.solver_ws.g2
@test results.g1[ws.backward_indices, 1:ws.backward_nbr] ≈ ws.solver_ws.g1
g1_1 = results.g1_1
g1_2 = results.g1_2
gg = g1_1[ws.forward_indices, :]
@test gg ≈ ws.solver_ws.g2
A2 = jacobian_orig[:, 15:19]
A1 = jacobian_orig[:, 6:14]
A0 = jacobian_orig[:, 1:5]
B = jacobian_orig[:, 20:22]

@test (A2*gg*g1_1[ws.backward_indices,:] +
       A1*g1_1[ws.current_indices,:] ≈ -A0 )
@test (A2*gg*g1_2[ws.backward_indices,:] +
       A1*g1_2[ws.current_indices,:] ≈ -B )

algo = "CR"
ws = LinearRationalExpectationsWs(algo,
                                  endogenous_nbr,
                                  exogenous_nbr,
                                  exogenous_deterministic_nbr,
                                  forward_indices,
                                  current_indices,
                                  backward_indices,
                                  both_indices,
                                  static_indices)
LinearRationalExpectations.get_abc!(ws, jacobian)
@test ws.a[:, [2, 3, 5, 7, 8]] == jacobian[3:end, 15:19]
@test ws.b[:, [2, 3, 4, 5, 7, 8]] == jacobian[3:end, [6, 7, 8, 9, 11, 13]]
@test ws.c[:, [1, 4, 6, 7, 8]] == jacobian[3:end, 1:5]


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

##########################################
# Timer2
##########################################
#
# This is a backported replacement for the new Timer code introduced in v0.4
#

type Timer2
    handle::Ptr{Void}
    cond::Condition
    isopen::Bool

    function Timer2(timeout::Real, repeat::Real=0.0)
        timeout ≥ 0 || throw(ArgumentError("timer cannot have negative timeout of $timeout seconds"))
        repeat ≥ 0 || throw(ArgumentError("timer cannot repeat $repeat times"))

        this = new(Libc.malloc(Base._sizeof_uv_timer), Condition(), true)
        err = ccall(:uv_timer_init,Cint,(Ptr{Void},Ptr{Void}),Base.eventloop(),this.handle)
        if err != 0
            #TODO: this codepath is currently not tested
            Libc.free(this.handle)
            this.handle = C_NULL
            throw(UVError("uv_make_timer",err))
        end

        Base.associate_julia_struct(this.handle, this)
        Base.preserve_handle(this)

        ccall(:uv_update_time, Void, (Ptr{Void},), Base.eventloop())
        @compat ccall(:uv_timer_start, Cint, (Ptr{Void},Ptr{Void},UInt64,UInt64),
              this.handle, Base.uv_jl_asynccb::Ptr{Void},
              UInt64(round(timeout*1000))+1, UInt64(round(repeat*1000)))
        this
    end
end

Base.wait(t::Timer2) = wait(t.cond)

Base.isopen(t::Timer2) = t.isopen

function Base.close(t::Timer2)
    if t.handle != C_NULL
        t.isopen = false
        ccall(:uv_timer_stop, Cint, (Ptr{Void},), t.handle)
        ccall(:jl_close_uv, Void, (Ptr{Void},), t.handle)
    end
end

function Base._uv_hook_close(t::Timer2)
    Base.unpreserve_handle(t)
    Base.disassociate_julia_struct(t)
    t.handle = C_NULL
    nothing
end

function Base._uv_hook_asynccb(t::Timer2)
    @compat if ccall(:uv_timer_get_repeat, UInt64, (Ptr{Void},), t.handle) == 0
        # timer is stopped now
        close(t)
    end
    notify(t.cond)
    nothing
end

# timer with repeated callback
function Base.Timer(cb::Function, timeout::Real, repeat::Real=0.0)
    t = Timer2(timeout, repeat)
    @schedule begin
        while isopen(t)
            wait(t)
            cb(t)
        end
    end
    t
end

Base.Timer(timeout::Real, repeat::Real=0.0) = Timer2(timeout, repeat)


text:
"""
    path(vertices, npts)

Creates a path in k-space. The path is formed by lines connecting the points provided in `vertices` as an array of arrays.
The return path has approx `npts` points.
The path is returned as a `KPath` object.
"""
function path(vertices, npts)

    for pt in vertices
        @assert length(pt)==length(vertices[1]) "all points must be of the same size"
    end

    v = [toSVector(u) for u in vertices]

    lengths = [norm(v[i+1] - v[i]) for i=1:length(v)-1]
    plengths = [sum(lengths[1:i]) for i=1:length(lengths)]

    totallength = plengths[end]

    dx = totallength/npts

    pNpts = floor.(lengths/dx)

    kpts = [v[1] + (v[2]-v[1])*(n-1)/pNpts[1] for n=1:pNpts[1]]
    arc = [norm((v[2]-v[1])*(n-1)/pNpts[1]) for n=1:pNpts[1]]

    for p=2:length(vertices)-1

        append!(kpts, [v[p] + (v[p+1]-v[p])*(n-1)/pNpts[p] for n=1:pNpts[p]])
        append!(arc,
                [
                arc[end] + norm((v[p]-v[p-1])/pNpts[p-1]) + norm((v[p+1]-v[p])*(n-1)/pNpts[p]) for n=1:pNpts[p]]
            )
    end

    append!(kpts, [v[end]])
    append!(arc, [arc[end] + norm((v[end]-v[end-1])/pNpts[end])])
    prepend!(plengths, [0.0])

    return KPath(kpts, arc, plengths, length(kpts))

end

"""
    mpmesh(B::AbstractMatrix, Ns::NTuple{D, Int}) where {D}

Creates a Ns[1]*Ns[2]*... Monkhorst-Pack mesh from the basis vectors B[:, 1], B[:, 2], ...
The mesh is returned as a `KMesh` object.
"""
function mpmesh(B::AbstractMatrix, Ns::NTuple{D, Int}) where {D}

    @assert size(B, 1) == size(B, 2) "number of basis vectors must be the same as a the dimension of the vectors"
    @assert D == size(B, 2) "dimensions of mesh must coincide with the number of basis vectors"

    itr = Iterators.product([1:Ni for Ni in Ns ]...)

    sB = SMatrix{size(B)...}(B)

    list = SVector[]
    for I in itr
        push!(list, sB*SVector(I...))
    end

    len = length(list)

    return KMesh(list, [1/len for k in list], len)
end


text:
module CTImages

export rescale!, rescale, rotate
export polar2cart
export AbstractCTImage, CTImage, CTSinogram, CTTomogram
export CTImageOrTomog, CTImageMat, CTSinogMat, CTTomogMat
export ctimage, ctsinogram, cttomogram

using IntervalSets, SimpleTraits
import ..Monads: mjoin
using ..Utils: linspace, half, _compute_radius, _compute_angle, _wrap_angle
import ..Utils: _atype
using ..Applicative, ..Monads, ..Interpolation
using ..Geometry:
    AbstractParallelBeamGeometry, AbstractFanBeamGeometry, ParallelBeamGeometry
import Base: size, convert, similar, eltype, getindex, setindex!, IndexStyle
using Base: @propagate_inbounds

function ctimage end
function ctsinogram end
function cttomogram end

abstract type AbstractCTImage{T<:Number} <: AbstractMatrix{T} end

@traitimpl Monad{AbstractCTImage}
mjoin(img::AbstractCTImage) = img.data

# struct CTImageStyle <: Broadcast.AbstractArrayStyle{2} end
# Base.BroadcastStyle(::Type{<:AbstractCTImage}) = CTImageStyle()
# CTImageStyle(::Val{2}) = CTImageStyle()
# CTImageStyle(::Val{N}) where N = Broadcast.DefaultArrayStyle{N}()

_atype(img::M) where {M <: AbstractCTImage} = _atype(M)
eltype(::Type{M}) where {M<:AbstractCTImage{T}} where T = T
size(img::AbstractCTImage) = mbind(size, img)
size(img::AbstractCTImage, dim::Int) = size(mjoin(img), dim)
@propagate_inbounds function getindex(img::AbstractCTImage, i::Int)
    @boundscheck checkbounds(img, i)
    @inbounds getindex(mjoin(img), i)
end
@propagate_inbounds function getindex(
    img::AbstractCTImage, I::Vararg{Int,2}
)
    @boundscheck checkbounds(img, I...)
    @inbounds getindex(mjoin(img), I...)
end
@propagate_inbounds function setindex!(img::AbstractCTImage, v, i::Int)
    @boundscheck checkbounds(img, i)
    @inbounds setindex!(mjoin(img), v, i)
end
@propagate_inbounds function setindex!(
    img::AbstractCTImage, v, I::Vararg{Int,2}
)
    @boundscheck checkbounds(img, I...)
    @inbounds setindex!(mjoin(img), v, I...)
end

const ctnames = (image = :CTImage, sinog = :CTSinogram, tomog = :CTTomogram)
const ctfn = (image = :ctimage, sinog = :ctsinogram, tomog = :cttomogram)

for nm in ctnames
    @eval begin
        struct $nm{T<:Number,M<:AbstractMatrix{T}} <: AbstractCTImage{T}
            data::M
        end
        @inline _atype(::Type{$nm{T,M}}) where {T,M} = M
        IndexStyle(::Type{T}) where T<:$nm = IndexStyle(_atype(T))
        @inline function $nm{T,M}(
            ::UndefInitializer,
            dims::Vararg{Union{Integer,AbstractUnitRange},2}
        ) where {T,M}
            $nm(M(undef, dims...))
        end
        @inline $nm{T,M}(img::$nm) where {T,M} = $nm{T,M}(mjoin(img))
        @inline $nm{T}(img::M) where {T,M<:AbstractMatrix{T}} = $nm{T,M}(img)
        @inline $nm(img::$nm) = img
        @inline convert(::Type{M}, img::AbstractMatrix) where {M<:$nm} =
            mreturn(M, img)
        @inline convert(::Type{M}, img::$nm) where {M<:$nm} =
            mreturn($nm, convert(_atype(M), mjoin(img)))
        @inline function similar(img::M) where {M<:$nm}
            mreturn(M, similar(mjoin(img)))
        end
        @inline function similar(
            img::M,
            ::Type{S},
            dims::Tuple{Vararg{Int,2}}
        ) where {S,M<:$nm}
            mreturn(M, similar(mjoin(img), S, dims))
        end
    end
end

CTImage(::Union{CTSinogram,CTTomogram}) =
    @error "Cannot construct a CT image from a sinogram or a tomogram"
CTSinogram(::Union{CTImage,CTTomogram}) =
    @error "Cannot construct a CT sinogram from an image or a tomogram"
CTTomogram(::Union{CTImage,CTSinogram}) =
    @error "Cannot construct a CT tomogram from an image or a sinogram"

const CTImageOrTomog{T,M} = Union{CTImage{T,M},CTTomogram{T,M}}
const CTImageMat{T} = CTImage{T,Matrix{T}}
const CTSinogMat{T} = CTSinogram{T,Matrix{T}}
const CTTomogMat{T} = CTTomogram{T,Matrix{T}}

"""
    rotate(mat::AbstractMatrix, α::Real; <keyword arguments>)

Rotate matrix `mat` about the center of angle `α` given in degrees.
If `rows` and `cols` are not given, the rotated matrix has the same
dimensions of the original matrix.

# Arguments
- `mat`: matrix to rotate.
- `α`: angle in degrees.
- `rows=nothing`: number of rows of the rotated matrix.
- `cols=nothing`: number of columns of the rotated matrix.
- `interpolation`: interpolation strategy. By default is
    `BilinearInterpolation`.
"""
function rotate(
    mat::AbstractMatrix{T},
    α::Real;
    rows::Optional{Integer} = nothing,
    cols::Optional{Integer} = nothing,
    background::Optional{Real} = nothing,
    interpolation::Optional{Interp} = nothing,
) where {T <: Number,Interp <: AbstractInterp2DOrNone}
    orows, ocols = size(mat)
    rows = maybe(orows, rows)
    cols = maybe(ocols, cols)
    sϕ, cϕ = sincos(deg2rad(α))
    interpolation = maybe(interpolate, interpolation)
    interp = interpolation(mat)
    x₀::T = (cols + 1) / 2
    y₀::T = (rows + 1) / 2
    x′₀::T = (ocols + 1) / 2
    y′₀::T = (orows + 1) / 2
    z::T = maybe(zero(T), background)
    rmat = similar(mat, rows, cols)
    fill!(rmat, z)
    for ix ∈ 1:rows, iy ∈ 1:cols
        x = T(ix) - x₀
        y = T(iy) - y₀
        x′ = x * cϕ + y * sϕ + x′₀
        y′ = y * cϕ - x * sϕ + y′₀
        if x′ ∈ 1..ocols && y′ ∈ 1..orows
            rmat[iy, ix] = interp(y′, x′)
        end
    end
    rmat
end

"""
    rescale!(x::AbstractArray;
        interval=nothing, calibration=nothing, window=nothing)

Rescale array x to the interval specified by `interval`.

If `calibration` is not `nothing`, then rescaling is done with
reference to the values given by `calibration`. In other words,
minimum and maximum are assumed to be the values specified by
`calibration`.

See also: [`rescale`](@ref)
"""
function rescale!(
    img::AbstractArray{T};
    interval::ClosedInterval = zero(T)..one(T),
    calibration::Optional{ClosedInterval{U}} = nothing,
    window::Optional{ClosedInterval{W}} = nothing,
) where {T <: Number,U <: Number,W <: Number}
    if isnothing(calibration)
        calibration = ClosedInterval(extrema(img)...)
    end
    a, b = T.(endpoints(interval))
    m, M = T.(endpoints(calibration))
    if m == M
        @warn "Cannot calibrate image as calibration values ($calibration) are equal, image unchanged"
        return img
    end
    f = (b - a) / (M - m)
    if m != zero(T)
        img .-= m
    end
    img .*= f
    if a != zero(T)
        img .+= a
    end
    isnothing(window) && return img
    a, b = T.(endpoints(window))
    map!(img, img) do x
        x < a && return a
        x > b && return b
        return x
    end
    img
end

"""
    rescale(x::AbstractArray, slope::Number, intercept::Number; window)

Linear rescaling of `x` as `x * slope + intercept`.

See also: [`rescale!`](@ref)
"""
function rescale(
    image::AbstractArray{T},
    slope::Number,
    intercept::Number;
    window::Optional{ClosedInterval{U}} = nothing,
) where {T <: Number,U <: Number}
    res = @. image * slope + intercept
    isnothing(window) && return res
    a, b = T.(endpoints(window))
    map!(res, res) do x
        x < a && return a
        x > b && return b
        return x
    end
end

"""
    rescale!(x::AbstractArray, slope::Number, intercept::Number)

In place linear rescaling of `x`.

See also: [`rescale`](@ref)
"""
function rescale!(
    img::AbstractArray{T},
    slope::Number,
    intercept::Number;
    window::Optional{ClosedInterval{U}} = nothing,
) where {T <: Number,U <: Number}
    @. img = slope * img + intercept
    isnothing(window) && return img
    a, b = T.(endpoints(window))
    map!(img, img) do x
        x < a && return a
        x > b && return b
        return x
    end
    img
end

"""
    rescale(x::AbstractArray;
        interval=nothing, calibration=nothing, window=nothing)

Rescale array x to the interval specified by `interval`.

If `calibration` is not `nothing`, then rescaling is done with
reference to the values given by `calibration`. In other words,
minimum and maximum are assumed to be the values specified by
`calibration`.

See also: [`rescale!`](@ref)
"""
function rescale(
    x::AbstractArray;
    interval = 0..1,
    calibration = nothing,
    window = nothing,
)
    rescale!(deepcopy(x); interval, calibration, window)
end

function polar2cart(
    mp::M,
    xs::AbstractVector,
    ys::AbstractVector;
    ν::Real = 1,
    scale::Optional{Real} = nothing,
    diag::Bool = false,
    background::Optional{Real} = nothing,
    transposed::Bool = false,
    interpolation::Optional{Interp} = nothing,
)::M where {
    M <: AbstractMatrix{T},
    Interp <: AbstractInterp2DOrNone,
} where {T <: Real}
    if transposed
        nr, nθ = size(mp)
        mp′ = similar(mp, nθ, nr)
        permutedims!(mp′, mp, (2,1))
        return polar2cart(mp′, xs, ys; background, interpolation)
    end
    ν = maybe(ν, scale)
    @assert ν > 0
    nθ, nr = size(mp)
    x₀, y₀ = half(xs), half(ys)
    κ::T = min(x₀, y₀)
    if diag
        ν *= hypot(x₀, y₀) / κ
    end
    δri::T = (nr-1) / (κ * ν)
    δθi::T = nθ / 2π
    xs′ = xs * δri
    ys′ = ys * δri
    nt1 = nθ + 1
    nthalf = nθ + 1//2
    interpolation = maybe(interpolate, interpolation)
    interp = interpolation(mp)
    rows, cols = length(ys), length(xs)
    z::T = maybe(zero(T), background)
    mc = similar(mp, rows, cols)
    fill!(mc, z)
    Threads.@threads for ix ∈ axes(mc, 2)
        @inbounds @fastmath @simd for iy ∈ axes(mc, 1)
            x::T = xs′[ix]
            y::T = ys′[iy]
            r::T = _compute_radius(x, y)
            θ::T = _compute_angle(x, y, δθi)
            if 1 <= r <= nr && 1 <= θ <= nt1
                if θ >= nthalf
                    θ = one(T)
                elseif θ > nθ
                    θ = nθ
                end
                mc[iy,ix] = interp(θ, r)
            end
        end
    end
    mc
end

@inline function polar2cart(
    mp::AbstractMatrix{T};
    rows::Optional{Integer} = nothing,
    cols::Optional{Integer} = nothing,
    transposed::Bool = false,
    kwargs...
) where {T <: Real,Interp <: AbstractInterp2DOrNone}
    mp = transposed ? permutedims(mp) : mp
    rows = isnothing(rows) ? maybe(2 * size(mp, 2), cols) : rows
    cols = maybe(rows, cols)
    y₀, x₀ = sincos(atan(rows, cols))
    xs = linspace(T, -x₀..x₀, cols)
    ys = linspace(T, -y₀..y₀, rows)
    polar2cart(mp, xs, ys; kwargs...)
end

function polar2cart(
    image::AbstractMatrix,
    geometry::AbstractParallelBeamGeometry;
    kwargs...
)
    polar2cart(image; geometry.rows, geometry.cols, kwargs...)
end

function polar2cart(
    image::AbstractMatrix,
    geometry::AbstractFanBeamGeometry;
    kwargs...
)
    polar2cart(image, ParallelBeamGeometry(geometry); kwargs...)
end

polar2cart(; kwargs...) = x -> polar2cart(x; kwargs...)
polar2cart(g::AbstractParallelBeamGeometry; kwargs...) =
    x -> polar2cart(x, g; kwargs...)

end # module


text:
module Patterns
import Base.&, Base.isequal, Base.>=, Base.>, Base.<=, Base.<, Base.==
using ..Intern

export Node, Predicate, Atom, Guard, never, always
export depsof
export Intension, intension, naught, anything
export encode, predsof, depsof, subs, resultof
export Pattern, suffix_bindings

# ---- Node -------------------------------------------------------------------

abstract Node{T}
typealias Predicate Node{Bool}

@interned type Atom{T} <: Node{T}
    value::T
end
Atom{T}(value::T) = Atom{T}(value)

depsof(node::Atom) = []

const never  = Atom(false)
const always = Atom(true)

@interned type Guard <: Node{Union{}}
    pred::Predicate
end
subs(d::Dict, node::Guard) = Guard(d[node.pred])

resultof(node::Node) = error("Undefined!")

# ---- Intension --------------------------------------------------------------

type Intension
    factors::Dict{Node,Predicate}
end

#const naught   = Intension((Node=>Predicate)[argnode => never])
const naught   = Intension(Dict(always => never))
const anything = Intension(Dict{Node,Predicate}())

predsof(x::Intension) = values(x.factors)

depsof(node::Guard) = [node.pred]
depsof(i::Intension,node::Node) = depsof(node)

function intension(factors::Predicate...)
    gs = Dict{Node,Predicate}()
    for g in factors
        if g === never; return naught; end
        if g === always; continue; end

        new_g = gs[g.arg] = get(gs, g.arg, always) & g
        if new_g === never; return naught; end
    end
    Intension(gs)
end

(&)(x::Intension, y::Intension) = intension(predsof(x)..., predsof(y)...)
isequal(x::Intension, y::Intension) = isequal(x.factors, y.factors)

# works around that type equivalence is weaker than isequal
# works together with (&)(::Isa, ::Isa)
>=(x::Intension, y::Intension) = isequal(x & y, y)     # (x & y) == y
>( x::Intension, y::Intension) = (x >= y) && !(y >= x) # (x >= y) && (x != y)
==(x::Intension, y::Intension) = (x >= y) && (y >= x)
<=(x::Intension, y::Intension) = y >= x
<( x::Intension, y::Intension) = y >  x

# ---- Pattern ----------------------------------------------------------------

type Pattern
    intent::Intension
    bindings::Dict{Symbol,Node}
    rev_bindings::Dict{Node,Symbol}

    function Pattern(intent::Intension, bindings::Dict{Symbol,Node})
        rev = (Node=>Symbol)[node => name for (name,node) in bindings]
        new(intent, bindings, rev)
    end
end
Pattern(intent::Intension) = Pattern(intent, Dict{Symbol,Node}())

function (&)(p::Pattern, q::Pattern)
    bindings = merge(p.bindings, q.bindings)
    @assert length(bindings) == length(p.bindings)+length(q.bindings)
    Pattern(p.intent & q.intent, bindings)
end

function suffix_bindings(p::Pattern, suffix::AbstractString)
    Pattern(p.intent, (Symbol=>Node)[symbol(string(name, suffix)) => node
                                     for (name,node) in p.bindings])
end

end # module


text:
@testset "sinus" begin
    A = rand(10)
    B = rand(10)
    p = rand(10)
    d = KernelFunctions.Sinus(p)
    @test Distances.parameters(d) == p
    @test evaluate(d, A, B) == sum(abs2.(sinpi.(A - B) ./ p))
    @test d(3.0, 2.0) == abs2(sinpi(3.0 - 2.0) / first(p))
end


text:
"""
	$(FUNCTIONNAME)(m::TriMatrix)

Wrap a [`TriMatrix`](@ref) in the appropriate `LinearAlgebra` view type for
optimized performance.
"""
function wraptri end

wraptri(m::TriUpperMatrix) = UpperTriangular(m)
wraptri(m::TriLowerMatrix) = LowerTriangular(m)
wraptri(m::TriSymmetricMatrix) = Symmetric(m)

Base.transpose(m::TriULMatrix) = TriMatrix(transpose_layout(TriLayout(m)), m.n, m.data; diag=m.diag)
Base.transpose(m::TriSymmetricMatrix) = m

LinearAlgebra.istriu(::TriUpperMatrix) = true
LinearAlgebra.istril(::TriLowerMatrix) = true
LinearAlgebra.issymmetric(::TriSymmetricMatrix) = true


text:
const romeo = unwrap # access unwrap function via alias romeo
const romeo! = unwrap!

"""
    mask_from_voxelquality(qmap::AbstractArray, threshold=0.5)

Creates a mask from a quality map. Another option is to use `robustmask(qmap)`

# Examples
```julia-repl
julia> qmap = romeovoxelquality(phase_3echo; TEs=[1,2,3]);
julia> mask = mask_from_voxelquality(qmap);
```

See also [`romeovoxelquality`](@ref), [`romeo`](@ref), [`robustmask`](@ref)
"""
function mask_from_voxelquality(qmap::AbstractArray, threshold=0.2, th2=0.8)
    qmap_bin = qmap .> threshold # NaN defaults to false (0)
    qmap_bin = fill_holes(qmap_bin)
    return gaussiansmooth3d(qmap_bin, (1,1,1)) .> th2
end

function fill_holes(mask; max_hole_size=length(mask) / 20)
    return .!imfill(.!mask, (1, max_hole_size)) # fills all holes up to max_hole_size (uses 6 connectivity as default for 3D)
end

function get_largest_connected_region(mask)
    labels = label_components(mask)
    return labels .== argmax(countmap(labels[labels .!= 0]))
end

function brain_mask(mask)
    bm = get_largest_connected_region(mask)
    return bm
end

function ROMEO.calculateweights(phase::AbstractArray{T,4}; TEs, template=2, p2ref=1, keyargs...) where T
    args = Dict{Symbol, Any}(keyargs)
    args[:phase2] = phase[:,:,:,p2ref]
    args[:TEs] = TEs[[template, p2ref]]
    if haskey(args, :mag) && size(args[:mag], 4) > 1
        args[:mag] = args[:mag][:,:,:,template]
    end
    return ROMEO.calculateweights(view(phase,:,:,:,template); args...)
end

"""
    calculateB0_unwrapped(unwrapped_phase, mag, TEs)

Calculates B0 in [Hz] from unwrapped phase.
TEs in [ms].
The phase offsets have to be removed prior.

See also [`mcpc3ds`](@ref) and [`romeo`](@ref)
"""
function calculateB0_unwrapped(unwrapped_phase, mag, TEs)
    dims = 4
    TEs = to_dim(TEs, 4)
    B0 = (1000 / 2π) * sum(unwrapped_phase .* mag .* mag .* TEs; dims) ./ sum(mag .* mag .* TEs.^2; dims) |> I -> dropdims(I; dims)
    B0[.!isfinite.(B0)] .= 0
    return B0
end

"""
    romeovoxelquality(phase::AbstractArray; keyargs...)

Calculates a quality for each voxel. The voxel quality can be used to create a mask.

# Examples
```julia-repl
julia> qmap = romeovoxelquality(phase_3echo; TEs=[1,2,3]);
julia> mask1 = mask_from_voxelquality(qmap);
julia> mask2 = robustmask(qmap);
```

Takes the same inputs as `romeo`/`unwrap`:
$(@doc unwrap)

See also [`mask_from_voxelquality`](@ref), [`romeo`](@ref), [`robustmask`](@ref)
"""
function romeovoxelquality(phase; keyargs...)
    weights = ROMEO.calculateweights(phase; type=Float32, rescale=x->x, keyargs...) # [0;1]
    qmap = dropdims(sum(weights; dims=1); dims=1)
    qmap[2:end,:,:] .+= weights[1,1:end-1,:,:]
    qmap[:,2:end,:] .+= weights[2,:,1:end-1,:]
    qmap[:,:,2:end] .+= weights[3,:,:,1:end-1]
    return qmap ./ 6 # [0;1]
end


text:
"""
    Modulus

Represents a non-negative integer modulus of up to 61 bits, e.g., for the plain modulus and the
coefficient modulus in instances of `EncryptionParameters`.

See also: [`EncryptionParameters`](@ref)
"""
mutable struct Modulus <: SEALObject
  handle::Ptr{Cvoid}

  function Modulus(value::Integer)
    handleref = Ref{Ptr{Cvoid}}(C_NULL)
    retval = ccall((:Modulus_Create1, libsealc), Clong,
                   (UInt8, Ref{Ptr{Cvoid}}),
                   value, handleref)
    @check_return_value retval
    return Modulus(handleref[])
  end

  function Modulus(handle::Ptr{Cvoid}; destroy_on_gc=true)
    object = new(handle)
    if destroy_on_gc
      finalizer(destroy!, object)
    end
    return object
  end
end

function destroy!(object::Modulus)
  if isallocated(object)
    @check_return_value ccall((:Modulus_Destroy, libsealc), Clong, (Ptr{Cvoid},), object)
    sethandle!(object, C_NULL)
  end

  return nothing
end

module SecLevelType
@enum SecLevelTypeEnum::Cint none=0 tc128=128 tc192=192 tc256=256
end

function bit_count(modulus::Modulus)
  bit_count = Ref{Cint}(0)
  retval = ccall((:Modulus_BitCount, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{Cint}),
                 modulus, bit_count)
  @check_return_value retval
  return Int(bit_count[])
end

function value(modulus::Modulus)
  value = Ref{UInt64}(0)
  retval = ccall((:Modulus_Value, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt64}),
                 modulus, value)
  @check_return_value retval
  return Int(value[])
end

function coeff_modulus_create(poly_modulus_degree, bit_sizes)
  modulusptrs = Vector{Ptr{Cvoid}}(undef, length(bit_sizes))
  retval = ccall((:CoeffModulus_Create, libsealc), Clong,
                 (UInt64, UInt64, Ref{Cint}, Ref{Ptr{Cvoid}}),
                 poly_modulus_degree, length(bit_sizes), collect(Cint, bit_sizes), modulusptrs)
  @check_return_value retval
  modulus = Modulus[Modulus(modulusptrs[i]) for i in 1:length(bit_sizes)]
  return modulus
end

function coeff_modulus_bfv_default(poly_modulus_degree, sec_level=SecLevelType.tc128)
  len = Ref{UInt64}(0)

  # First call to obtain length (modulus result pointer is null)
  retval = ccall((:CoeffModulus_BFVDefault, libsealc), Clong,
                 (UInt64, Cint, Ref{UInt64}, Ptr{Ptr{Cvoid}}),
                 poly_modulus_degree, sec_level, len, C_NULL)
  @check_return_value retval

  # Second call to obtain modulus
  modulusptrs = Vector{Ptr{Cvoid}}(undef, len[])
  retval = ccall((:CoeffModulus_BFVDefault, libsealc), Clong,
                 (UInt64, Cint, Ref{UInt64}, Ref{Ptr{Cvoid}}),
                 poly_modulus_degree, sec_level, len, modulusptrs)
  @check_return_value retval

  modulus = Modulus[Modulus(ptr) for ptr in modulusptrs]
  return modulus
end

function plain_modulus_batching(poly_modulus_degree, bit_size)
  return coeff_modulus_create(poly_modulus_degree, [bit_size])[1]
end


text:
export sacred_load

function absolute_path(path::String)
    if first(path) == '/'
        return path
    else
        return joinpath(pwd(), path)
    end
end

function sacred_load(depot_path)

    depot_path = absolute_path(depot_path)
    depot_entries = readdir(depot_path, join=true)

    df = DataFrame(:name=>Vector{Union{String,Missing}}(),
                   :id=>Vector{Union{Missing, Int}}())

    for entry_path in depot_entries
        first(basename(entry_path)) == '_' && continue
        first(basename(entry_path)) == '.' && continue

        #@info "parsing " id=basename(entry_path)

        el_files = readdir(entry_path)

        @info "loading entry $entry_path"
        data = read(joinpath(entry_path, "run.json"), String)
        el_json = JSON3.read(data)
        config_json = JSON3.read(read(joinpath(entry_path, "config.json"), String))
        experiment = el_json[:experiment]

        push!(df)
        row = nrow(df)

        # Experiment
        setcol!(df, row, :name, experiment[:name])
        setcol!(df, row, :id, parse(Int, basename(entry_path)))
        setcol!(df, row, :status, el_json.status)
        if el_json.status == "RUNNING"
            t = el_json.heartbeat
            cut = findlast('.', t)
            dt = now() - DateTime(t[begin:cut-1], "y-m-dTH:M:S")
            if dt >= Day(5)
                setcol!(df, row, :status, "FAILED $(floor(dt, Day(1)))")
            elseif dt >= Hour(8)
                setcol!(df, row, :status, "RUNNING/FAILED $(floor(dt, Hour(1)))")
            else
                setcol!(df, row, :status, "RUNNING")
            end
        end
        setcol!(df, row, :path, entry_path)

        ks = keys(config_json)
        for k in ks
            setcol!(df, row, k, config_json[k])
        end

        # Sources
        sources_list = experiment[:sources]
        sources = Vector{String}()
        for (source_orig, source_stored) in sources_list
            push!(sources, joinpath(depot_path, source_stored))
        end
        setcol!(df, row, :sources, sources)

        # Metrics stuff
        metrics = MetricsData(joinpath(entry_path, "metrics.json"), MetricsDataType())
        setcol!(df, row, :metrics, metrics)

        # Artifacts stuff
        artifacts = el_json[:artifacts]
        artifacts_paths = Vector{Any}()
        for artifact in artifacts
            path = joinpath(entry_path, artifact)
            push!(artifacts_paths, path)

            if endswith(path, ".log")
                # convert to jsonlog netket
                nklog = MetricsData(path, NetketJsonDataType())
                setcol!(df, row, :nklog, nklog)
            end
        end
        setcol!(df, row, :artifacts, artifacts_paths)
    end

    return df
end


text:
# parameters
struct Params
    λ :: Float64
    ψ :: Float64
    ν :: Float64
    α :: Float64
    ϕ :: Float64
    θˡ:: Float64
    θʰ:: Float64
    θᴱ:: Float64
    ε :: Float64
    ρ :: Float64
    γ :: Float64
    γᴱ:: Float64
    σ :: Float64
    Lˢ:: Float64
    ω :: Float64

    # transformed
    Π :: Float64
    𝛎 :: Array{Float64,1}
    𝛂 :: Array{Float64,1}
    𝚯 :: Array{Float64,1}
    #
    𝚯ⁿˢ :: Array{Float64,1}
    θᴱⁿˢ:: Float64
    ϕⁿˢ :: Float64

    # policy variables
    sⁱ:: Float64
    sᶠ:: Float64
    sᴱ:: Float64
end

function Params(λ,ψ,ν, α, ϕ, θˡ, θʰ, θᴱ, ε, ρ, γ, γᴱ, σ, Lˢ, ω; sⁱ = 0.0, sᶠ = 0.0, sᴱ = 0.0)
    # with policy
    Π = (1.0/(ε-1.0))*((ε-1.0)/ε)^ε
    𝛎 = [0.0, ν]
    𝛂 = [1.0 - α, α]
    𝚯 = [θˡ, θʰ]

    𝚯ⁿˢ  = 𝚯
    θᴱⁿˢ = θᴱ
    ϕⁿˢ  = ϕ

    # R&D subsidy
    𝚯 = ((1.0/(1.0 - sⁱ))^((1.0 - γ)/γ)).*𝚯

    # fixed cost subsidy (exit tax)
    ϕ = (1.0 - sᶠ)*ϕ;

    # entry subsidy
    θᴱ= ((1.0/(1.0-sᴱ))^((1.0 - γᴱ)/γᴱ))*θᴱ

    Params(λ, ψ, ν, α, ϕ, θˡ, θʰ, θᴱ, ε, ρ, γ, γᴱ, σ, Lˢ, ω, Π, 𝛎, 𝛂, 𝚯, 𝚯ⁿˢ, θᴱⁿˢ, ϕⁿˢ, sⁱ, sᶠ, sᴱ)
end

# easy modification of a parametrization
function Params(p::Params, modifiedParams::NamedTuple)
    # copy p
    for name in fieldnames(Params)
        eval(Meta.parse("$name = getproperty($p,:$name)"))
    end
    # replace with modifiedParams
    for name in fieldnames(typeof(modifiedParams))
        eval(Meta.parse("$name = getproperty($modifiedParams,:$name)"))
    end

    Params(λ, ψ, ν, α, ϕ, θˡ, θʰ, θᴱ, ε, ρ, γ, γᴱ, σ, Lˢ, ω; sⁱ = sⁱ, sᶠ = sᶠ, sᴱ = sᴱ)
end

# equilibirum object
mutable struct EqObj
             ws:: Float64
         cactiv:: Array{Float64,1}
            eyq:: Array{Float64,1}
           qbar:: Float64
      cactivtot:: Float64
          cinac:: Float64
              x:: Array{Float64,1}
           xout:: Float64
         optval:: Array{Float64,1}
             τs:: Array{Float64,1}
              τ:: Float64
           qmin:: Array{Float64,1}
              g:: Float64
        solFAll:: ODESolution
          PhiHG:: Float64
          PhiLG:: Float64
       solFRest:: ODESolution
        PhiHExo:: Float64
        PhiLExo:: Float64
       PhiHObso:: Float64
       PhiLObso:: Float64
        qbarAct:: Float64
              r:: Float64
           cfix:: Float64
             wu:: Float64
             cx:: Array{Float64,1}
           crnd:: Float64
           cout:: Float64
    skilled_lab:: Float64
        FAllend:: Float64
        qMinAll:: Float64
        qMaxAll:: Float64
        node   :: Array{Float64,1}
        weight :: Array{Float64,1}
        qAllDist::Array{Float64,1}

    EqObj() = new()
end


text: println(isinteractive())

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
#!/usr/bin/julia

#
## Sidef translated to Julia
#

import Base.-,
       Base.+,
       Base.*,
       Base./,
       Base.//,
       Base.^,
       Base.$,
       Base.&,
       Base.|,
       Base.<=,
       Base.>=,
       Base.<,
       Base.>,
       Base.==,
       Base.!=,
       Base.!,
       Base.sqrt,
       Base.print,
       Base.println,
       Base.setindex!,
       Base.getindex;

abstract Sidef_Object
abstract Sidef_Types_Nil_Nil
abstract Sidef_Types_Bool_True
abstract Sidef_Types_Bool_False

immutable Sidef_Types_Bool_Bool <: Sidef_Object
    value::Bool
end

immutable Sidef_Types_String_String <: Sidef_Object
    value::AbstractString
end

immutable Sidef_Types_Number_Number <: Sidef_Object
    value::Number
end

immutable Sidef_Types_Block_Block <: Sidef_Object
    value::Function
end

immutable Sidef_Types_Array_Array <: Sidef_Object
    value::Array{Any}
end

immutable Sidef_Types_Hash_Hash <: Sidef_Object
    value::Dict{Any,Any}
end

const NIL = Sidef_Types_Nil_Nil
const TRUE = Sidef_Types_Bool_True
const FALSE = Sidef_Types_Bool_False

#
## Object methods
#
function say(s::Sidef_Object)
    println(s.value)
    TRUE
end

function print(s::Sidef_Object)
    print(s.value)
    TRUE
end

#
## String methods
#
function +(a::Sidef_Types_String_String, b::Sidef_Types_String_String)
    Sidef_Types_String_String(a.value * b.value)
end

function *(a::Sidef_Types_String_String, b::Sidef_Types_Number_Number)
    Sidef_Types_String_String(repeat(a.value, Int(b.value)))
end

#
## Number methods
#
for sym in Symbol[:+, :-, :*, :/, ://, :|, :&, :^, :$]
    @eval function $sym(a::Sidef_Types_Number_Number, b::Sidef_Types_Number_Number)
        Sidef_Types_Number_Number($sym(a.value, b.value))
    end
end

for sym in Symbol[:<=, :>=, :<, :>]
    @eval function $sym(a::Sidef_Types_Number_Number, b::Sidef_Types_Number_Number)
        $sym(a.value, b.value) ? TRUE : FALSE
    end
end

function ==(a::Sidef_Types_Number_Number, b::Sidef_Types_Number_Number)
    (a.value == b.value) ? TRUE : FALSE
end

function !=(a::Sidef_Types_Number_Number, b::Sidef_Types_Number_Number)
    (a.value == b.value) ? FALSE : TRUE
end

function sqrt(a::Sidef_Types_Number_Number)
    Sidef_Types_Number_Number(sqrt(a.value))
end

for p in Dict{Symbol, Symbol}(
    :+ => :add,
    :- => :sub,
    :* => :mul,
    :/ => :div,
    :^ => :pow,
)
    @eval function $(p.second)(a::Sidef_Types_Number_Number...)
        Sidef_Types_Number_Number(mapreduce((x) -> x.value, $(p.first), a))
    end
end

function !(a::Sidef_Types_Number_Number)
    Sidef_Types_Number_Number(factorial(a.value))
end

#
## Array methods
#

# TODO: Change "Sidef_Types_Array_Array" to "Type{Sidef_Types_Array_Array}"
function call(h::Sidef_Types_Array_Array, v::Any...)
    Sidef_Types_Array_Array(Any[v...])
end

function getindex(a::Sidef_Types_Array_Array, i::Sidef_Types_Number_Number)
    (a.value)[i.value+1]
end

function setindex!(a::Sidef_Types_Array_Array, v::Any, i::Sidef_Types_Number_Number)
    #setindex!(a.value, v, i.value)
    (a.value)[i.value+1] = v
end

#
## Hash methods
#

function call(h::Type{Sidef_Types_Hash_Hash}, s::Any...)
    d = Dict{Any,Any}()
    for i in 1:2:length(s)
        d[s[i]] = s[i+1]
    end
    Sidef_Types_Hash_Hash(d)
end

function getindex(h::Sidef_Types_Hash_Hash, k::Any)
    (h.value)[k]
end

function setindex!(h::Sidef_Types_Hash_Hash, v::Any, k::Any)
    (h.value)[k] = v
end

#
## Block methods
#
function call(b::Sidef_Types_Block_Block, args...)
    (b.value)(args...)
end

function *(b::Sidef_Types_Block_Block, n::Sidef_Types_Number_Number)
    for i = 1:Int(n.value)
        (b.value)(Sidef_Types_Number_Number(i))
    end
end

const Number345121681 = Sidef_Types_Number_Number(big"1");
const Number388550481 = Sidef_Types_Number_Number(big"2");
const Number387990241 = Sidef_Types_Number_Number(big"12");
const Number387972001 = Sidef_Types_Number_Number(big"0");
const Number387602161 = Sidef_Types_Number_Number(big"5");
const Number387588241 = Sidef_Types_Number_Number(big"3");
const Number387245121 = Sidef_Types_Number_Number(big"4");
const Number387231441 = Sidef_Types_Number_Number(big"42");
const Number387016081 = Sidef_Types_Number_Number(big"32");
const Number387002401 = Sidef_Types_Number_Number(big"1"//big"2");
const Number386991841 = Sidef_Types_Number_Number(big"10");
const String386054321 = Sidef_Types_String_String("a");
const String386056001 = Sidef_Types_String_String("b");
const String431397441 = Sidef_Types_String_String("");
const Number429418641 = Sidef_Types_Number_Number(big"123");
const Number417909201 = Sidef_Types_Number_Number(big"99");
const Number428219841 = Sidef_Types_Number_Number(big"7");
const Hash428257680 = Sidef_Types_Hash_Hash
const String433191281 = Sidef_Types_String_String("c");
const String438638641 = Sidef_Types_String_String(" ");
const String438653681 = Sidef_Types_String_String("+");
const String438637441 = Sidef_Types_String_String("-");
const String438667361 = Sidef_Types_String_String("|");
const String438672401 = Sidef_Types_String_String("/");
const Number455838401 = Sidef_Types_Number_Number(big"6");
const true316266241 = Sidef_Types_Bool_Bool(1);
const false316266001 = Sidef_Types_Bool_Bool(0);

fib37669136 = Function
fac38796504 = Function
cuboid43862520 = Function

(x38722496, z38720984) = (Any,Any)
(x38628888) = (Any)
(arr43140128) = (Any)
(h42824928) = (Any)
(a43321864) = (Any)

fib37669136 = Sidef_Types_Block_Block(function(_38260232::Any...)
_anys38260232 = Any[]
for i in 1:(1 - length(_38260232)) push!(_anys38260232, NIL); end
_38260232 = (_38260232..., _anys38260232...)
    n38260952, = _38260232
    (<=((n38260952),(Number345121681)) == TRUE ?begin n38260952 end:begin +((call((fib37669136),-((n38260952),(Number345121681)))),call((fib37669136),-((n38260952),(Number388550481)))) end)end);
say(call((fib37669136),(Number387990241)));
fac38796504 = Sidef_Types_Block_Block(function(_37719184::Any...)
_anys37719184 = Any[]
for i in 1:(1 - length(_37719184)) push!(_anys37719184, NIL); end
_37719184 = (_37719184..., _anys37719184...)
    n38760864, = _37719184
    (==((n38760864),(Number387972001)) == TRUE ?begin (Number345121681) end:begin *((n38760864),call((fac38796504),-((n38760864),(Number345121681)))) end)end);
say(call((fac38796504),(Number387602161)));
say(^(((Number387588241)),(Number387245121)));
say(((Number387588241))$(Number387602161));
(x38722496, z38720984)=((Number387231441), (Number387016081));
say(-((x38722496),(Number388550481)));
say(z38720984);
say(sqrt(((Number388550481)),));
say((Number387002401));
say(/(((Number386991841)),(Number387588241)));
say(add(((Number345121681)),(Number388550481),(Number387588241)));
say(mul(((Number388550481)),(Number387588241),(Number387245121)));
say(pow(((Number388550481)),(Number387588241)));
say(!(((Number387602161)),));
(x38628888)=((Number387231441));
say(x38628888);
(x38628888) += ((Number388550481));
say(x38628888);
(x38628888) //= ((Number387602161));
say(x38628888);
print((String386054321));
print((String386056001));
say((String431397441));
(x38628888) ^= ((Number388550481));
say(x38628888);
(arr43140128)=(Sidef_Types_Array_Array(Any[(Number429418641), (Number417909201), (Number387231441)]));
say(arr43140128[(Number388550481)]);
((arr43140128[(Number388550481)])=((Number428219841)));
say(arr43140128[(Number345121681)]);
say(arr43140128[(Number388550481)]);
(h42824928)=(call(((Hash428257680)),(String386054321),(Number345121681),(String386056001),(Number388550481)));
say(h42824928[(String386054321)]);
say(h42824928[(String386056001)]);
((h42824928[(String433191281)])=((Number387231441)));
say(h42824928[(String433191281)]);
(a43321864)=(call((Sidef_Types_Array_Array(Any[])),(Number345121681),(Number388550481),(Number387588241)));
say(a43321864);
say(a43321864[(Number345121681)]);
((a43321864[(Number388550481)])=((Number387231441)));
cuboid43862520 = Sidef_Types_Block_Block(function(_43864648::Any...)
_anys43864648 = Any[]
for i in 1:(8 - length(_43864648)) push!(_anys43864648, NIL); end
_43864648 = (_43864648..., _anys43864648...)
    x43864056, y43865848, z43866376, s43863288, c43864936, h38760576, v43866664, d43867168, = _43864648
    (x43864056 == NIL) && (x43864056 = (Number345121681));
    (y43865848 == NIL) && (y43865848 = (Number345121681));
    (z43866376 == NIL) && (z43866376 = (Number345121681));
    (s43863288 == NIL) && (s43863288 = (String438638641));
    (c43864936 == NIL) && (c43864936 = (String438653681));
    (h38760576 == NIL) && (h38760576 = (String438637441));
    (v43866664 == NIL) && (v43866664 = (String438667361));
    (d43867168 == NIL) && (d43867168 = (String438672401));
    say((+((+((+((*(((String438638641)),+((z43866376),(Number345121681)))),c43864936)),*((h38760576),x43864056))),c43864936)),);
    *((Sidef_Types_Block_Block(    function(_45086384::Any...)
_anys45086384 = Any[]
for i in 1:(1 - length(_45086384)) push!(_anys45086384, NIL); end
_45086384 = (_45086384..., _anys45086384...)
        i45086528, = _45086384
        say((+((+((+((+((+((*(((String438638641)),+((-((z43866376),i45086528)),(Number345121681)))),d43867168)),*((s43863288),x43864056))),d43867168)),*((s43863288),-((i45086528),(>((i45086528),y43865848) == TRUE ?begin -((i45086528),y43865848) end:begin (Number345121681) end))))),(==((-((i45086528),(Number345121681))),y43865848) == TRUE ?begin c43864936 end:begin (>((i45086528),y43865848) == TRUE ?begin d43867168 end:begin v43866664 end) end))),)    end)),z43866376);
    say((+((+((+((c43864936),*((h38760576),x43864056))),c43864936)),+((*((s43863288),(<((z43866376),y43865848) == TRUE ?begin z43866376 end:begin y43865848 end))),(<((z43866376),y43865848) == TRUE ?begin v43866664 end:begin (==((z43866376),y43865848) == TRUE ?begin c43864936 end:begin d43867168 end) end)))),);
    *((Sidef_Types_Block_Block(    function(_45523240::Any...)
_anys45523240 = Any[]
for i in 1:(1 - length(_45523240)) push!(_anys45523240, NIL); end
_45523240 = (_45523240..., _anys45523240...)
        i45523384, = _45523240
        say((+((+((+((v43866664),*((s43863288),x43864056))),v43866664)),(>((z43866376),y43865848) == TRUE ?begin (>=((i45523384),z43866376) == TRUE ?begin +((*((s43863288),x43864056)),c43864936) end:begin +((*((s43863288),-((y43865848),i45523384))),d43867168) end) end:begin (>((-((y43865848),i45523384)),z43866376) == TRUE ?begin +((*((s43863288),z43866376)),v43866664) end:begin +((*((s43863288),-((y43865848),i45523384))),(==((-((y43865848),i45523384)),z43866376) == TRUE ?begin c43864936 end:begin d43867168 end)) end) end))),)    end)),y43865848);
    say((+((+((c43864936),*((h38760576),x43864056))),c43864936)),)end);
call((cuboid43862520),(Number388550481),(Number387588241),(Number387245121));
call((cuboid43862520),(Number345121681),(Number345121681),(Number345121681));
call((cuboid43862520),(Number455838401),(Number388550481),(Number345121681));
call((cuboid43862520),(Number388550481),(Number387245121),(Number345121681));
say((true316266241));
say((false316266001));

max_stars_repo_name: trizen/Sidef-Deparse-Julia

max_issues_repo_name: trizen/Sidef-Deparse-Julia

max_forks_repo_name: trizen/Sidef-Deparse-Julia


text:
@testset "GameViewGtk" begin

    @testset "get_world_to_screen_coeffs Tests" begin
        @test string(sb.get_world_to_screen_coeffs(900, 1500, [9.0, 6.0], [1.0, 1.0])) == string((81.81818181818181, [0.0, 422.727]))
    end;

    @testset "world_to_screen Tests" begin
        @test sb.world_to_screen([-4.5, -3.0], [4.5, 3.0], [1.0, 1.0], 75.0, [0.0, 450.0]) == [75.0, 525.0]
    end;

    @testset "screen_to_world Tests" begin
        @test sb.screen_to_world([112.5, 525.0], [4.5, 3.0], [1.0, 1.0], 75.0, [0.0, 450.0]) == [-4.0, -3.0]
    end;

end;