text:
using StatsBase: autocov

"""
    check_cache_type(terminator, prob, sol)
"""
check_cache_type(::T, ::P, ::S) where {T, P, S} =
    error("Terminator type is incompatible with the ",
          "problem and solution types.\n",
          "terminator: $T\n",
          "problem: $P\n",
          "solution: $S")

should_terminate!(::T, ::S) where {T, S} =
    error("`should_terminate!` not defined for $T and $S")

struct NullTerminator <: Terminator end
check_cache_type(::NullTerminator, ::Any, ::LESolution) = true
should_terminate!(::NullTerminator, ::LESolution) = false

@enum ConvErrorType UnstableConvError=1 StableConvError=2

abstract type ConvDetail end
struct UnstableConvDetail <: ConvDetail
    var::Float64
    tail_cov::Float64
    tail_ok::Bool
end

struct ConvergenceHistory
    orth::Vector{Int}
    kinds::Vector{ConvErrorType}
    errors::Vector{Vector{Float64}}
    thresholds::Vector{Vector{Float64}}
    details::Vector{Vector{ConvDetail}}
end

ConvergenceHistory(dim_lyap::Int) =
    ConvergenceHistory(
        Vector{Int}(),
        Vector{ConvErrorType}(),
        [Vector{Float64}() for _ in 1:dim_lyap],
        [Vector{Float64}() for _ in 1:dim_lyap],
        [Vector{ConvDetail}() for _ in 1:dim_lyap],
    )

function new_error!(history::ConvergenceHistory, n::Int, kind::ConvErrorType)
    push!(history.orth, n)
    push!(history.kinds, kind)
end

function record_error!(history::ConvergenceHistory,
                       ::Type{Val{UnstableConvError}},
                       i::Int, err, th, args...)
    push!(history.errors[i], err)
    push!(history.thresholds[i], th)
    push!(history.details[i], UnstableConvDetail(args...))
end

function record_error!(history::ConvergenceHistory,
                       i::Int, err, th, detail)
    push!(history.errors[i], err)
    push!(history.thresholds[i], th)
    push!(history.details[i], detail)
end

function assert_consistent(history::ConvergenceHistory)
    @assert [length(history.orth)] ==
        [length(history.kinds)] ==
        unique(length.(history.errors)) ==
        unique(length.(history.thresholds)) ==
        unique(length.(history.details))
end

absolute_error(history::ConvergenceHistory, i::Int) = history.errors[i]

relative_error(history::ConvergenceHistory, i::Int) =
    history.errors[i] ./ history.thresholds[i]

mutable struct AutoCovTerminator <: Terminator
    rtol::Float64
    atol::Float64
    chunks::Int
    max_mle::Float64
    max_tail_corr::Float64
    tail_ratio::Float64
    max_fp_cv::Float64
    min_periodic_corr::Float64
    skip_ratio::Float64
    next_check::Int
end
# TODO: separate state (next_check) and settings.
#
# TODO: Separate terminator for unstable (chaotic) and stable cases
# and add a terminator combinator for composing them.  The stable case
# can be decomposed into more sub-components (see stable_le_error).
#
# TODO: rtol and atol are used for unstable (chaotic) and stable cases
# but they have different semantics in each case.  Find a way to
# translate one to the other and/or use different option.  Using
# different option would be trivial once the combinator is
# implemented.

function AutoCovTerminator(;
        rtol = 0.1,
        atol = 0.1,
        chunks = 2,
        max_mle = 0.01,
        max_tail_corr = 0.05,
        tail_ratio = 0.1,
        skip_ratio = 0.1,
        max_fp_cv = 0.001,
        min_periodic_corr = 0.8,
        first_check = 400,
        )
    @assert chunks > 1
    AutoCovTerminator(rtol, atol,
                      chunks, max_mle,
                      max_tail_corr,
                      tail_ratio,
                      skip_ratio,
                      max_fp_cv,
                      min_periodic_corr,
                      first_check)
end

check_cache_type(::AutoCovTerminator, ::Any, ::LESolution{true, true}) = true
# Need RecOS=true and RecFTLE=true.

function should_terminate!(tmnr::AutoCovTerminator, sol::LESolRecFTLE)
    atol = tmnr.atol
    rtol = tmnr.rtol
    n = sol.num_orth
    if n < tmnr.next_check
        return false
    end

    les = lyapunov_exponents(sol)
    if les[1] < tmnr.max_mle
        return should_terminate_stable!(tmnr, sol)
    end

    ftle = @view sol.ftle_history[1:n]
    dim_lyap = size(ftle[1], 1)

    history = sol.convergence
    new_error!(history, n, UnstableConvError)
    rerr = 1.0
    ok = true
    for i in 1:dim_lyap
        xs = [x[i] for x in ftle]
        m = mean(xs)
        err, c = correlated_sem_with_cov(xs)
        th = atol + abs(m) * rtol
        rerr = max(rerr, err / th)
        tail_width = ceil(Int, length(c) * tmnr.tail_ratio)
        tail_cov = mean(abs, @view c[end-tail_width:end])
        tail_ok = tail_cov < abs(c[1]) * tmnr.max_tail_corr
        ok = ok && err < th
        ok = ok && tail_ok
        record_error!(history, Val{UnstableConvError}, i, err, th,
                      c[1], tail_cov, tail_ok)
    end
    assert_consistent(history)

    if !ok
        # aiming for err < th/2 (so that's why "4 *")
        tmnr.next_check = ceil(Int, 4 * rerr^2 * n)
        tmnr.next_check = min(tmnr.next_check, length(sol.ftle_history))
    end
    sol.converged = ok
    return ok
end

function should_terminate_stable!(tmnr::AutoCovTerminator, sol::LESolRecFTLE)
    atol = tmnr.atol
    rtol = tmnr.rtol
    n = sol.num_orth
    skip = max(1, floor(Int, n * tmnr.skip_ratio))
    dim_lyap = size(sol.ftle_history[1], 1)

    history = sol.convergence
    new_error!(history, n, StableConvError)
    rerr = 1.0
    ok = true
    for i in 1:dim_lyap
        ftle = @view ftle_history(sol, i)[skip:end]
        m = mean(ftle)
        th = atol + abs(m) * rtol
        err, detail = stable_le_error(ftle)
        if isfinite(err)
            rerr = max(rerr, err / th)
        end
        ok = ok && err < th
        record_error!(history, i, err, th, detail)
    end
    assert_consistent(history)

    if !ok
        # aiming for err < th/2 (so that's why "2 *")
        tmnr.next_check = ceil(Int, 2 * rerr * n)
        tmnr.next_check = min(tmnr.next_check, length(sol.ftle_history))
    end
    sol.converged = ok
    return ok
end

"""
    correlated_sem_with_cov(xs [, cutoff ≈ √length(xs)])

Standard error of mean of (possibly) correlated data.
https://stats.stackexchange.com/a/274672
"""
function correlated_sem_with_cov(xs::AbstractVector,
                                 cutoff = ceil(Int, sqrt(length(xs))))
    c = autocov(xs, 0:cutoff)
    σ = c[1]
    g = @view c[2:end]
    n = length(xs)
    k = 1:cutoff
    sem² = (σ + 2sum((n - k) ./ n .* g)) / n
    sem = sqrt(max(σ / n, sem²))    # TODO: should I use max?
    return (sem, c)
end

struct FixedPointConvDetail <: ConvDetail
    var::Float64
end

struct NonNegativeAutoCovConvDetail <: ConvDetail
    var::Float64
    min_corr::Float64
end

struct NonPeriodicConvDetail <: ConvDetail
    var::Float64
    max_corr::Float64
end

struct PeriodicConvDetail <: ConvDetail
    period::Int
end

stable_le_error(ftle, opts) =
    stable_le_error(ftle;
                    max_fp_cv = opts.max_fp_cv,
                    min_periodic_corr = opts.min_periodic_corr,
                    )

function stable_le_error(ftle::AbstractVector;
                         max_fp_cv = 0.001,
                         min_periodic_corr = 0.8,
                         cutoff = length(ftle) ÷ 2)

    v = var(ftle, corrected=false)  # = autocov(ftle, 0:0)[1]
    m = mean(ftle)

    # Handle non-fluctuating LE first.  Likely a fixed point.
    if sqrt(v) / abs(m) < max_fp_cv
        # TODO: detect fixed-point directly in the phase space.
        return 0.0, FixedPointConvDetail(v)
    end

    # Detect "period" in ftle.  Here, "period" is refered to the index
    # with maximum covariance *after* it crosses the negative value.
    g = autocov(ftle, 1:cutoff)
    i = findfirst(x -> x < 0, g)
    if i == 0
        # Dynamics assumed to be periodic but autocovariance does not
        # cross 0.  Probably length(ftle) is too small.
        return (Inf, NonNegativeAutoCovConvDetail(v, minimum(g) / v))
    end
    min_cov = min_periodic_corr * v
    j = findfirst(x -> x > min_cov, @view g[i + 1:end])
    if j in (0, length(g) - 1)
        max_cov = maximum(@view g[i + 1:end])
        return (Inf, NonPeriodicConvDetail(v, max_cov / v))
    end
    period = i + j

    err = periodic_mean_error(ftle, period)
    return err, PeriodicConvDetail(period)
end

function periodic_mean_error(xs::AbstractVector, period::Int)
    meanstat = OnlineStats.Mean()
    s0 = OnlineStats.Series(meanstat)
    extstat = OnlineStats.Extrema()
    s1 = OnlineStats.Series(extstat)

    OnlineStats.fit!(s0, @view xs[1:end - period])
    for x in @view xs[end - period + 1:end]
        OnlineStats.fit!(s1, mean(meanstat))
        OnlineStats.fit!(s0, x)
    end
    OnlineStats.fit!(s1, mean(meanstat))

    min, max = OnlineStats.value(extstat)
    return max - min
end

"""Choose appropriate terminator."""
function Terminator(flag::Bool, prob, sol; terminator_options...)
    if flag
        tmnr = AutoCovTerminator(; terminator_options...)
    else
        @assert isempty(terminator_options)
        tmnr = NullTerminator()
    end
    check_cache_type(tmnr, prob, sol)
    return tmnr
end

function Terminator(tmnr::Terminator, prob, sol)
    check_cache_type(tmnr, prob, sol)
    return tmnr
end

has_convergence_history(::Any, i::Int) = false
has_convergence_history(solver::LESolver, i) =
    has_convergence_history(solver.sol, i)
has_convergence_history(sol::LESolution, i) =
    has_convergence_history(sol.convergence, i)
has_convergence_history(history::ConvergenceHistory, i::Int) =
    1 <= i <= length(history.errors)

convergence_history(solver::LESolver) = convergence_history(solver.sol)
convergence_history(sol::LESolution) = sol.convergence


text:
"""
    u(t::Float64)::Float64

Unit Step Function

``\\mathrm{u}(t)= \\begin{cases} 1, &   t \\geq 0 \\\\ 0, &   t < 0 \\end{cases}``

# Examples
```@example
using ISA
t = 2.5
u(t)
```
"""
function u(t::Float64)::Float64 #unit step function
    return ifelse( t < 0, zero(t), one(t) )
end

"""
    𝒩ᵤ(x::Float64; μ::Float64, σ::Float64)::Float64

Unnormalized Gaussian Function

``\\mathcal{N}_\\mathrm{u}(x;\\mu,\\sigma) = \\exp\\left(-\\frac{1}{2}\\left(\\frac{x-\\mu}{\\sigma}\\right)^2\\right).``

# Examples
```@example
using ISA
x = 2.0
μ = 1.0
σ = 0.1
𝒩ᵤ(x; μ, σ)
```
"""
function 𝒩ᵤ(x::Float64; μ::Float64, σ::Float64)::Float64
    return exp(-0.5 * ((x - μ) / σ)^2)
end

"""
    𝒩(x::Float64; μ::Float64, σ::Float64)::Float64

Unit Area Gaussian Function

``\\mathcal{N}(x;\\mu,\\sigma) = \\frac{1}{\\sqrt{2\\pi}\\sigma}\\exp\\left(-\\frac{1}{2}\\left(\\frac{x-\\mu}{\\sigma}\\right)^2\\right).``

# Examples
```@example
using ISA
x = 2.0
μ = 1.0
σ = 0.1
𝒩(x; μ, σ)
```
"""
function 𝒩(x::Float64; μ::Float64, σ::Float64)::Float64
    return 1 / (σ * sqrt(2 * π)) * exp(-0.5 * ((x - μ) / σ)^2)
end

"""
    δn(t::Float64,σ::Float64=1.0e-11)::Float64

Normalized Dirac Delta Approximation

``δ_n(t;\\sigma) = \\exp\\left(-\\frac{1}{2}\\left(\\frac{t^2}{\\sigma^2}\\right)\\right). ``

# Examples
```@example
using ISA
t = 1.0e-10
σ = 1.0e-11
δn(t; σ)
```
"""
function δn(t::Float64,σ::Float64=1.0e-11)::Float64
   return 𝒩ᵤ(t, μ=0.0, σ=σ)
end

"""
    δ(t::Float64,σ::Float64=1.0e-11)::Float64

Dirac Delta Approximation

``δ(t;\\sigma) = \\frac{1}{\\sqrt{2\\pi}\\sigma}\\exp\\left(-\\frac{1}{2}\\left(\\frac{t^2}{\\sigma^2}\\right)\\right). ``

# Examples
```@example
using ISA
t = 1.0e-10
σ = 1.0e-11
δ(t; σ)
```
"""
function δ(t::Float64,σ::Float64=1.0e-11)::Float64
   return 𝒩(t, μ=0.0, σ=σ)
end

ext: jl


text:
include("../required.jl")

#import data
df = CSV.read("WoodEtal_Append1_v2.csv")

# Separate based on WebScale
scale = [df[df.WebScale .== u,:] for u in unique(df.WebScale)]
scale2 = [scale[5][scale[5].WebID .== j,:] for j in unique(scale[5].WebID)]

#notes which scale[x] = scale
#scale[1] = T
#scale[2] = Q
#scale[3] = A
#scale[4] = L
#scale[5] = S

# select for edge interaction columns only
edgelist_T = scale[1][:, [:PredTSN, :PreyTSN]]

# view all entries that contain "san" in PredTSN
pred_san = edgelist_T[startswith.(edgelist_T.PredTSN, "san"),:]
# remove "san"
pred_san.PredTSN = replace.(pred_san.PredTSN, "san" .=> "")
pred_san
#view unique predID without san
unique(pred_san,1)
edgelist_T.PredTSN == "102" #manually check if in original df

# # alternative
# backup = copy(edgelist_T)
# old_ids = unique(backup.PredTSN)
# # keep the indexes before the change
# pred_indexes = startswith.(edgelist_T.PredTSN, "san")
#
# # remove the "san" from the dataframe itself
 edgelist_T.PredTSN = replace.(edgelist_T.PredTSN, "san" .=> "")
#
# #view unique preyIDs without san (make sure they are not repeated in original)
# check_uniques = [in(new_id, old_ids) for new_id in edgelist_T.PredTSN[pred_indexes]]
# filter(isone, check_uniques)

### Do the same for PreyTSN
backup2 = copy(edgelist_T)
old_ids2 = unique(backup2.PreyTSN)
# keep the indexes before the change
prey_indexes = startswith.(edgelist_T.PreyTSN, "san")
# remove the "san" from the dataframe itself
edgelist_T.PreyTSN = replace.(edgelist_T.PreyTSN, "san" .=> "")
edgelist_T
#check unique preyIDs without san
check_uniques2 = [in(new_id2, old_ids2) for new_id2 in edgelist_T.PreyTSN[prey_indexes]]
filter(isone, check_uniques2)

edgelist_T

# convert "string" to "integer"
edgelist_T.PredTSN = parse.(Int64, edgelist_T.PredTSN)
edgelist_T.PreyTSN = parse.(Int64, edgelist_T.PreyTSN)

max_stars_repo_name: kiristern/reseaux

max_issues_repo_name: kiristern/reseaux

max_forks_repo_name: kiristern/reseaux


text:
# #1...	Dissimilarity changes
# ##########################################
include("fctSR.jl")
using PyCall
using PyPlot
mutable struct plotSR
    modus::Int64
    file::Int64
    # SEG::SuperResolution.SegmentSR
    # EVAL::SuperResolution.EvalSR

    function plotSR(SEG::SuperResolution.SegmentSR,EVAL::SuperResolution.EvalSR,fileINPUT,modINPUT=5)
        obj=new()
        obj.modus=modINPUT  # 1..original method, 2..AIC ?
        obj.file=fileINPUT
        select_mode=[   :K_KL_estimated_by_mean,                    #   modINPUT=1
                        :K_KL_estimated_by_median,                  #   modINPUT=2
                        :K_KL_estimated_by_max,                     #   modINPUT=3
                        :K_KL_estimated_by_minDelta,                #   modINPUT=4
                        :K_KL_estimated_by_95quantile]              #   modINPUT=5

        print_mode=[    "K_KL_estimated_by_mean",                    #1
                        "K_KL_estimated_by_median",                  #2
                        "K_KL_estimated_by_max",                     #3
                        "K_KL_estimated_by_minDelta",                #4
                        "K_KL_estimated_by_95quantile"]              #5

        H=heatSR(SEG.loc)
        #getfield(EVAL.wModel[1],fieldnames(Main.SuperResolution.ModelSR)[3])
        for IDX in EVAL._selectIDX
            IDX_norm=cartIDX(EVAL._selectIDX,IDX)
            tmpidx=collect(1:length(IDX_norm))[IDX_norm]
            #collect(fieldnames(SuperResolution.ModelSR) .== Symbol(:dim))
            print("*** plotting IDX = $(IDX) = $(tmpidx) ***\n")
            estimatedlocs=EVAL.wModel[IDX_norm][1].loc[Int64.(length.(EVAL.wModel[IDX_norm][1].loc)./2).==getproperty(EVAL, select_mode[modINPUT])[IDX_norm]][1]
            estimatedmolid=EVAL.wModel[IDX_norm][1].member[Int64.(length.(EVAL.wModel[IDX_norm][1].loc)./2).==getproperty(EVAL, select_mode[modINPUT])[IDX_norm]][1]
            estimatedmolid=relabel(estimatedmolid)
            ################################
            size_geom=exp.(exp.(exp.(getproperty.(EVAL.wModel[IDX_norm][1].geom[Int64.(length.(EVAL.wModel[IDX_norm][1].loc)./2).==getproperty(EVAL, select_mode[modINPUT])[IDX_norm]][1], :p))))
            ################################
            plt.figure()
            plt.title("fileINPUT$(obj.file)_modINPUT$(obj.modus)_part($(IDX[1])-$(IDX[2]))_standW_obs$(SEG.N[CartesianIndex(IDX)])_est$(getproperty(EVAL, select_mode[modINPUT])[IDX_norm])_GrndTrth$(EVAL.K_ground_truth[IDX_norm])")
            plt.scatter(SEG.loc[CartesianIndex(IDX)][:,1],SEG.loc[CartesianIndex(IDX)][:,2],label="data",marker="o",
                            s=log.(SEG.loc[CartesianIndex(IDX)][:,3].-minimum(SEG.loc[CartesianIndex(IDX)][:,3])),
                            c=estimatedmolid,alpha=.2)
            plt.scatter(SEG.MU[CartesianIndex(IDX)][:,1],SEG.MU[CartesianIndex(IDX)][:,2],label="means[GrndTrth]",color="green",marker="x",s=50.0,alpha=0.50)
            plt.scatter(SEG.ground_truth[CartesianIndex(IDX)][:,1],SEG.ground_truth[CartesianIndex(IDX)][:,2],label="loc[GrndTrth]",color="orange",marker="x",s=50.0,alpha=.50)
            plt.scatter(estimatedlocs[:,1],estimatedlocs[:,2],label="loc[est]",color="red",marker=".",s=10 .* size_geom,alpha=0.50)
            plt.legend(loc="best")
            aD=plt.gca()
            aD.set_xlabel("x [nm]")
            aD.set_ylabel("y [nm]")
            PyPlot.savefig(joinpath(@__DIR__,string("SR_Plots/",DATUM,"fileINPUT$(obj.file)_modINPUT$(obj.modus)_part($(IDX[1])-$(IDX[2]))_standW_obs$(SEG.N[CartesianIndex(IDX)])_est$(getproperty(EVAL, select_mode[modINPUT])[IDX_norm])_GrndTrth$(EVAL.K_ground_truth[IDX_norm]).pdf")),format="pdf",dpi=300)
        end
    end
end
# #using Plots
# X=Main.SuperResolution.SEGMENTtest.loc
# H=heatSR(X)
#
# Hhist1=map(i->sum(H[i,:]),1:size(H,1))
# Hhist2=map(i->sum(H[:,i]),1:size(H,1))
# n=size(H,1)
# plt[:figure]()
# plt[:title]("part$(1)-$(n)_obs$(sum(H))_est$(0)-sim$(0)")
# plt[:bar](collect(1:n),1 .* Hhist1[:],label="colSum",color="blue",alpha=0.5)#[end-tmpSelect:end],label="minKLD",color="blue",alpha=0.5)
# plt[:bar](collect(1:n),1 .* Hhist2[:],label="rowSum",color="green",alpha=0.5)#[end-tmpSelect:end],label="Δ minKLD",color="green",alpha=0.5)
# plt[:legend](loc="best")
# scl=Int64(ceil(sqrt(n))/2)
# lb=map(i->collect(1:n)[scl*i],1:Int64(floor(n/scl)))
# plt[:xticks](lb,fontsize=3)
# aD=plt[:gca]()
# aD[:set_xlabel]("col/row")
# aD[:set_ylabel]("sum")
# savefig("/Users/Patrick/Documents/GitHub/SuperResolution/src/SR_Plots/ORI_KLC_part$(1)-$(n)_obs$(sum(H))_est$(0)-sim$(0).pdf")
# # ##########################################
# plt[:figure]()
# plt[:title]("part$(partN)-$(partN2)_obs$(N)_est$(solvNumClust)-sim$(size(unique(𝕄sim),1))")
# plt[:scatter](𝕐sim[:,1],𝕐sim[:,2],label="molID[est]",marker="o",s=30.0,alpha=.2, c=Acoll[solvIDX])
# plt[:scatter](𝛍sim[:,1],𝛍sim[:,2],label="means[sim]",color="green",marker="x",s=50.0,alpha=0.50)
# plt[:scatter](ℝsim[:,1],ℝsim[:,2],label="loc[sim]",color="orange",marker="x",s=50.0,alpha=.50)
# plt[:scatter](𝕐simclAll[solvIDX][:,1],𝕐simclAll[solvIDX][:,2],label="loc[est]",color="red",marker=".",s=10.0,alpha=1.0)
# plt[:legend](loc="best")
# map(i->annotate(AfinalString[i],
# xy=[𝕐simclAll[solvIDX][i,1];𝕐simclAll[solvIDX][i,2]],
# xytext=[𝕐simclAll[solvIDX][i,1];𝕐simclAll[solvIDX][i,2]],fontsize=2.0),1:uniqueLength)
# aD=plt[:gca]()
# aD[:set_xlabel]("x [nm]")
# aD[:set_ylabel]("y [nm]")
# savefig("/Users/Patrick/Documents/GitHub/SuperResolution/src/SR_Plots/ORI_KLC_part$(partN)-$(partN2)_obs$(N)_est$(solvNumClust)-sim$(size(unique(𝕄sim),1)).pdf",format="pdf",dpi=300)
#
#
# # using pyplot() make the above syntax plt[:figure]() unfunctional..
# pyplot()
# HM=heatmap(H)
# savefig(HM,"/Users/Patrick/Documents/GitHub/SuperResolution/src/SR_Plots/01_ORI_KLC_heatmap.pdf")

# #2...	2D Model
# ##########################################

#
#
# using PyCall
# using PyPlot
#
# function plotSR(SEG::SuperResolution.SegmentSR,RESTR_SEG::Array{Tuple{Int64,Int64},1})
# #1...Dissimilarity changes
# ##########################################
# #heatmap for observations in over segmented space
# X=SEG.loc
# H=heatSR(X)
# pyplot()
# heatmap(H)
# savefig(HM,"/Users/Patrick/Documents/GitHub/SuperResolution/src/SR_Plots/01_ORI_KLC_heatmap.pdf",format="pdf",dpi=300)
#
# Hhist1=map(i->sum(H[i,:]),1:size(H,1))
# Hhist2=map(i->sum(H[:,i]),1:size(H,1))
#
# n=size(H,1)
#
# pyplot()
# heatmap(H)
# savefig(HM,"/Users/Patrick/Documents/GitHub/SuperResolution/src/SR_Plots/ORI_KLC_heatmap.pdf",format="pdf",dpi=300)
#
# plt[:figure]()
# plt[:title]("part$(1)-$(n)_obs$(sum(H))_est$(0)-sim$(0)")
# #plt[:bar](xRangeSel,1 .* maxKLDeval[end-tmpSelect:end],label="maxKLD",color="red",alpha=0.5)
# #plt[:bar](xRangeSel,1 .* meanKLDeval[end-tmpSelect:end],label="meanKLD",color="grey",alpha=0.5)
# plt[:bar](collect(1:n),1 .* Hhist1[:],label="colSum",color="blue",alpha=0.5)#[end-tmpSelect:end],label="minKLD",color="blue",alpha=0.5)
# plt[:bar](collect(1:n),1 .* Hhist2[:],label="rowSum",color="green",alpha=0.5)#[end-tmpSelect:end],label="Δ minKLD",color="green",alpha=0.5)
# plt[:legend](loc="best")
# #lb=collect(1:n)
# scl=Int64(ceil(sqrt(n))/2)
# lb=map(i->collect(1:n)[scl*i],1:Int64(floor(n/scl)))
# plt[:xticks](lb,fontsize=3)
# #lb=map(i->"$(Knum[i])",collect(1:(ι))[1:end])#[end-tmpSelect:end])
# #plt[:xticks](xRangeSel,lb,fontsize=10.0/log(length(xRangeSel)))
# aD=plt[:gca]()
# aD[:set_xlabel]("col/row")
# aD[:set_ylabel]("sum")
# savefig("/Users/Patrick/Documents/GitHub/SuperResolution/src/SR_Plots/ORI_KLC_part$(1)-$(n)_obs$(sum(H))_est$(0)-sim$(0).pdf",format="pdf",dpi=300)
# ##########################################
# 	plt[:figure]()
# 	plt[:title]("part$(partN)-$(partN2)_obs$(N)_est$(solvNumClust)-sim$(size(unique(𝕄sim),1))")
# 	#plt[:bar](xRangeSel,1 .* maxKLDeval[end-tmpSelect:end],label="maxKLD",color="red",alpha=0.5)
# 	#plt[:bar](xRangeSel,1 .* meanKLDeval[end-tmpSelect:end],label="meanKLD",color="grey",alpha=0.5)
# 	plt[:bar](xRangeSel,1 .* minKLDeval[:],label="minKLD",color="blue",alpha=0.5)#[end-tmpSelect:end],label="minKLD",color="blue",alpha=0.5)
# 	plt[:bar](xRangeSel,1 .* abs_minKLDeval[:],label="Δ minKLD",color="green",alpha=0.5)#[end-tmpSelect:end],label="Δ minKLD",color="green",alpha=0.5)
# 	plt[:legend](loc="best")
# 	lb=map(i->"$(Knum[i])",collect(1:(ι))[1:end])#[end-tmpSelect:end])
# 	lb=lb[1:Int64(ceil(log(length(lb)))):end]
# 	plt[:xticks](xRangeSel[1:Int64(ceil(log(length(xRangeSel)))):end],lb,fontsize=10.0/log(length(xRangeSel)))
# 	#lb=map(i->"$(Knum[i])",collect(1:(ι))[1:end])#[end-tmpSelect:end])
# 	#plt[:xticks](xRangeSel,lb,fontsize=10.0/log(length(xRangeSel)))
# 	aD=plt[:gca]()
# 	aD[:set_xlabel]("κ")
# 	aD[:set_ylabel]("KLDvalue")
# 	savefig("/Users/Patrick/Documents/GitHub/SuperResolution/src/SR_Plots/ORI_KLC_part$(partN)-$(partN2)_obs$(N)_est$(solvNumClust)-sim$(size(unique(𝕄sim),1))_EVAL.pdf",format="pdf",dpi=300)
#
# 	#EVALUATE ERRORS FOR EACH CLUSTER
# 	#detection of molecules per Tetramere
# 	#pairwise evaluation sim data vs. computed results
#
# 	#generate string of cluster IDs for plotting
# 	uniqueLength=length(unique(Acoll[solvIDX-1]))
# 	AfinalString=map(i->string(Aall[solvIDX-1][i]),1:uniqueLength)
# #2...Solution 2D Plot with corresponding cluster coloring
# 	plt[:figure]()
# 	plt[:title]("part$(partN)-$(partN2)_obs$(N)_est$(solvNumClust)-sim$(size(unique(𝕄sim),1))")
# 	plt[:scatter](𝕐sim[:,1],𝕐sim[:,2],label="molID[est]",marker="o",s=30.0,alpha=.2, c=Acoll[solvIDX])
# 	plt[:scatter](𝛍sim[:,1],𝛍sim[:,2],label="means[sim]",color="green",marker="x",s=50.0,alpha=0.50)
# 	plt[:scatter](ℝsim[:,1],ℝsim[:,2],label="loc[sim]",color="orange",marker="x",s=50.0,alpha=.50)
# 	plt[:scatter](𝕐simclAll[solvIDX][:,1],𝕐simclAll[solvIDX][:,2],label="loc[est]",color="red",marker=".",s=10.0,alpha=1.0)
# 	plt[:legend](loc="best")
# 	map(i->annotate(AfinalString[i],
# 	xy=[𝕐simclAll[solvIDX][i,1];𝕐simclAll[solvIDX][i,2]],
# 	xytext=[𝕐simclAll[solvIDX][i,1];𝕐simclAll[solvIDX][i,2]],fontsize=2.0),1:uniqueLength)
# 	aD=plt[:gca]()
# 	aD[:set_xlabel]("x [nm]")
# 	aD[:set_ylabel]("y [nm]")
# 	savefig("/Users/Patrick/Documents/GitHub/SuperResolution/src/SR_Plots/ORI_KLC_part$(partN)-$(partN2)_obs$(N)_est$(solvNumClust)-sim$(size(unique(𝕄sim),1)).pdf",format="pdf",dpi=300)
#
# 	return DataSim, RealLoc
# end
#
# export importSR

export plotSR

max_stars_repo_name: PGSch/SuperResolution

max_issues_repo_name: PGSch/SuperResolution

max_forks_repo_name: PGSch/SuperResolution


text:
############################################################################################
#Markov Kernel
markov_MV_data = [randn(_rng, 2) for _ in 1:N]
param_HMM_MV = (;
    μ = Param([[-2., -2.], [0., 0.], [2., 2.]],
        [MvNormal( [i, i], [1., 1.] ) for i in -1.0:1:1.]),
    σ = Param([[ 2.2 0.2 ; 0.2 2.2], [ 1.1 0.1 ; 0.1 1.1], [ 2.2 0.2 ; 0.2 2.2] ],
        [InverseWishart(10, [0.8 0.5 ; 0.5 .8] ) for i in 1:3]),
    p = Param([[.6, .2, .2],  [.1, .1, .8], [.1, .1, .8]], [Dirichlet(3,3) for i in 1:3]),
    latent = Param(markov_latent, Fixed()),
)
struct HMM_MV <: ModelName end
hmm_MV = ModelWrapper(HMM_MV(), param_HMM_MV)
markov_MV_objective = Objective(hmm_MV, markov_MV_data, :latent)

################################################################################
function get_dynamics(model::ModelWrapper{<:HMM_MV}, θ)
    @unpack μ, σ, p = θ
    dynamicsᵉ = [ MvNormal(μ[iter], σ[iter]) for iter in eachindex(μ) ]
    dynamicsˢ = [ Categorical( p[iter] ) for iter in eachindex(μ) ]
    return dynamicsᵉ, dynamicsˢ
end

function BaytesFilters.dynamics(objective::Objective{<:ModelWrapper{<:HMM_MV}})
    @unpack model, data = objective
    dynamicsᵉ, dynamicsˢ = get_dynamics(model, model.val)

    initialˢ = Categorical(fill(1 / length(dynamicsᵉ), length(dynamicsᵉ)))
    transition(particles, iter) = dynamicsˢ[particles[iter - 1]]
    evidence(particles, iter) = dynamicsᵉ[particles[iter]]

    return Markov(initialˢ, transition, evidence)
end
#dynamics(markov_MV_objective)

############################################################################################
"Forward Filter HMM"
function filter_forward(objective::Objective{<:ModelWrapper{HMM_MV}})
    @unpack model, data = objective
## Map Parameter to observation and state probabilities
    @unpack p = model.val
    dynamicsᵉ, _ = get_dynamics(model, model.val)
#!NOTE: Working straight with parameter instead of distributions here as easier to implement
    dynamicsˢ           = transpose( reduce(hcat, p ) ) #[ Categorical(p[iter]) for iter in eachindex(μ) ]
    initialˢ            = get_stationary( dynamicsˢ )
    structure = BaytesCore.ByRows()
## Assign Log likelihood
    ℓℒᵈᵃᵗᵃ = zeros(Float64, size(data,1), size(dynamicsˢ, 1) )
    Base.Threads.@threads for state in Base.OneTo( size(ℓℒᵈᵃᵗᵃ, 2) )
    for iter in Base.OneTo( size(ℓℒᵈᵃᵗᵃ, 1) )
            ℓℒᵈᵃᵗᵃ[iter, state] += logpdf( dynamicsᵉ[state], grab(data, iter, structure) )
        end
    end
## Initialize
    Nstates = size(dynamicsˢ, 1)
    α = zeros( size(ℓℒᵈᵃᵗᵃ) ) # α represents P( sₜ | e₁:ₜ ) for each t, which is numerically more stable than classical forward probabilities p(sₜ, e₁:ₜ)
    c = vec_maximum( @view(ℓℒᵈᵃᵗᵃ[1,:]) ) # Used for stable logsum() calculations for incremental log likelihood addition, i.e.: log( exp(x) + exp(y) ) == x + log(1 + exp(y -x ) ) for y >= x.
    ℓℒ = 0.0 #Log likelihood container
## Calculate initial probabilities
    for stateₜ in Base.OneTo(Nstates)
         α[1,stateₜ] += initialˢ[stateₜ] * exp(ℓℒᵈᵃᵗᵃ[1,stateₜ]-c) #Calculate initial p(s₁, e₁) ∝ P(s₁) * p(e₁ | s₁)
    end
## Normalize P( s₁, e₁ ) and return normalizing constant P(e₁), which can be used to calculate ℓℒ = P(e₁:ₜ) = p(e₁) ∏_k p(eₖ | e₁:ₖ₋₁)
    norm = sum( @view(α[1,:]) )
    α[1,:] ./= norm
    ℓℒ += log(norm) + c# log(norm) = p(eₜ | e₁:ₜ₋₁), c is constant that is added back after removing on top for numerical stability
## Loop through sequence
    for t = 2:size(α, 1)
        c = vec_maximum( @view(ℓℒᵈᵃᵗᵃ[t,:]) ) # Used for stable logsum() calculations, i.e.: log( exp(x) + exp(y) ) == x + log(1 + exp(y -x ) ) for y >= x.
        ## Calculate ∑_sₜ₋₁ P( sₜ | sₜ₋₁) * P(sₜ₋₁ | e₁:ₜ₋₁) - we sum over sₜ₋₁ - states are per row
        for stateₜ in Base.OneTo(Nstates)
            for stateₜ₋₁ in Base.OneTo(Nstates)
                α[t,stateₜ] += dynamicsˢ[stateₜ₋₁, stateₜ] * α[t-1,stateₜ₋₁] # * for both log and non-log version as inside log( sum(probability(...)))
            end
            ## Then multiply with ℒᵈᵃᵗᵃ corrector P( eₜ | sₜ )
            α[t,stateₜ] *= exp(ℓℒᵈᵃᵗᵃ[t,stateₜ]-c) # - c for higher numerical stability, will be added back to likelihood increment below
        end
        ## Normalize α and obtain P( eₜ | e₁:ₜ₋₁)
        norm = sum( @view(α[t,:]) )
        α[t,:] ./= norm
        ## Add normalizing constant p(eₜ | e₁:ₜ₋₁) to likelihood term
        ℓℒ += log(norm)+c
    end
    return (α, ℓℒ)
end
filter_forward(markov_MV_objective)

############################################################################################
function ModelWrappers.simulate(rng::Random.AbstractRNG, model::ModelWrapper{HMM_MV}; Nsamples = 1000)
    dynamicsᵉ, dynamicsˢ = get_dynamics(model, model.val)
    latentⁿᵉʷ = Vector{eltype(model.val.latent)}(undef, Nsamples)
    observedⁿᵉʷ = [zeros(2) for _ in 1:Nsamples]

    stateₜ = rand(rng, dynamicsˢ[ rand(1:length(dynamicsˢ) ) ] )
    latentⁿᵉʷ[1] = stateₜ
    observedⁿᵉʷ[1] = rand(rng, dynamicsᵉ[stateₜ])

    for iter in 2:size(observedⁿᵉʷ,1)
            stateₜ = rand(rng, dynamicsˢ[ stateₜ ]) #stateₜ for t-1 overwritten
            latentⁿᵉʷ[iter] = stateₜ
            observedⁿᵉʷ[iter] = rand(rng, dynamicsᵉ[stateₜ])
    end
    return observedⁿᵉʷ, latentⁿᵉʷ
end


text:
@testset "Unconstrained" begin
  @testset "Basic" begin
    # These basic tests are simple ways of checking that your solver isn't breaking.
    # nlp = ADNLPModel(
    #   # x -> (x[1] - 1)^2 + (x[2] - 2)^2 / 4,
    #   # zeros(2)
    #   x -> (x[1] - 1.0)^2 + 100 *(x[2] - x[1]^2)^2,
    # [-1.2; 1.0]
    # )
    nlp = CUTEstModel("ROSENBR")
    # nlp = CUTEstModel("OSBORNEA")
    # nlp = CUTEstModel("DMN37142LS")

    # output = with_logger(NullLogger()) do
    #   uncsolver(nlp)

    # output = bfgs_bl(nlp)
    # println("bfgs-bl")
    # print(output)

    # output = bfgs_rc(nlp)
    # println("bfgs-rc")
    # print(output)

    # output = l_bfgs_rcst(nlp)
    # println("l-bfgs-rc")
    # print(output)

    output = newtoncombusca(nlp)
    println("newtoncombusca")
    print(output)

    #  @test isapprox(output.solution, [1.0; 1.0], rtol=1e-2)
    #  @test output.objective < 1e-3
    #  @test output.dual_feas < 1e-3
    #  @test output.status == :first_order

    # output = with_logger(NullLogger()) do
    #   gradiente(nlp)
    # end
    # #@test isapprox(output.solution, [1.0; 2.0], rtol=1e-4)
    # @test output.objective < 1e-4
    # @test output.dual_feas < 1e-4
    # @test output.status == :first_order

  end

  # @testset "Failures" begin
  #   # This checks that your solver is failing accordingly
  #   nlp = ADNLPModel(
  #     x -> (x[1] - 1.0)^2 + 100 *(x[2] - x[1]^2)^2,
  #   [-0.1; 1.0]
  #   )
  #   output = with_logger(NullLogger()) do
  #     bfgs_rc(nlp, max_eval = 1)
  #   end
  #   @test output.status == :max_eval
  #   output = with_logger(NullLogger()) do
  #     bfgs_rc(nlp, max_iter = 1)
  #   end
  #   @test output.status == :max_iter
  #   output = with_logger(NullLogger()) do
  #     bfgs_rc(nlp, max_time = 1e-8)
  #   end
  #   @test output.status == :max_time
  # end
end


text:
# set voltage limit
# SU
export psp_setvlimit

"""
    psp_setvlimit(io, voltage_limit)

Set the voltage limit.
"""
function psp_setvlimit(io_psp::IO,vlim::Real)
  @assert vlim<100 "voltage limit must be < 100V"
  vlimstring = @sprintf("%02d",vlim)
  set_voltage_limit_string = "SU "*vlimstring*"\r"
  write(io_psp, set_voltage_limit_string)
  return nothing
end


text:
using SparseArrays

A = [2 3; 1 0; 0 6]

@test @inferred sparsify(A, 0.4, Inf, 200, true) ≈ sparse([2.0 3.0; 1.0 0; 0 6.0])
@test @inferred sparsify(A', 0.4, Inf, 200, true) ≈ sparse([2.0 3.0; 1.0 0; 0 6.0]')
A_id = @inferred sparsity_pattern(A, 0.4, Inf)
A_id_t = @inferred sparsity_pattern(A', 0.4, Inf)
@test A_id ≈ sparse([1 1; 1 0; 0 1])
@test A_id_t ≈ sparse([1 1; 1 0; 0 1]')
@test @inferred bin_sparse_matrix!(A, A_id, 200) ≈ sparse([1 3; 2 0; 0 4])
@test @inferred bin_sparse_matrix!(A', A_id_t, 200) ≈ sparse([1 3 0; 2 0 4])


text:
module ProbeBase

using Preferences

abstract type AbstractPayload end

function null_payload(x...) end

struct ProbeSpec
    lineno::LineNumberNode
    kind::Symbol
    argtypes::Vector{DataType}
    payload::Base.RefValue{Ptr{Cvoid}}
    semaphore::Threads.Atomic{Int}
end
ProbeSpec(lineno::LineNumberNode, kind::Symbol, argtypes::Vector, payload::Ptr{Cvoid}) =
    ProbeSpec(lineno, kind, argtypes,
              Ref(payload), Threads.Atomic{Int}(PROBES_ENABLED))
ProbeSpec(lineno::LineNumberNode, kind::Symbol, argtypes::Vector) =
    ProbeSpec(lineno, kind, argtypes,
              # null_payload is varargs, and doesn't use arguments anyway
              Ref(@cfunction(null_payload, Cvoid, ())),
              Threads.Atomic{Int}(PROBES_ENABLED))

const PROBES = Dict{Module,Dict{Symbol,ProbeSpec}}()
const PROBES_FUNC = Dict{Module,Dict{Symbol,Any}}()
const PROBES_ENABLED = parse(Bool, @load_preference("default_enabled", "false"))
const PROBES_LOCK = Threads.ReentrantLock()

default_enabled(val::Bool) = @set_preferences!("default_enabled" => repr(val))

function set_probe_payload(mod::Module, category::Symbol, payload::Ptr{Cvoid})
    lock(PROBES_LOCK) do
        PROBES[mod][category].payload[] = payload
    end
end
function set_probe_payload(f, mod::Module, category::Symbol)
    spec = nothing
    lock(PROBES_LOCK) do
        spec = PROBES[mod][category]
        PROBES_FUNC[mod][category] = f
    end
    ptr = eval(:(@cfunction($f, Cvoid, ($(spec.argtypes...),))))
    set_probe_payload(mod, category, ptr)
end
function enable_probe(mod::Module, category::Symbol)
    spec = lock(PROBES_LOCK) do
        PROBES[mod][category]
    end
    Threads.atomic_add!(spec.semaphore, 1)
end
function disable_probe(mod::Module, category::Symbol)
    spec = lock(PROBES_LOCK) do
        PROBES[mod][category]
    end
    Threads.atomic_sub!(spec.semaphore, 1)
end
function enable_probes(mod::Module)
    lock(PROBES_LOCK) do
        for spec in values(PROBES[mod])
            Threads.atomic_add!(spec.semaphore, 1)
        end
    end
end
function disable_probes(mod::Module)
    lock(PROBES_LOCK) do
        for spec in values(PROBES[mod])
            Threads.atomic_sub!(spec.semaphore, 1)
        end
    end
end

probe_enabled(spec::ProbeSpec) = spec.semaphore[] > 0
@generated function probe_trigger(spec::ProbeSpec, args...)
    Targs = Expr(:tuple, map(nameof, args)...)
    ex = Expr(:call, :ccall, :(spec.payload[]), :Cvoid, Targs)
    for (idx, arg) in enumerate(args)
        push!(ex.args, :(args[$idx]))
    end
    ex
end
function probe_maybe_trigger(spec::ProbeSpec, kind::Symbol, args...)
    if probe_enabled(spec)
        probe_trigger(spec, args...)
    end
end

function parse_args(args)
    argnames = Symbol[]
    argtypes = DataType[]
    argvalues = Expr[]
    for arg in args.args
        if Meta.isexpr(arg, :(=))
            push!(argnames, arg.args[1].args[1])
            push!(argtypes, eval(arg.args[1].args[2]))
            push!(argvalues, esc(arg.args[2]))
        elseif Meta.isexpr(arg, :(::))
            push!(argnames, arg.args[1])
            push!(argtypes, eval(arg.args[2]))
            push!(argvalues, esc(arg.args[1]))
        else
            throw(ArgumentError("Cannot handle expr with head: $(arg.head)"))
        end
    end
    (;argnames, argtypes, argvalues)
end
function register_probe(mod::Module, category::Symbol, spec::ProbeSpec; enable=false)
    if enable
        spec.semaphore[] = 1
    end
    lock(PROBES_LOCK) do
        module_specs = get!(PROBES, mod) do
            Dict{Symbol,ProbeSpec}()
        end
        module_specs[category] = spec
        (get!(PROBES_FUNC, mod) do
            Dict{Symbol,Any}()
        end)[category] = nothing
    end
end
macro probe(kind::QuoteNode, category::QuoteNode, args=Expr(:tuple))
    argspec = parse_args(args)
    spec = ProbeSpec(__source__, kind.value, argspec.argtypes)
    register_probe(__module__, category.value, spec)
    :(probe_maybe_trigger($spec, $kind, $(argspec.argvalues...)))
end
#= FIXME
macro probe_start(category::Symbol)
    :(probe($spec, :start))
end
macro probe_stop(category::Symbol)
    :(probe($spec, :stop))
end
macro probe_event(category::Symbol)
    :(probe($spec, :event))
end
=#

end # module


text:
#!/usr/bin/env julia

max_threads = 12

using Distributed
if Distributed.nprocs() <= 2
    threads = min(trunc(Int, Sys.CPU_THREADS*0.75), max_threads)
    println("threads: $(threads)/$(Sys.CPU_THREADS)")

    proc_ids = Distributed.addprocs(threads)
    println("process ids: $(proc_ids)")
end

processes = Distributed.nprocs()-1 # save one for the master process
println("processes ids: $(processes)")

@everywhere function test(pid)
    println("pid: $(pid) - $(gethostname()) - $(pwd())")
end

solution2_files = pmap(test, 1:processes)

max_stars_repo_name: ccoffrin/julia-examples

max_issues_repo_name: ccoffrin/julia-examples

max_forks_repo_name: ccoffrin/julia-examples


text:
using LinearAlgebra
using CSV
using DataFrames

# selecting the index combination from given range
function getIndexes(i,j,d)
    if i == j
        return [i]
    elseif i > j && i > 2
        return 0
    else
        niz1 = []
        niz2 = collect(i:j-1)
        for k in 1:d
            if !(k in niz2)
                push!(niz1,k)
            end
        end
        return length(niz1) == 0 ? 0 : niz1
    end
end

#selecting columns from indexes as matrix
function getSelectedColumns(selectedColumns,X)
    if selectedColumns != 0
        localX = ones(size(X)[1],1)
        for i in selectedColumns
            localX = [localX X[:,i]]
        end
        return localX
    else
        return 0
    end
end

#load data
Diabete_Data = CSV.read("diabetes.data",DataFrame)

n,d = size(Diabete_Data)
d = d-1

# compute the optimal theta

X = zeros(n,d)
y = zeros(n,1)

for i=1:d
    X[:,i] = Diabete_Data[:,i]
end

y = Diabete_Data[:,end]

# for storing MSEs
savedMSE = []

#for  storing top contender colums
savedIndexes = []

#  age
#  sex
#  bmi
#  map	   mean arterial pressure
#  s1  tc :   total cholesterol
#  s2  ldl    low density lipoprotein
#  s3  hdl    high density lipoprotein
#  s4  tch
#  s5  ltg
#  s6  glu

# string representation of said colums
nameByIndex = ["AGE" "SEX" "BMI" "MAP" "TC" "LDL" "HDL" "TCH" "LTG" "GLU"]

#Calculating all possible solutions
for i in 1:d-1
    for j in 1:d
        localX = []
        temp = getIndexes(i,j,d)
        selectedColumns = temp == 0 ? continue : temp
        localX = getSelectedColumns(selectedColumns,X)
        localX = localX == 0 ? continue : localX
        theta_opt = localX\y
        # achieved loss
        MSE = norm(localX*theta_opt-y,2)^2/n
        push!(savedMSE,MSE)
        push!(savedIndexes,selectedColumns)
    end
end

# Number of top results
numOfBests = 10

#Finding top results
for i in 1:numOfBests
   bestIndex = argmin(savedMSE)
    bestMSE = savedMSE[bestIndex]
    leastColums = savedIndexes[bestIndex]
    println("\n [$i] BEST SOLUTION")
    println("Best MSE of: $(bestMSE) with ")
    for j in leastColums
        println("$(nameByIndex[j]) ")
    end
    println("parameters as most impactful")
    deleteat!(savedMSE,bestIndex)
    deleteat!(savedIndexes,bestIndex)
end

max_stars_repo_path: Primer Predavanje/analyse.jl

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_path: Primer Predavanje/analyse.jl

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_path: Primer Predavanje/analyse.jl

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
"""
    abstract type AbstractSingularityHandler{R}

Used for handling localized integrand singularities in `R`.
"""
abstract type AbstractSingularityHandler{T} end

"""
    struct IMT{A,P} <: AbstractSingularityHandler{ReferenceLine}

One-dimensional change of variables mapping `[0,1] -> [0,1]` with the property that
all derivatives vanish at the point `x=0`.

See [Davis and Rabinowitz](https://www.elsevier.com/books/methods-of-numerical-integration/davis/978-0-12-206360-2).
"""
struct IMT{A,P} <: AbstractSingularityHandler{ReferenceLine}
end
IMT(;a=1,p=1) = IMT{a,p}()

domain(::IMT) = ReferenceLine()
image(::IMT)  = ReferenceLine()

function (f::IMT{A,P})(x) where {A,P}
    exp(A * (1 - 1 / x[1]^P))
end

derivative(f::IMT{A,P},x) where {A,P} = f(x) * A * P * 1 / x[1]^(P + 1)

jacobian(f::IMT,x) = derivative(f, x) |> SMatrix{1,1}

"""
    struct Kress{P} <: AbstractSingularityHandler{ReferenceLine}

Change of variables mapping `[0,1]` to `[0,1]` with the property that the first
`P-1` derivatives of the transformation vanish at `x=0`.
"""
struct Kress{P} <: AbstractSingularityHandler{ReferenceLine}
end
Kress(;order=5) = Kress{order}()

domain(k::Kress) = ReferenceLine()
image(k::Kress)  = ReferenceLine()

# NOTE: fastmath is needed here to allow for various algebraic simplifications
# which are not exact in floating arithmetic. Maybe reorder the operations *by
# hand* to avoid having to use fastmath? In any case, benchmark first.
function (f::Kress{P})(y) where {P}
    x = y[1]
    v = (x) -> (1 / P - 1 / 2) * ((1 - x))^3 + 1 / P * ((x - 1)) + 1 / 2
    return 2v(x)^P / (v(x)^P + v(2 - x)^P)
end

function derivative(f::Kress{P}, y) where {P}
    x = y[1]
    v = (x) -> (1 / P - 1 / 2) * ((1 - x))^3 + 1 / P * ((x - 1)) + 1 / 2
    vp = (x) -> -3 * (1 / P - 1 / 2) * ((1 - x))^2 + 1 / P
    return 2 * (P * v(x)^(P - 1) * vp(x) * (v(x)^P + v(2 - x)^P) - (P * v(x)^(P - 1) * vp(x) - P * v(2 - x)^(P - 1) * vp(2 - x) ) * v(x)^P ) /
        (v(x)^P + v(2 - x)^P)^2
end

jacobian(f::Kress,x) = derivative(f, x) |> SMatrix{1,1}

struct KressR{P} <: AbstractSingularityHandler{ReferenceLine}
end
KressR(;order=5) = KressR{order}()

domain(k::KressR) = ReferenceLine()
image(k::KressR)  = ReferenceLine()

# NOTE: fastmath is needed here to allow for various algebraic simplifications
# which are not exact in floating arithmetic. Maybe reorder the operations *by
# hand* to avoid having to use fastmath? In any case, benchmark first.
function (f::KressR{P})(y) where {P}
    x = 1-y[1]
    v = (x) -> (1 / P - 1 / 2) * ((1 - x))^3 + 1 / P * ((x - 1)) + 1 / 2
    return 1 - 2v(x)^P / (v(x)^P + v(2 - x)^P)
end

function derivative(f::KressR{P}, y) where {P}
    x = 1-y[1]
    v = (x) -> (1 / P - 1 / 2) * ((1 - x))^3 + 1 / P * ((x - 1)) + 1 / 2
    vp = (x) -> -3 * (1 / P - 1 / 2) * ((1 - x))^2 + 1 / P
    return 2 * (P * v(x)^(P - 1) * vp(x) * (v(x)^P + v(2 - x)^P) - (P * v(x)^(P - 1) * vp(x) - P * v(2 - x)^(P - 1) * vp(2 - x) ) * v(x)^P ) /
        (v(x)^P + v(2 - x)^P)^2
end

jacobian(f::KressR,x) = derivative(f, x) |> SMatrix{1,1}

"""
    struct KressP{P} <: AbstractSingularityHandler{ReferenceLine}

Like [`Kress`](@ref), this change of variables maps the interval `[0,1]` onto
itself, but the first `P` derivatives of the transformation vanish at **both**
endpoints.

This change of variables can be used to *periodize* integrals in the following
sense. Suppose we wish to compute the integral of `f(x)` from `0` to `1` where
`f is not a `1-periodic function. If `ϕ` is an object of type `KressP`, then
using it as a change of variables in the integration yields a similar integral
from `0` to `1` (the interval `0≤0≤1` is mappend onto itself), but with
integrand given by `g(x) = f(ϕ(x))ϕ'(x)`. Since `ϕ'` vanishes (together with `P`
of its derivatives), the function `g(x)` is now periodic (up to derivatives of
order up to `P`) at the endpoints. Thus quadrature rules designed for periodic
functions like the [`TrapezoidalOpen`](@ref) can be used to obtain high order
convergence of `g`, which in turn yields a modified quadrature rule when viewed
as a quadrature rule for `f`.
"""
struct KressP{P} <: AbstractSingularityHandler{ReferenceLine}
end
KressP(;order=5) = KressP{order}()

domain(k::KressP) = ReferenceLine()
image(k::KressP)  = ReferenceLine()

@fastmath function (f::KressP{P})(y) where {P}
    x = y[1]
    v = (x) -> (1 / P - 1 / 2) * ((1 - 2x))^3 + 1 / P * ((2x - 1)) + 1 / 2
    return v(x)^P / (v(x)^P + v(1 - x)^P)
end

@fastmath function derivative(f::KressP{P}, y) where {P}
    x = y[1]
    v =  (x) -> (1 / P - 1 / 2) * ((1 - 2x))^3 + 1 / P * ((2x - 1)) + 1 / 2
    vp = (x) -> -6 * (1 / P - 1 / 2) * ((1 - 2x))^2 + 2 / P
    return (P * v(x)^(P - 1) * vp(x) * (v(x)^P + v(1 - x)^P) - (P * v(x)^(P - 1) * vp(x) - P * v(1 - x)^(P - 1) * vp(1 - x) ) * v(x)^P ) /
        (v(x)^P + v(1 - x)^P)^2
end

jacobian(f::KressP,x) = derivative(f, x) |> SMatrix{1,1}

"""
    struct Duffy <: AbstractSingularityHandler{RefereceTriangle}

Change of variables mapping the `ReferenceSquare` to the `RefereceTriangle` with
the property that the jacobian vanishes at the `(1,0)` vertex of the triangle.

Useful for integrating functions with a singularity on the `(1,0)` edge of the
reference triangle.
"""
struct Duffy <: AbstractSingularityHandler{ReferenceTriangle} end

domain(::Duffy) = ReferenceSquare()
image(::Duffy)  = ReferenceTriangle()

function (::Duffy)(u)
    SVector(u[1], (1 - u[1]) * u[2])
end

function jacobian(::Duffy, u)
    SMatrix{2,2,Float64}(1, 0, -u[2][1], (1 - u[1][1]))
end

# TODO: generalize to `N` dimensions
"""
    struct TensorProductSingularityHandler{S} <: AbstractSingularityHandler{ReferenceSquare}

A tensor product of two one-dimensional `AbstractSingularityHandler`s for performing
integration over the `ReferenceSquare`.
"""
struct TensorProductSingularityHandler{S} <: AbstractSingularityHandler{ReferenceSquare}
    shandler::S
end

domain(::TensorProductSingularityHandler) = ReferenceSquare()
image(::TensorProductSingularityHandler) = ReferenceSquare()

function TensorProductSingularityHandler(q...)
    TensorProductSingularityHandler(q)
end

TensorProductSingularityHandler{Tuple{P,Q}}() where {P,Q} = TensorProductSingularityHandler(P(), Q())

function (f::TensorProductSingularityHandler)(x)
    shandler = f.shandler
    @assert length(shandler) == length(x)
    N = length(shandler)
    svector(i -> shandler[i](x[i]), N)
end

function jacobian(f::TensorProductSingularityHandler, x)
    shandler = f.shandler
    @assert length(shandler) == length(x)
    N = length(shandler)
    if N == 2
        jx = derivative(shandler[1], x[1])
        jy = derivative(shandler[2], x[2])
        return SMatrix{2,2}(jx, 0, 0, jy)
    else
        notimplemented()
    end
end

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
using LargeColumns
using Base.Test

using LargeColumns:
    MmappedColumns, SinkColumns, meta_path,
    # internals
    fixed_Tuple_types, representative_value, ensure_proper_subpath,
    write_layout, read_layout, get_column_subset, get_columns

@testset "utilities" begin
    @test fixed_Tuple_types(Tuple{Int, Int}) ≡ (Int, Int)
    @test fixed_Tuple_types(Tuple{Int, Float64, Char}) ≡ (Int, Float64, Char)
    @test_throws ArgumentError fixed_Tuple_types(Tuple{Int, Vararg{Int}})

    @test representative_value(Int) isa Int
    @test representative_value(Float64) isa Float64
    @test representative_value(Date) isa Date
    @test representative_value(Tuple{Date,Int}) isa Tuple{Date,Int}
    @test_throws ArgumentError representative_value(Vector{Int})
end

@testset "meta subpath calculation" begin
    @test ensure_proper_subpath("/tmp", "foo") == "/tmp/foo"
    @test ensure_proper_subpath("/tmp/dir/", "foo") == "/tmp/dir/foo"
    @test_throws ArgumentError ensure_proper_subpath("/tmp/dir/", "../foo")
    @test ensure_proper_subpath("/tmp/", "dir/../foo") == "/tmp/foo"
    @test_throws ArgumentError ensure_proper_subpath("/tmp/", "dir/../../foo")
    @test_throws ArgumentError ensure_proper_subpath("/tmp", "/root")
end

@testset "layout information" begin
    dir = mktempdir()
    N = rand(1:10_000)
    S = Tuple{Date,Int}
    write_layout(dir, N, S)
    @test meta_path(dir, ".") == joinpath(dir, "meta/")
    @test meta_path(dir, "foo.jld2") == joinpath(dir, "meta/foo.jld2")
    @test isdir(meta_path(dir, ".")) # test that meta was created
    @test read_layout(dir) ≡ (N, S)
end

@testset "meta path" begin
    @test meta_path("/tmp/", "test") == "/tmp/meta/test"
    @test_throws ArgumentError meta_path("/tmp/", "../foo")
end

@testset "write values, get back as mmapped" begin
    dir = mktempdir()

    # write
    sink = SinkColumns(dir, Tuple{Int, Float64})
    for i in 1:9
        push!(sink, (i, Float64(i)))
    end
    push!(sink, (10, 10))       # test conversion
    @test length(sink) == 10
    @test eltype(sink) == Tuple{Int, Float64}
    flush(sink)     # NOTE calling both `flush` and `close` is not a strong test
    close(sink)

    # test meta path
    @test meta_path(sink, "foo") == meta_path(dir, "foo")

    # append
    sink = SinkColumns(dir, true)
    for i in 11:15
        push!(sink, (i, Float64(i)))
    end
    @test length(sink) == 15
    @test eltype(sink) == Tuple{Int, Float64}
    close(sink)

    # mmap
    cols = MmappedColumns(dir)
    @test eltype(cols) == Tuple{Int, Float64}
    @test length(cols) == 15
    @test cols[3] ≡ (3, 3.0)                                        # getindex Int
    @test cols[(end-2):end] == [(13, 13.0), (14, 14.0), (15, 15.0)] # getindex to_indices
    @test cols == [(i, Float64(i)) for i in 1:15]
end

@testset "mmap standalone, opened multiple times" begin
    dir = mktempdir()
    N = 39
    # create
    cols = MmappedColumns(dir, N, Tuple{Int}) # create
    col = cols.columns[1]
    col .= randperm(N)          # random permutation
    Mmap.sync!(cols)
    # reopen and sort
    cols = MmappedColumns(dir)
    sort!(cols.columns[1])
    Mmap.sync!(cols)
    # reopen and test
    cols = MmappedColumns(dir)
    @test cols == [(i,) for i in 1:N]
    # test meta path
    @test meta_path(cols, "foo") == meta_path(dir, "foo")
end

@testset "mmap getindex and setindex" begin
    dir = mktempdir()
    N = 10
    # setindex! — fill with values
    cols = MmappedColumns(dir, N, Tuple{Int, Char})
    for i in 1:N
        cols[i] = i, Char(i + 'a')
    end
    Mmap.sync!(cols)
    # getindex — reopen and check
    for i in 1:N
        @test cols[i] ≡ (i, Char(i + 'a'))
    end
    A = [(N+i, Char(i + 'a')) for i in 2:7]
    cols[3:8] = A
    @test cols[3:8] == A
end

@testset "path creation tests" begin
    dir = tempname()
    @test !isdir(dir)           # verify that it does not exist
    cols = MmappedColumns(dir, 10, Tuple{Float64, Date})
    # test that directories are created
    @test isdir(dir)
    @test isdir(joinpath(dir, "meta"))
    # test layout and data files
    @test isfile(joinpath(dir, "layout.jld2"))
    @test isfile(joinpath(dir, "1.bin"))
    @test isfile(joinpath(dir, "2.bin"))
    # test that no other files are created
    @test sort(readdir(dir)) == ["1.bin", "2.bin", "layout.jld2", "meta"]
    @test readdir(joinpath(dir, "meta")) == String[]
end

@testset "mmapped column accessors" begin
    dir = mktempdir()
    N = 10
    cols = MmappedColumns(dir, N, Tuple{Int, Float64, Float32})
    for i in 1:N
        cols[i] = i, Float64(i), Float32(i)
    end
    Mmap.sync!(cols)
    # subset
    c12 = get_column_subset(cols, 1:2)
    @test eltype(c12) ≡ Tuple{Int, Float64}
    for i in 1:N
        c12[i] ≡ (i, Float64(i))
    end
    # invalid subset
    @test_throws ArgumentError get_column_subset(cols, 1)
    # columns
    cc = get_columns(cols)[2:3]
    @test eltype.(cc) == (Float64, Float32)
    @test cc == (Float64.(1:N), Float32.(1:N))
    @test cols.columns == get_columns(cols)
end


text:
"""
    Action

Defines what is drawn in a defined frame range.

# Fields
- `frames::Frames`: A range of frames for which the `Action` is called
- `id::Union{Nothing, Symbol}`: An id which can be used to save the result of `func`
- `func::Function`: The drawing function which draws something on the canvas.
    It gets called with the arguments `video, action, frame`
- `anim::Animation`: defines the interpolation function for the transitions
- `transitions::Vector{Transition}` a list of transitions
    which can be performed before the function gets called.
- `internal_transitions::Vector{InternalTransition}`:
    Similar to `transitions` but holds the concrete information whereas `Transition` can
    hold links to other actions which need to be computed first.
    See [`compute_transition!`](@ref)
- `opts::Any` can hold any options defined by the user
"""
mutable struct Action <: AbstractAction
    frames::Frames
    id::Union{Nothing,Symbol}
    func::Function
    anim::Animation
    transitions::Vector{Transition}
    internal_transitions::Vector{InternalTransition}
    subactions::Vector{SubAction}
    current_setting::ActionSetting
    opts::Dict{Symbol,Any}
end

"""
    CURRENT_ACTION

holds the current action in an array to be declared as a constant
The current action can be accessed using CURRENT_ACTION[1]
"""
const CURRENT_ACTION = Array{Action,1}()

"""
    Action(frames, func::Function, args...)

The most simple form of an action (if there are no `args`/`kwargs`) just calls
`func(video, action, frame)` for each of the frames it is defined for.
`args` are defined it the next function definition and can be seen in action
    in this example [`javis`](@ref)
"""
Action(frames, func::Function, args...; kwargs...) =
    Action(frames, nothing, func, args...; kwargs...)

"""
    Action(frames_or_id::Symbol, func::Function, args...)

This function decides whether you wrote `Action(frames_symbol, ...)`,
    or `Action(id_symbol, ...)`
If the symbol `frames_or_id` is not a `FRAMES_SYMBOL` then it is used as an id_symbol.
"""
function Action(frames_or_id::Symbol, func::Function, args...; kwargs...)
    if frames_or_id in FRAMES_SYMBOL
        Action(frames_or_id, nothing, func, args...; kwargs...)
    else
        Action(:same, frames_or_id, func, args...; kwargs...)
    end
end

"""
    Action(func::Function, args...)

Similar to the above but uses the same frames as the action above.
"""
Action(func::Function, args...; kwargs...) =
    Action(:same, nothing, func, args...; kwargs...)

"""
    Action(frames, id::Union{Nothing,Symbol}, func::Function,
           transitions::Transition...; kwargs...)

Fallback constructor for an Action which doesn't define an animation.
A linear animation is assumed.
"""
function Action(
    frames,
    id::Union{Nothing,Symbol},
    func::Function,
    transitions::Transition...;
    kwargs...,
)

    Action(frames, id, func, easing_to_animation(linear()), transitions...; kwargs...)
end

"""
    Action(frames, id::Union{Nothing,Symbol}, func::Function, easing::Union{ReversedEasing, Easing},
           args...; kwargs...)

Fallback constructor for an Action which does define an animation using an easing function.

# Example
```
javis(
    video, [
        BackgroundAction(1:100, ground),
        Action((args...)->t(), sineio(), Translation(250, 0))
    ]
)
```
"""
function Action(
    frames,
    id::Union{Nothing,Symbol},
    func::Function,
    easing::Union{ReversedEasing,Easing},
    args...;
    kwargs...,
)

    Action(frames, id, func, easing_to_animation(easing), args...; kwargs...)
end

"""
    Action(frames, id::Union{Nothing,Symbol}, func::Function,
           transitions::Transition...; kwargs...)

# Arguments
- `frames`: defines for which frames this action is called
- `id::Symbol`: Is used if the `func` returns something which
    shall be accessible by other actions later
- `func::Function` the function that is called after the `transitions` are performed
- `transitions::Transition...` a list of transitions that are performed before
    the function `func` itself is called

The keywords arguments will be saved inside `.opts` as a `Dict{Symbol, Any}`
"""
function Action(
    frames,
    id::Union{Nothing,Symbol},
    func::Function,
    anim::Animation,
    transitions::Transition...;
    kwargs...,
)
    if isempty(CURRENT_VIDEO)
        throw(ErrorException("A `Video` must be defined before an `Action`"))
    end
    CURRENT_VIDEO[1].defs[:last_frames] = frames
    opts = Dict(kwargs...)
    subactions = SubAction[]
    if haskey(opts, :subactions)
        subactions = opts[:subactions]
        delete!(opts, :subactions)
    end
    Action(
        frames,
        id,
        func,
        anim,
        collect(transitions),
        [],
        subactions,
        ActionSetting(),
        opts,
    )
end

"""
    BackgroundAction(frames, func::Function, args...; kwargs...)

Create an Action where `in_global_layer` is set to true such that
i.e the specified color in the background is applied globally (basically a new default)
"""
function BackgroundAction(frames, func::Function, args...; kwargs...)
    Action(frames, nothing, func, args...; in_global_layer = true, kwargs...)
end

"""
    BackgroundAction(frames, id::Symbol, func::Function, args...; kwargs...)

Create an Action where `in_global_layer` is set to true and saves the return into `id`.
"""
function BackgroundAction(frames, id::Symbol, func::Function, args...; kwargs...)
    Action(frames, id, func, args...; in_global_layer = true, kwargs...)
end


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

JLLWrappers.@generate_main_file_header("Qt5QuickControls2")
JLLWrappers.@generate_main_file("Qt5QuickControls2", UUID("bf3ac11c-603e-589e-b4b7-e696ac65aa4a"))
end  # module Qt5QuickControls2_jll


text:
# Use
#    @warnpcfail precompile(args...)
# if you want to be warned when a precompile directive fails
macro warnpcfail(ex::Expr)
    modl = __module__
    file = __source__.file === nothing ? "?" : String(__source__.file)
    line = __source__.line
    quote
        $(esc(ex)) || @warn """precompile directive
     $($(Expr(:quote, ex)))
 failed. Please report an issue in $($modl) (after checking for duplicates) or remove this directive.""" _file=$file _line=$line
    end
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:name, :study_id),Tuple{Nothing,Nothing}},Type{CTScanner{FBPInfo{G,A} where A<:AbstractFBP,D} where D where G},FanBeamGeometry{Float32,DefaultTomograph},CTImage{Float32,Array{Float32,2}}})
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:name, :study_id),Tuple{Nothing,Nothing}},Type{CTScanner{FBPInfo{G,A} where A<:AbstractFBP,D} where D where G},FanBeamGeometry{Float64,DefaultTomograph},CTImage{Float64,Array{Float64,2}}})
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:name, :study_id),Tuple{Nothing,Nothing}},Type{CTScanner{FBPInfo{G,A} where A<:AbstractFBP,D} where D where G},ParallelBeamGeometry{Float32,DefaultTomograph},CTImage{Float32,Array{Float32,2}}})
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:name, :study_id),Tuple{Nothing,Nothing}},Type{CTScanner{FBPInfo{G,A} where A<:AbstractFBP,D} where D where G},ParallelBeamGeometry{Float64,DefaultTomograph},CTImage{Float64,Array{Float64,2}}})
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:name, :study_id),Tuple{Symbol,String}},Type{CTScanner},FBPInfo{ParallelBeamGeometry{Float32,DefaultTomograph},FBP{MartaCT.RadonAlgorithm.Filters.RamLak}},MartaCT.CTScan.CTImageData.CTData{Float32,Array{Float32,2}}})
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:name, :study_id),Tuple{Symbol,String}},Type{CTScanner},FBPInfo{ParallelBeamGeometry{Float64,DefaultTomograph},FBP{MartaCT.RadonAlgorithm.Filters.RamLak}},MartaCT.CTScan.CTImageData.CTData{Float64,Array{Float64,2}}})
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:nϕ,),Tuple{Int64}},Type{FanBeamGeometry},CircleImage{Float32}})
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:nϕ,),Tuple{Int64}},Type{FanBeamGeometry},CircleImage{Float64}})
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:nϕ,),Tuple{Int64}},Type{ParallelBeamGeometry},CircleImage{Float32}})
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:nϕ,),Tuple{Int64}},Type{ParallelBeamGeometry},CircleImage{Float64}})
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:width,),Tuple{Int64}},Type{CircleImage},Type{Float32}})
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:width,),Tuple{Int64}},Type{CircleImage},Type{Float64}})
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:width,),Tuple{Int64}},Type{GrayScaleLine},Type{Float32}})
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:width,),Tuple{Int64}},Type{GrayScaleLine},Type{Float64}})
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:width,),Tuple{Int64}},Type{GrayScalePyramid},Type{Float32}})
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:width,),Tuple{Int64}},Type{GrayScalePyramid},Type{Float64}})
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:width,),Tuple{Int64}},Type{WhiteRect},Type{Float32}})
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:width,),Tuple{Int64}},Type{WhiteRect},Type{Float64}})
    Base.precompile(Tuple{Core.kwftype(typeof(circle_image)),NamedTuple{(:radius, :calibration_value, :background, :rows),Tuple{Int64,Float32,Float32,Int64}},typeof(circle_image),Type{Float32}})
    Base.precompile(Tuple{Core.kwftype(typeof(circle_image)),NamedTuple{(:radius, :calibration_value, :background, :rows),Tuple{Int64,Float64,Float64,Int64}},typeof(circle_image),Type{Float64}})
    Base.precompile(Tuple{Core.kwftype(typeof(gray_scale_image)),NamedTuple{(:swidth, :sheight, :gray_scale, :background),Tuple{Int64,Int64,IntervalSets.Interval{:closed,:closed,Float32},Float32}},typeof(gray_scale_image),Type{Float32}})
    Base.precompile(Tuple{Core.kwftype(typeof(gray_scale_image)),NamedTuple{(:swidth, :sheight, :gray_scale, :background),Tuple{Int64,Int64,IntervalSets.Interval{:closed,:closed,Float64},Float64}},typeof(gray_scale_image),Type{Float64}})
    Base.precompile(Tuple{Type{CTScanner{FBPInfo{G,A} where A<:AbstractFBP,D} where D where G},CircleImage{Float32}})
    Base.precompile(Tuple{Type{CTScanner{FBPInfo{G,A} where A<:AbstractFBP,D} where D where G},CircleImage{Float64}})
    Base.precompile(Tuple{Type{CTScanner{FBPInfo{G,A} where A<:AbstractFBP,D} where D where G},FanBeamGeometry{Float64,DefaultTomograph},CircleImage{Float64}})
    Base.precompile(Tuple{Type{CTScanner{FBPInfo{G,A} where A<:AbstractFBP,D} where D where G},GrayScalePyramid{Float32}})
    Base.precompile(Tuple{Type{CTScanner{FBPInfo{G,A} where A<:AbstractFBP,D} where D where G},GrayScalePyramid{Float64}})
    Base.precompile(Tuple{Type{CTScanner{FBPInfo{G,A} where A<:AbstractFBP,D} where D where G},WhiteRect{Float32}})
    Base.precompile(Tuple{typeof(calibrate_tomogram),CTScanner{FBPInfo{ParallelBeamGeometry{Float32,DefaultTomograph},FBP{MartaCT.RadonAlgorithm.Filters.RamLak}},MartaCT.CTScan.CTImageData.CTData{Float32,Array{Float32,2}}},CircleImage{Float32}})
    Base.precompile(Tuple{typeof(calibrate_tomogram),CTScanner{FBPInfo{ParallelBeamGeometry{Float32,DefaultTomograph},FBP{MartaCT.RadonAlgorithm.Filters.RamLak}},MartaCT.CTScan.CTImageData.CTData{Float32,Array{Float32,2}}},CircleParams{Float32}})
    Base.precompile(Tuple{typeof(calibrate_tomogram),CTScanner{FBPInfo{ParallelBeamGeometry{Float32,DefaultTomograph},FBP{MartaCT.RadonAlgorithm.Filters.RamLak}},MartaCT.CTScan.CTImageData.CTData{Float32,Array{Float32,2}}},ImageParams{Float32}})
    Base.precompile(Tuple{typeof(calibrate_tomogram),CTScanner{FBPInfo{ParallelBeamGeometry{Float64,DefaultTomograph},FBP{MartaCT.RadonAlgorithm.Filters.RamLak}},MartaCT.CTScan.CTImageData.CTData{Float64,Array{Float64,2}}},CircleImage{Float64}})
    Base.precompile(Tuple{typeof(calibrate_tomogram),CTScanner{FBPInfo{ParallelBeamGeometry{Float64,DefaultTomograph},FBP{MartaCT.RadonAlgorithm.Filters.RamLak}},MartaCT.CTScan.CTImageData.CTData{Float64,Array{Float64,2}}},CircleParams{Float64}})
    Base.precompile(Tuple{typeof(calibrate_tomogram),CTScanner{FBPInfo{ParallelBeamGeometry{Float64,DefaultTomograph},FBP{MartaCT.RadonAlgorithm.Filters.RamLak}},MartaCT.CTScan.CTImageData.CTData{Float64,Array{Float64,2}}},ImageParams{Float64}})
    Base.precompile(Tuple{typeof(project_image),Some{CTImage{Float32,Array{Float32,2}}},RadonInfo{FanBeamGeometry{Float32,DefaultTomograph},Radon}})
    Base.precompile(Tuple{typeof(project_image),Some{CTImage{Float32,Array{Float32,2}}},RadonInfo{ParallelBeamGeometry{Float32,DefaultTomograph},Radon}})
    Base.precompile(Tuple{typeof(project_image),Some{CTImage{Float64,Array{Float64,2}}},RadonInfo{FanBeamGeometry{Float64,DefaultTomograph},Radon}})
    Base.precompile(Tuple{typeof(project_image),Some{CTImage{Float64,Array{Float64,2}}},RadonInfo{ParallelBeamGeometry{Float64,DefaultTomograph},Radon}})
end


text:
module API

import LLVM.API: LLVMValueRef, LLVMModuleRef, LLVMTypeRef, LLVMContextRef
using Enzyme_jll
using Libdl
using CEnum

struct EnzymeAAResultsRef
    a::Ptr{Cvoid}
    b::Ptr{Cvoid}
    c::Ptr{Cvoid}
end
const EnzymeTypeAnalysisRef = Ptr{Cvoid}
const EnzymeAugmentedReturnPtr = Ptr{Cvoid}

struct IntList
    data::Ptr{Int64}
    size::Csize_t
end
IntList() = IntList(Ptr{Int64}(0),0)

@cenum(CConcreteType,
  DT_Anything = 0,
  DT_Integer = 1,
  DT_Pointer = 2,
  DT_Half = 3,
  DT_Float = 4,
  DT_Double = 5,
  DT_Unknown = 6
)

struct EnzymeTypeTree end
const CTypeTreeRef = Ptr{EnzymeTypeTree}

EnzymeNewTypeTree() = ccall((:EnzymeNewTypeTree, libEnzyme), CTypeTreeRef, ())
EnzymeNewTypeTreeCT(T, ctx) = ccall((:EnzymeNewTypeTreeCT, libEnzyme), CTypeTreeRef, (CConcreteType, LLVMContextRef), T, ctx)
EnzymeNewTypeTreeTR(tt) = ccall((:EnzymeNewTypeTreeTR, libEnzyme), CTypeTreeRef, (CTypeTreeRef,), tt)

EnzymeFreeTypeTree(tt) = ccall((:EnzymeFreeTypeTree, libEnzyme), Cvoid, (CTypeTreeRef,), tt)
EnzymeSetTypeTree(dst, src) = ccall((:EnzymeSetTypeTree, libEnzyme), Cvoid, (CTypeTreeRef, CTypeTreeRef), dst, src)
EnzymeMergeTypeTree(dst, src) = ccall((:EnzymeMergeTypeTree, libEnzyme), Cvoid, (CTypeTreeRef, CTypeTreeRef), dst, src)
EnzymeTypeTreeOnlyEq(dst, x) = ccall((:EnzymeTypeTreeOnlyEq, libEnzyme), Cvoid, (CTypeTreeRef, Int64), dst, x)
EnzymeTypeTreeShiftIndiciesEq(dst, dl, offset, maxSize, addOffset) =
    ccall((:EnzymeTypeTreeShiftIndiciesEq, libEnzyme), Cvoid, (CTypeTreeRef, Cstring, Int64, Int64, UInt64),
        dst, dl, offset, maxSize, addOffset)

struct CFnTypeInfo
    arguments::Ptr{CTypeTreeRef}
    ret::CTypeTreeRef

    known_values::Ptr{IntList}
end

@cenum(CDIFFE_TYPE,
  DFT_OUT_DIFF = 0,  # add differential to an output struct
  DFT_DUP_ARG = 1,   # duplicate the argument and store differential inside
  DFT_CONSTANT = 2,  # no differential
  DFT_DUP_NONEED = 3 # duplicate this argument and store differential inside,
                     # but don't need the forward
)

function EnzymeGetGlobalAA(mod)
    ccall((:EnzymeGetGlobalAA, libEnzyme), EnzymeAAResultsRef, (LLVMModuleRef,), mod)
end

function EnzymeFreeGlobalAA(aa)
    ccall((:EnzymeFreeGlobalAA, libEnzyme), Cvoid, (EnzymeAAResultsRef,), aa)
end

# Create the derivative function itself.
#  \p todiff is the function to differentiate
#  \p retType is the activity info of the return
#  \p constant_args is the activity info of the arguments
#  \p returnValue is whether the primal's return should also be returned
#  \p dretUsed is whether the shadow return value should also be returned
#  \p additionalArg is the type (or null) of an additional type in the signature
#  to hold the tape.
#  \p typeInfo is the type info information about the calling context
#  \p _uncacheable_args marks whether an argument may be rewritten before loads in
#  the generated function (and thus cannot be cached).
#  \p augmented is the data structure created by prior call to an augmented forward
#  pass
#  \p AtomicAdd is whether to perform all adjoint updates to memory in an atomic way
#  \p PostOpt is whether to perform basic optimization of the function after synthesis
function EnzymeCreatePrimalAndGradient(todiff, retType, constant_args, TA, global_AA,
                                       returnValue, dretUsed, topLevel, additionalArg, typeInfo,
                                       uncacheable_args, augmented, atomicAdd, postOpt)
    ccall((:EnzymeCreatePrimalAndGradient, libEnzyme), LLVMValueRef,
        (LLVMValueRef, CDIFFE_TYPE, Ptr{CDIFFE_TYPE}, Csize_t, EnzymeTypeAnalysisRef,
         EnzymeAAResultsRef, UInt8, UInt8, UInt8, LLVMTypeRef, CFnTypeInfo,
         Ptr{UInt8}, Csize_t, EnzymeAugmentedReturnPtr, UInt8, UInt8),
        todiff, retType, constant_args, length(constant_args), TA, global_AA, returnValue,
        dretUsed, topLevel, additionalArg, typeInfo, uncacheable_args, length(uncacheable_args),
        augmented, atomicAdd, postOpt)
end

# Create an augmented forward pass.
#  \p todiff is the function to differentiate
#  \p retType is the activity info of the return
#  \p constant_args is the activity info of the arguments
#  \p returnUsed is whether the primal's return should also be returned
#  \p typeInfo is the type info information about the calling context
#  \p _uncacheable_args marks whether an argument may be rewritten before loads in
#  the generated function (and thus cannot be cached).
#  \p forceAnonymousTape forces the tape to be an i8* rather than the true tape structure
#  \p AtomicAdd is whether to perform all adjoint updates to memory in an atomic way
#  \p PostOpt is whether to perform basic optimization of the function after synthesis
function EnzymeCreateAugmentedPrimal(todiff, retType, constant_args, TA, global_AA, returnUsed,
                                     typeInfo, uncacheable_args, forceAnonymousTape, atomicAdd, postOpt)
    ccall((:EnzymeCreateAugmentedPrimal, libEnzyme), EnzymeAugmentedReturnPtr,
        (LLVMValueRef, CDIFFE_TYPE, Ptr{CDIFFE_TYPE}, Csize_t,
         EnzymeTypeAnalysisRef, EnzymeAAResultsRef, UInt8,
         CFnTypeInfo, Ptr{UInt8}, Csize_t, UInt8, UInt8, UInt8),
        todiff, retType, constant_args, length(constant_args), TA, global_AA, returnUsed,
        typeInfo, uncacheable_args, length(uncacheable_args), forceAnonymousTape, atomicAdd, postOpt)
end

# typedef bool (*CustomRuleType)(int /*direction*/, CTypeTree * /*return*/,
#                                CTypeTree * /*args*/, size_t /*numArgs*/,
#                                LLVMValueRef);
const CustomRuleType = Ptr{Cvoid}

function CreateTypeAnalysis(triple, rulenames, rules)
    @assert length(rulenames) == length(rules)
    ccall((:CreateTypeAnalysis, libEnzyme), EnzymeTypeAnalysisRef, (Cstring, Ptr{Cstring}, Ptr{CustomRuleType}, Csize_t), triple, rulenames, rules, length(rules))
end

function FreeTypeAnalysis(ta)
    ccall((:FreeTypeAnalysis, libEnzyme), Cvoid, (EnzymeAAResultsRef,), ta)
end

function EnzymeExtractReturnInfo(ret, data, existed)
    @assert length(data) == length(existed)
    ccall((:EnzymeExtractReturnInfo, libEnzyme),
           Cvoid, (EnzymeAugmentedReturnPtr, Ptr{Int64}, Ptr{UInt8}, Csize_t),
           ret, data, existed, length(data))
end

function EnzymeExtractFunctionFromAugmentation(ret)
    ccall((:EnzymeExtractFunctionFromAugmentation, libEnzyme), LLVMValueRef, (EnzymeAugmentedReturnPtr,), ret)
end

function EnzymeExtractTapeTypeFromAugmentation(ret)
    ccall((:EnzymeExtractTapeTypeFromAugmentation, libEnzyme), LLVMTypeRef, (EnzymeAugmentedReturnPtr,), ret)
end

end


text:
export
  SparseDenseMatrix,
  SparseDenseVector,
  Vertex,
  VertexSet,
  vertexsetsize,
  subspace,
  vlist,
  length

import Base: ==, hash, getindex, length

macro argumentcheck(ex, msgs...)
    msg = isempty(msgs) ? ex : msgs[1]
    if isa(msg, AbstractString)
        msg = msg # pass-through
    elseif !isempty(msgs) && (isa(msg, Expr) || isa(msg, Symbol))
        #message is an expression needing evaluating
        msg = :(Main.Base.string($(esc(msg))))
    #elseif isdefined(Main, :Base) && isdefined(Main.Base, :string) && applicable(Main.Base.string, msg)
        #msg = Main.Base.string(msg)
    #else
        # string() might not be defined during bootstrap
        #msg = :(Main.Base.string($(Expr(:quote, msg))))
    end
    return :($(esc(ex)) ? $(nothing) : throw(Main.Base.ArgumentError($msg)))
end

"""
    type Vertex

Type consisting of list of `Int`, describing the labels of vectors from the
canonical basis corresponding to the `Vertex`. To get the vector label one can
use `Vertex()` function, or `Vertex[i]` for a unique label.

See [1] for the more information and usage exmaples.

[1] K. Domino, A. Glos, M. Ostaszewski, Superdiffusive quantum stochastic walk
definable on arbitrary directed graph, Quantum Information & Computation,
Vol.17 No.11&12, pp. 0973-0986, arXiv:1701.04624.
"""
struct Vertex
  subspace::Vector{Int}

  Vertex(v::Vector{Int}) = all(v.>0) ? new(v) : throw(ArgumentError("Vector should consist of positive elments"))
end

"""
    subspace(v)

Returns the subspace connected to vertex `v`.

```jldoctest; setup = :(using QSWalk)
julia> v = Vertex([1,2,3])
Vertex([1, 2, 3])

julia> subspace(v)
3-element Array{Int64,1}:
 1
 2
 3
```
"""
subspace(v::Vertex) = v.subspace

==(v::Vertex, w::Vertex) = subspace(v) == subspace(w)
hash(v::Vertex) = hash(subspace(v))

==(v::Tuple{Vertex, Vertex}, w::Tuple{Vertex, Vertex}) = [subspace(v[1]), subspace(v[2])] == [subspace(w[1]), subspace(w[2])]
hash(v::Tuple{Vertex, Vertex}) = hash([subspace(v[1]), subspace(v[2])])
getindex(v::Vertex, i::Int) = v.subspace[i]

length(v::Vertex) = length(subspace(v))

function checkvertexset(partition::Vector{Vector{Int}})
  joined = Int[]
  for lin = partition
    append!(joined, lin)
  end
  length(joined) == length(Set(joined))
end

checkvertexset(vset::Vector{Vertex}) = checkvertexset([subspace(v) for v = vset])

"""
    type VertexSet

Type consisting of a list of `Vertex` objects. It describes the partition of the
linear subspace. Object of this type should be constructed using
`make_vertex_set` or by `nm_lind` functions. In order to get a
list of the vertices from an object of type `vertexset`, one should use
`vlist` function, or, for a specific `Vertex`, an getindex function
`vertexset[i]`.
"""
struct VertexSet
  vertices::Vector{Vertex}

  VertexSet(vset::Vector{Vertex}) = checkvertexset(vset) ? new(vset) : throw(ArgumentError("Vertices should correspond to orthogonal linear spaces"))
end

VertexSet(vset::Vector{Vector{Int}}) = VertexSet([Vertex(v) for v = vset])

"""
    vlist(vset)

Returns the list of vertices for given `vset` of type `VertexSet`.

```jldoctest; setup = :(using QSWalk)
julia> vset = VertexSet([[1], [2, 3, 4], [5, 6, 7], [8, 9]])
VertexSet(Vertex[Vertex([1]), Vertex([2, 3, 4]), Vertex([5, 6, 7]), Vertex([8, 9])])

julia> vlist(vset)
4-element Array{Vertex,1}:
 Vertex([1])
 Vertex([2, 3, 4])
 Vertex([5, 6, 7])
 Vertex([8, 9])
```
"""
vlist(vset::VertexSet) = vset.vertices

==(v::VertexSet, w::VertexSet) = v.vertices == w.vertices

getindex(vset::VertexSet, i::Int) = vlist(vset)[i]
getindex(vset::VertexSet, veci::Vector{Int}) = vlist(vset)[veci]
length(vset::VertexSet) = length(vlist(vset))

"""
    vertexsetsize(vertexset)

Return the dimension of the linearspace corresponding to given `vertexset`.

# Examples

```jldoctest; setup = :(using QSWalk)
julia> vertexsetsize(VertexSet([[1, 2, 3], [4, 5]]))
5
```
"""
function vertexsetsize(vset::VertexSet)
  sum(length.(vlist(vset)))
end

"""

    incidence_list(A[; epsilon])

Return list of indices of non-zero elements of the matrix `A`. The `i`-th element of the
result list is the vector of indices `j` for which `abs(A[j, i]) >=  epsilon`.

# Examples

```jldoctest; setup = :(using QSWalk)
julia> A = [1 2 3; 0 3. 4.; 0 0 5.]
3×3 Array{Float64,2}:
 1.0  2.0  3.0
 0.0  3.0  4.0
 0.0  0.0  5.0

julia> QSWalk.incidence_list(A)
3-element Array{Array{Int64,1},1}:
 [1]
 [1, 2]
 [1, 2, 3]

julia> QSWalk.incidence_list(A, epsilon = 2.5)
3-element Array{Array{Int64,1},1}:
 []
 [2]
 [1, 2, 3]

```
"""
function incidence_list(A::AbstractMatrix{<:Number};
                        epsilon::Real = eps())
  @argumentcheck epsilon >= 0 "epsilon needs to be nonnegative"
  @argumentcheck size(A, 1) == size(A, 2) "A matrix must be square"
  [findall(x -> abs(x)>= epsilon, A[:, i]) for i = 1:size(A, 1)]
end

"""

    reversed_incidence_list(A[; epsilon])

For given matrix `A` the function returns list of indices. The `i`-th element of
result list is the vector of indices for which `abs(A[i, j]) >=  epsilon`.

# Examples

```jldoctest; setup = :(using QSWalk)
julia> A = [1 2 3; 0 3. 4.; 0 0 5.]
3×3 Array{Float64,2}:
 1.0  2.0  3.0
 0.0  3.0  4.0
 0.0  0.0  5.0

julia> QSWalk.reversed_incidence_list(A)
3-element Array{Array{Int64,1},1}:
 [1, 2, 3]
 [2, 3]
 [3]

julia> QSWalk.reversed_incidence_list(A, epsilon = 2.5)
3-element Array{Array{Int64,1},1}:
 [3]
 [2, 3]
 [3]
```
"""
function reversed_incidence_list(A::AbstractMatrix; epsilon::Real = eps())
  incidence_list(transpose(A), epsilon = epsilon)
end

"""

    revinc_to_vertexset(revincidence_list)

Return `vertexset` of type `VertexSet` corresponding to given list of indices.
The function map to consecutive element list orthogonal subspaces. The
dimensions of the subspaces equal to size of each element of revincidence_list.
The exception is empty list, for which onedimensional subspace is attached.

# Examples

```jldoctest; setup = :(using QSWalk)
julia> vset = [Int64[], [2], [1, 2, 3]]
3-element Array{Array{Int64,1},1}:
 []
 [2]
 [1, 2, 3]

julia> vlist(QSWalk.revinc_to_vertexset(vset))
3-element Array{Vertex,1}:
 Vertex([1])
 Vertex([2])
 Vertex([3, 4, 5])

```
"""
function revinc_to_vertexset(revincidence_list::Vector{Vector{Int}})
  vertexset = Vector{Int}[]
  start = 1
  for i = revincidence_list
    if length(i)!= 0
      push!(vertexset, collect(start:(start+length(i)-1)))
      start+= length(i)
    else
      push!(vertexset, [start])
      start +=  1
    end
  end
  VertexSet(vertexset)
end

"""

    make_vertex_set(A[, epsilon])

Creates object of type `VertexSet` which represents how vertices are located in
matrix. Should be used only in the nondefault use of `evolve_generator` function.
It is always equal to the second element if output of `evolve_generator` function.

# Examples

```jldoctest; setup = :(using QSWalk)
julia> A = [1 2 3; 0 3. 4.; 0 0 5.]
3×3 Array{Float64,2}:
 1.0  2.0  3.0
 0.0  3.0  4.0
 0.0  0.0  5.0

julia> vlist(make_vertex_set(A))
3-element Array{Vertex,1}:
 Vertex([1, 2, 3])
 Vertex([4, 5])
 Vertex([6])

julia> vlist(make_vertex_set(A, epsilon = 2.5))
3-element Array{Vertex,1}:
 Vertex([1])
 Vertex([2, 3])
 Vertex([4])
```
"""
function make_vertex_set(A::AbstractMatrix; epsilon::Real = eps())
  @argumentcheck epsilon >= 0 "epsilon needs to be nonnegative"
  @argumentcheck size(A, 1) == size(A, 2) "A matrix must be square"
  revinc_to_vertexset(reversed_incidence_list(A, epsilon = epsilon))
end


text:
# -------------------------------------------------------------------------------------------------
# Interagency technical report "Global and Regional Sea Level Rise Scenarios for the United States:
# Updated Mean Projections and Extreme Water Level Probabilities Along U.S. Coastlines"
#
# This script calls all the routines to compute and analyze the projections, trajectories and
# observations. All results are stored in netCDF files. This script also stores the results for
# some of the figures in the report as text files, which can be read by the GMT scripts used
# to make these figures.
# -------------------------------------------------------------------------------------------------

dir_code = homedir()*"/Projects/2021_NOAA_TR/Code/"
include(dir_code*"ConvertNCA5ToGrid.jl")
include(dir_code*"Masks.jl")
include(dir_code*"GlobalProjections.jl")
include(dir_code*"RegionalProjections.jl")
include(dir_code*"LocalProjections.jl")
include(dir_code*"GriddedProjections.jl")
include(dir_code*"SaveFigureData.jl")

function main()
    # Define some projection settings
    settings = DefSettings()

    # Routines to do some preprocessing. Not needed to run: output has been provided
    # Currently commented out:
    # Masks.CreateMask(settings)                       # Create a region mask for each of the regions
    ConvertNCA5ToGrid.RunConvertNCA5ToGrid(settings) # Covert the NCA5 projections to single GMSL, tide-gauge and gridded projections

    # Create the projections
    GlobalProjections.RunGlobalProjections(settings)
    RegionalProjections.RunRegionalProjections(settings)
    LocalProjections.RunLocalProjections(settings)
    GriddedProjections.RunGriddedProjections(settings)

    # Save the data for the report figures
    SaveFigureData.RunSaveFigureData(settings)
    return nothing
end

function DefSettings()
    settings = Dict()
    settings["dir_project"] = homedir()*"/Projects/2021_NOAA_TR/"

    # File names
    settings["fn_region_mask"] = settings["dir_project"] * "Data/region_mask.nc"
    settings["fn_basin_codes"] = settings["dir_project"] * "Data/basin_codes.nc"
    settings["fn_bathymetry"]  = settings["dir_project"] * "Data/GEBCO_bathymetry_05.nc"

    # Observed sea level
    settings["fn_tg_data"]             = settings["dir_project"] * "Data/US_tg_monthly.xlsx"
    settings["fn_gmsl_20c_ensembles"]  = settings["dir_project"] * "Data/GMSL_ensembles.nc"
    settings["fn_gmsl_20c_mean"]       = settings["dir_project"] * "Data/global_timeseries_measures.nc"
    settings["fn_gmsl_altimetry_GSFC"] = settings["dir_project"] * "Data/GMSL_TPJAOS_5.0_199209_202106.txt"
    settings["fn_GIA"]                 = settings["dir_project"] * "Data/GIA_Caron_stats_05.nc"
    settings["fn_altimetry"]           = settings["dir_project"] * "Data/CDS_monthly_1993_2020.nc"
    settings["fn_GRD"]                 = settings["dir_project"] * "Data/grd_1992_2020.nc"

    # Climate indices
    settings["fn_gmsl_ENSO"] = settings["dir_project"]*"Data/ClimIdx/enso_correction.mat"
    settings["fn_NAO"] = settings["dir_project"]*"Data/ClimIdx/nao_pc_monthly.txt"
    settings["fn_PDO"] = settings["dir_project"]*"Data/ClimIdx/pdo.txt"
    settings["fn_MEI_1"] = settings["dir_project"]*"Data/ClimIdx/MEIext.data"
    settings["fn_MEI_2"] = settings["dir_project"]*"Data/ClimIdx/meiv2.data"

    # NCA5 projections
    settings["dir_NCA5_raw"] = homedir() * "/Data/NCA5/proj_2021_09_14/"
    settings["dir_NCA5"] =settings["dir_project"]*"Data/NCA5_projections/"

    # Files to save the Technical Report projections
    settings["fn_proj_glb"] = settings["dir_project"]*"Results/TR_global_projections.nc"
    settings["fn_proj_reg"] = settings["dir_project"]*"Results/TR_regional_projections.nc"
    settings["fn_proj_lcl"] = settings["dir_project"]*"Results/TR_local_projections.nc"
    settings["fn_proj_gri"] = settings["dir_project"]*"Results/TR_gridded_projections.nc"

    # Scenarios
    settings["NCA5_scenarios"]  = ["Low","IntLow","Int","IntHigh","High"]
    settings["regions"]         = ["USA","EC", "SE", "GCE",  "GCW",  "SWC","NWC","PAC","CAR","ALN","ALS"] # V1.0 names
    settings["regions"]         = ["USA","NE", "SE", "EGOM", "WGOM", "SW", "NW", "PAC","CAR","NAL","SAL"] # V1.1 names
    settings["processes"]       = ["AIS","GIS","glaciers","landwaterstorage","oceandynamics","total","verticallandmotion"]
    settings["processes_clim"]  = ["AIS","GIS","glaciers","landwaterstorage","oceandynamics"]

    # Trajectory years
    settings["years"] = [1900:2150...]
    settings["years_trajectory_global"] = [1970:2100...]
    settings["years_trajectory"] = [1970:2050...]
    settings["years_tg"] = [1920:2020...]
    settings["years_baseline"] = [2000]

    # Figure directories
    settings["dir_gmt"] = settings["dir_project"]*"GMT/"
    settings["dir_fig_1_gmsl_usa"] = settings["dir_gmt"]*"fig_1_gmsl_usa/"
    settings["dir_fig_2_regional"] = settings["dir_gmt"]*"fig_2_regional/"
    settings["dir_fig_3_map"] = settings["dir_gmt"]*"fig_3_map/"
    settings["dir_fig_4_gmsl"] = settings["dir_gmt"]*"fig_4_gmsl/"
    settings["dir_fig_5_divergence"] = settings["dir_gmt"]*"fig_5_divergence/"
    return settings
end

main()


text:
"""
    params(object)
    params(model::AbstractModel)

Returns a tuple of all `Param`s in the model or arbitrary object.
"""
function params end

"""
    printparams(object)
    printparams(io::IO, object)

Prints a table of all `Param`s in the object, similar to what
is printed in the repl for `AbstractModel`.
"""
function printparams end

"""
    stripparams(object)

Strips all `AbstractParam` from an object, replacing them with the `val`
field, or a combination of `val` and `units` if a `units` field exists.
"""
function stripparams end

"""
    update!(m::MutableModel, table)

Update the model in-place from an object that implements the Tables.jl interface.

Note: the parent object can be immutable, it will be completely rebuilt.
But the wrapper `AbstractModel` is mutable, such as `Model` or `InteractModel`.
"""
function update! end

"""
    update(m::AbstractModel, table)

Update the model from an object that implements the Tables.jl interface,
returning a new, updated object.
"""
function update end

"""
    withunits(object, [fieldname])
    withunits(model::AbstractModel, [fieldname])
    withunits(param::AbstractParam, [fieldname])

Returns the field specifed by `fieldname` (by default `:val`) for a single `Param`,
or a tuple of the `Param`s in a `Model` or arbitrary object.

If there is a `units` field the returned value will be a combination of the specied field
and the `units` fields.

If there is no units field or a specific `Param`s `units` fields contains `nothing`,
the field value is returned unchanged.
"""
function withunits end

"""
    stripunits(model::AbstractModel, xs)
    stripunits(param::AbstractParam, x)

Returns the `x` or `xs` divided by their corresponding units field, if it exists.

It there is no units field, and x has units, it will be returned with units! It
you want to simply remove all units, using Unitful.ustrip.
"""
function stripunits end

# Low-level, non-exported interface

"""
    setparent!(model::MutableModel, x)

Internal interface method to define for custom `AbstractModel` with a different field
for `parent`.

Set the parent object. Must be defined if the parent field of an `AbstractModel`
is not `:parent`.
"""
function setparent! end

"""
    setparent(model::AbstractModel, x)

Internal interface method to define for custom `AbstractModel` with a different field
for `parent`.

Set the parent object and return the rebuilt model. Must be defined
if the parent field of an `AbstractModel` is not `:parent`.
"""
function setparent end


text:
mutable struct Greedy
    num_lever::Int64
    num_play::Array{}
end


text:
using Flux, Statistics, DelimitedFiles
using Flux: onehotbatch, onecold, crossentropy, throttle
using Base.Iterators: repeated

# covtype.jl
println("covtype.jl")

# simple dense network with fluxml to reproduce covtype dataset results
# https://archive.ics.uci.edu/ml/datasets/covertype
# based on mlp.jl and housing.jl examples from fluxml/model-zoo

# ADAM hyperparamters, cranked up for fast/loose learning
adam_η = 0.05
adam_β = (0.8, 0.999)

# dense network with softmax output
m = Chain(
  Dense(54, 14, relu),
  Dense(14, 7),
  softmax)

# loss / accuracy / optimizer fns
loss(x, y) = crossentropy(m(x), y)
accuracy(x, y) = mean(onecold(m(x)) .== onecold(y))
optimizer = ADAM(adam_η, adam_β)

# load data
println("loading data")

# training and validation counts from original study
training_count = 11340
validation_count = 3780

# training batch size/count
batch_size = 20 # repeated dataset per training
batch_count = 10 # number of trainings

# read csv data
rawdata = readdlm("covtype-min.data", ',')
data = rawdata[1:(training_count + validation_count),:]

# in-line normalize datasets to 0.0-1.0, assume 0-intercept
normalize(x) = x ./ maximum(x)

# split column sets to normalize seperately
elevation_meters = normalize(data[:, 1:1])
aspect_azimuth = normalize(data[:, 2:2])
slope_degrees = normalize(data[:, 3:3])
horizontal_distances = normalize(hcat(data[:, 4:6], data[:, 10:10]))
hillshades = normalize(data[:, 7:9])
binary_columns = normalize(data[:, 11:54])

# input data
x = hcat(elevation_meters,aspect_azimuth,slope_degrees,horizontal_distances,
    hillshades,binary_columns)'
training_set_x = x[:,1:training_count]
test_set_x = x[:,(1 + training_count):(training_count + validation_count)]

println("  training input", size(training_set_x), "; test input", size(test_set_x))

# last column is target Cover_type 1-7
training_set_y = vec(data[1:training_count, 55:55]')
test_set_y = vec(data[(1 + training_count):(training_count + validation_count), 55:55]')
training_onehot_y = onehotbatch(training_set_y, 1:7)
test_onehot_y = onehotbatch(test_set_y, 1:7)

# train
println("training")

# train batch_count times
for i in 1:batch_count
  # entire dataset repeated batch_size times
  dataset = repeated((training_set_x, training_onehot_y), batch_size)
  Flux.train!(loss, params(m), dataset, optimizer)
  println("  ", accuracy(training_set_x, training_onehot_y), " ($i)")
end

# validate
println("validation")
println("  ", accuracy(test_set_x, test_onehot_y))

# sample of actual/expected output
println("covtype", floor.(Int, test_set_y[1:20]))
println("guess  ", onecold(m(test_set_x[:,1:20])))

max_stars_repo_name: masonblier/fluxml-covtype

max_issues_repo_name: masonblier/fluxml-covtype

max_forks_repo_name: masonblier/fluxml-covtype


text:
#
# Jammalamadaka circular correlation
#
# See Jammalamadaka, S. R., & Sengupta, A. (2001). Topics in Circular
# Statistics. World Scientific, p. 176
immutable JammalamadakaR{Normalized} <: NormalizedPairwiseStatistic{Normalized} end
JammalamadakaR() = JammalamadakaR{false}()
Base.eltype{T<:Real}(::JammalamadakaR, X::AbstractArray{Complex{T}}) = T

function sinmeanphasediff!(out, X)
    for j = 1:size(X, 2)
        # Compute mean phase
        m = zero(eltype(X))
        @simd for i = 1:size(X, 1)
            @inbounds m += X[i, j]
        end

        # Normalize to unit length
        m /= abs(m)

        # Compute difference from mean phase
        @simd for i = 1:size(X, 1)
          @inbounds out[i, j] = imag(X[i, j]*conj(m))
      end
    end
    out
end

allocwork{T<:Real}(t::JammalamadakaR{true}, X::AbstractVecOrMat{Complex{T}}) =
    Array(T, size(X, 1), size(X, 2))
function computestat!{T<:Real}(t::JammalamadakaR{true}, out::AbstractMatrix{T},
                               work::Matrix{T}, X::AbstractVecOrMat{Complex{T}})
    chkinput(out, X)

    # Sins minus mean phases
    workX = sinmeanphasediff!(work, X)

    # Products of phase differences
    At_mul_A!(out, workX)
    cov2coh!(out, out, Base.AbsFun())
end

immutable JammalamadakaRWorkXY{T<:Real}
    workX::Matrix{T}
    workY::Matrix{T}
    sumworkX::Matrix{T}
    sumworkY::Matrix{T}
end
function allocwork{T<:Real}(t::JammalamadakaR{true}, X::AbstractVecOrMat{Complex{T}}, Y::AbstractVecOrMat{Complex{T}})
    JammalamadakaRWorkXY{T}(Array(T, size(X, 1), size(X, 2)),
                            Array(T, size(Y, 1), size(Y, 2)),
                            Array(T, 1, nchannels(X)),
                            Array(T, 1, nchannels(Y)))
end
function computestat!{T<:Real}(t::JammalamadakaR{true}, out::AbstractMatrix{T}, work::JammalamadakaRWorkXY{T},
                               X::AbstractVecOrMat{Complex{T}}, Y::AbstractVecOrMat{Complex{T}})
    chkinput(out, X, Y)

    # Sins minus mean phases
    workX = sinmeanphasediff!(work.workX, X)
    workY = sinmeanphasediff!(work.workY, Y)

    # Products of phase differences
    At_mul_B!(out, workX, workY)
    cov2coh!(out, workX, workY, work.sumworkX, work.sumworkY, out, Base.AbsFun())
end


text:
using LocalitySensitive
using Test
using SparseArrays
using TextAnalysis

@testset "SimHash" begin
    @testset "Estimate cosine similarity" begin
        sh = SimHash(1)
        features_weights =  [("ab", 1), ("bc",1)]
        f = fingerprint(sh, features_weights)
        @test 1 == estimate_cosine(f, f)
        f2 = fingerprint(sh, [("fg", 1), ("de",1)])
        @test estimate_cosine(f, f2) < 0.2
    end

    @testset "Fingerprint sparse vector" begin
        sh = SimHash(1)
        f = fingerprint(sh, sparse([1, 0]))
        @test 1 == estimate_cosine(f, f)
        f2 = fingerprint(sh, sparse([0, 1]))
        @test estimate_cosine(f, f2) < 0.1
    end

    @testset "Find similar" begin
        sh = SimHash(1)
        shind = SimHashIndex(0.2)
        f = fingerprint(sh, sparse([1, 0]))
        push!(shind, f)
        @test [] == find_similar(shind, fingerprint(sh, sparse([0, 1])))
        @test [1] == find_similar(shind, fingerprint(sh, sparse([1, 0])))
        @test [1] == find_similar(shind, fingerprint(sh, sparse([0.5, 0])))
    end

    @testset  "Find all pair indices" begin
        sh = SimHash(1)
        shind = SimHashIndex(0.3)
        @test [] == similar_pairs(shind)
        docs = [
            "this is an english sentence",
            "this is an english sentence this is an english sentence",
            "this too, is an english sentence",
            "no resemblance whatsoever"
        ]
        docs = Document.(docs)
        crps = Corpus(docs)
        update_lexicon!(crps)
        m = DocumentTermMatrix(crps)
        features = tf_idf(m)
        for i in 1:length(docs)
            push!(shind, fingerprint(sh, features[i,:]))
        end
        @test [(1, 2), (1, 3), (2, 3)] == sort(similar_pairs(shind))
    end
end


text:
# ============================================================================ #
abstract type AbstractObjective end
abstract type SSE <: AbstractObjective end
abstract type WeightedSSE <: AbstractObjective end

abstract type TuningFunction{T<:AbstractObjective} end
# ============================================================================ #
mutable struct TuningData
    ts::Vector{Float64}
    evt::Vector{Float64}
    lab::Vector{Float64}
    dur::Float64
    tf::Float64
end
# ============================================================================ #
mutable struct FitData
    x::Vector{Float64}
    y::Vector{Float64}
    err::Vector{Float64}
end
# ---------------------------------------------------------------------------- #
FitData(x::Vector{Float64}, y::Vector{Float64}) = FitData(x, y, zeros(length(x)))
# ============================================================================ #
mutable struct TuningResult{T}
    res::Dict{Symbol, Float64}
    p::Vector{Float64}
    min::Float64
    flag::Symbol
    data::FitData
end
# ---------------------------------------------------------------------------- #
function TuningResult(::Type{T}, fd::FitData, p::Vector{Float64},
    mn::Float64, flag::Symbol) where {T<:TuningFunction}

    return TuningResult{T}(Dict(k=>v for (k,v) in zip(labels(T), p)), p, mn, flag, fd)
end
# ---------------------------------------------------------------------------- #
Base.getindex(tr::TuningResult, s::Symbol) = tr.res[s]
# ============================================================================ #
struct Contrast{T} <: TuningFunction{T} end
labels(::Type{Contrast{T}}) where T = [:max, :n, :c50, :baseline]
function parameters(::Type{Contrast{T}}, fd::FitData) where T
    xmn, xmx = extrema(fd.x)
    ymn, ymx = extrema(fd.y)

    p0 = [ymx, 1.5, xmn + ((xmx - xmn)*0.3), ymn]

    lb = [ymn, 0.6, xmn, 0.0]
    ub = [ymx*1.2, 6.0, xmx, ymn * 1.2]

    return p0, lb, ub
end
function generate(::Type{Contrast{T}}, x::Vector{Float64}, p::Vector) where T
    return p[1] .* (x.^p[2] ./ (x.^p[2] .+ p[3].^p[2])) .+ p[4]
end
# ============================================================================ #
struct LogisticContrast{T} <: TuningFunction{T} end
labels(::Type{LogisticContrast{T}}) where T = [:max, :k, :c50, :baseline]
function parameters(::Type{LogisticContrast{T}}, fd::FitData) where T
    p0, lb, ub = parameters(Contrast{SSE}, fd)
    p0[2] = 1.0
    # ub[2] = 15.0
    return p0, lb, ub
end
function generate(::Type{LogisticContrast{T}}, x::Vector{Float64}, p::Vector) where T
    return ((p[1] - p[4]) ./ (1.0 .+ exp.(-p[2].*(x .- p[3])))) .+ p[4]
end
# ============================================================================ #
struct Area{T} <: TuningFunction{T} end
labels(::Type{Area{T}}) where T = [:ke, :ki, :se, :si]
function parameters(::Type{Area{T}}, fd::FitData) where T
    mx = maximum(fd.y)
    p0 = [maximum(fd.y), 1.0, 1.0, 4.0]
    lb = [0.0, 0.0, fd.x[1], fd.x[1]]
    ub = [Inf, Inf, fd.x[end], fd.x[end]]
    return p0, lb, ub
end
function generate(::Type{Area{T}}, x::Vector{Float64}, p::Vector) where T
    step = 0.01
    ge = zeros(length(x))
    gi = copy(ge)
    for k = 1:length(x)
        x_sub = 2.0 .* range(-x[k]/2, x[k]/2, step=step)
        ge[k] = trapz(exp.(-1 .*((x_sub) ./ p[3]).^2))
        gi[k] = trapz(exp.(-1 .*((x_sub) ./ p[4]).^2))
    end
    return p[1].*ge.*step .- p[2].*gi.*step
end
# ============================================================================ #
trapz(y::AbstractVector{<:Real}) = 0.5 * y[1] + sum(y[2:end-1]) + 0.5 * y[end]
# ============================================================================ #


text:
b = ones(Float32, 2, 3)
ccall((:make_twice_array, "libexample"), Cvoid, (Any,), b)
println(b)


text:
isincident(v::VRepElement, h::HRepElement) = v in hyperplane(h)
isincident(h::HRepElement, v::VRepElement) = v in hyperplane(h)
isincident(p::Polyhedron, v::Index, h::Index) = isincident(get(p, v), get(p, h))

abstract type Incident{T, ElemT<:RepElement{T}, PT<:Polyhedron{T}, IdxT<:Index{T}} end

function Base.get(p::Polyhedron{T}, inc::Incident{T, ElemT}) where {T, ElemT}
    el = get(p, inc.idx)
    incs = _inctype(inc)[]
    for idx in Indices{T, ElemT}(p)
        eli = get(p, idx)
        if isincident(el, eli)
            push!(incs, _incel(inc, idx, eli))
        end
    end
    incs
end

struct IncidentElements{T, ElemT, PT, IdxT} <: Incident{T, ElemT, PT, IdxT}
    p::PT
    idx::IdxT
end
IncidentElements{T, ElemT}(p, idx) where {T, ElemT<:RepElement{T}} = IncidentElements{T, ElemT, typeof(p), typeof(idx)}(p, idx)
_inctype(inc::IncidentElements{T, ElemT}) where {T, ElemT} = ElemT
_incel(inc::IncidentElements, idx, eli) = eli

struct IncidentIndices{T, ElemT, PT, IdxT} <: Incident{T, ElemT, PT, IdxT}
    p::PT
    idx::IdxT
end
IncidentIndices{T, ElemT}(p, idx) where {T, ElemT<:RepElement{T}} = IncidentIndices{T, ElemT, typeof(p), typeof(idx)}(p, idx)
_inctype(inc::IncidentIndices{T, ElemT}) where {T, ElemT} = Index{T, ElemT}
_incel(inc::IncidentIndices, idx, eli) = idx

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra