text:
function getSeatingArea(file)
    lines = map(collect, readlines(file))
    rows = length(lines)
    cols = length(lines[1])
    seatingArea::Matrix{Char} = Matrix{Char}(undef, rows, cols)
    for i in 1:rows
        line = lines[i]
        for j in 1:cols
            seat = line[j]
            seatingArea[i,j] = seat
        end
    end
    return seatingArea
end

function isWithin(move, bounds)
    rows = bounds[1]
    cols = bounds[2]
    within = (1 <= move[1] <= rows) && (1 <= move[2] <= cols)
    #println(move, " ", bounds, " ", within)
    return within
end

function isOccupiedSeat(move, seatingArea)
    isOccupied = seatingArea[move[1],move[2]] == '#'
    return isOccupied
end

function adjecentOccupiedSeats(row, col, seatingArea)
    rows = size(seatingArea, 1)
    cols = size(seatingArea, 2)
    bounds = (rows, cols)
    adjSeats::Vector{Int} = []

    NW = (row-1, col-1)
    N =  (row-1, col)
    NE = (row-1, col+1)
    W =  (row,   col-1)
    E =  (row,   col+1)
    SW = (row+1, col-1)
    S =  (row+1, col)
    SE = (row+1, col+1)
    moves = [NW,N,NE,W,E,SW,S,SE]

    for move in moves
        if isWithin(move, bounds) && isOccupiedSeat(move, seatingArea)
            push!(adjSeats,1)
        end
    end
    #printAdjSeats(adjSeats)
    return adjSeats
end

function printAdjSeats(a::Vector{Int})
    if length(a) > 0
        println(a[1],a[2],a[3])
        println(a[4],'*',a[5])
        println(a[6],a[7],a[8])
    else
        println("empty adjecentSeats vector")
    end
end

function applyRulesOnSeat(row, col, seatingArea)
    adjSeats::Vector{Int} = adjecentOccupiedSeats(row, col, seatingArea)
    seatToCheck = seatingArea[row,col]
    if seatToCheck == '#' && reduce(+, adjSeats) >= 4
        return 'L'
    elseif seatToCheck == 'L' && reduce(+, adjSeats) == 0
        return '#'
    else
        return seatToCheck
    end
end

function countOccuppied(seatingArea)
    rows = size(seatingArea, 1)
    cols = size(seatingArea, 2)
    occurrencies = 0
    for i in 1:rows
        for j in 1:cols
            if seatingArea[i,j] == '#'
                occurrencies += 1
            end
        end
    end
    return occurrencies
end

function runRound(seatingArea)
    newSA = copy(seatingArea)
    rows = size(seatingArea, 1)
    cols = size(seatingArea, 2)
    for i in 1:rows
        for j in 1:cols
            newSA[i,j] = applyRulesOnSeat(i, j, seatingArea)
        end
    end
    #printSA(newSA)
    return newSA
end

function printSA(seatingArea)
    rows = size(seatingArea, 1)
    cols = size(seatingArea, 2)
    for i in 1:rows
        for j in 1:cols
            print(seatingArea[i,j])
        end
        print("\n")
    end
    println("")
end

function solveProblemA()
    seatingArea = getSeatingArea("data.txt")
    #printSA(seatingArea)
    newSeatingArea = runRound(seatingArea)
    while(seatingArea != newSeatingArea)
        seatingArea = newSeatingArea
        newSeatingArea = runRound(newSeatingArea)
    end
    occupiedSeats = countOccuppied(newSeatingArea)
    println("succes, result=", occupiedSeats)
end

solveProblemA()


text:
#################### Phylogenetic No-U-Turn Sampler ####################

#################### Types and Constructors ####################
mutable struct PNUTSTune <: SamplerTune
    logf::Union{Function,Missing}
    stepsizeadapter::NUTSstepadapter
    adapt::Bool
    epsilon::Float64
    delta::Float64
    moves::Vector{Int}
    tree_depth::Int
    tree_depth_trace::Vector{Int}
    acc_p_r::Vector{Int}

    PNUTSTune() = new()

    function PNUTSTune(
        x::Vector{T},
        epsilon::Float64,
        logfgrad::Union{Function,Missing};
        target::Real = 0.6,
        tree_depth::Int = 10,
        targetNNI::Float64 = 0.5,
        delta::Float64 = 0.003,
    ) where {T<:GeneralNode}

        new(
            logfgrad,
            NUTSstepadapter(
                0,
                0,
                0,
                NUTS_StepParams(0.5, target, 0.05, 0.75, 10, targetNNI),
            ),
            false,
            epsilon,
            delta,
            Int[],
            tree_depth,
            Int[],
            Int[],
        )
    end
end

PNUTSTune(
    x::Vector{T},
    logfgrad::Function,
    ::NullFunction,
    delta::Float64 = 0.003,
    target::Real = 0.6;
    args...,
) where {T<:GeneralNode} =
    PNUTSTune(x, nutsepsilon(x[1], logfgrad, delta, target), logfgrad; args...)

PNUTSTune(
    x::Vector{T},
    logfgrad::Function,
    delta::Float64,
    target::Real;
    args...,
) where {T<:GeneralNode} =
    PNUTSTune(x, nutsepsilon(x[1], logfgrad, delta, target), logfgrad; args...)

PNUTSTune(x::Vector; epsilon::Real, args...) = PNUTSTune(x, epsilon, missing, args...)

const PNUTSVariate = Sampler{PNUTSTune,T} where {T<:GeneralNode}

#################### Sampler Constructor ####################

"""
    PNUTS(params::ElementOrVector{Symbol}; args...)

Construct a `Sampler` object for PNUTS sampling. The Parameter is assumed to be
a tree.

Returns a `Sampler{PNUTSTune}` type object.

* params: stochastic node to be updated with the sampler.

* args...: additional keyword arguments to be passed to the PNUTSVariate constructor.
"""
function PNUTS(
    params::ElementOrVector{Symbol};
    delta::Float64 = 0.003,
    target::Float64 = 0.6,
    epsilon::Float64 = -Inf,
    args...,
)
    tune = PNUTSTune(
        GeneralNode[],
        epsilon,
        logpdfgrad!;
        delta = delta,
        target = target,
        args...,
    )
    Sampler(params, tune, Symbol[], false)
end

#################### Sampling Functions ####################

#################### Sampling Functions ####################

sample!(v::PNUTSVariate; args...) = sample!(v, v.tune.logfgrad; args...)

function sample!(
    v::PNUTSVariate{<:Vector{<:GeneralNode}},
    logfgrad::Function;
    adapt::Bool = false,
    args...,
)
    tune = v.tune
    adapter = tune.stepsizeadapter
    const_params = tune.stepsizeadapter.params

    if adapter.m == 0 && isinf(tune.epsilon)
        tune.epsilon = nutsepsilon(v.value[1], logfgrad, tune.delta, const_params.δ)
    end
    setadapt!(v, adapt)
    if tune.adapt
        adapter.m += 1

        nuts_sub!(v, tune.epsilon, logfgrad)

        adaptstat = adapter.metro_acc_prob > 1 ? 1 : adapter.metro_acc_prob

        HT = const_params.δ - adaptstat

        HT2 = const_params.τ - adapter.avg_nni

        HT -= abs_adapter(HT2)
        HT /= 2

        η = 1.0 / (adapter.m + const_params.t0)

        adapter.s_bar = (1.0 - η) * adapter.s_bar + η * HT
        x = const_params.μ - adapter.s_bar * sqrt(adapter.m) / const_params.γ

        x_η = adapter.m^-const_params.κ
        adapter.x_bar = (1.0 - x_η) * adapter.x_bar + x_η * x
        tune.epsilon = exp(x)

    else
        if (adapter.m > 0)
            tune.epsilon = exp(adapter.x_bar)
        end

        nuts_sub!(v, tune.epsilon, logfgrad)
    end
    v
end

function setadapt!(v::PNUTSVariate{<:Vector{<:GeneralNode}}, adapt::Bool)
    tune = v.tune
    if adapt && !tune.adapt
        tune.stepsizeadapter.m = 0
        tune.stepsizeadapter.params =
            update_step(tune.stepsizeadapter.params, log(10.0 * tune.epsilon))
    end
    tune.adapt = adapt
    v
end

function nuts_sub!(
    v::PNUTSVariate{<:AbstractArray{<:GeneralNode}},
    epsilon::Float64,
    logfgrad::Function,
)
    x = deepcopy(v.value[1])
    delta = v.tune.delta
    blv = get_branchlength_vector(x)
    r = randn(length(blv))

    blv = get_branchlength_vector(x)
    set_branchlength_vector!(x, molifier.(blv, delta))
    logf, grad = logfgrad(x)
    xminus = Tree_HMC_State(deepcopy(x), r, grad, logf)
    xplus = Tree_HMC_State(deepcopy(x), r, grad, logf)

    lu = log(rand())
    logp0 = hamiltonian(xminus)

    nni = 0
    tnni = 0
    j = 0
    n = 1

    meta = NUTSMeta()

    acc_p_r = 0
    while j < v.tune.tree_depth
        pm = 2 * (rand() > 0.5) - 1
        meta.alpha = 0
        meta.nalpha = 0
        meta.accnni = 0
        meta.nni = 0
        if pm == -1

            xminus, _, xprime, nprime, sprime =
                buildtree(xminus, pm, j, epsilon, logfgrad, logp0, lu, delta, meta)

        else

            _, xplus, xprime, nprime, sprime =
                buildtree(xplus, pm, j, epsilon, logfgrad, logp0, lu, delta, meta)

        end#if pm
        v.tune.stepsizeadapter.metro_acc_prob = meta.alpha / meta.nalpha

        tnni += meta.nni
        if !sprime
            break
        end
        # sprime is true so checking is not necessary

        if rand() < nprime / n
            acc_p_r += 1
            v.value[1] = xprime.x
        end

        n += nprime
        nni += meta.nni

        j += 1

        s = nouturn(xminus, xplus, epsilon, logfgrad, delta)

        if !s
            break
        end
    end
    v.tune.stepsizeadapter.avg_nni = tnni == 0 ? 0.0 : nni / tnni
    push!(v.tune.moves, nni)
    push!(v.tune.tree_depth_trace, j)
    push!(v.tune.acc_p_r, acc_p_r)
    v
end

function buildtree(
    x::Tree_HMC_State{T},
    pm::Int64,
    j::Integer,
    epsilon::Float64,
    logfgrad::Function,
    logp0::Real,
    lu::Real,
    delta::Float64,
    meta::NUTSMeta,
) where {T<:GeneralNode}

    if j == 0
        xprime = transfer(x)
        nni = 0.0
        if !xprime.extended
            nni = refraction!(xprime, pm * epsilon, logfgrad, delta)
        else
            nni = xprime.nni
            xprime.extended = false
        end

        logpprime = hamiltonian(xprime)

        nprime = Int((logp0 + lu) < logpprime)
        sprime = (logp0 + lu) < logpprime + 1000.0

        meta.nni = nni

        xminus = transfer(xprime)
        xplus = transfer(xprime)
        meta.alpha = min(1.0, exp(logpprime - logp0))
        if rand() < meta.alpha
            meta.accnni += nni
        end
        meta.nalpha = 1

    else
        xminus, xplus, xprime, nprime, sprime =
            buildtree(x, pm, j - 1, epsilon, logfgrad, logp0, lu, delta, meta)
        if sprime
            meta1 = NUTSMeta()
            if pm == -1
                xminus, _, xprime2, nprime2, sprime2 =
                    buildtree(xminus, pm, j - 1, epsilon, logfgrad, logp0, lu, delta, meta1)
            else

                _, xplus, xprime2, nprime2, sprime2 =
                    buildtree(xplus, pm, j - 1, epsilon, logfgrad, logp0, lu, delta, meta1)
            end # if pm
            update!(meta, meta1)
            if rand() < nprime2 / (nprime + nprime2)
                transfer!(xprime, xprime2)
            end
            nprime += nprime2

            xminus_bar = transfer(xminus)
            xplus_bar = transfer(xplus)
            sprime = sprime2 && nouturn(xminus_bar, xplus_bar, epsilon, logfgrad, delta)
            if pm == -1
                transfer!(xminus, xminus_bar)
            else
                transfer!(xplus, xplus_bar)
            end

        end #if sprime
    end #if j

    xminus, xplus, xprime, nprime, sprime
end

function nouturn(
    xminus::Tree_HMC_State,
    xplus::Tree_HMC_State,
    epsilon::Float64,
    logfgrad::Function,
    delta::Float64,
)
    curr_l, curr_h = BHV_bounds(xminus.x, xplus.x)

    nni_m = refraction!(xminus, -epsilon, logfgrad, delta)
    nni_p = refraction!(xplus, epsilon, logfgrad, delta)
    xminus.extended = true
    xplus.extended = true
    xminus.nni = nni_m
    xplus.nni = nni_p
    curr_t_l, curr_t_h = BHV_bounds(xminus.x, xplus.x)
    return curr_h < curr_t_l
end

ext: jl

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
module SbmlInterface

using ModelingToolkit
export ODESystem, ODEProblem, ReactionSystem
export @variables, @parameters

using OrdinaryDiffEq
export solve, FunctionMap, Euler, Heun, Ralston, Midpoint, RK4, ExplicitRK, OwrenZen3, OwrenZen4, OwrenZen5,
       BS3, BS5, RK46NL, DP5, Tsit5, DP8, Vern6, Vern7, Vern8, TanYam7, TsitPap8,
       Vern9, Feagin10, Feagin12, Feagin14, CompositeAlgorithm, Anas5, RKO65, FRK65, PFRK87, RKM

ENV["PYTHON"] = Sys.which("python")
using PyCall

const libsbml = PyNULL()
function __init__()
    copy!(libsbml, pyimport("libsbml"))
end

include("sbml2modelingtoolkit.jl")
export simulatesbml
export sbml2odeproblem,sbml2odesystem
export getmodel,getreactions,getparameters,getinitialconditions

end


text:
circle_verts = [(1.41*sin(2π*n/20), 1.41*cos(2π*n/20)) for n=-10:10]
lrnn_verts = [(1.2*sin(2π*n/20), 1.2*(1+cos(2π*n/20))) for n=-10:10]
lstm_verts = vcat([(1.0*sin(2π*n/20), 1.2 + 1.0*cos(2π*n/20)) for n=-10:10],
                  [(1.4*sin(2π*n/20), 1.2 + 1.4*cos(2π*n/20)) for n=-10:10])
lgru_verts = vcat([(0.5*sin(2π*n/20), 1.2 + 1.4*cos(2π*n/20)) for n=-10:10],
                  [(1.0*sin(2π*n/20), 1.2 + 1.0*cos(2π*n/20)) for n=-10:10])

"""
    layerplotattributes()

Retrive plot attributes for each specific type of layer.
"""
layerplotattributes(::Any) = (mrkrsize = 12, mrkrshape =  [:circle], mrkrcolor = [:gray])
layerplotattributes(::Flux.Dense) = (mrkrsize = 12, mrkrshape =  [:circle], mrkrcolor = [:lightgreen])
layerplotattributes(::Flux.RNNCell) = (mrkrsize = 12, mrkrshape =  [Main.Plots.Shape(lrnn_verts), :circle], mrkrcolor = [false, :lightskyblue1])
layerplotattributes(::Flux.LSTMCell) = (mrkrsize = 12, mrkrshape =  [Main.Plots.Shape(lstm_verts), :circle], mrkrcolor = [false, :skyblue2])
layerplotattributes(::Flux.GRUCell) = (mrkrsize = 12, mrkrshape =  [Main.Plots.Shape(lgru_verts), :circle], mrkrcolor = [false, :skyblue3])
layerplotattributes(r::Flux.Recur) = layerplotattributes(r.cell)
layerplotattributes(::Flux.Conv) = (mrkrsize = 10, mrkrshape = [:square], mrkrcolor = [:plum])
function layerplotattributes(s::Symbol)
    if s == :input_layer
        return (mrkrsize = 12, mrkrshape = [Main.Plots.Shape(circle_verts), :rtriangle], mrkrcolor = [false, :yellow])
    elseif s == :output_layer
        return (mrkrsize = 12, mrkrshape = [:rtriangle], mrkrcolor = [:orange])
    else
        return (mrkrsize = 0, mrkrshape = [:none], mrkrcolor = [:none])
    end
end

"""
    get_layer_type(m::Flux.Chain, i::Int)

Return `:input_layer` if "layer number" `i` is zero and the layer `m[i]` itself,
otherwise.

For use with [`layerplotattributes`](@ref), for properly retrieving the
corresponding plot attributes.
"""
get_layer_type(m::Flux.Chain, i::Int) = i == 0 ? :input_layer : m[i]

"""
    plot(m::Flux.Chain, input_data::Union{Nothing,Array} = nothing)

Plot a Flux.Chain neural network according to recipe.
"""
@recipe function plot(m::Flux.Chain, input_data::Union{Nothing,Array} = nothing)
    m = f32(m)
    if input_data !== nothing
        input_data = convert.(Float32, input_data)
        chain_dimensions = get_dimensions(m, input_data)
    else
        chain_dimensions = get_dimensions(m)
    end
    # chain_dimensions = get_dimensions(m, input_data)
    max_width, = maximum(chain_dimensions)
    mg = chaingraph(m, input_data)
    connections = neuron_connections(m, input_data)

    axis --> false
    xrotation --> 60
    xticks --> begin
        ll = 0:length(m)
        (ll,
            vcat(
                ["input \nlayer "],
                ["hidden \n   layer $n" for n in ll[1:end-2]],
                ["output \nlayer   "]
            )
        )
    end
    xlims --> length(m).*(-0.2,1.2)
    yticks --> false
    ylims --> (-0.1,1.2)
    legend --> false

    # draw connections
    @series begin
        seriescolor --> :gray
        dataseries = [
            ([get_prop(mg, e.src, :loc_x), get_prop(mg, e.dst, :loc_x)],
             [get_prop(mg, e.src, :loc_y), get_prop(mg, e.dst, :loc_y)]) for e in edges(mg)]
        return dataseries
    end

    # draw neurons
    @series begin
        scale_sz(sz, max_width) = min(sz, 7.5*sz/max_width)
        seriestype --> :scatter
        markersize --> vcat(
            [scale_sz(layerplotattributes(get_layer_type(m, get_prop(mg, v, :layer_number))).mrkrsize, max_width)
                for v in vertices(mg)],
            [scale_sz(layerplotattributes(get_layer_type(m, get_prop(mg, v, :layer_number))).mrkrsize, max_width)
                for v in vertices(mg)
                    if length(layerplotattributes(get_layer_type(m, get_prop(mg, v, :layer_number))).mrkrshape)>1],
            [scale_sz(layerplotattributes(get_layer_type(m, get_prop(mg, v, :layer_number))).mrkrsize, max_width)
                for v in vertices(mg) if get_prop(mg, v, :layer_number) == length(m)]
        )
        markershape --> vcat(
            [layerplotattributes(get_layer_type(m, get_prop(mg, v, :layer_number))).mrkrshape[1]
                for v in vertices(mg)],
            [layerplotattributes(get_layer_type(m, get_prop(mg, v, :layer_number))).mrkrshape[end]
                for v in vertices(mg)
                    if length(layerplotattributes(get_layer_type(m, get_prop(mg, v, :layer_number))).mrkrshape)>1],
            [layerplotattributes(:output_layer).mrkrshape[end] for v in vertices(mg) if get_prop(mg, v, :layer_number) == length(m)]
        )
        markercolor --> vcat(
            [layerplotattributes(get_layer_type(m, get_prop(mg, v, :layer_number))).mrkrcolor[1]
                for v in vertices(mg)],
            [layerplotattributes(get_layer_type(m, get_prop(mg, v, :layer_number))).mrkrcolor[end]
                for v in vertices(mg)
                    if length(layerplotattributes(get_layer_type(m, get_prop(mg, v, :layer_number))).mrkrshape)>1],
            [layerplotattributes(:output_layer).mrkrcolor[end]
                for v in vertices(mg)
                    if get_prop(mg, v, :layer_number) == length(m)]
        )
        dataseries = vcat(
            [(get_prop(mg, v, :loc_x), get_prop(mg, v, :loc_y)) for v in vertices(mg)],
            [(get_prop(mg, v, :loc_x), get_prop(mg, v, :loc_y)) for v in vertices(mg)
                if length(layerplotattributes(get_layer_type(m, get_prop(mg, v, :layer_number))).mrkrshape)>1],
            [(get_prop(mg, v, :loc_x), get_prop(mg, v, :loc_y)) for v in vertices(mg)
                if get_prop(mg, v, :layer_number) == length(m)]
        )
        return dataseries
    end

    # display layer type
    for ln in 0:length(m)
        @series begin
            nj = chain_dimensions[ln+1]
            series_annotations --> Main.Plots.series_annotations(
                [ln == 0 ? "input" : string(m[ln])], Main.Plots.font("Sans", 8, rotation=20)
            )
            return [ln], [(nj[1]/2 + 1 + max_width/2)/(max_width+1)]
        end
    end
end

@recipe function plot(d::Union{Flux.Dense, Flux.Recur, Flux.RNNCell, Flux.LSTMCell, Flux.GRUCell})
    Flux.Chain(d)
end

"""
    plot((m,s)::Tuple{Flux.Chain,Tuple{Int}})

Plot a Flux.Chain neural network `m` with input dimensions `s`.

Useful when the neural network starts with a variable-input layer, such as a convolutional layer.

YET TO BE IMPLEMENTED.
"""
@recipe function plot((m,s)::Tuple{Flux.Chain,Tuple{Int}})
    nothing
end


text:
"""
==================== SciTools.jl ====================

Author: Damyn Chipman

[insert description here]

===== HEADER FILE =====

Contents:
 : LineSegment.jl

"""
module SciTools

# ===== Imports =====
using LinearAlgebra

# ===== Global Variables =====
const DEG2RAD = pi/180
const RAD2DEG = 180/pi

# ===== Files =====
files = ["LineSegment", "Body2D"]
for header_name in files
    include("SciTools_"*header_name*".jl")
end

# ==================== END OF SciTools.jl ====================
end


text:
using MyApp

push!(ARGS, "arg")
MyApp.julia_main()


text:
# Define species concentrations
C_ER                # μM        ER-located collagen
C_H                 # μM        Collagen-HSP47 complex in ERGIC
C_G                 # μM        Golgi-located collagen
C_PG                # μM        Post-Golgi compartments-located collagen
C_PM                # μM        Plasma membrane-located collagen
C_E                 # μM        Extracellular collagen
H_ER                # μM        ER-located HSP47
H_G                 # μM        Golgi-located HSP47
K                   # μM        Pka
S                   # μM        SEC61A2
T                   # μM        TANGO1
P                   # μM        PDE4D
V                   # μM        VPS33B
M                   # μM        MMP14
D                   # μM        CTSK

# Define mass action kinetics parameters
μ_typ   = 0.0342    # μM        Median typical concentration
k₁      = 0.00822   # μM        SEC61A2 mean concentration
k₂      = 0.0105    # μM        TANGO1 mean concentration
k₃      = 0.0342    # μM        PDE4D mean concentration
k₆      = 0.00335   # μM        VPS33B mean concentration
k₇      = 0.00457   # μM        MMP14 mean concentration
k₁₅     = 0.0342    # μM        CTSK mean concentration
t̅       = 24.0      # h         Circadian period
k₄_a    = 4.84      # h⁻¹       HSP47 protein synthesis rate
k₄_b    = 0.00413   # h⁻¹       ER HSP47 degradation rate
k₄_c    = 0.00413   # h⁻¹       Golgi HSP47 degradation rate
k₅_a    = 1.89      # h⁻¹       PKA protein synthesis rate
k₅_b    = 0.0069    # h⁻¹       PKA protein degradation rate
k₈      = 0.0693    # h⁻¹       Col-I synthesis rate
k₁₀     = 3600.0    # h⁻¹       ERGIC to Golgi Col-I transition rate
k₁₂     = 0.0417    # h⁻¹       Post-golgi compartments to plasma membrane Col-I transition rate
k₅_c    = 29.2      # μM⁻¹      Rate constant for repression of PKA by PDE4D
k₋₉     = 1.22      # μM⁻¹h⁻¹   Col-I-HSP47 dissociation rate
k₁₁     = 1.22      # μM⁻¹h⁻¹   Golgi to post-golgi compartments Col-I transition rate
k₁₃     = 12.2      # μM⁻¹h⁻¹   Plasma membrane to extracellular Col-I transition rate
k₁₄     = 1.22      # μM⁻¹h⁻¹   CTSK-dependent extracellular Col-I degradation rate
k₁₆     = 1.22      # μM⁻¹h⁻¹   Golgi to ER HSP47 transition rate
k₉      = 0.356     # μM⁻²h⁻¹   Col-I-HSP47 complex association rate
G̅                   # Unitless  Mean pulse function value determined by mean value theorem
t₀                  # h         Phase

max_stars_repo_name: chris-revell/SecretoryPathwayGUI

max_issues_repo_name: chris-revell/SecretoryPathwayGUI

max_forks_repo_name: chris-revell/SecretoryPathwayGUI


text:
using MarketBasket
using Base.Test

function _transactions()
      Vector{Symbol}[[:b,:d], [:a,:b,:d], [:b,:c,:d], [:d,:f,:h], [:b,:g,:i],
                    [:b,:i], [:a,:b,:i], [:b,:e,:i], [:b,:c,:e,:i], [:c,:i],
                    [:d,:h], [:d,:h,:f], [:b,:c,:d,:f], [:h], [:d,:g,:h]]
end

function _clusters()
    c1 = Cluster([1,2,3,4,5])
    c2 = Cluster([6,7,8,9,10])
    c3 = Cluster([11,12,13,14,15])
    [c1,c2,c3]
end

import Base: ==
==(c1::Cluster, c2::Cluster) = sort(c1.inds) == sort(c2.inds)

@testset "Clustering" begin
    n_transactions = length(_transactions())
    transactions = Transactions(collect(1:n_transactions), _transactions())
    @test length(transactions) == 15

    @test transactions[6] == [:b,:i]

    min_support = 0.6
    max_ceiling = 0.3
    slr_threshold = 1.5
    clusters = _clusters()
    refine!(clusters, max_ceiling, min_support, transactions, slr_threshold)
    @test clusters[1] == Cluster([1,2,3,13])
    @test clusters[2] == Cluster([6,7,8,9,10,5])
    @test clusters[3] == Cluster([11,12,15,14,4])
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Elfutils_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Elfutils")
JLLWrappers.@generate_main_file("Elfutils", UUID("ab5a07f8-06af-567f-a878-e8bb879eba5a"))
end  # module Elfutils_jll


text:
module JLBoostTrees

import Base: show, *, print, println, +, getproperty
import Base: vcat

export JLBoostTree, JLBoostTreeModel, show, trees
export WeightedJLBoostTree, *, AbstractJLBoostTree, print, println, vcat
export getproperty

# this needs to be defined before everything else
include("tree-and-tree-models.jl")

include("children.jl")
include("filter_tbl_by_splits.jl")
include("get_leaf_nodes.jl")
include("traits.jl")
include("treedepth.jl")

end # end module


text:
module Utils

import ..Types: PkgMirror

const CURRENT = Ref{Union{PkgMirror, Nothing}}(nothing)

current() = CURRENT[]

const CACHEPATH = joinpath(@__DIR__, "../cache")

cachefile(name::AbstractString...) = joinpath(CACHEPATH, name...)

function setcache(value, name::AbstractString...)
    filename = cachefile(name...)
    path = dirname(filename)
    isdir(path) || mkpath(path)
    try
        open(filename, "w") do fo
            println(fo, value)
        end
    catch
        return false
    end
    true
end

function getcache(name::AbstractString...)
    filename = cachefile(name...)
    isfile(filename) || return nothing
    value = ""
    open(filename) do fi
        value = readline(fi)
    end
    value
end

function delcache(name::AbstractString...)
    filename = cachefile(name...)
    isfile(filename) && Base.rm(filename)
end

hascache(name::AbstractString...) = isfile(cachefile(name...))

function download_cache(url::AbstractString, localfile::AbstractString)
    download(url, cachefile(localfile))
end

function replace_in_file(filename::AbstractString, old_new::Pair...)
    x = read(filename, String)
    x = replace(x, old_new...)
    write(filename, x)
end

function __init__()
    isdir(CACHEPATH) || mkdir(CACHEPATH)
end

end


text:
# ============================================================
# Solver
# ============================================================

export SolverSet,
       set_setup,
       set_geometry,
       set_velocity,
       set_property,
       set_ib,
       initialize,
       init_fvm,
       init_ptc!,
       solve!,
       timestep,
       reconstruct!,
       evolve!,
       update!,
       update_boundary!

include("solver_set.jl")
include("solver_reconstruction.jl")
include("solver_evolution.jl")
include("solver_step.jl")
include("solver_update.jl")
include("evolve_boundary.jl")
include("update_boundary.jl")
include("solver_main.jl")
include("solver_init.jl")
include("solver_particle.jl")


text:
a=parse(Int,readuntil(stdin,' '))
b=parse(Int,readline())
print("answer:",a+b)

max_stars_repo_name: Rratic/JuliaGuide

max_issues_repo_name: Rratic/JuliaGuide

max_forks_repo_name: Rratic/JuliaGuide


text:
using Documenter, NbodyGradient

makedocs(sitename="NbodyGradient",
    modules = [NbodyGradient],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    pages = [
        "Index" => "index.md",
        "Tutorials" => ["basic.md", "gradients.md"],
        "API" => "api.md"
    ]
)

deploydocs(
    repo = "github.com/ericagol/NbodyGradient.jl.git",
    devbranch="master"
)


text:
module Items

import SearchLight: AbstractModel, DbId
import Base: @kwdef

using StructTypes

export Item

StructTypes.StructType(::Type{DbId}) = StructTypes.Struct()

@kwdef mutable struct Item <: AbstractModel
  id::DbId = DbId()
  a::String = ""
  b::Int = 0
end

StructTypes.StructType(::Type{Item}) = StructTypes.Mutable()

end

max_stars_repo_name: jaantollander/genie-webapp-db

max_issues_repo_name: jaantollander/genie-webapp-db

max_forks_repo_name: jaantollander/genie-webapp-db


text:
"""
    CircularThermometer{T <: Real} <: AbstractEncoder

Circular thermometer encoder, used to represent real values as binary patterns.
"""
struct CircularThermometer{T <: Real} <: AbstractEncoder
    min::Union{T,AbstractVector{T}}
    max::Union{T,AbstractVector{T}}
    resolution::Int # TODO: This could be unsigned
    block_len::Int

    """
        CircularThermometer{T}(min, max, resolution::Int) where {T <: Real}

    Constructs a circular thermometer encoder. `min` and `max` denotes the
    admissible numerical limits of encoded values and `resolution` specifies the
    number of bits dedicated to the encoding of a scalar. Values to be encoded
    must be of type T.
    """
    function CircularThermometer{T}(min, max, resolution::Int) where {T <: Real}
        max <= min && throw(ArgumentError(
            "the maximum of the circular thermometer may not be lesser or equal to its minimum"
        ))

        resolution <= 1 && throw(DomainError(
            resolution,
            "circular thermometer's resolution must be greater then one"
        ))

        new(min, max, resolution, floor(Int, resolution / 2))
    end

    """
    CircularThermometer{T}(min::AbstractVector{T}, max::AbstractVector{T}, resolution::Int) where {T <: Real}

    Constructs a circular thermometer encoder. `min` and `max` denotes the
    admissible numerical limits of encoded values and `resolution` specifies the
    number of bits dedicated to the encoding of a scalar. Values to be encoded
    are expected to have the same shape as `min` and `max`. Values to be encoded
    must be of type T.
    """
    function CircularThermometer{T}(min::AbstractVector{T}, max::AbstractVector{T}, resolution::Int) where {T <: Real}
        length(min) != length(max) && throw(DimensionMismatch(
            "min and max vectors must have the same length"
        ))

        !all(min .≤ max) && throw(ArgumentError(
            "the maximum of the circular thermometer may not be lesser or equal to its minimum"
        ))

        resolution <= 1 && throw(DomainError(
            resolution,
            "circular thermometer's resolution must be greater then one"
        ))

        new(min, max, resolution, floor(Int, resolution / 2))
    end
end

"""
    CircularThermometer(min::T, max::T, resolution::Int) where {T <: Real}

Constructs a circular thermometer encoder. `min` and `max` denotes the admissible
numerical limits of encoded values and `resolution` specifies the number of
bits dedicated to the encoding of a scalar. Values to be encoded must be of type
T.
"""
CircularThermometer(min::T, max::T, resolution::Int) where {T <: Real} = CircularThermometer{T}(min, max, resolution)

"""
    CircularThermometer(min::Real, max::Real, resolution::Int)

Constructs a circular thermometer encoder. `min` and `max` denotes the admissible
numerical limits of encoded values and `resolution` specifies the number of
bits dedicated to the encoding of a scalar. Values to be encoded must be same as
the promotion of `min` and `max` types.
"""
CircularThermometer(min::Real, max::Real, resolution::Int) = CircularThermometer(promote(min, max)..., resolution)

"""
    CircularThermometer(::Type{T}, min::Real, max::Real, resolution) where {T <: Real}

Constructs a circular thermometer encoder. `min` and `max` denotes the admissible
numerical limits of encoded values and `resolution` specifies the number of
bits dedicated to the encoding of a scalar. First argument determines the type
of values to be encoded.
"""
CircularThermometer(::Type{T}, min::Real, max::Real, resolution) where {T <: Real} = CircularThermometer(convert.(T, (min, max))..., resolution)

"""
    CircularThermometer(::Type{T}, min::Real, max::Real, resolution) where {T <: Real}

Constructs a circular thermometer encoder. `min` and `max` denotes the admissible
numerical limits of encoded values and `resolution` specifies the number of
bits dedicated to the encoding of a scalar. First argument determines the type
of values to be encoded.
"""
CircularThermometer(min::AbstractVector{T}, max::AbstractVector{T}, resolution::Int) where {T <: Real} = CircularThermometer{T}(min, max, resolution)

"""
    CircularThermometer(min::AbstractVector{<:Real}, max::AbstractVector{<:Real}, resolution::Int)

Constructs a circular thermometer encoder. `min` and `max` denotes the admissible
numerical limits of encoded values and `resolution` specifies the number of bits
dedicated to the encoding of a scalar. Values to be encoded are expected to have
the same shape as `min` and `max`. Values to be encoded must be same as the
promotion of `min` and `max` types.
"""
CircularThermometer(min::AbstractVector{<:Real}, max::AbstractVector{<:Real}, resolution::Int) = CircularThermometer(promote(min, max)..., resolution)

"""
    CircularThermometer(::Type{T}, min::AbstractVector{<:Real}, max::AbstractVector{<:Real}, resolution)

Constructs a circular thermometer encoder. `min` and `max` denotes the admissible
numerical limits of encoded values and `resolution` specifies the number of bits
dedicated to the encoding of a scalar. Values to be encoded are expected to have
the same shape as `min` and `max`. First argument determines the type of values
to be encoded.
"""
CircularThermometer(::Type{T}, min::AbstractVector{<:Real}, max::AbstractVector{<:Real}, resolution) where {T <: Real} = CircularThermometer(convert.(Vector{T}, (min, max))..., resolution)

function encode!(encoder::CircularThermometer{T}, min::T, max::T, x::T, pattern::AbstractVector{Bool}) where {T <: Real}
    if x < min
        x = min
    elseif x > max
        x = max
    end

    shifts = floor(Int, (x - min) / (max - min) * encoder.resolution)
    l = shifts % encoder.resolution
    u = (shifts + encoder.block_len + 1) % encoder.resolution
    u = (u == 0) ? encoder.resolution : u

    for i in eachindex(pattern)
        if l < u
            pattern[i] = (l < i < u) ? 1 : 0
        else
            pattern[i] = (1 ≤ i < u || l < i ≤ encoder.resolution) ? 1 : 0
        end
    end

    pattern
end

function encode!(encoder::CircularThermometer{T}, x::T, pattern::AbstractVector{Bool}) where {T <: Real}
    validate_input_output(encoder, x, pattern)

    encode!(encoder, encoder.min[1], encoder.max[1], x, pattern)
end

max_stars_repo_name: RafaelFK/ramnet

max_issues_repo_name: RafaelFK/ramnet

max_forks_repo_name: RafaelFK/ramnet


text:
import YAML

function import_blueprints(fileName::AbstractString)
	# because loading a json is so much faster than loading a yaml for some reason
	blueprints = YAML.load_file(joinpath(dataDir, "blueprints.yaml"))
	json_write(fileName, blueprints; indent=0)
end

#=
import ODBC

const dbName = "esb_DATADUMP"

function import_decryptors(fileName::AbstractString)
	ODBC.connect(dbName)
	res = ODBC.query(
		"""
		select ta.typeID, typeName, attributeName, valueFloat
		from dbo.dgmTypeAttributes ta, dbo.dgmAttributeTypes at, dbo.invTypes it
		where ta.typeID in (select typeID
						    from dbo.invTypes it, dbo.invMarketGroups mg
						    where it.marketGroupID=mg.marketGroupID
							      and mg.marketGroupName='Decryptors')
			  and ta.attributeID=at.attributeID
			  and ta.typeID = it.typeID
		order by typeID, ta.attributeID
		"""
	)
	ODBC.disconnect()
	# ODBC.DataFrames.writetable("data/Decryptors.csv", res)
	decID = -1
	decryptors = Dict{Any, Any}()
	for row = 1:size(res, 1)
		if decID != res[row, :typeID]
			decID = res[row, :typeID]
			dec = Dict{AbstractString, Any}()
			dec["id"] = decID
			dec["name"] = res[row, :typeName]
			dec["attributes"] = Dict{AbstractString, Float64}()
			push!(decryptors, dec)
		end
		dec["attributes"][res[row, :attributeName]] = res[row, :valueFloat]
	end
	json_write(fileName, decryptors)
end

function import_station_assembly_lines(fileName::AbstractString)
	ODBC.connect(dbName)
	res = ODBC.query(
		"""
		select als.stationID, stationName, s.solarSystemID, constellationID, s.regionID, baseTimeMultiplier, baseMaterialMultiplier, baseCostMultiplier, minCostPerHour, ast.activityID, activityName
		from dbo.ramAssemblyLineStations als, dbo.ramAssemblyLineTypes ast, dbo.staStations s, dbo.ramActivities a
		where als.assemblyLineTypeID = ast.assemblyLineTypeID and als.stationID = s.stationID and ast.activityID = a.activityID
		"""
	)
	#ODBC.DataFrames.writetable("data/AssemblyLines.csv", res)
	ODBC.disconnect()
	json_write(fileName, make_array(res))
end

=#


text:
export CnoidalWaveSerreGreenNaghdi, CnoidalSGN
using Elliptic

"""
    CnoidalWaveSerreGreenNaghdi(param; P=1)

Compute the Serre-Green-Naghdi cnoidal wave with prescribed `h₀<h₁<h₂`.
`h₁` is the minimum, `h₂` is the maximum of the wave.
As `h₀ -> h₁`, the cnoidal wave converges towards the solitary wave.
See for instance [Gavrilyuk, Nkonga, Shyue and Truskinovsky](https://doi.org/10.1088/1361-6544/ab95ac).

# Arguments
- `param :: NamedTuple`: parameters of the problem containing `h₀<h₁<h₂` and dimensionless parameters `ϵ` and `μ`, and number of collocation points `N`.
- `P :: Int`: (keyword, optional, default = 1) the number of periods of the cnoidal wave in the constructed mesh.

# Return values
`(η,u,v,mesh,param)` with
- `η :: Vector{Float64}`: surface deformation;
- `u :: Vector{Float64}`: layer-averaged velocity;
- `v :: Vector{Float64}`: derivative of the trace of the velocity potential at the surface;
- `mesh :: Mesh`: mesh collocation points;
- `param :: NamedTuple`: useful parameters
"""
function CnoidalWaveSerreGreenNaghdi(
                param :: NamedTuple;
                P = 1 :: Int)

        ϵ = param.ϵ
        μ = param.μ

        h₀=param.h₀
        h₁=param.h₁
        h₂=param.h₂
        c = sqrt(h₀*h₁*h₂)
        m = sqrt((h₂-h₁)/(h₂-h₀))
        κ = sqrt(3*(h₂-h₀))/(2*c)/sqrt(μ)
        λ = Elliptic.K(m^2)/κ
        mesh = Mesh((L=P*λ,N=param.N))
        formula = h₁ .-1 .+ (h₂-h₁)*(Jacobi.cn.(κ*mesh.x,m^2).^2)

        a₀ = h₀
        a₁ = h₂ - h₀
        formula2 = a₀ .-1 .+ a₁*(Jacobi.dn.(κ*mesh.x,m^2).^2)
        H₀= a₀+a₁*Elliptic.E(m^2)/Elliptic.K(m^2)
        u2 = c*(1 ./ H₀ .- 1 ./(1 .+ formula2))
        param = (h₀=h₀,h₁=h₁,h₂=h₂,a₀=a₀,a₁=a₁,H₀=H₀,c=c,λ=λ,m=m,κ=κ)

        η = formula/ϵ
        h = 1 .+ ϵ*η
        u = c*η./h
		k = mesh.k
        Dx=  1im * k
		F₀ = sqrt(μ)*Dx
        DxF(v) = real.(ifft(F₀ .* fft(v)))
		v = u - 1/3 ./h .* (DxF(h.^3 .*DxF(u)))
        #
        # h2 = 1 .+ ϵ*formula2
        # v2 = u2 - 1/3 ./h2 .* (DxF(h2.^3 .*DxF(u2)))

        return (η,u,v,mesh,param)

end

"""
    CnoidalSGN(param; P=1)

Build the initial data associated with `CnoidalWaveSerreGreenNaghdi(param; P=1)`, of type `InitialData`,
to be used in initial-value problems `Problem(model, initial::InitialData, param)`.
"""
struct CnoidalSGN <: InitialData

	η
	v
	label :: String
	info  :: String

	function CnoidalSGN(param; P=1)
		(η,u,v,mesh,para)=CnoidalWaveSerreGreenNaghdi(param; P)
		init = Init(mesh,η,v)
		label = "Green-Naghdi cnoidal wave"
		info = "Cnoidal travelling wave for the Serre-Green-Naghdi model.\n\
		├─velocity c = $(para.c)\n├───period P = $(2(para.λ))\n\
		├─maximum h₂ = $(para.h₂) (from bottom)\n└─minimum h₁ = $(para.h₁) (from bottom)."

		new( init.η,init.v,label,info  )
	end
end


text:
module Connect

using Sockets

using ..Client: Core, Connection, Version
using ..Reader: read_msg
using ..Requests: startApi

import ..from_ibtime

"""
    connect([host, ]port, clientId, connectOptions="", optionalCapabilities="")

Connect to host `host` on port `port` and set client ID `clientId`.

Return a [`Connection`](@ref) instance.
"""
function connect(host, port::Int, clientId::Int, connectOptions::String="", optionalCapabilities::String="")

  s = Sockets.connect(host, port)

  @assert isopen(s)

  # Init string
  m, M = Int(typemin(Version)), Int(typemax(Version))

  buf = IOBuffer()
  print(buf, m==M ? "v$m" : "v$m..$M")
  isempty(connectOptions) || print(buf, " ", connectOptions)

  # Handshake
  Core.write_one(s, buf, true)

  res = read_msg(s)

  @assert length(res) == 2

  @info "Connected" V=res[1] T=res[2]

  v = parse(Int, res[1])
  m ≤ v ≤ M || error("Unsupported version.")

  ib = Connection(s, clientId, connectOptions, Version(v), from_ibtime(res[2]))

  startApi(ib, clientId, optionalCapabilities)

  ib
end

function connect(port::Int, clientId::Int, connectOptions::String="", optionalCapabilities::String="")

  localip = getalladdrinfo("localhost")

  @assert !isempty(localip)

  connect(localip[1], port, clientId, connectOptions, optionalCapabilities)
end

"""
    disconnect(ib)

Close the socket connection.
"""
disconnect(ib::Connection) = close(ib.socket)

end


text:
"""
	module Boundaries
"""
module Boundaries

export Boundary,
noBC,
DirichletBC,
NeumannBC,
FloquetBC,
MatchedBC,
AbstractBC,
AbstractLocalBC,
AbstractRealBC,
AbstractComplexBC,
PML,
cPML,
noBL,
AbstractBL,
AbstractRealBL,
AbstractRealBC

using ...Defaults
using ..Shapes
using RecipesBase

abstract type AbstractBC end
abstract type AbstractRealBC <: AbstractBC end
abstract type AbstractComplexBC <: AbstractBC end
abstract type AbstractLocalBC <: AbstractRealBC end

struct noBC <: AbstractBC end
struct DirichletBC <: AbstractLocalBC end
struct NeumannBC <: AbstractLocalBC end
struct FloquetBC <: AbstractComplexBC end
struct MatchedBC <: AbstractComplexBC end

# because displaying struct with no fields is broken
Base.show(io::IO,bc::noBC) = print(io,"noBC")
Base.show(io::IO,bc::DirichletBC) = print(io,"Dirichlet")
Base.show(io::IO,bc::NeumannBC) = print(io,"Neumann")
Base.show(io::IO,bc::FloquetBC) = print(io,"Floquet")
Base.show(io::IO,bc::MatchedBC) = print(io,"Matched")

abstract type AbstractBL end
abstract type AbstractRealBL <: AbstractBL end
abstract type AbstractComplexBL <: AbstractBL end

"""
	struct PML(depth) -> pml

absorbing PML with depth `depth`. The location of the PML is later determined
by associating it with a shape in a call to `Boundary`.

see also: [`Boundary`](@ref), [`cPML`](@ref), [`noBL`](@ref)

---

	(::PML)(x,y) -> σ

evaluate the PML conductivity `σ` at the point `x,y` in local frame.
"""
struct PML <: AbstractComplexBL
	depth::Float64 # depth of PML layer
    startx::Float64
    stopx::Float64
	starty::Float64
    stopy::Float64

	function PML(depth,startx,stopx,starty,stopy)
		@assert depth>0 "cannot have PML with 0 depth"
		return new(depth,startx,stopx,starty,stopy)
	end
    PML(depth::Real) = PML(depth,NaN,NaN,NaN,NaN)
    Base.show(io::IO,pml::PML) = print(io,"PML")
end
"""
	struct cPML(depth) -> pml

amplifying conjugate PML with depth `depth`. The location of the cPML is later determined
by associating it with a shape in a call to `Boundary`.

see also: [`Boundary`](@ref), [`PML`](@ref), [`noBL`](@ref)

---

	(::cPML)(x,y) -> σ

evaluate the PML conductivity `σ` at the point `x,y` in local frame
"""
struct cPML <: AbstractComplexBL
	depth::Float64 # depth of PML layer
    startx::Float64
    stopx::Float64
	starty::Float64
    stopy::Float64

    function cPML(depth,startx,stopx,starty,stopy)
        @assert depth>0 "cannot have cPML with 0 depth"
        return new(depth,startx,stopx,starty,stopy)
    end
	cPML(depth::Real) = cPML(depth,NaN,NaN,NaN,NaN)
    Base.show(io::IO,cpml::cPML) = print(io,"cPML")
end
"""
	struct noBL(depth) -> nobl

lack of boundary layer. `depth` is not used, but it there for consistency.
The location is later determined by associating it with a shape in a call to `Boundary`.

see also: [`Boundary`](@ref), [`PML`](@ref), [`cPML`](@ref)

---

	(::noBL)(x,y) -> σ

returns `σ` = 0. This is here for consistency.
"""
struct noBL <: AbstractRealBL
	depth::Float64 # depth of PML layer
    startx::Float64
    stopx::Float64
	starty::Float64
    stopy::Float64

    noBL(args...) = new(0,NaN,NaN,NaN,NaN)
    Base.show(io::IO,nobl::noBL) = print(io,"noBL")
end

TBL = (PML,cPML,noBL)
@eval for tbl ∈ TBL
	function (bl::tbl)(;depth=bl.depth, startx=bl.startx, stopx=bl.stopx, starty=bl.starty, stopy=bl.stopy)
		return tbl(depth,startx,stopx,starty,stopy)
	end
	function (bl::tbl)(x::Number,y::Number)
		if tbl!==noBL
			@assert  isfinite(bl.stopx-bl.startx) || isfinite(bl.stopy-bl.starty) "boundary layer not well defined. check startx and stopx fields (and y)"
			z = isfinite(bl.stopx-bl.startx) ? x : y
			start = isfinite(bl.stopx-bl.startx) ? bl.startx : bl.starty
			stop = isfinite(bl.stopx-bl.startx) ? bl.stopx : bl.stopy
			if tbl<:PML
				return conductivity_profile(z,start,stop)
			elseif tbl<:cPML
				return conj_conductivity_profile(z,start,stop)
			end
		else
			return complex(0.0)
		end
	end
end

"""
	conj(::PML) -> cpml::cPML
	conj(::cPML) -> pml::PML
	conj(::noBL) -> nobl::noBL

conjugate a boundary layer
"""
Base.conj(pml::PML) = cPML(pml.depth,pml.startx,pml.stopx,pml.starty,pml.stopy)
Base.conj(cpml::cPML) = PML(cpml.depth,cpml.startx,cpml.stopx,cpml.starty,cpml.stopy)
Base.conj(nbl::noBL) = noBL()

# used inside PML. This is the place to modify the PML profile
function conductivity_profile(z::Real,start::Real,stop::Real)
	if (start≤stop && start≤z) || (start≥stop && start≥z)
		depth = abs(stop-start)
	    β = 2*sqrt(depth)
	    α = -(1/4)*exp(complex(0,SCALING_ANGLE))*log(EXTINCTION)/depth/((2exp(β)-β)/(2β*expm1(β))-1/β^2)
        u = abs(z-start)/depth
        σ = α*u*expm1(β*u)/expm1(β)
    else
        σ = complex(0.0)
    end
    return σ
end

# used inside cPML
conj_conductivity_profile(x::Real,start::Real,stop::Real) = -conj(conductivity_profile(x,start,stop))

function integrated_conductivity_profile(z::Real,start::Real,stop::Real)
	if (start≤stop && start≤z) || (start≥stop && start≥z)
		depth = abs(stop-start)
		β = 2*sqrt(depth)
		α = -(1/4)*exp(complex(0,SCALING_ANGLE))*log(EXTINCTION)/depth/((2exp(β)-β)/(2β*expm1(β))-1/β^2)
		u = abs(z-start)/depth
		Σ = (α*depth/expm1(β))*(u*exp(β*u)-(1+1/β)*expm1(β*u)/β)
	else
		Σ = complex(0.0)
	end
	return Σ
end

conj_integrated_conductivity_profile(x::Real,start::Real,stop::Real) = -conj(integrated_conductivity_profile(x,start,stop))

"""
	struct Boundary

	Boundary(shape::Shape,bcs,bls) -> bnd

`bcs` is either tuple of boundary conditions, or it is a single boundary condition which applies to all sides.
`bls` is either tuple of boundary layers, or it is a single boundary layer which applies to all sides.

if fewer boundary conditions or layers are supplied than the number of sides of
`shape`, it assumes `noBC` and `noBL`
"""
struct Boundary{TS,TBC,TBL}
    shape::TS
    bcs::TBC
    bls::TBL

	function Boundary(shape::AbstractShape{N},bcs::AbstractBC=noBC(),bls::AbstractBL=noBL()) where N
		bcs = fill(bcs,N)
		bls = fill(bls,N)
		return Boundary(shape,(bcs...,),(bls...,))
	end
	function Boundary(shape::AbstractShape{N},bcs::AbstractBC,bls) where N
		bcs = fill(bcs,N)
		return Boundary(shape,(bcs...,),bls)
	end
	function Boundary(shape::AbstractShape{N},bcs,bls::AbstractBL=noBL()) where N
		bls = fill(bls,N)
		return Boundary(shape,bcs,(bls...,))
	end
    function Boundary(shape::AbstractShape{N},
        bcs::Tuple,
        bls::Tuple
        ) where {N}

		if length(bcs)<N
			@warn "fewer boundary conditions provided ($(length(bcs))) than shape requires ($N). rest assumed to be noBC"
		elseif length(bcs)>N
			@warn "more boundary conditions provided ($(length(bcs))) than shape requires ($N). ignoring all bc's after $N"
			bcs = bcs[1:N]
		end
		if length(bls)>N
			@warn "more boundary layers provided ($(length(bls))) than shape requires ($N). ignoring all bl's after $N"
			bls = bls[1:N]
		end

		bcs_type = promote_type(typeof.(bcs)...)
		length(bcs)≤N ? bcs_type=promote_type(bcs_type,noBC) : nothing
		BCS = bcs_type[bcs...]
		for i ∈ 1:(N-length(bcs))
			append!(BCS,[noBC()])
		end

        bls_type = promote_type(typeof.(bls)...)
        length(bls)≤N ? bls_type=promote_type(bls_type,noBL) : nothing
        BLS = bls_type[bls...]
		for i ∈ 1:(N-length(bls))
	        append!(BLS,[noBL()])
        end

		fix_bls_by_shape!(BLS,shape)

		BLS = (BLS...,)
        BCS = (BCS...,)

        return new{typeof(shape),typeof(BCS),typeof(BLS)}(shape,BCS,BLS)
    end

	function Base.show(io::IO,bnd::Boundary)
		N = get_number_of_sides(bnd.shape)
		print(io,"Boundary: \n")
		print(io,"\tshape: $(bnd.shape)")
		for i ∈ 1:N
			println(io)
			print(io,"\tside ",i,": ")
			print(io,bnd.bcs[i],"/",bnd.bls[i])
		end
	end
end

"""
	conj(::Boundary) -> bnd

Make a new `Boundary` object with conjugated boundary layers.
"""
Base.conj(bnd::Boundary) = Boundary(bnd.shape,bnd.bcs,conj.(bnd.bls))

# set location of boundary layer. Valid only for squares, rectangles, circle,
# since for now these are the only shapes supporting PMLs/cPMLs.
fix_bls_by_shape!(bls,s::Square) = fix_bls_by_shape!(bls,Rectangle(s))
function fix_bls_by_shape!(bls,shape::Rectangle)
	stop = shape.x0-shape.a/2
	start = stop + bls[1].depth
	bls[1] = bls[1](startx=start, stopx=stop)

	stop = shape.x0+shape.a/2
	start = stop - bls[2].depth
	bls[2] = bls[2](startx=start, stopx=stop)

	stop = shape.y0-shape.b/2
	start = stop + bls[3].depth
	bls[3] = bls[3](starty=start, stopy=stop)

	stop = shape.y0+shape.b/2
	start = stop - bls[4].depth
	bls[4] = bls[4](starty=start, stopy=stop)
end
function fix_bls_by_shape!(bls,c::Circle)
	stop=c.R
	start=stop-bls[1].depth
	bls[1] = bls[1](startx=start,stopx=stop)
end
fix_bls_by_shape!(bls,shape::AbstractShape) = nothing

"""
	plot(::Boundary)

PML's plotted in blue patches (b/c absorbing, so cold), cPML's in red (b/c emitting, so hot)

DirichletBC are solid black lines, NeumannBC dashed, FloquetBC dashdot, MatchedBC/noBC nothing
"""
@recipe function f(bnd::Boundary)
	aspect_ratio --> 1
	legend --> false

	@series begin
		fillcolor --> BOUNDARY_COLOR
		color --> BOUNDARY_COLOR
		markershape := :none
		bnd.shape
	end
	@series (bnd,1)
	@series (bnd,1,1)
end
# plot boundary layers
@recipe function f(bnd::Boundary,_1::Int)
	shape = bnd.shape
	cosθ = shape.cosθ
	sinθ = shape.sinθ
	typeof(shape)<:Square ? shape = Rectangle(shape) : nothing
	for i ∈ eachindex(bnd.bls)
		bls = bnd.bls[i]
		depth = bnd.bls[i].depth
		if typeof(shape)<:Rectangle
			if i==1
				x = -(shape.a-depth)/2
				y = 0
				bl_shape = Rectangle(depth,shape.b,shape.x0+cosθ*x-sinθ*y,shape.y0+sinθ*x+cosθ*y,shape.θ)
			elseif i==4
				x = 0
				y = (shape.b-depth)/2
				bl_shape = Rectangle(shape.a,depth,shape.x0+cosθ*x-sinθ*y,shape.y0+sinθ*x+cosθ*y,shape.θ)
			elseif i==2
				x = +(shape.a-depth)/2
				y = 0
				bl_shape = Rectangle(depth,shape.b,shape.x0+cosθ*x-sinθ*y,shape.y0+sinθ*x+cosθ*y,shape.θ)
			else
				x = 0
				y = -(shape.b-depth)/2
				bl_shape = Rectangle(shape.a,depth,shape.x0+cosθ*x-sinθ*y,shape.y0+sinθ*x+cosθ*y,shape.θ)
			end
		elseif typeof(shape)<:Circle
			bls = bnd.bls[1]
			depth = bls.depth
			bl_shape = Annulus(shape.R-depth,shape.R,shape.x0,shape.y0)
		end
		if typeof(bls)<:AbstractComplexBL
			@series begin
				alpha --> 0
				typeof(bls)<:PML ? fillcolor --> BOUNDARY_PML_COLOR : fillcolor --> BOUNDARY_CPML_COLOR
				typeof(bls)<:PML ? color --> BOUNDARY_PML_COLOR : color --> BOUNDARY_CPML_COLOR
				bl_shape
			end
		end
	end
end
# plot boundary conditions
@recipe function f(bnd::Boundary,_1::Int,_2::Int)
	shape = bnd.shape
	for i ∈ eachindex(bnd.bcs)
		bcs = bnd.bcs[i]
		if typeof(shape)<:Rectangle
			if i==1
				x = [-1,-1]*shape.a/2
				y = [-1,+1]*shape.b/2
			elseif i==2
				x = [-1,+1]*shape.a/2
				y = [+1,+1]*shape.b/2
			elseif i==3
				x = [+1,+1]*shape.a/2
				y = [-1,+1]*shape.b/2
			else
				x = [-1,+1]*shape.a/2
				y = [-1,-1]*shape.b/2
			end
			if typeof(bcs)<:AbstractLocalBC
				@series begin
					if typeof(bcs)<:DirichletBC
						ls --> BOUNDARY_DIRICHLET_LINETYPE
					elseif typeof(bcs)<:NeumannBC
						ls --> BOUNDARY_NEUMANN_LINETYPE
					else
						ls --> BOUNDARY_FLOQUET_MATCHED_LINETYPE
					end
					seriestype --> :line
					markershape := :none
					color --> BOUNDARY_BC_COLOR
					shape.x0 .+ shape.cosθ*x-shape.sinθ*y, shape.y0 .+ shape.sinθ*x+shape.cosθ*y
				end
			end
		else
			if typeof(bcs)<:Union{AbstractLocalBC,FloquetBC}
				@series begin
					if typeof(bcs)<:DirichletBC
						ls --> BOUNDARY_DIRICHLET_LINETYPE
					elseif typeof(bcs)<:NeumannBC
						ls --> BOUNDARY_NEUMANN_LINETYPE
					else
						ls --> BOUNDARY_FLOQUET_MATCHED_LINETYPE
					end
					color --> BOUNDARY_BC_COLOR
					markershape := :none
					fillalpha --> 0
					alpha --> 1
					shape
				end
			end
		end
	end
end

end # module


text:
using SPNetworks
using Documenter

makedocs(;
    modules=[SPNetworks],
    authors="Denis Maua <denis.maua@gmail.com>",
    repo="https://github.com/denismaua/SPNetworks.jl/blob/{commit}{path}#L{line}",
    sitename="SPNetworks.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://denismaua.github.io/SPNetworks.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/denismaua/SPNetworks.jl",
)


text:
using Printf
using HTTP

function get_destination(i, create)
  parent = (
    joinpath(
      joinpath(
        dirname(dirname(Base.source_path())), "data"
      ), @sprintf("case_%05d", i)
    )
  )
  file = joinpath(parent, "imaging.nii.gz")

  if create
    mkpath(parent)
  end

  return file
end

function download(i)
  imaging_url = @sprintf("https://kits19.sfo2.digitaloceanspaces.com/master_%05d.nii.gz", i)
  tmp = tempdir()
  HTTP.download(imaging_url, tmp)
  mv(joinpath(tmp, @sprintf("master_%05d.nii.gz", i)), get_destination(i, true))
end

n_to_download = 0
left_to_download = Int64[]
for i = 0:299
  dst = get_destination(i, false)
  if !isfile(dst)
    append!(left_to_download, i)
    global n_to_download = n_to_download + 1
  end
end

println(@sprintf("%d cases to download...", n_to_download))
for i = 1:n_to_download
  println(@sprintf("%d/%d...", i, n_to_download))
  download(left_to_download[i])
end


text:
# -*- coding: utf-8 -*-
# Wx.jl

module Wx

export WxApp

function WxApp()
  true
end

end


text:
"""
    RiemannProblem{Model<:AbstractBalanceLaw{1}, U, T<:Real}

A Riemann problem in one space dimension for `model` with left and right values
`uₗ`,`uᵣ` at `x₀`,`t₀`.
"""
struct RiemannProblem{Model<:AbstractBalanceLaw{1}, U, T<:Real} <: AbstractRiemannProblem
  model::Model
  uₗ::U
  uᵣ::U
  x₀::T
  t₀::T
end

"""
    RiemannProblem(model::AbstractBalanceLaw{1}, uₗ, uᵣ, x₀::Real=0, t₀::Real=0)

Create the `RiemannProblem` for `model` with left and right values  `uₗ`,`uᵣ` at
`x₀`,`t₀`.
"""
function RiemannProblem(model::AbstractBalanceLaw{1}, uₗ, uᵣ, x₀::Real=0, t₀::Real=0)
  @assert(typeof(uₗ) == typeof(uᵣ))
  x₀, t₀ == promote(x₀, t₀)
  RiemannProblem{typeof(model),typeof(uₗ),typeof(x₀)}(model, uₗ, uᵣ, x₀, t₀)
end

function Base.convert(::Type{RiemannProblem{Model,U,T1}}, prob::RiemannProblem{Model,U,T2}) where {Model,U,T1,T2}
  RiemannProblem{Model,U,T1}(prob.model, prob.uₗ, prob.uᵣ, convert(T1,prob.x₀), convert(T1,prob.t₀))
end

function Base.promote_rule(::Type{RiemannProblem{Model,U,T1}}, ::Type{RiemannProblem{Model,U,T2}}) where {Model,U,T1,T2}
  RiemannProblem{Model,U,promote_type(T1,T2)}
end

function (prob::RiemannProblem)(x::Real)
  @unpack x₀, uₗ, uᵣ = prob

  x < x₀ ? uₗ : uᵣ
end

function Base.show(io::IO, prob::RiemannProblem)
  print(io,
    "Riemann problem for ", prob.model, "\n  ",
    "with uₗ = ", prob.uₗ, " and uᵣ = ", prob.uᵣ, "\n  ",
    "at (x₀,t₀) = (", prob.x₀, ", ", prob.t₀, ")")
end

################################################################################

"""
An abstract type representing the solution of a scalar Riemann problem.
"""
abstract type ScalarRiemannSolution <: AbstractRiemannSolution end

function Base.show(io::IO, sol::AbstractRiemannSolution)
  σ⁻, σ⁺ =  minmax_speeds(sol)
  print(io, "Solution  of the ", sol.prob,
        "\n  speeds: σ⁻ = ", σ⁻, ", σ⁺ = ", σ⁺ )
end

@recipe function f(sol::AbstractRiemannSolution)
  xguide --> L"\xi"
  label  --> L"u"
  legend --> false

  σ⁻, σ⁺ =  minmax_speeds(sol)
  if σ⁻ ≈ σ⁺
    σ = one(σ⁻)
  else
    σ = min(abs(σ⁻ - σ⁺), max(abs(σ⁻), abs(σ⁺))/10)
  end
  σ⁻ = σ⁻ - σ
  σ⁺ = σ⁺ + σ

  ξ = range(σ⁻, σ⁺, length=10^3)
  u = sol.(ξ)

  ((ξ, u, sol.prob.model),)
end

@recipe function f(ξumodel::Tuple{Ξ,U,AbstractBalanceLaw{1}}) where {Ξ,U}
    ξ, u, model = ξumodel
    ((ξ, u), )
end

################################################################################

"""
    RiemannProblemTuple{N, Model<:AbstractBalanceLaw{1}, U, T<:Real}

An NTuple of consecutive Riemann problems.
"""
struct RiemannProblemTuple{N, Model<:AbstractBalanceLaw{1}, U, T<:Real}
  tup::NTuple{N, RiemannProblem{Model,U,T}}

  function RiemannProblemTuple{N,Model,U,T}(tup::NTuple{N, RiemannProblem{Model,U,T}}) where{N,Model,U,T}
    for i in 1:N-1
      tup[i].t₀ == tup[i+1].t₀ ||
        error("The Riemann problems must have the same initial time `t₀`.")
    end
    for i in 1:N-1
      tup[i].x₀ < tup[i+1].x₀ ||
        error("The Riemann problems must be sorted in increasing order.")
    end
    for i in 1:N-1
      tup[i].uᵣ == tup[i+1].uₗ ||
        error("The Riemann problems must be compatible.")
    end
    new(tup)
  end
end

function RiemannProblemTuple(tup::NTuple{N,RiemannProblem{Model,U,T}}) where {N,Model,U,T}
  RiemannProblemTuple{N,Model,U,T}(tup)
end

Base.length(tup::RiemannProblemTuple) = length(tup.tup)
Base.iterate(tup::RiemannProblemTuple) = iterate(tup.tup)
Base.iterate(tup::RiemannProblemTuple, state) = iterate(tup.tup, state)

function *(prob1::RiemannProblem{Model,U}, prob2::RiemannProblem{Model,U}) where {Model,U}
  RiemannProblemTuple(promote(prob1, prob2))
end

function *(tup::RiemannProblemTuple{N,Model,U}, prob::RiemannProblem{Model,U}) where {N,Model,U}
  RiemannProblemTuple(promote(tup..., prob))
end

function *(prob::RiemannProblem{Model,U}, tup::RiemannProblemTuple{N,Model,U}) where {N,Model,U}
  RiemannProblemTuple(promote(prob, tup...))
end

function *(tup1::RiemannProblemTuple{N,Model,U}, tup2::RiemannProblemTuple{M,Model,U}) where {N,M,Model,U}
  RiemannProblemTuple(promote(tup1..., tup2...))
end

function Base.show(io::IO, prob::RiemannProblemTuple{N,Model,U,T}) where {N,Model,U,T}
  println(io, "Tuple of ", N, " consecutive Riemann problems:")
  for i in 1:N
    print(io, prob.tup[i], "\n")
  end
end

"""
An NTuple of consecutive Riemann solutions.
"""
struct RiemannSolutionTuple{N,Sol<:AbstractRiemannSolution}
  tup::NTuple{N, Sol}

  function RiemannSolutionTuple{N,Sol}(tup::NTuple{N, Sol}) where {N,Sol}
    for i in 1:N-1
      tup[i].prob.x₀ < tup[i+1].prob.x₀ ||
        error("The Riemann problems must have the same initial time `t₀`.")
    end
    for i in 1:N-1
      tup[i].prob.x₀ < tup[i+1].prob.x₀ ||
        error("The Riemann problems must be sorted in increasing order.")
    end
    new(tup)
  end
end

function RiemannSolutionTuple(tup::NTuple{N,Sol}) where {N,Sol<:AbstractRiemannSolution}
  RiemannSolutionTuple{N,Sol}(tup)
end

function Base.show(io::IO, prob::RiemannSolutionTuple{N,Sol}) where {N,Sol}
  println(io, "Tuple of ", N, " solutions to consecutive Riemann problems:")
  for i in 1:N
    println(io, prob.tup[i])
  end
end

# for broadcasting; treat as scalar
Base.broadcastable(sol::RiemannSolutionTuple) = Ref(sol)

function solve(prob::RiemannProblemTuple)
  RiemannSolutionTuple(map(solve, prob.tup))
end

@inline get_t₀(sol::RiemannSolutionTuple) = sol.tup[1].prob.t₀

"""
Tha maximal time before any of the solutions interact.
"""
function compute_tmax(sol::RiemannSolutionTuple{N,Sol}) where {N,Sol}
  tup = sol.tup
  t₀ = get_t₀(sol)

  tmax = typemax(tup[1].prob.t₀)
  for i in 1:N-1
    _, σₗ⁺ = minmax_speeds(tup[i])
    σᵣ⁻, _ = minmax_speeds(tup[i+1])
    if σₗ⁺ > σᵣ⁻
      tmax = min(tmax, t₀+(tup[i+1].prob.x₀-tup[i].prob.x₀)/(σₗ⁺-σᵣ⁻))
    end
  end
  tmax
end

function (sol::RiemannSolutionTuple{N,Sol})(t::Real, x::Real) where {N,Sol}
  tup = sol.tup
  t₀ = sol.tup[1].prob.t₀

  t > compute_tmax(sol) && error("The solutions of the Riemann problems interact.")

  for i in 1:N
    σ⁻, σ⁺ = minmax_speeds(tup[i])
    if (x-tup[i].prob.x₀)/(t-t₀) <= σ⁺
      return tup[i](t,x)
    end
  end
  tup[end](t,x)
end

@recipe function f(sol::RiemannSolutionTuple)
  xguide --> L"x"
  yguide --> L"u"

  t₀ = get_t₀(sol)
  tmax = compute_tmax(sol)
  t = t₀ + (tmax-t₀)/2

  label --> "\$ u($t) \$"
  legend --> false

  σ⁻, _ = minmax_speeds(sol.tup[1])
  xmin = sol.tup[1].prob.x₀ + (t-t₀)*σ⁻
  _, σ⁺ = minmax_speeds(sol.tup[end])
  xmax = sol.tup[end].prob.x₀ + (t-t₀)*σ⁺

  if xmin ≈ xmax
    Δx = one(xmin)
  else
    Δx = min(abs(xmin - xmax), max(abs(xmin), abs(xmax))/10)
  end
  xmin = xmin - Δx
  xmax = xmax + Δx

  x = range(xmin, xmax, length=10^3)
  u = sol.(t,x)

  ((x, u, sol.tup[1].prob.model),)
end


text:
module VortexCollisions

    using Distributed
    using LinearAlgebra
    using SharedArrays

    export Grid2d, get_field, frequencies

    include("grid.jl")

    export FourierTransform, get_trans, prfft!, frfft!, irfft!

    include("fourier_transform.jl")

    export apply_operator!, get_gradient, get_inverse_laplacian

    include("fourier_operators.jl")

    export fourier_quadrature

    include("fourier_quadrature.jl")

    export trapezoidal_quadrature

    include("trapezoidal_quadrature.jl")

    export CollisionOperator

    include("collision_operator.jl")

    export FokkerPlanckOperator, collision_operator!

    include("fokker_planck_operator.jl")

    export timestep!

    include("time_integration.jl")

    export initialize_workers, run_simulation

    include("simulation.jl")

end


text:
using DataFrames
using Test
include("../src/Plots.jl")

dict = Dict("I1" => [1,2,3,4], "I2" => [3,3,3,3])
expected_out = DataFrame(I1 = [1,2,3,4], I2 = [3,3,3,3])

@test DictToDf(dict) == expected_out

input = DataFrame(Genes = ["G1", "G1"], Isoforms = ["I1", "I3"],
Cell1 = [5,0],Cell2 = [1,1])

@test MeanIsoPerGenePerCell(input) == 1.5

genes_dict = Dict("G1" => dict)
expected_out = [2]

@test MeanNumIsoformsPerGenePerCellReal(genes_dict)[1] == expected_out

genes_dict = Dict("G1" => dict, "G2" => Dict("I1" => [0,0,1,1]))

@test MeanNumIsoformsPerGenePerCellReal(genes_dict)[1] == [2, 0.5]
@test MeanNumIsoformsPerGenePerCellReal(genes_dict)[2] == Dict("G1" => 2, "G2" => 0.5)

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
#-----------------------------------------------------------------------------------
# Tuning the constant stepsize for SGD and SMD
#-----------------------------------------------------------------------------------

using Random
using PyPlot
using LaTeXStrings, Printf

include("solve_cov_est.jl");

#-------------------------------------
#   Parameters
#-------------------------------------
maxIter = 100;
r = 2;  # rank
d = 100;
stoch_err = 0.01;  # standard deviation of errors in stochastic measurements b
init_radius = 1.0;  # 2-norm measure of relative deviation between Xtrue and Xinit

subgrad_stepsizes = [1.25, 1.5, 1.75] * 1e-4;
mirror_stepsizes = [0.4, 0.5, 0.6];
step_colors = ["#1f78b4"] #, "#33a02c", "red"]; #, "blue"];

Random.seed!(321);  # for reproducibility

#-------------------------------------
#   Initialize plots
#-------------------------------------
subgrad_plot = figure(figsize=[10,6]);
xlabel(L"Iteration $k$");
ylabel(L"Relative distance to solution set $\min_U \; \frac{\| X_k U - X_{true} \|_2^2 }{\| X_{true} \|_2^2}$");
title(@sprintf("SGD for covariance estimation (r=%i, d=%i)", r,d));
xlim(0,maxIter)
ylim(1e-4, 1e1)

mirror_plot = figure(figsize=[10,6]);
xlabel(L"Iteration $k$");
ylabel(L"Relative distance to solution set $\min_U \; \frac{\| X_k U - X_{true} \|_2^2 }{\| X_{true} \|_2^2}$");
title(@sprintf("SMD for covariance estimation (r=%i, d=%i)", r,d));
xlim(0,maxIter)
ylim(1e-4, 1e1)

#-----------------------------------------------------------------------------------------
#   Run stochastic subgradient method for each dimension
#-----------------------------------------------------------------------------------------
for i=1:length(step_colors)
    ηsubgrad = subgrad_stepsizes[i];
    ηmirror = mirror_stepsizes[i];

    # Generate true matrix we are searching for
    Xtrue = randn(d,r);

    # Generate initial point
    pert = randn(d,r);
    Xinit = Xtrue + init_radius * (norm(Xtrue, 2) / norm(pert, 2)) * pert;

    # Run the two methods
    (subgrad_err_hist, ~) = solve_cov_est(Xinit, Xtrue, stoch_err,  maxIter,
                                                   fill(ηsubgrad, maxIter), method="subgradient", clipped=true)
    (mirror_err_hist, ~) = solve_cov_est(Xinit, Xtrue, stoch_err,  maxIter,
                                                    fill(ηmirror, maxIter), method="mirror", clipped=true)

    # Add errors to plots
    plt[:figure](subgrad_plot[:number])
    semilogy(subgrad_err_hist, color=step_colors[i], label=@sprintf("step=%1.2e", ηsubgrad));

    plt[:figure](mirror_plot[:number])
    semilogy(mirror_err_hist, color=step_colors[i], label=@sprintf("step=%1.2e", ηmirror));
end

# Add legends to plots and save final plots
plt[:figure](subgrad_plot[:number])
legend(loc="upper right")
savefig("plots/sgd_stepsize_test.pdf");

plt[:figure](mirror_plot[:number])
legend(loc="upper right")
savefig("plots/smd_stepsize_test.pdf");

max_stars_repo_name: kelliemac/stochastic-model-based

max_issues_repo_name: kelliemac/stochastic-model-based

max_forks_repo_name: kelliemac/stochastic-model-based


text:
abstract type AbstractHypercubeTransform end

"""
    `$(FUNCTIONNAME)(c::AbstractHypercubeTransform)`
Returns the dimension of the hypercube transform.
"""
function dimension end

"""
    `$(FUNCTIONNAME)(c)`
Constructs the object that contains the necessary information to move from
the unit hypercube to the distribution space. This is the usual function to use
when construct the transformation.

There are a few different behaviors depending on the type of the object.

 - If `c::Distribution` then this will store the distributions.
 - If `c::Tuple{AbstractHypercubeTransform}` then this will store the tuple

# Examples
```julia
ascube(Normal())
ascube(MultivariateNormal())
ascube((Normal(), Normal(2.0)))
ascube( (α = Uniform(), β = Normal()) )
```
"""
function ascube end

"""
    `$(FUNCTIONNAME)(c::AbstractHypercubeTransform, p)`
Transforms from the hypercube with coordinates `p`, to the parameter space
defined by the transformation `c`.

The behavior of this function depends on the nature of `c`.
 - If `c` is a <: Distributions.Distributions and has a quantile method
this will just call the quantile function. If no quantile function is defined
then a custom transformation depending on the type of `c` will be called. If
no custom transformation exists then an error will be raised.
 - If `c` is a Tuple of transformations then transform will iterate through the
 tuple using a similar method to the
 [TransformVariables.jl](https://github.com/tpapp/TransformVariables.jl) method.
"""
function transform end

struct ScalarHC{D} <: AbstractHypercubeTransform
    dist::D
    ScalarHC(d) = new{typeof(d)}(d)
end

# Trait to decide whether a quantile function is defined. This should be compile time right?
struct HasQuant end
struct NoQuant end
has_quantile(::ScalarHC{D}) where {D} = static_hasmethod(quantile, Tuple{D,Float64}) ? HasQuant() : NoQuant()

dimension(::ScalarHC) = 1
dist(d::ScalarHC) = d.dist
ascube(d::Dists.UnivariateDistribution) = ScalarHC(d)

"""
    $(SIGNATURES)
Computes the transformation from the unit hypercube to the distribution space.
"""
function transform(c::AbstractHypercubeTransform, x)
    return _transform(has_quantile(c), c, x)
end

function _transform(::HasQuant, c::AbstractHypercubeTransform, x)
    return quantile(dist(c), x)
end

function _transform(::NoQuant, c::AbstractHypercubeTransform, x)
    throw("No quantile for distribution $(c.dist), implement transform manually")
end

function _step_transform(c::ScalarHC, x::AbstractVector, index)
    transform(c, x[index]), index+1
end

abstract type VectorHC <: AbstractHypercubeTransform end

function transform(c::VectorHC, x::AbstractVector)
    return first(_step_transform(c, x, firstindex(x)))
end

struct ArrayHC{T,M} <: VectorHC
    dist::T
    dims::NTuple{M,Int}
end

dimension(c::ArrayHC) = prod(c.dims)

function ArrayHC(d)
    return ArrayHC{typeof(d),length(size(d))}(d, size(d))
end

ascube(d::Union{Dists.MultivariateDistribution,Dists.Matrixvariate}) = ArrayHC(d)
dist(d::ArrayHC) = d.dist

function _step_transform(h::ArrayHC{<:Dists.Product, M}, p::AbstractVector, index) where {M}
    out = Vector{eltype(p)}(undef, dimension(h))
    for i in 1:dimension(h)
        out[i] = first(_step_transform(ascube(dist(h).v[i]), p, index-1+i))
    end
    return out, index+dimension(h)
end

dimension(c::ArrayHC{<:Dists.Dirichlet,M}) where {M} = prod(c.dims)-1
function _step_transform(h::ArrayHC{<:Dists.Dirichlet, M}, p::AbstractVector, index) where {M}
    d = dist(h)
    α = d.alpha
    T = promote_type(eltype(p), eltype(α))
    #println(T)
    out  = zeros(T, dimension(h)+1)
    #println(eltype(out))
    dstart = Dists.Beta(T(α[1]), sum(@view(α[2:end])))
    #println(typeof(dstart))
    out[1] = quantile(Dists.Beta(T(α[1]), T(sum(@view(α[2:end])))), p[index])
    for i in 2:dimension(h)
        ϕ = quantile(Dists.Beta(T(α[i]),T(sum(@view(α[i+1:end])))), p[index-1+i])
        out[i] = T((1-sum(out))*ϕ)
    end
    out[end] = 1-sum(out)
    return out, index+dimension(h)
end

ascube(d::MT.ProductMeasure) = ArrayHC(d)

function _step_transform(h::ArrayHC{<:MT.ProductMeasure,M}, p::AbstractVector, index) where {M}
    out = Vector{eltype(p)}(undef, dimension(h))
    m = MT.marginals(dist(h))
    for (i,mi) in enumerate(m)
        out[i] = first(_step_transform(ascube(mi), p, index-1+i))
    end
    return out, index+dimension(h)
end


text:
using GZip

# get the format of a file
function getformat(filename::AbstractString)
	format_extensions = Dict(
			"fastq"=>["fastq", "fq", "fastq.gz", "fq.gz", "fastq.zip", "fq.zip"],
			"fasta"=>["fasta", "fa", "fasta.gz", "fa.gz", "fasta.zip", "fa.zip"],
			"sam"=>["sam", "sam.gz", "sam.zip"],
			"bam"=>["bam"],
			"bed"=>["bed", "bed.gz", "bed.zip"]
		)
	for (format, extensions) in format_extensions
		for ext in extensions
			if endswith(filename, "." * ext)
				return format
			end
		end
	end
	#no format match, return unknown
	return "unknown"
end

function iszipped(filename::AbstractString)
	return endswith(filename, ".gz") || endswith(filename, ".zip")
end

# detect the stream type by its filename and mode
function streamtype(filename::AbstractString, mode::AbstractString)
	if iszipped(filename)
		return GZipStream
	else
		return IOStream
	end
end

function verify_openmode(mode::AbstractString)
	if contains(mode, "r") && contains(mode, "w")
		error("Open mode cannot be both read(r) and write(w)")
		return false
	elseif !contains(mode, "r") && !contains(mode, "w")
		error("No write(w) or read(r) mode specified")
		return false
	else
		return true
	end
end

"the common open function to be used by other open files like fastq_open, bed_open"
function opengene_open(filename::AbstractString, mode::AbstractString="r")
	if !verify_openmode(mode)
		return false
	end
	# WAR to fix gz file reading issue of Libz
	# https://github.com/BioJulia/Libz.jl/issues/9
	if GZipStream == streamtype(filename, mode)
        return GZip.open(filename, mode)
	elseif contains(mode, "w")
        return open(filename, mode)
    else
        return open(filename, mode) |> streamtype(filename, mode)
	end
end

# write a dataframe without quote
function write_dataframe(stream, data::DataFrame)
    # we can use printtable function from DataFrames package if no quotemark printing is supported
    # like: printtable(stream,data,header=false, separator='\t', quotemark='\0')
    for i in 1:nrow(data)
        for j in 1:ncol(data)
            print(stream, data[i, j])
            if j<ncol(data)
                print(stream, "\t")
            else
                print(stream, "\n")
            end
        end
    end
end


text:
### A Pluto.jl notebook ###
# v0.12.4

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ d2bf60aa-1a49-11eb-3706-5f7507fc511c
begin
	using Images
	using PlutoUI
	using ImageFiltering
	using Statistics
	using LinearAlgebra
end

# ╔═╡ b2a8e540-19af-11eb-0214-9d5f04ac2486
md"""
# Seam Carving

1. We use convolution with Sobel filters for "edge detection".
2. We use that to write an algorithm that removes "uninteresting" bits of an image in order to shrink it.
"""

# ╔═╡ 455ad8f4-1a4f-11eb-38aa-b9da92a25c17
md"""
## Utility functions
"""

# ╔═╡ 988f2d30-1a50-11eb-0a3d-a15827ed7248
function hbox(x, y, gap=16; sy=size(y), sx=size(x))
	w,h = (max(sx[1], sy[1]),
		   gap + sx[2] + sy[2])

	slate = fill(RGB(1,1,1), w,h)
	slate[1:size(x,1), 1:size(x,2)] .= RGB.(x)
	slate[1:size(y,1), size(x,2) + gap .+ (1:size(y,2))] .= RGB.(y)
	slate
end

# ╔═╡ 14fc8df0-1a64-11eb-13af-cdd3cfc72f2e
vbox(x,y, gap=16) = hbox(x', y')'

# ╔═╡ d5fa6566-1a4f-11eb-0528-a35576836b80
function shrink_image(image, ratio=5)
	height, width = size(image)
	new_height = height ÷ ratio - 1
	new_width = width ÷ ratio - 1
	list = [
		mean(image[
				ratio * i:ratio * (i + 1),
				ratio * j:ratio * (j + 1)
		])
		for i in 1:new_height, j in 1:new_width
	]
	#reshape(list, new_height, new_width)
end

# ╔═╡ ba4af658-1a4f-11eb-217a-7f9165d30dbd
function convolve(M, kernel)
	height, width = size(kernel)

	half_height = height ÷ 2
	half_width = width ÷ 2

	new_image = similar(M)

	# (i, j) loop over the original image
	m, n = size(M)
	@inbounds for i in 1:m
		for j in 1:n
			# (k, l) loop over the neighbouring pixels
			accumulator = 0 * M[1, 1]
			for k in -half_height:-half_height + height - 1
				for l in -half_width:-half_width + width - 1
					Mi = i - k
					Mj = j - l
					# First index into M
					if Mi < 1
						Mi = 1
					elseif Mi > m
						Mi = m
					end
					# Second index into M
					if Mj < 1
						Mj = 1
					elseif Mj > n
						Mj = n
					end

					accumulator += kernel[k, l] * M[Mi, Mj]
				end
			end
			new_image[i, j] = accumulator
		end
	end

	return new_image
end

# ╔═╡ 571690f6-1a4f-11eb-2cb7-8d75d2dfc5eb
function show_colored_array(array)
	pos_color = RGB(0.36, 0.82, 0.8)
	neg_color = RGB(0.99, 0.18, 0.13)
	to_rgb(x) = max(x, 0) * pos_color + max(-x, 0) * neg_color
	to_rgb.(array) / maximum(abs.(array))
end

# ╔═╡ 53bf61d0-19b4-11eb-2212-d91065dae28f
@bind image_url Select([
"https://cdn.shortpixel.ai/spai/w_1086+q_lossy+ret_img+to_webp/https://wisetoast.com/wp-content/uploads/2015/10/The-Persistence-of-Memory-salvador-deli-painting.jpg",
"https://upload.wikimedia.org/wikipedia/commons/thumb/1/17/Gustave_Caillebotte_-_Paris_Street%3B_Rainy_Day_-_Google_Art_Project.jpg/1014px-Gustave_Caillebotte_-_Paris_Street%3B_Rainy_Day_-_Google_Art_Project.jpg",
"https://web.mit.edu/facilities/photos/construction/Projects/stata/1_large.jpg",
])

# ╔═╡ 53d99a46-19b4-11eb-21fa-859486e6289b
img = load(download(image_url))

# ╔═╡ 53f647ea-19b4-11eb-2356-755781371a3a
# arbitrarily choose the brightness of a pixel as a mean of rgb
brightness(c::AbstractRGB) = 0.5 * c.r + 0.59 * c.g + 0.11 * c.b

# ╔═╡ 65932dac-1a67-11eb-2b1d-a74cd795c111
function edgeness(img)
	Sy, Sx = Kernel.sobel()
	b = brightness.(img)
	∇y = convolve(b, Sy)
	∇x = convolve(b, Sx)
	sqrt.(∇x.^2 + ∇y.^2)
end

# ╔═╡ 5413928c-19b4-11eb-1890-230a8077f4ae
Gray.(img)

# ╔═╡ 542f1244-19b4-11eb-207c-131b5699ad24
Gray.(brightness.(img))

# ╔═╡ ec7a7c5e-1a52-11eb-0b3f-49a1c0c2ab4b
shrink_image(img)

# ╔═╡ 2b5e3d8c-1a64-11eb-3a64-cb4e1c314db2
show_colored_array(brightness.(img))

# ╔═╡ ad7c6f50-1a4b-11eb-35fe-a7331d599b45
md"""
## Edge detection filter

We use the Sobel edge detection filter.

```math
\begin{align}

G_x &= \begin{bmatrix}
1 & 0 & -1 \\
2 & 0 & -2 \\
1 & 0 & -1 \\
\end{bmatrix}*A\\
G_y &= \begin{bmatrix}
1 & 2 & 1 \\
0 & 0 & 0 \\
-1 & -2 & -1 \\
\end{bmatrix}*A
\end{align}
```

Where: $A$ is the target image. We can think of these as derivates in the $x$ and $y$ directions.

Then we combine them by finding the magnitude of the **gradient** (in the sense of the multivatiate calculus) be defining.

$$G_\text{total} = \sqrt{G_x^2 + G_y^2}.$$
"""

# ╔═╡ ad6289f0-1a4b-11eb-1c1e-b1158e269574
Sy, Sx = Kernel.sobel();

# ╔═╡ ad4684c6-1a4b-11eb-3c63-5f888a1d834b
[show_colored_array(Sx) show_colored_array(Sy)]

# ╔═╡ ad2c6d34-1a4b-11eb-3f16-7332941d79fa
begin
	img_brightness = brightness.(img)
	∇x = convolve(img_brightness, Sx)
	∇y = convolve(img_brightness, Sy)
	hbox(show_colored_array(∇x), show_colored_array(∇y))
end

# ╔═╡ ad0d9044-1a4b-11eb-1a23-b7160b7793ed
vbox(
	hbox(img[200:end, 50:250], img[200:end, 50:250]),
	hbox(show_colored_array.((∇x[200:end, 50:250], ∇y[200:end, 50:250]))...)
)

# ╔═╡ ac812a78-1a4b-11eb-216d-91a51c6b4742
begin
	edged = edgeness(img)
	hbox(img, Gray.(edged) / maximum(abs.(edged)))
end

# ╔═╡ ac652cc4-1a4b-11eb-2fee-dd2887e80796
md"""
## Seam carving idea

The idea of seam carving is to find the path from the top of the image to the bottom of the imaga where the path minimized the edgness.

In other words, this path minimized the number of edges it crosses.

"""

# ╔═╡ ac48d830-1a4b-11eb-3296-4d12c4153138

# ╔═╡ ac28a5f6-1a4b-11eb-16e8-e18a949846cb

# ╔═╡ ac002a2c-1a4b-11eb-3c4e-cb541f081ff1

# ╔═╡ abe1468e-1a4b-11eb-00e2-a394e9ae1504

# ╔═╡ abbef57a-1a4b-11eb-2c48-b12aa928410f

# ╔═╡ Cell order:
# ╠═d2bf60aa-1a49-11eb-3706-5f7507fc511c
# ╟─b2a8e540-19af-11eb-0214-9d5f04ac2486
# ╟─455ad8f4-1a4f-11eb-38aa-b9da92a25c17
# ╟─988f2d30-1a50-11eb-0a3d-a15827ed7248
# ╟─14fc8df0-1a64-11eb-13af-cdd3cfc72f2e
# ╟─d5fa6566-1a4f-11eb-0528-a35576836b80
# ╟─ba4af658-1a4f-11eb-217a-7f9165d30dbd
# ╟─571690f6-1a4f-11eb-2cb7-8d75d2dfc5eb
# ╟─65932dac-1a67-11eb-2b1d-a74cd795c111
# ╟─53bf61d0-19b4-11eb-2212-d91065dae28f
# ╠═53d99a46-19b4-11eb-21fa-859486e6289b
# ╠═53f647ea-19b4-11eb-2356-755781371a3a
# ╠═5413928c-19b4-11eb-1890-230a8077f4ae
# ╠═542f1244-19b4-11eb-207c-131b5699ad24
# ╠═ec7a7c5e-1a52-11eb-0b3f-49a1c0c2ab4b
# ╠═2b5e3d8c-1a64-11eb-3a64-cb4e1c314db2
# ╟─ad7c6f50-1a4b-11eb-35fe-a7331d599b45
# ╠═ad6289f0-1a4b-11eb-1c1e-b1158e269574
# ╠═ad4684c6-1a4b-11eb-3c63-5f888a1d834b
# ╟─ad2c6d34-1a4b-11eb-3f16-7332941d79fa
# ╟─ad0d9044-1a4b-11eb-1a23-b7160b7793ed
# ╟─ac812a78-1a4b-11eb-216d-91a51c6b4742
# ╠═ac652cc4-1a4b-11eb-2fee-dd2887e80796
# ╠═ac48d830-1a4b-11eb-3296-4d12c4153138
# ╠═ac28a5f6-1a4b-11eb-16e8-e18a949846cb
# ╠═ac002a2c-1a4b-11eb-3c4e-cb541f081ff1
# ╠═abe1468e-1a4b-11eb-00e2-a394e9ae1504
# ╠═abbef57a-1a4b-11eb-2c48-b12aa928410f

max_stars_repo_name: martin-fabbri/colab-notebooks

max_issues_repo_name: martin-fabbri/colab-notebooks

max_forks_repo_name: martin-fabbri/colab-notebooks


text:
function cable_transform(y, z)
    v1 = [0.0, 0.0, 1.0]
    # if norm(z) > norm(y)
    #     v2 = y[1:3,1] - z[1:3,1]
    # else
    #     v2 = z[1:3,1] - y[1:3,1]
    # end
    v2 = y[1:3,1] - z[1:3,1]
    normalize!(v2)
    ax = cross(v1, v2)
    ang = acos(v1'*v2)
    R = AngleAxis(ang, ax...)

    if any(isnan.(R))
        R = I
    else
        nothing
    end

    compose(Translation(z), LinearMap(R))
end

function cast3d(v::AbstractVector{T}) where {T}
	l = length(v)
	if l == 3
		return v
	elseif l == 2
		return [v[1], 0.0, v[2]]
	else
		error("Wrong vector length")
	end
end

function bivector_rotation(a::AbstractVector{T}, b::AbstractVector{T}) where {T}
	na = norm(a)
	nb = norm(b)
	v = cross(a,b)/(na*nb)
	s = norm(v)
	c = (a'*b)/(na*nb)
	vx = skew(v)
	α = (1 - c)/s^2
	R = I + vx + vx*vx*α
	return R
end

function default_background!(vis; grid::Bool=false, axes::Bool=false, zoom=1.0)
    setvisible!(vis["/Background"], true)
    setprop!(vis["/Background"], "top_color", RGBA(1.0, 1.0, 1.0, 1.0))
    setprop!(vis["/Background"], "bottom_color", RGBA(1.0, 1.0, 1.0, 1.0))
	setvisible!(vis["/Axes"], axes)
    # setvisible!(vis["/Grid"], grid)
	# setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", zoom)
end

function get_line_material(size::Real; α=1.0)
    orange_col = [1,153/255,51/255]
    blue_col = [51/255,1,1]
    black_col = [0,0,0]
    orange_mat = LineBasicMaterial(color=color=RGBA(orange_col...,α), linewidth=size)
    blue_mat = LineBasicMaterial(color=color=RGBA(blue_col...,α), linewidth=size)
    black_mat = LineBasicMaterial(color=color=RGBA(black_col...,α), linewidth=size)
    return orange_mat, blue_mat, black_mat
end

function get_material(;α=1.0)
    orange_col = [1,153/255,51/255]
    blue_col = [51/255,1,1]
    black_col = [0,0,0]
    orange_mat = MeshPhongMaterial(color=RGBA(orange_col...,α))
    blue_mat = MeshPhongMaterial(color=RGBA(blue_col...,α))
    black_mat = MeshPhongMaterial(color=RGBA(black_col...,α))
    return orange_mat, blue_mat, black_mat
end


text:
using ArgMinimist
using Base.Test

function checkPositional()
  pos, kw = parse_args(["hello","world"])
  return pos == ["hello","world"] && length(kw)==0
end

function checkKeyword()

  expKW = Dict{Symbol,Any}()
  expKW[:jam] = false
  expKW[:honey] = true
  expKW[:cafe] = true

  expKW[:a] = true
  expKW[:b] = true
  expKW[:c] = true

  pos, kw = parse_args(["hello","world",
                          "--no-jam",
                          "--with-honey",
                          "-abc",
                          "--cafe"])

  return pos == ["hello","world"] && kw == expKW
end

function checkNamed()

  expKW = Dict{Symbol,Any}()
  expKW[:jam] = "yes"
  expKW[:honey] = "no"
  expKW[:cafe] = "ok"

  args = ["yes","no","--cafe=ok"]

  pos, kw = parse_named_args(["jam","honey","cafe"], args)
  return length(pos)==0 && kw == expKW
end
# write your own tests here
@test checkPositional()
@test checkKeyword()
@test checkNamed()


text:
"""
```
compute_Ψ(Γ₁, Γ₂, Γ₃, Γ₄, Γ₅, Γ₆, JV = []; schur_fnct = schur!)
compute_Ψ(m::RALLinearizedSystem; zero_entropy_jacobian = false, schur_fnct = schur!)
compute_Ψ(m::RiskAdjustedLinearization; zero_entropy_jacobian = false, schur_fnct = schur!)
```

solves via QZ decomposition for ``\\Psi_n`` in the quadratic matrix equation

``math
\\begin{aligned}
0 =  JV + \\Gamma_3 + \\Gamma_4\\Psi_n + (\\Gamma_5 + \\Gamma_6\\Psi_n) (\\Gamma_2 + \\Gamma_1 \\Psi_n).
\\end{aligned}
``

See the documentation of `RiskAdjustedLinearizations.jl` for details about what these matrices are.

### Inputs
For the first method, all the required inputs must have type `AbstractMatrix{<: Number}`. The `JV` term is empty by default,
in which case `qzdecomp` assumes that `JV` is the zero matrix, which corresponds to the
case of the deterministic steady state.

The second and third methods are wrappers for the first method.

Internally, there are in-place versions (`compute_Ψ!`) to avoid allocations when
`compute_Ψ!` is called repeatedly by one of the numerical algorithms (e.g. `relaxation!`).

### Keywords
- `schur_fnct::Function`: specifies which Generalized Schur algorithm is desired. By default,
    the implementation from BLAS is called, but the user may want to use the Generalized Schur algorithm
    from packages like `GenericLinearAlgebra.jl` and `GenericSchur.jl`.
- `zero_entropy_jacobian::Bool`: if true, then we assume the Jacobian of the entropy is all zeros (i.e. in the deterministic steady state).
    This keyword allows the user to avoid adding zeros unnecessarily and/or let the cached entropy Jacobian remain undefined.
"""
function compute_Ψ(Γ₁::AbstractMatrix{S}, Γ₂::AbstractMatrix{S}, Γ₃::AbstractMatrix{S}, Γ₄::AbstractMatrix{S},
                   Γ₅::AbstractMatrix{S}, Γ₆::AbstractMatrix{S}, JV::AbstractMatrix{S} = Matrix{S}(undef, 0, 0);
                   schur_fnct::Function = schur!) where {S <: Number}

    Nzy = sum(size(Γ₅))
    AA  = Matrix{Complex{S}}(undef, Nzy, Nzy)
    BB  = similar(AA)

    return compute_Ψ!(AA, BB, Γ₁, Γ₂, Γ₃, Γ₄, Γ₅, Γ₆, JV; schur_fnct = schur_fnct)
end

function compute_Ψ(m::RALLinearizedSystem; zero_entropy_jacobian::Bool = false, schur_fnct::Function = schur!) where {S <: Number}
    if zero_entropy_jacobian
        return compute_Ψ(m[:Γ₁], m[:Γ₂], m[:Γ₃], m[:Γ₄], m[:Γ₅], m[:Γ₆]; schur_fnct = schur_fnct)
    else
        return compute_Ψ(m[:Γ₁], m[:Γ₂], m[:Γ₃], m[:Γ₄], m[:Γ₅], m[:Γ₆], m[:JV]; schur_fnct = schur_fnct)
    end
end

@inline function compute_Ψ(m::RiskAdjustedLinearization; zero_entropy_jacobian::Bool = false, schur_fnct::Function = schur!) where {S <: Number}
    return compute_Ψ(m.linearization; zero_entropy_jacobian = zero_entropy_jacobian, schur_fnct = schur_fnct)
end

function compute_Ψ!(AA::AbstractMatrix{Complex{S}}, BB::AbstractMatrix{Complex{S}},
                    Γ₁::AbstractMatrix{S}, Γ₂::AbstractMatrix{S}, Γ₃::AbstractMatrix{S}, Γ₄::AbstractMatrix{S},
                    Γ₅::AbstractMatrix{S}, Γ₆::AbstractMatrix{S}, JV::AbstractMatrix{S} = Matrix{S}(undef, 0, 0);
                    schur_fnct::Function = schur!) where {S <: Number}

    if !isempty(JV)
        Γ₃ += JV
    end
    Ny, Nz = size(Γ₅)

    # Populate AA
    AA[1:Ny, 1:Nz]                 = Γ₅
    AA[1:Ny, (Nz + 1):end]         = Γ₆
    AA[(Ny + 1):end, 1:Nz]         = Matrix{Complex{S}}(I, Nz, Nz) # faster but makes allocations, unlike Diagonal(Ones{Complex{S}}(Nz))
    AA[(Ny + 1):end, (Nz + 1):end] = Zeros{Complex{S}}(Nz, Ny)

    # Populate BB
    BB[1:Ny, 1:Nz]                 = -Γ₃
    BB[1:Ny, (Nz + 1):end]         = -Γ₄
    BB[(Ny + 1):end, 1:Nz]         =  Γ₁
    BB[(Ny + 1):end, (Nz + 1):end] =  Γ₂

    # Compute QZ and back out Ψ
    schurfact = schur_fnct(AA, BB)
	ordschur!(schurfact, [abs(αᵢ) >= abs(βᵢ) for (αᵢ, βᵢ) in
                          zip(schurfact.α, schurfact.β)]) # eigenvalues = schurfact.β / schurfact.α

    # Note: We could calculate stability conditions here but, typically,
    # compute_Ψ is called enough within the relaxation algorithm
    # that calculating the stability conditions repeatedly takes more time
    # than just calling eigen once

    return real(schurfact.Z[Nz + 1:end, 1:Nz] / schurfact.Z[1:Nz, 1:Nz])
end

function compute_Ψ!(AA::AbstractMatrix{Complex{S}}, BB::AbstractMatrix{Complex{S}},
                    m::RALLinearizedSystem; zero_entropy_jacobian::Bool = false, schur_fnct::Function = schur!) where {S <: Number}
    if zero_entropy_jacobian
        return compute_Ψ!(AA, BB, m[:Γ₁], m[:Γ₂], m[:Γ₃], m[:Γ₄], m[:Γ₅], m[:Γ₆]; schur_fnct = schur_fnct)
    else
        return compute_Ψ!(AA, BB, m[:Γ₁], m[:Γ₂], m[:Γ₃], m[:Γ₄], m[:Γ₅], m[:Γ₆], m[:JV]; schur_fnct = schur_fnct)
    end
end

@inline function compute_Ψ!(AA::AbstractMatrix{Complex{S}}, BB::AbstractMatrix{Complex{S}},
                            m::RiskAdjustedLinearization; zero_entropy_jacobian::Bool = false, schur_fnct::Function = schur!) where {S <: Number}
    return compute_Ψ!(AA, BB, m.linearization; zero_entropy_jacobian = zero_entropy_jacobian, schur_fnct = schur_fnct)
end


text:
module Misc

using Statistics, Clustering
using JSON
import DataStructures: PriorityQueue, peek
import SpecialFunctions: lbeta, lgamma
import Base: sum, maximum, minimum
import Statistics: mean
import LinearAlgebra: norm

export nanmean, nanmax, nansum, nanmin, sum, mean, maximum, minimum
export logmeanexp, logvarexp, logsumexp, lbeta
export KL, rmse, mae, hoyer
export discretize, categorize, cooccurrences
export finucan, SNR
export save_json, load_json

@inline sum(X::AbstractArray,dims) = dropdims(sum(X;dims=dims); dims=dims)
@inline mean(X::AbstractArray,dims) = dropdims(mean(X;dims=dims); dims=dims)
@inline maximum(X::AbstractArray,dims) = dropdims(maximum(X;dims=dims); dims=dims)
@inline minimum(X::AbstractArray,dims) = dropdims(minimum(X;dims=dims); dims=dims)

@inline _nanmean(X::AbstractArray) = mean(filter(!isnan,X))
nanmean(X::AbstractArray; dims=nothing) = dims == nothing ? _nanmean(X) : mapslices(_nanmean, X; dims=dims)
nanmean(X::AbstractArray, dims) = dropdims(nanmean(X; dims=dims),dims=dims)

@inline _nansum(X::AbstractArray) = sum(filter(!isnan,X))
nansum(X::AbstractArray; dims=nothing) = dims == nothing ? _nansum(X) : mapslices(_nansum, X; dims=dims)
nansum(X::AbstractArray, dims) = dropdims(nansum(X; dims=dims),dims=dims)

@inline _nanmax(X::AbstractArray) = maximum(filter(!isnan,X))
nanmax(X::AbstractArray; dims=nothing) = dims == nothing ? _nanmax(X) : mapslices(_nanmax, X; dims=dims)
nanmax(X::AbstractArray, dims) = dropdims(nanmax(X; dims=dims),dims=dims)

@inline _nanmin(X::AbstractArray) = minimum(filter(!isnan,X))
nanmin(X::AbstractArray; dims=nothing) = dims == nothing ? _nanmin(X) : mapslices(_nanmin, X; dims=dims)
nanmin(X::AbstractArray, dims) = dropdims(nanmin(X; dims=dims),dims=dims)

@inline _logsumexp(X::AbstractArray) = begin x⁺ = maximum(X); log(sum(exp, X .- x⁺)) + x⁺ end
logsumexp(X::AbstractArray; dims=nothing) = dims == nothing ? _logsumexp(X) : mapslices(_logsumexp, X; dims=dims)
logsumexp(X::AbstractArray, dims) = dropdims(logsumexp(X; dims=dims),dims=dims)

@inline _logmeanexp(X::AbstractArray) = _logsumexp(X) - log(length(X))
logmeanexp(X::AbstractArray; dims=nothing) = dims == nothing ? _logmeanexp(X) : mapslices(_logmeanexp, X; dims=dims)
logmeanexp(X::AbstractArray, dims) = dropdims(logmeanexp(X; dims=dims),dims=dims)

@inline _logvarexp(X::AbstractArray; μ=_logmeanexp(X)) = begin EX = _logsumexp(X); log(var(exp.(X .- EX); mean=exp(μ - EX))) + 2.0 * EX end
logvarexp(X::AbstractArray; dims=nothing) = dims == nothing ? _logvarexp(X) : mapslices(_logvarexp, X; dims=dims)
logvarexp(X::AbstractArray, dims) = dropdims(logvarexp(X; dims=dims),dims=dims)

@inline _lbeta(X::AbstractArray) =  sum(lgamma, X) - lgamma(sum(X))
lbeta(X::AbstractArray; dims=nothing) = dims == nothing ? _lbeta(X) : mapslices(_lbeta, X; dims=dims)
lbeta(X::AbstractArray, dims) = dropdims(lbeta(X; dims=dims),dims=dims)

@inline _lbeta(γ::Number,I::Integer) = I*lgamma(γ) - lgamma(I*γ)
lbeta(γ::Number, sz::Vararg{Integer}; dims=nothing) = dims == nothing ? _lbeta(γ,prod(sz)) : fill(_lbeta(γ,prod(sz[dims])), map(t -> t[1] in dims ? 1 : t[2], enumerate(sz))...)

function rmse(X::AbstractArray, Xᵖ::AbstractArray)
    return sqrt(nanmean((X .- Xᵖ) .^ 2 ))
end

function mae(X::AbstractArray, Xᵖ::AbstractArray)
    return nanmean(abs.(X .- Xᵖ))
end

function KL(X::AbstractArray, Xᵖ::AbstractArray)
    KL_div = 0.0
    for (x,xᵖ) in zip(X,Xᵖ)
        if !isnan(x) && !isnan(xᵖ)
            KL_div += (x ≈ 0.0) ? 0.0 : x * (log(x) - log(xᵖ)) + xᵖ - x
        end
    end
    return KL_div
end

function SNR(X::AbstractArray, Xᵖ::AbstractArray; base=10.0)
    return 2.0*base*(log(base,norm(X)) - log(base,norm(X .- Xᵖ)))
end

function hoyer(X::AbstractArray)
    N = length(X)
    return (sqrt(N) - sum(abs,X)/sqrt(sum(X .* X))) / (sqrt(N) - 1)
end

function discretize(X::AbstractArray{ℜ}, L::Ƶ; ϵ=1e-10) where {ℜ <: Real, Ƶ <: Integer}
    Xₙ = X .- minimum(X) .+ ϵ
    Xₙ ./= maximum(Xₙ) + ϵ
    return Ƶ.(floor.(L*Xₙ)) .+ 1
end

function categorize(X::AbstractArray, L::Ƶ=2; cluster::Bool=true) where {Ƶ <: Integer}
    if eltype(X) <: AbstractFloat
    	if cluster
        	Mₓ = kmeans(reshape(X,1,:),L)
        	Categories = Dict(zip(sortperm(Mₓ.centers[:]),1:L))
        	return map(x -> Categories[x], reshape(Mₓ.assignments, size(X)))
        else
        	return discretize(X,L)
        end
    elseif eltype(X) <: Integer
            return X .- minimum(X) .+ 1
    else
        Labels = sort(unique(X))
        Categories = Dict(zip(Labels,1:length(Labels)))
        return map(x -> Categories[x], X)
    end
end

function cooccurrences(Tab::AbstractArray{Ƶ,2},dims) where {Ƶ <: Integer}
    T, N = size(Tab)
    X = zeros(Ƶ,dims...)
    for t in 1:T
        X[Tab[t,:]...] += 1
    end
    return X
end

function cooccurrences(Tab::AbstractArray{Ƶ,2}) where {Ƶ <: Integer}
	I = nanmax(Tab;dims=1)[1,:]
	return cooccurrences(Tab,I)
end

function finucan(N::ℜ, p::Array{ℜ}; ϵ::ℜ=1e-16) where {ℜ<:Real}
    K, k = length(p), 0
    p ./= sum(p) + ϵ
    Sᴿ = (N + K/2.0).*(p .+ ϵ)
    S = floor.(Sᴿ)
    N₀ = sum(S)

    F = Array{ℜ}(undef,K)

    if N₀ < N
        F .= Sᴿ .- S
        Q = PriorityQueue{Int,ℜ}([k => ((1.0 - F[k])./(S[k] + 1.0)) for k=1:K])
        while N₀ < N
            k, Qₖ = peek(Q)
            S[k] += 1.0
            F[k] -= 1.0
            Q[k] = (1 - F[k])/(S[k] + 1.0)
            N₀ += 1.0
        end
    elseif N₀ > N
        F = Sᴿ .- S
        Q = PriorityQueue{Int,ℜ}([k => (F[k]/(S[k] + ϵ)) for k=1:K])
        while N₀ > N
            k, Qₖ = peek(Q)
            S[k] -= 1.0
            F[k] += 1.0
            Q[k] = F[k]/(S[k] + ϵ)
            N₀ -= 1.0
        end
    end
    return S
end

function save_json(filename; variables...)
    json_str = JSON.json(Dict(variables));
    open(filename, "w") do f
        write(f, json_str)
    end
    return json_str
end

function load_json(filename)
    json_str = nothing
    open(filename, "r") do f
        json_str = read(f, String)
    end
    return JSON.parse(json_str)
end

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
struct IslandIteratorState
    pos
end

struct IslandIterator
    read_start_pos
    read_length
end

function Base.iterate(iter::IslandIterator, state = nothing)

    pos = isnothing(state) ? 1 : state.pos
    no_reads = 1

    pos == length(iter.read_start_pos) && return ((1,iter.read_length), IslandIteratorState(pos+1))
    pos>length(iter.read_start_pos) && return nothing

    while iter.read_start_pos[pos]+iter.read_length>=iter.read_start_pos[pos+1]
        pos += 1
        no_reads+=1
        pos == length(iter.read_start_pos) && return  ((no_reads,iter.read_length+iter.read_start_pos[pos]-iter.read_start_pos[pos-no_reads+1]), IslandIteratorState(pos + 1))
    end

    return ((no_reads,iter.read_length+iter.read_start_pos[pos]-iter.read_start_pos[pos-no_reads+1]), IslandIteratorState(pos + 1))
end

Base.IteratorSize(::Type{IslandIterator}) = Base.SizeUnknown()


text:
using Flux, Flux.Tracker
using NeuralProcesses, NeuralProcesses.Experiment
using Distributions
using Plots
using Random
using PyCall

@pyimport pickle

function save(path, obj)
    out = open(path, "w")
    pickle.dump(obj, out)
    close(out)
end

convcnp = NeuralProcesses.untrack(best_model("models/convcnp/loglik/sawtooth.bson"))
convnp = NeuralProcesses.untrack(best_model("models/convnp/loglik/sawtooth.bson"))

p = Sawtooth()

xt = collect(range(-2, 2, length=600))
yt = reshape(rand(p(xt, 1e-10)), :)
inds = sort(randperm(length(xt))[1:7])
xc = xt[inds]
yc = yt[inds]

save("sawtooth/data.pickle", Dict(
    "xc" => xc,
    "yc" => yc,
    "xt" => xt,
    "yt" => yt
))

xc = reshape(xc, :, 1, 1)
yc = reshape(yc, :, 1, 1)
xt = reshape(xt, :, 1, 1)
yt = reshape(yt, :, 1, 1)

# Embedding:
xz = convcnp.encoder.disc(xc, xt, margin=16f0)
xz, z = code(convcnp.encoder.coder[1], xc, yc, xz)

inds = reshape((xz .>= -2) .& (xz .<= 2), :)
plot(xz[inds, 1, 1], z[inds, 1, 1])
save("sawtooth/embedding.pickle", Dict(
    "xz" => xz[inds, 1, 1],
    "z1" => z[inds, 1, 1],
    "z2" => z[inds, 2, 1]
))

xz, z = code(convcnp.encoder.coder[2], xz, z, xz)
z = NeuralProcesses.sample(z, num_samples=1)

# CNN:
xz, z = code(convcnp.decoder[1], xz, z, xt)
save("sawtooth/embedding_cnn.pickle", Dict(
    "xz" => xz[inds, 1, 1],
    "z1" => z[inds, 1, 1],
    "z2" => z[inds, 2, 1]
))

xz, z = code(convcnp.decoder[2], xz, z, xt)
xz, z = code(convcnp.decoder[3], xz, z, xt)

xc = reshape(xc, :)
yc = reshape(yc, :)
xt = reshape(xt, :)
yt = reshape(yt, :)

μ = reshape(mean(z), :)
σ = reshape(std(z), :)

# Prediction:
save("sawtooth/pred.pickle", Dict(
    "xt" => xt,
    "mu" => μ,
    "std" => σ
))

# Also predict with ConvNP:
_, _, _, samples = predict(convnp, xc, yc, xt; margin=16f0)
save("sawtooth/pred_convnp.pickle", Dict(
    "xt" => xt,
    "samples" => samples
))


text:
using Random
import ParticleInCell.F90

@testset "push" begin

    dimx, dimy = 1.0, 1.0
    nx, ny = 100, 100
    mesh = TwoDGrid(dimx, nx, dimy, ny)
    e = 1.0
    b = 2.0
    nbpart = 1
    pj = zeros(7, 1)
    dt = 0.1
    x0, y0 = 0.1, 0.1

    pj = ParticleGroup{2,2}(nbpart, charge = 1.0, mass = 1.0, n_weights = 1)
    pf = ParticleGroup{2,2}(nbpart, charge = 1.0, mass = 1.0, n_weights = 1)

    pj.array[1, 1] = x0
    pj.array[2, 1] = y0
    pj.array[5, 1] = 1.0

    ex = ones(nx+1, ny+1)
    ey = zeros(nx+1, ny+1)
    bz = ones(nx+1, ny+1)

    pf.array .= copy(pj.array)

    time = 0

    kernel = CloudInCell()

    for istep = 1:50

        push_v!(pj, kernel, mesh, ex, ey, bz, dt)
        push_x!(pj, mesh, dt)

        F90.push_v!(pf, kernel, mesh, ex, ey, bz, dt)
        F90.push_x!(pf, mesh, dt)

        @test all(pj.array .== pf.array)
    end

end


text:
import SpecialFunctions: erfc
using CLArrays

function blackscholes(sptprice, strike, rate, volatility, time)
    logterm = log( sptprice / strike)
    powterm = .5f0 * volatility * volatility
    den = volatility * sqrt(time)
    d1 = (((rate + powterm) * time) + logterm) / den
    d2 = d1 - den
    NofXd1 = cndf2(d1)
    NofXd2 = cndf2(d2)
    futureValue = strike * exp(- rate * time)
    c1 = futureValue * NofXd2
    call_ = sptprice * NofXd1 - c1
    put  = call_ - futureValue + sptprice
    return put
end

function cndf2(x)
    0.5f0 + 0.5f0 * erfc(0.707106781f0 * x)
end

using BenchmarkTools

T = Float32
N = 10^7
sptprice   = T[42.0 for i = 1:N]
initStrike = T[40.0 + (i / N) for i = 1:N]
rate       = T[0.5 for i = 1:N]
volatility = T[0.2 for i = 1:N]
time       = T[0.5 for i = 1:N]
result     = similar(time)
comparison = blackscholes.(sptprice, initStrike, rate, volatility, time)
_sptprice = CLArray(sptprice)
_initStrike = CLArray(initStrike)
_rate = CLArray(rate)
_volatility = CLArray(volatility)
_time = CLArray(time)
_result = CLArray(result)

@time begin
    _result .= blackscholes.(_sptprice, _initStrike, _rate, _volatility, _time)
    (GPUArrays.synchronize)(_result)
end