text:
mutable struct ShowSettings
    tableformat::PrettyTableFormat
    floatformat::String
    border::Crayon
    header::Crayon
    subheader::Crayon
    fixed::Crayon
    error::Crayon
    highlighted::Crayon
    section::Crayon
    fixedpars::Bool
    ShowSettings() = new(unicode_rounded, "%10.4g",
                         crayon"light_blue", crayon"light_blue negative bold",
                         crayon"dark_gray bold", crayon"dark_gray",
                         crayon"light_red blink", crayon"negative", crayon"green bold",
                         true)
end

const showsettings = ShowSettings()

function printtable(args...; kw...)
    pretty_table(args...; tf=showsettings.tableformat, kw...,
                 border_crayon=showsettings.border,
                 header_crayon=showsettings.header,
                 subheader_crayon=showsettings.subheader, crop=:horizontal)
end

section(io, args...) = println(io, showsettings.section, args...)

function left(s::String, maxlen::Int)
    (length(s) <= maxlen)  &&  (return s)
    return s[1:maxlen]
end

function show(io::IO, dom::AbstractCartesianDomain)
    section(io, "Cartesian domain (ndims: ", ndims(dom), ", length: ", length(dom), ")")
    table = Matrix{Union{Int,Float64}}(undef, ndims(dom), 6)
    for i in 1:ndims(dom)
        a = dom[i]
        b = 0
        if length(a) > 1
            b = a .- circshift(a, 1)
            b = b[2:end]
        end
        table[i, 1] = i
        table[i, 2] = length(a)
        table[i, 3:6] = [minimum(a), maximum(a), minimum(b), maximum(b)]
    end
    printtable(io, table, ["Dim", "Size", "Min val", "Max val", "Min step", "Max step"],
               formatters=ft_printf(showsettings.floatformat, 3:6))
end

function show(io::IO, dom::AbstractLinearDomain)
    section(io, "Linear domain (ndims: ", ndims(dom), ", length: ", length(dom), ")")
    table = Matrix{Union{Int,Float64}}(undef, ndims(dom), 3)
    for i in 1:ndims(dom)
        table[i, 1] = i
        table[i, 2] = getaxismin(dom, i)
        table[i, 3] = getaxismax(dom, i)
    end
    printtable(io, table, ["Dim", "Min val", "Max val"],
               formatters=ft_printf(showsettings.floatformat, 2:3))
end

# Special case for Domain_1D: treat it as a Cartesian domain, despite it is a Linear one.
function show(io::IO, dom::Domain_1D)
    section(io, "Domain (ndims: ", ndims(dom), ", length: ", length(dom), ")")
    table = Matrix{Union{Int,Float64}}(undef, ndims(dom), 6)
    for i in 1:ndims(dom)
        a = dom[i]
        b = 0
        if length(a) > 1
            b = a .- circshift(a, 1)
            b = b[2:end]
        end
        table[i, 1] = i
        table[i, 2] = length(a)
        table[i, 3:6] = [minimum(a), maximum(a), minimum(b), maximum(b)]
    end
    printtable(io, table, ["Dim", "Size", "Min val", "Max val", "Min step", "Max step"],
               formatters=ft_printf(showsettings.floatformat, 3:6))
end

function show(io::IO, data::AbstractData)
    section(io, typeof(data), "   length: ", (length(data.val)))
    table = Matrix{Union{String,Float64}}(undef, 0, 7)

    names = fieldnames(typeof(data))
    error = Vector{Bool}()
    for name in names
        a = getfield(data, name)
        nan = length(findall(isnan.(a))) + length(findall(isinf.(a)))
        a = a[findall(isfinite.(a))]
        push!(error, nan > 0)
        table = vcat(table, [string(name) minimum(a) maximum(a) mean(a) median(a) std(a) (nan > 0  ?  string(nan)  :  "") ])
    end
    printtable(io, table, ["", "Min", "Max", "Mean", "Median", "Std. dev.", "Nan/Inf"],
               formatters=ft_printf(showsettings.floatformat, 2:6),
               highlighters=(Highlighter((data,i,j) -> error[i], showsettings.error)))
end

function show(io::IO, par::Parameter)
    if par.free
        println(io, "Value : ", par.val, "  [", par.low , " : ", par.high, "]")
    else
        println(io, "Value : ", par.val, "   (FIXED)")
    end
end

function preparetable(comp::AbstractComponent, cname="")
    table = Matrix{Union{String,Float64}}(undef, 0, 5)
    fixed = Vector{Bool}()
    error = Vector{Bool}()

    ctype = split(string(typeof(comp)), ".")
    (ctype[1] == "DataFitting")  &&   (ctype = ctype[2:end])
    ctype = join(ctype, ".")

    for (pname, param) in getparams(comp)
        parname = string(pname[1])
        if pname[2] >= 1
            parname *= "[" * string(pname[2]) * "]"
        end
        parname *= (param.free  ?  ""  :  " (FIXED)")
        (!showsettings.fixedpars)  &&  (!param.free)  &&  continue
        range = strip(@sprintf("%7.2g:%-7.2g", param.low, param.high))
        (range == "-Inf:Inf")  &&  (range = "")
        table = vcat(table, [cname ctype parname param.val range])
        push!(fixed, !param.free)
        push!(error, !(param.low <= param.val <= param.high))
        cname = ""
        ctype = ""
    end
    return (table, fixed, error)
end

show(io::IO, comp::AbstractComponent) =
    show(io, OrderedDict(Symbol("?") => comp), OrderedDict(Symbol("?") => true))

function show(io::IO, dict::OrderedDict{Symbol, T}, cfree::OrderedDict{Symbol, Bool}) where T <: AbstractComponent
    (length(dict) > 0)  ||  (return nothing)
    table = Matrix{Union{String,Float64}}(undef, 0, 5)
    fixed = Vector{Bool}()
    error = Vector{Bool}()
    hrule = Vector{Int}()
    push!(hrule, 0, 1)
    for (cname, comp) in dict
        (t, f, e) = preparetable(comp, string(cname) .* (cfree[cname]  ?  ""  :  " (FIXED)"))
        table = vcat(table, t)
        append!(fixed, f .| (.!cfree[cname]))
        append!(error, e)
        push!(hrule, length(error)+1)
    end
    printtable(io, table , ["Component" "Type" "Param." "Value" "Range"], alignment=:l,
               hlines=hrule, formatters=ft_printf(showsettings.floatformat, [4]),
               highlighters=(Highlighter((data,i,j) -> fixed[i], showsettings.fixed),
                             Highlighter((data,i,j) -> (error[i] &&  (j in (3,4))), showsettings.error)))
end

show(io::IO, mime::MIME"text/plain", model::Model) = show(io, model)
function show(io::IO, model::Model)
    section(io, "Components:")
    length(model.comps) != 0  || (return nothing)
    show(io, model.comps, model.cfree)

    for i in 1:length(model.preds)
        println(io)
        section(io, "Prediction #$i ")
        show(io, model.preds[i])
    end
end

function show(io::IO, pred::Prediction)
    (length(pred.cevals) == 0)  &&  (return nothing)
    table = Matrix{Union{String,Int,Float64}}(undef, length(pred.cevals) + 1, 6)
    error = Vector{Bool}()

    i = 1
    for (cname, ceval) in pred.cevals
        result = ceval.eval
        v = view(result, findall(isfinite.(result)))
        (length(v) == 0)  &&  (v = [NaN])
        nan = length(findall(isnan.(result)))
        inf = length(findall(isinf.(result)))
        table[i, 1] = string(cname)
        table[i, 2] = ceval.counter
        table[i, 3:5] = [minimum(v), maximum(v), mean(v)]
        table[i, 6] = (nan+inf > 0  ?  string(nan)  :  "")
        push!(error, (nan+inf > 0))
        i += 1
    end

    result = pred.eval
    v = view(result, findall(isfinite.(result)))
    nan = length(findall(isnan.(result)))
    inf = length(findall(isinf.(result)))
    table[i, 1] = "Reduced"
    table[i, 2] = pred.counter
    table[i, 3:5] = [minimum(v), maximum(v), mean(v)]
    table[i, 6] = (nan+inf > 0  ?  string(nan)  :  "")
    push!(error, (nan+inf > 0))

    printtable(io, table, ["Label", "Counter", "Min", "Max", "Mean", "NaN/Inf"], alignment=:l,
               hlines=[0,1,length(pred.cevals)+1,length(pred.cevals)+2],
               formatters=ft_printf(showsettings.floatformat, 3:5),
               highlighters=(Highlighter((data,i,j) -> (error[i] && j==5), showsettings.error)))
end

show(io::IO, par::BestFitPar) = println(io, par.val, " ± ", par.unc,
                                        (par.val == par.calc  ?  ""  :
                                         " (calculated value: " * string(par.calc) * ")"))

function preparetable(comp::BestFitComp)
    table = Matrix{Union{String,Float64}}(undef, 0, 5)
    fixed = Vector{Bool}()
    error = Vector{Bool}()
    watch = Vector{Bool}()

    for (pname, param) in comp
        if isa(param, Vector{BestFitPar})
            for ii in 1:length(param)
                par = param[ii]
                (!showsettings.fixedpars)  &&  (!par.free)  &&  continue
                spname = string(pname) * "[" * string(ii) * "]"
                table = vcat(table, ["" spname par.val par.unc par.calc])
                push!(fixed, !par.free)
                push!(error, !isfinite(par.unc))
                push!(watch, par.val != par.calc)
            end
        else
            par = param
            (!showsettings.fixedpars)  &&  (!par.free)  &&  continue
            spname = string(pname)
            table = vcat(table, ["" spname par.val par.unc par.calc])
            push!(fixed, !par.free)
            push!(error, !isfinite(par.unc))
            push!(watch, par.val != par.calc)
        end
    end
    return (table, fixed, error, watch)
end

function show(io::IO, comp::BestFitComp)
    (table, fixed, error, watch) = preparetable(comp)
    (length(table) == 0)  &&  return
    printtable(io, table , ["Component" "Param." "Value" "Uncert." "Calculated"], alignment=:l,
               hlines=[0,1,size(table)[1]+1], formatters=ft_printf(showsettings.floatformat, [3,4,5]),
               highlighters=(Highlighter((data,i,j) -> (fixed[i]  &&  (j in (2,3,4))), showsettings.fixed),
                             Highlighter((data,i,j) -> (watch[i]  &&  (j==5)), showsettings.highlighted),
                             Highlighter((data,i,j) -> (error[i]  &&  (!fixed[i])  &&  (j==4)), showsettings.error)))
end

function show(io::IO, res::BestFitResult)
    section(io, "Best Fit results:")

    table = Matrix{Union{String,Float64}}(undef, 0, 5)
    fixed = Vector{Bool}()
    error = Vector{Bool}()
    watch = Vector{Bool}()
    hrule = Vector{Int}()
    push!(hrule, 0, 1)
    for (cname, comp) in res.comps
        if length(comp) > 0
            (t, f, e, w) = preparetable(comp)
            (length(t) > 0)  ||  continue
            t[1,1] = string(cname)
            table = vcat(table, t)
            append!(fixed, f)
            append!(error, e)
            append!(watch, w)
            push!(hrule, length(error)+1)
        end
    end
    printtable(io, table , ["Component" "Param." "Value" "Uncert." "Calculated"], alignment=:l,
               hlines=hrule, formatters=ft_printf(showsettings.floatformat, [3,4,5]),
               highlighters=(Highlighter((data,i,j) -> (fixed[i]  &&  (j in (2,3,4))), showsettings.fixed),
                             Highlighter((data,i,j) -> (watch[i]  &&  (j==5)), showsettings.highlighted),
                             Highlighter((data,i,j) -> (error[i]  &&  (!fixed[i])  &&  (j==4)), showsettings.error)))

    println(io)
    println(io, @sprintf("    #Data  : %10d              Cost : %-10.5g", res.ndata, res.cost))
    println(io, @sprintf("    #Param : %10d              Red. : %-10.4g", res.ndata-res.dof, res.cost / res.dof))
    print(  io, @sprintf("    DOF    : %10d              ", res.dof))
    if res.log10testprob < -3
        println(io, @sprintf("Prob.: 10^%-10.4g", res.log10testprob))
    else
        println(io, @sprintf("Prob.: %10.4g", 10^res.log10testprob))
    end
    printstyled(io, "    Status :  ")
    if res.status == :Optimal
        printstyled(color=:green, io, @sprintf("%-15s", "Optimal"))
    elseif res.status == :NonOptimal
        printstyled(color=printcolorerr(), io, @sprintf("%-15s", "Non Optimal"))
    elseif res.status == :Warn
         printstyled(color=printcolorerr(), io, @sprintf("%-15s", "Warning"))
    elseif res.status == :Error
        printstyled(color=printcolorerr(), io, @sprintf("%-15s", "Error"))
    else
        printstyled(color=printcolorerr(), io, @sprintf("%-15s", "Unknown (" * string(res.status) * "), see fitter output"))
    end
    println(io, @sprintf("        Elapsed: %-10.4g s", res.elapsed))
end


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

JLLWrappers.@generate_main_file_header("libLLVM")
JLLWrappers.@generate_main_file("libLLVM", UUID("8f36deef-c2a5-5394-99ed-8e07531fb29a"))
end  # module libLLVM_jll

ext: jl


text:
@everywhere begin
  using Pkg
  Pkg.activate("juliaenv")
  using Distributions
  using SharedArrays
end

"""
    onesimulation(N, n_rare_cat, mu)

Simulate data from a multinomial distribution with 3 categories
and `N` units total.
- `n_rare_cat` = 1 or 2 rare categories, with very small probabilities
- `mu`: expected (or mean) number of observations in each rare category.
  The probability that each observation is in a given rare category is
  `p=mu/N`.

With 2 rare categories, the probabilities are [p,p,1-2p].
With 1 rare category, the probabilities are [p,(1-p)/2,(1-p)/2].

The data (counts in each categories) are then analyzed for a
goodness-of-fit test:
1. first using the Pearson test statistics: sum of (oi - ei)^2 / ei
   where ei = expected count in category i, oi = observed count
2. second using Qlog statistics: 2 * sum of L(ni,ei) * (oi - ei)/ei
   where L(x,y) is the logarithmic mean of x and y:
   L(x,y) = (x-y) / (log(x)-log(y)).
   L(x,y) is always between x and y, as expected for a mean.
   See for example Lorenzen (1995) for a justification for using Qlog.

Both test statistics are compared to a chi-square distribution
with 2 degrees of freedom to obtain p-values.

output: (pearson_stat, Qlog_stat, pearson_pvalue, Qlog_pvalue)

examples:

```julia
onesimulation( 50, 1, 0.1)
onesimulation(500, 2, 2.0)
```
"""

@everywhere function onesimulation(N,nrarecat,mu)
  p = mu/N
  @assert p>=0.0 && p<=0.5 "the rare category/ies cannot have a negative or >.5 probability"
  if nrarecat==1
    probs = [p,(1.0-p)/2,(1.0-p)/2]
    @assert p<=1.0
  else
    probs = [p, p, (1.0-2p)]
  end
  df = 2 # because 3 categories
  observed = rand(Multinomial(N,probs), 1) # 3x1 array
  phat = observed./N # estimated probabilities of the 3 categories
  pearson = N * sum((phat - probs).^2 ./probs)
  Qlog = 0.0
  for i in 1:3
    if observed[i] == 0 continue; end # no need to do xxx / log(0) = 0.0
    ratio = phat[i]/probs[i]
    ratio ≉ 1.0 || continue
    dev = ratio - 1.0 # deviation on ratio scale
    # if ratio ≈ 1 then dev/log(ratio) ≈ 1, but NaN if ratio = 1 exactly
    Qlog += probs[i] * dev^2 / log(ratio)
  end
  Qlog *= 2N
  # cdf = cumulative distribution function: P(random value <= x)
  # ccdf = complementary cdf:               P(random value >  x)
  p_pearson = ccdf(Chisq(df), pearson)
  p_Qlog    = ccdf(Chisq(df), Qlog)
  return pearson,Qlog, p_pearson,p_Qlog
end

nreps = 200 # start with nreps=2 only
samplesizes = [30, 100, 1000]
nrarecats = [1, 2]
mus = [0.1, 1., 2.]

Nssize    = length(samplesizes)
Nnrarecat = length(nrarecats)
Nmu       = length(mus)

ntotal = nreps * Nssize * Nnrarecat * Nmu

samplesize = SharedArray(repeat(samplesizes, inner = nreps, outer = Nnrarecat * Nmu))
nrarecat = SharedArray(repeat(nrarecats,     inner = nreps * Nssize, outer = Nmu))
mu = SharedArray(repeat(mus,                 inner = nreps * Nssize * Nnrarecat))
# hcat(samplesize, nrarecat, mu) # to double-check, when nreps=2

pearson   = SharedArray{Float64}(ntotal)
qlog      = SharedArray{Float64}(ntotal)
p_pearson = SharedArray{Float64}(ntotal)
p_qlog    = SharedArray{Float64}(ntotal)

@sync @distributed for i in 1:ntotal
  res = onesimulation(samplesize[i], nrarecat[i], mu[i])
  #@show samplesize[i], nrarecat[i], mu[i] # when nreps = 2
  pearson[i], qlog[i], p_pearson[i], p_qlog[i] = res
end

# save results, to analyze later interactively
using DataFrames
using CSV
df = DataFrame(
  samplesize = samplesize,
  numrarecategories = nrarecat,
  numrareexpected = mu,
  pearson = pearson,
  qlog = qlog,
  pval_pearson = p_pearson,
  pval_qlog = p_qlog
)
CSV.write("runsimulations_$(nreps)replicates.csv", df)

max_stars_repo_name: morrowcj/computingtools

max_issues_repo_name: morrowcj/computingtools

max_forks_repo_name: morrowcj/computingtools


text:
@testset "Apply Vulkan prefixes" begin
    @test prefix_vk(c1_novk) == c1
    @test prefix_vk(c2_novk) == c2
    @test prefix_vk(s2_novk) == s2
    @test prefix_vk(c1_novk.value) == c1.value
    @test prefix_vk(c2_novk.value) == c2.value
    @test prefix_vk(s2_novk.value) == s2.value
end

@testset "Remove Vulkan prefixes" begin
    @test remove_prefix(s2) == s2_novk
    @test remove_prefix(c1) == c1_novk
    @test remove_prefix(c2) == c2_novk
end


text:
"""
Area scaling factor for a function f with partial derivatives fx, fy
    * fx::Function
    * fy::Function

Returns:
    * [Function] : Area scaling factor at x.
"""
function af(fx, fy)
    return x -> norm(fx(x)×fy(x))
end

function af(f)
    # Find partial derivatives
    fx(x) = ForwardDiff.derivative(z -> f([z, x[2]]), x[1])
    fy(x) = ForwardDiff.derivative(z -> f([x[1], z]), x[2])

    # Compute area scaling factor
    return x -> norm(fx(x)×fy(x))
end

"""
Derivative of area scaling factor for a function f with derivatives
    * fx::Function
    * fy::Function
    * fxx::Function
    * fxy::Function
    * fyy::Function

Returns:
    * [Function] : Derivative of area scaling factor at x.
"""
function Daf(af, fx, fy, fxx, fxy, fyy)
    return x -> [
        (fxx(x) × fy(x) - fxy(x) × fx(x)) ⋅ (fx(x) × fy(x)) / af(x)
        (fxy(x) × fy(x) - fyy(x) × fx(x)) ⋅ (fx(x) × fy(x)) / af(x)
    ]
end

"""
Get q-map and derivative of a given function
    * f::Function

Returns:
    * q::Function, Dq::Function
"""
function Qmap(f)
    # Find partial derivatives
    fx(x) = ForwardDiff.derivative(z -> f([z, x[2]]), x[1])
    fy(x) = ForwardDiff.derivative(z -> f([x[1], z]), x[2])

    # Compute area scaling factor
    a = af(fx, fy)

    # Return q-map and derivative
    q(x) = √a(x) * f(x)
    Dq(x) = ForwardDiff.jacobian(q, x)
    return q, Dq
end

"""
Get q-map and derivative of a given function with derivatives
    * f::Function
    * fx::Function
    * fy::Function

Returns:
    * q::Function, Dq::Function
"""
function Qmap(f, fx, fy)
    a = af(fx, fy)

    q(x) = √a(x) * f(x)
    Dq(x) = ForwardDiff.jacobian(q, x)
    return q, Dq
end

"""
Get q-map and derivative of a given function with derivatives
    * f::Function
    * fx::Function
    * fy::Function
    * fxx::Function
    * fxy::Function
    * fyy::Function

Returns:
    * q::Function, Dq::Function
"""
function Qmap(f, fx, fy, fxx, fxy, fyy)
    a = af(fx, fy)
    Da = Daf(a, fx, fy, fxx, fxy, fyy)

    Df(x) = [fx(x) fy(x)]
    q(x) = √a(x) * f(x)
    Dq(x) = 1. / (2. * √a(x)) * f(x) * Da(x)' + √a(x) * Df(x)
    return q, Dq
end

"""
Compute the centroid of the parametric surface.
    * f::Function
    * af:: Function

Returns:
    * [Array{Float}] : Vector value of surface centre
"""
function surface_centroid(f, af)
    return hcubature(x -> f(x)*af(x), [0, 0], [1, 1])[1] / hcubature(af, [0, 0], [1, 1])[1]
end

"""
Translate function so that centroid coincides with the origin.
    * f::Function
    * af:: Function

Returns:
    * [Function] : Translated Surface
"""
function centered(f, af)
    c = surface_centroid(f, af)
    return x -> f(x) - c
end

"""
Compute the parametric surface area.
    * f::Function
    * af:: Function

Returns:
    * [Float] : Surface Area
"""
function surface_area(f, af)
    return hcubature(af, [0, 0], [1, 1])[1]
end

"""
Scale surface to unit area.
    * f::Function
    * af:: Function

Returns:
    * [Function] : Scaled Surface
"""
function scale_surface(f, af)
    s = surface_area(f, af)
    return x -> f(x) / s
end

"""
Translate and scale surface
    * f::Function

Returns:
    * [Function] : Scaled Surface
"""
function normalize_surface(f)
    # Find partial derivatives
    fx(x) = ForwardDiff.derivative(z -> f([z, x[2]]), x[1])
    fy(x) = ForwardDiff.derivative(z -> f([x[1], z]), x[2])

    # Compute area scaling factor
    a = af(fx, fy)

    c = surface_centroid(f, a)
    s = surface_area(f, a)
    return x -> (f(x) - c) / s
end

"""
Translate and scale surface
    * f::Function
    * af:: Function

Returns:
    * [Function] : Scaled Surface
"""
function normalize_surface(f, af)
    c = surface_centroid(f, af)
    s = surface_area(f, af)
    return x -> (f(x) - c) / s
end


text:
plural(::Type{Nanosecond}) = Nanoseconds
plural(::Type{Microsecond}) = Microseconds
plural(::Type{Millisecond}) = Milliseconds
plural(::Type{Second}) = Seconds
plural(::Type{Minute}) = Minutes
plural(::Type{Hour}) = Hours
plural(::Type{Day}) = Days
plural(::Type{Week}) = Weeks
plural(::Type{Month}) = Months
plural(::Type{Year}) = Years

Nanoseconds(x::Period) = Nanosecond(x)
Nanoseconds(x::CompoundPeriod) = sum(map(Nanosecond, x.periods))
Nanoseconds(x::ReverseCompoundPeriod) = sum(map(Nanosecond, x.cperiod.periods))

nanoseconds(x::T) where {T<:Union{Period,Compounds}} = Nanoseconds(x).value

TimeUnits(x::CompoundPeriod) = sum(map(mintype(x), x.periods))
TimeUnits(x::ReverseCompoundPeriod) = sum(map(mintype(x), x.cperiod.periods))
TimeUnits(x::Period) = x

Microseconds(x::T) where {T<:Nanosecond} = floor(x, Microsecond)
Microseconds(x::T) where {T<:Period} = Microsecond(x)
Microseconds(x::T) where {T<:Compounds} = Microseconds(TimeUnits(x))

microseconds(x::T) where {T<:Union{Period,Compounds}} = Microseconds(x).value

Milliseconds(x::T) where {T<:Union{Nanosecond, Microsecond}} = floor(x, Millisecond)
Milliseconds(x::T) where {T<:Period} = Millisecond(x)
Milliseconds(x::T) where {T<:Compounds} = Milliseconds(TimeUnits(x))

milliseconds(x::T) where {T<:Union{Period,Compounds}} = Milliseconds(x).value

Seconds(x::T) where {T<:Union{Nanosecond, Microsecond, Millisecond}} = floor(x, Second)
Seconds(x::T) where {T<:Period} = Second(x)
Seconds(x::T) where {T<:Compounds} = Seconds(TimeUnits(x))

seconds(x::T) where {T<:Union{Period,Compounds}} = Seconds(x).value

Minutes(x::T) where {T<:Union{Nanosecond, Microsecond, Millisecond, Second}} = floor(x, Minute)
Minutes(x::T) where {T<:Period} = Minute(x)
Minutes(x::T) where {T<:Compounds} = Minutes(TimeUnits(x))

minutes(x::T) where {T<:Union{Period,Compounds}} = Minutes(x).value

Hours(x::T) where {T<:Union{Nanosecond, Microsecond, Millisecond, Second, Minute}} = floor(x, Hour)
Hours(x::T) where {T<:Period} = Hour(x)
Hours(x::T) where {T<:Compounds} = Hours(TimeUnits(x))

hours(x::T) where {T<:Union{Period,Compounds}} = Hours(x).value

Days(x::T) where {T<:Union{Nanosecond, Microsecond, Millisecond, Second, Minute, Hour}} = floor(x, Day)
Days(x::T) where {T<:Period} = Day(x)
Days(x::T) where {T<:Compounds} = Days(TimeUnits(x))

days(x::T) where {T<:Union{Period,Compounds}} = Days(x).value

Weeks(x::T) where {T<:Union{Nanosecond, Microsecond, Millisecond, Second, Minute, Hour, Day}} = floor(x, Week)
Weeks(x::T) where {T<:Period} = Week(x)
Weeks(x::T) where {T<:Compounds} = Weeks(TimeUnits(x))

weeks(x::T) where {T<:Union{Period,Compounds}} = Weeks(x).value

Months(x::T) where {T<:Union{Nanosecond, Microsecond, Millisecond, Second, Minute, Hour, Day}} = Month(0)
Months(x::T) where {T<:Year} = floor(x, Month)
Months(x::T) where {T<:Period} = Month(x)
Months(x::CompoundPeriod) = sum(map(Month, x.periods))
Months(x::ReverseCompoundPeriod) = sum(map(Month, x.cperiod.periods))

months(x::T) where {T<:Union{Period,Compounds}} = Months(x).value

Years(x::T) where {T<:Union{Nanosecond, Microsecond, Millisecond, Second, Minute, Hour, Day}} = Year(0)
Years(x::T) where {T<:Month} = floor(x, Year)
Years(x::T) where {T<:Period} = Year(x)
Years(x::CompoundPeriod) = sum(map(Year, x.periods))
Years(x::ReverseCompoundPeriod) = sum(map(Year, x.cperiod.periods))

years(x::T) where {T<:Union{Period,Compounds}} = Years(x).value

for (P, Q) in ((:Years,:years), (:Months,:months), (:Weeks,:weeks), (:Days,:days),
               (:Hours,:hours), (:Minutes,:minutes), (:Seconds,:seconds),
               (:Milliseconds,:milliseconds), (:Microseconds,:microseconds),
               (:Nanoseconds,:nanoseconds))
    @eval begin
        $P(x::DateTime) = $P(CompoundPeriod(x))
        $Q(x::DateTime) = $Q(CompoundPeriod(x))
    end
end

for (P, Q) in ((:Years,:years), (:Months,:months), (:Days, :days))
    @eval begin
        $P(x::Date) = $P(CompoundPeriod(x))
        $Q(x::Date) = $Q(CompoundPeriod(x))
    end
end

for (P, Q) in ((:Hours,:hours), (:Minutes,:minutes), (:Seconds,:seconds),
               (:Milliseconds,:milliseconds), (:Microseconds,:microseconds),
               (:Nanoseconds,:nanoseconds))
    @eval begin
        $P(x::Time) = $P(CompoundPeriod(x))
        $Q(x::Time) = $Q(CompoundPeriod(x))
    end
end


text:
using SparseArrays, LinearAlgebra

include("./limited_dijkstra.jl")

function sortperm_unique_zerofirst(values)
    p = sortperm(values, rev=true)
    sorted = values[p]

    unique_vals = []
    unique_perm = []
    ctr = 0
    for i in 1:length(sorted)
        if i == 1 || sorted[i] != sorted[i-1] # new value
            push!(unique_perm, [])
            ctr += 1
            push!(unique_perm[ctr], p[i])
            push!(unique_vals, sorted[i])
        else sorted[i] == sorted[i-1] # existing value
            push!(unique_perm[ctr], p[i])
        end
    end

    unique_vals_zerofirst = []
    unique_perm_zerofirst = []

    for i in 1:length(unique_vals)
        v = unique_vals[i]
        if v > 0
            push!(unique_vals_zerofirst, unique_vals[i])
            push!(unique_perm_zerofirst, unique_perm[i])
        elseif v==0
            pushfirst!(unique_vals_zerofirst, unique_vals[i])
            pushfirst!(unique_perm_zerofirst, unique_perm[i])
        end
    end

    return unique_vals_zerofirst, unique_perm_zerofirst
end

function _limited_edge_betweenness_reducer(
        a::Tuple{Array{Array{Float64,1},1}, Array{SparseMatrixCSC{Float64,Int64},1}},
        b::Tuple{Array{Array{Float64,1},1}, Array{SparseMatrixCSC{Float64,Int64},1}}
    )::Tuple{Array{Array{Float64,1},1}, Array{SparseMatrixCSC{Float64,Int64},1}}
    return (a[1]+b[1], a[2]+b[2])
end

function limited_edge_betweenness_centrality(g::AbstractGraph,
    vs::AbstractVector=vertices(g),
    distmx::AbstractMatrix{T}=weights(g);
    normalize=true,
    endpoints=false,
    limits=[zero(T)]) where T <: Real

    n_v = nv(g)
    k = length(vs)
    isdir = is_directed(g)

    betweenness = zeros(n_v)
    edge_betweenness = spzeros(n_v, n_v)

    #sort
    (limits_sorted, limits_perm) = sortperm_unique_zerofirst(limits)

    (betweenness, edge_betweenness) = @distributed _limited_edge_betweenness_reducer for s in vs
    #for s in vs
        temp_betweenness = [zeros(n_v) for i in 1:length(limits_sorted)]
        temp_edge_betweenness = [spzeros(n_v, n_v) for i in 1:length(limits_sorted)]

        if degree(g, s) > 0  # this might be 1?
            state = limited_dijkstra_shortest_paths(g, s, distmx; allpaths=true, trackvertices=true, limit=limits_sorted[1])

            for i in 1:length(limits_sorted)
                trim_dijkstra_state!(state, limits_sorted[i])
                if endpoints
                    _edge_accumulate_endpoints!(temp_betweenness[i], temp_edge_betweenness[i], state, g, s)
                else
                    _edge_accumulate_basic!(temp_betweenness[i], temp_edge_betweenness[i], state, g, s)
                end
            end

        end
        (temp_betweenness, temp_edge_betweenness)
    end

    if !isdir
        edge_betweenness = [SparseMatrixCSC(UpperTriangular(i + i')) for i in edge_betweenness]
    end

    # rescale
    for i in 1:length(betweenness)
        LightGraphs._rescale!(betweenness[i],
        n_v,
        normalize,
        isdir,
        k)
    end

    for i in 1:length(edge_betweenness)
        _edge_rescale!(edge_betweenness[i],
        n_v,
        normalize,
        isdir,
        k)
    end

    #unsort
    betweenness_unsorted::Vector{Any} = [nothing for limit in limits]
    edge_betweenness_unsorted::Vector{Any} = [nothing for limit in limits]
    for (p, i) in zip(limits_perm, 1:length(limits_perm))
        for perm in p
            betweenness_unsorted[perm] = betweenness[i]
            edge_betweenness_unsorted[perm] = edge_betweenness[i]
        end
    end

    return (betweenness_unsorted, edge_betweenness_unsorted)
end

function limited_betweenness_centrality(g::AbstractGraph,
    vs::AbstractVector=vertices(g),
    distmx::AbstractMatrix{T}=weights(g);
    normalize=true,
    endpoints=false,
    limits=[zero(T)]) where T <: Real

    n_v = nv(g)
    k = length(vs)
    isdir = is_directed(g)

    betweenness = zeros(n_v)
    edge_betweenness = spzeros(n_v, n_v)

    #sort
    (limits_sorted, limits_perm) = sortperm_unique_zerofirst(limits)

    betweenness = @distributed (+) for s in vs
    #for s in vs
        temp_betweenness = [zeros(n_v) for i in 1:length(limits_sorted)]

        if degree(g, s) > 0  # this might be 1?
            state = limited_dijkstra_shortest_paths(g, s, distmx; allpaths=true, trackvertices=true, limit=limits_sorted[1])

            for i in 1:length(limits_sorted)
                trim_dijkstra_state!(state, limits_sorted[i])
                if endpoints
                    LightGraphs._accumulate_endpoints!(temp_betweenness[i], state, g, s)
                else
                    LightGraphs._accumulate_basic!(temp_betweenness[i], state, g, s)
                end
            end
        end
        temp_betweenness
    end

    # rescale
    for i in 1:length(betweenness)
        LightGraphs._rescale!(betweenness[i],
        n_v,
        normalize,
        isdir,
        k)
    end

    #unsort
    betweenness_unsorted::Vector{Any} = [nothing for limit in limits]
    for (p, i) in zip(limits_perm, 1:length(limits_perm))
        for perm in p
            betweenness_unsorted[perm] = betweenness[i]
        end
    end

    return betweenness_unsorted
end

# Accumulation

function _edge_accumulate_basic!(betweenness::Vector{Float64},
    edge_betweenness::AbstractMatrix{Float64},
    state::LightGraphs.DijkstraState,
    g::AbstractGraph,
    si::Integer)

    n_v = length(state.parents) # this is the ttl number of vertices
    δ = zeros(n_v)
    σ = state.pathcounts
    P = state.predecessors

    # make sure the source index has no parents.
    P[si] = []
    # we need to order the source vertices by decreasing distance for this to work.
    S = reverse(state.closest_vertices) #Replaced sortperm with this
    for w in S
        coeff = (1.0 + δ[w]) / σ[w]
        for v in P[w]
            if v > 0
                δ[v] += (σ[v] * coeff)
                edge_betweenness[v, w] += (σ[v] * coeff)
            end
        end
        if w != si
            betweenness[w] += δ[w]
        end
    end
    return nothing
end

function _edge_accumulate_endpoints!(betweenness::Vector{Float64},
    edge_betweenness::AbstractMatrix{Float64},
    state::LightGraphs.DijkstraState,
    g::AbstractGraph,
    si::Integer)

    n_v = length(state.parents) # this is the ttl number of vertices
    δ = zeros(n_v)
    σ = state.pathcounts
    P = state.predecessors

    # make sure the source index has no parents.
    P[si] = []
    # we need to order the source vertices by decreasing distance for this to work.
    S = reverse(state.closest_vertices) #Replaced sortperm with this

    betweenness[s] += length(S) - 1    # 289

    for w in S
        coeff = (1.0 + δ[w]) / σ[w]
        for v in P[w]
            δ[v] += σ[v] * coeff
            edge_betweenness[v, w] += (σ[v] * coeff)
        end
        if w != si
            betweenness[w] += (δ[w] + 1)
        end
    end
    return nothing

end

# rescale

function _edge_rescale!(edge_betweenness::AbstractMatrix{Float64}, n::Integer, normalize::Bool, directed::Bool, k::Integer)
    if normalize
        if n <= 2
            do_scale = false
        else
            do_scale = true
            scale = 1.0 / (n * (n - 1))
        end
    else
        if !directed
            do_scale = true
            scale = 1.0 / 2.0
        else
            do_scale = false
        end
    end
    if do_scale
        if k > 0
            scale = scale * n / k
        end
        edge_betweenness .*= scale

    end
    return nothing
end

max_stars_repo_name: kymok/local-betweenness

max_issues_repo_name: kymok/local-betweenness

max_forks_repo_name: kymok/local-betweenness


text:
# INTERFACE
sample_action(m, ds, idx, A, rng::AbstractRNG) = rand(rng, action_distribution(m, ds, idx, A))

observe(m, ds, new_context, dt) = new_context
observe(m::Void, ds, new_context, dt) = nothing

# statetype(t::Type, c::Type) = throw(MethodError(statetype, (t,c)))
# statetype(x::Any, c) = statetype(typeof(x))

# NO CRASH DRIVER

@with_kw struct NoCrashMLDriver{M1, M2, S}
    long::M1  = NORMAL_IDM
    lat::M2   = nothing
    safety::S = nothing
end
NoCrashMLDriver(dt::Float64) = NoCrashMLDriver(NORMAL_IDM, nothing, MaxSafeAccel(dt))

initial_state(d::NoCrashMLDriver, c::Context) = MLDriverState(initial_state(d.long), initial_state(d.lat))

struct MLDriverState{S1, S2, H<:HighwayDrivingContext}
    long::S1
    lat::S2
    c::H
end

# function statetype(::Type{NoCrashMLDriver{M1,M2,S}}) where {M1, M2, S}
#     return MLDriverState{statetype(M1), statetype(M2), }
# end

function observe(m::NoCrashMLDriver, ds, new_context, dt)
    long = observe(m.long, ds.long, new_context, dt)
    lat = observe(m.lat, ds.lat, new_context, dt)
    return MLDriverState(long, lat, new_context)
end

function action_distribution(m::NoCrashMLDriver, state::MLDriverState, idx, A::Type{AccelSpeed})
    c = state.c
    max_safe = max_safe_acc(m.safety, state.long, idx)
    a::LongAccel = action(m.long, state.long, idx, LongAccel)
    mid = min(a.a, max_safe)
    lower = min(-1e-5, max(-max_acc(m.long)/2, -brake_limit(vehicle(c, idx))-mid))
    upper = min(max_acc(m.long)/2, max(max_safe-mid, 1e-5))
    long_dist = TriangularDist(mid+lower, mid+upper, mid)

    lat_dist = DeltaDist(0.0)

    return AccelSpeedDist(long_dist, lat_dist)
end

struct AccelSpeedDist{D1, D2}
    long::D1
    lat::D2
end

rand(rng::AbstractRNG, d::AccelSpeedDist) = AccelSpeed(rand(rng, d.long), rand(rng, d.lat))
function pdf(d::AccelSpeedDist, val::AbstractVector{Float64})
    @assert length(val) == 2
    return pdf(d.long, val[1])*pdf(d.lat, val[2])
end
sampletype(::Type{AccelSpeedDist}) = AccelSpeed

struct DeltaDist val::Float64 end

rand(rng::AbstractRNG, d::DeltaDist) = d.val
pdf(d::DeltaDist, val::Float64) = convert(Float64, val == d.val)
sampletype(::Type{DeltaDist}) = Float64


text:
# Downloads Requirements
using Pkg; Pkg.activate("../"); Pkg.instantiate()

using AutoGrad, Knet, KnetLayers, Morse

Morse.main(ARGS)


text:
module Mixin

using Base.Meta:
    isexpr

macro interface(ex)
    @assert isexpr(ex, :struct, 3) "struct expression should be given"

    _, typedecl, body = ex.args
    T = get_type_sym(typedecl)
    flddecls = get_fielddecls(body)

    getter = create_getter(T, flddecls)
    actualtypedecl = :(abstract type $(typedecl) end)
    create_mixin_asserter(__module__, actualtypedecl, flddecls)

    return quote
        $(actualtypedecl) # already evaled so actually not needed to be returned, but just for explicit macro expansion
        $(getter) # define getter in user code
    end
end

get_type_sym(typedecl) = isexpr(typedecl, :(<:)) ? first(typedecl.args) : typedecl
get_fielddecls(body) = return filter(body.args) do x
    isexpr(x, :(::)) || x isa Symbol
end

get_field_sym(flddecl) = isexpr(flddecl, :(::)) ? first(flddecl.args) : flddecl
get_type(mod, flddecl) = isexpr(flddecl, :(::)) ? Core.eval(mod, last(flddecl.args)) : Any

function create_getter(T, flddecls)
    body = Expr(:block)

    for flddecl in flddecls
        isexpr(flddecl, :(::)) || continue # untyped

        sig, ty = flddecl.args
        get_typed_ex = :(if sym === $(QuoteNode(sig))
            return getfield(obj, sym)::$(esc(ty))
        end)
        push!(body.args, get_typed_ex)
    end

    fallback = :(return getfield(obj, sym))
    push!(body.args, fallback)

    return :(function Base.getproperty(obj::$(esc(T)), sym::Symbol)
        $(body)
    end)
end

function assert_mixin end # will be overloaded on each interface declaration

function create_mixin_asserter(defmod, actualtypedecl, flddecls)
    # get actual types
    Core.eval(defmod, actualtypedecl)
    T = get_type_sym(first(actualtypedecl.args))
    Ty = Core.eval(defmod, T)
    fldsym2ty_mixin = Dict(
        get_field_sym(flddecl) => get_type(defmod, flddecl)
        for flddecl in flddecls
    )

    asserter = :(function assert_mixin(mod, interface::$(Type{Ty}), body)
        flddecls = $(get_fielddecls)(body)
        fldsym2ty = Dict(
            $(get_field_sym)(flddecl) => $(get_type)(mod, flddecl)
            for flddecl in flddecls
        )

        foreach($(fldsym2ty_mixin)) do (fldsym_mixin, ty_mixin)
            @assert haskey(fldsym2ty, fldsym_mixin) string("interface ", $(Ty), " requires ", fldsym_mixin, "::", ty_mixin)
            ty = fldsym2ty[fldsym_mixin]
            @assert ty <: ty_mixin string("mixin field ", $(Ty), '.', fldsym_mixin, " is expected to be declared as subtype of ", ty_mixin, " but declared as ", fldsym_mixin, "::", ty)
        end
    end)

    Core.eval(@__MODULE__, asserter)
end

macro mixin(ex)
    @assert isexpr(ex, :struct, 3) "struct expression should be given"
    _, typedecl, body = ex.args
    @assert isexpr(typedecl, :(<:)) "mixed interface should be given as supertype"
    T = last(typedecl.args)

    Ty = Core.eval(__module__, T)
    assert_mixin(__module__, Ty, body)

    return esc(ex)
end

export
    @interface, @mixin

end


text:
using Documenter, PandaModels

makedocs(;
    modules=[PandaModels],
    authors="e2nIEE",
    repo="https://github.com/e2nIEE/PandaModels.jl/blob/{commit}{path}#L{line}",
    sitename="PandaModels.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://e2nIEE.github.io/PandaModels.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Manual" => [
        "Getting Started" => "quickguide.md",
        ],
        "Tutorials" => [
        # "Power Flow" => "pm_pf.md",
        "Optimal Power Flow" => "pm_opf.md",
        "Optimal Transmission Switching" => "pm_ots.md",
        # "Timeseries and Multinetwork" => "pm_ts_mn.md",
        "Optimal MultiNetwork Storage" => "pm_omns.md",
        "Transmission Network Expansion Planning" => "pm_tnep.md",
        # "Optimal Voltage Deviation" => "pm_vd.md",
        # "Radial Distribution Network" => "pm_rds.md",
        ],
        "Developer" => [
                "Develop Mode" => "develop.md",
                # "Model Guidlines" => "newmodel.md",
                # "Call Model in pandapower" => "newmodelpp.md",
                "Add Test" => "test.md",
        ],
    ],

    # repo = "https://github.com/e2nIEE/PandaModels.jl/blob/{commit}{path}#L{line}",
    # sitename = "PandaModels.jl",
    doctest = true,
    linkcheck = true,
    # format = Documenter.HTML(
    #     # See https://github.com/JuliaDocs/Documenter.jl/issues/868
    #     prettyurls = get(ENV, "CI", nothing) == "true",
    #     analytics = "UA-178297470-1",
    #     collapselevel = 1,
    #     )
)

deploydocs(;
    repo="github.com/e2nIEE/PandaModels.jl",
    # push_preview = true
)


text:
using Plexifications
using Base.Test: @test, @test_throws

@test_throws ErrorException begin
    inv(BiComplex(0, 0.0))
end

@test begin
    iszerodivisor(BiComplex(0,1,1,0))
end

@test begin
    !iszerodivisor(BiComplex(1,2,3,4))
end

@test begin
    a = Plexifications.randomBigFloat()
    isreal(BiComplex(a))
end

@test begin
    a = 1
    !isreal(BiComplex(a, a))
end

@test begin
    z = random(BiComplex{Int})
    z == +(z)
end

@test begin
    a = Plexifications.randomBigFloat()
    real(BiComplex(a)) == a
end

@test begin
    io = IOBuffer()
    show(io, BiComplex(1, 2, 3, 4))
    l = String(take!(io))
    r = "[1: 1, H: 2, G: 3, GH: 4]"
    l == r
end

@test begin
    x = random(BiComplex{BigFloat})
    y = random(BiComplex{BigFloat})
    x != y
end

@test begin
    a = Plexifications.randomBigInt()
    z = random(BiComplex{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    z = random(BiComplex{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(Cmplex{BigInt})
    y = random(BiComplex{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    l = x * y
    r = y * x
    l == r
end

@test begin
    x = random(Cmplex{BigInt})
    y = random(BiComplex{BigInt})
    l = x * y
    r = y * x
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(Cmplex{BigInt})
    y = random(BiComplex{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    z = random(BiComplex{BigInt})
    l = (x + y) + z
    r = x + (y + z)
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    z = random(BiComplex{BigInt})
    l = (x * y) * z
    r = x * (y * z)
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    barH(barH(x)) == x
end

@test begin
    x = random(BiComplex{BigInt})
    barG(barG(x)) == x
end

@test begin
    x = random(BiComplex{BigInt})
    barG(barH(x)) == barH(barG(x))
end

@test begin
    x = random(BiComplex{Rational{BigInt}})
    inv(inv(x)) == x
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    l = barH(x * y)
    r = barH(x) * barH(y)
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    l = barH(x + y)
    r = barH(x) + barH(y)
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    l = barH(x - y)
    r = barH(x) - barH(y)
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    l = barG(x * y)
    r = barG(x) * barG(y)
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    l = barG(x + y)
    r = barG(x) + barG(y)
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    l = barG(x - y)
    r = barG(x) - barG(y)
    l == r
end

@test begin
    x = random(BiComplex{Rational{BigInt}})
    y = random(BiComplex{Rational{BigInt}})
    l = inv(x * y)
    r = inv(y) * inv(x)
    l == r
end

@test begin
    x = random(BiComplex{Rational{BigInt}})
    y = random(BiComplex{Rational{BigInt}})
    l = inv(x / y)
    r = y / x
    l == r
end

@test begin
    x = random(BiComplex{Rational{BigInt}})
    y = random(BiComplex{Rational{BigInt}})
    l = inv(x \ y)
    r = y \ x
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    z = random(BiComplex{BigInt})
    l = x * (y + z)
    r = (x * y) + (x * z)
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    l = a * (x + y)
    r = (a * x) + (a * y)
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    x = random(BiComplex{BigInt})
    l = a * x
    r = x * a
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    z = random(BiComplex{BigInt})
    l = (x + y) * z
    r = (x * z) + (y * z)
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    z = random(BiComplex{BigInt})
    l = x * (y - z)
    r = (x * y) - (x * z)
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    z = random(BiComplex{BigInt})
    l = (x - y) * z
    r = (x * z) - (y * z)
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    norm(x) > big(0)
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    l = norm(x * y)
    r = norm(x) * norm(y)
    l == r
end

@test_throws ErrorException begin
    inv(BiComplex(0, 0.0))
end

@test_throws ErrorException begin
    random(BiComplex{Int}) / 0
end

@test_throws ErrorException begin
    0 \ random(BiComplex{Int})
end

@test begin
    a = Plexifications.randomBigRational()
    z = random(BiComplex{Rational{BigInt}})
    l = z / a
    r = a \ z
    l == r
end

@test begin
    a = Plexifications.randomBigRational()
    z = random(BiComplex{Rational{BigInt}})
    l = a / z
    r = z \ a
    l == r
end

@test begin
    z = random(BiComplex{BigInt})
    z == z + zero(BiComplex{BigInt})
end

@test begin
    z = random(BiComplex{BigInt})
    z == z * one(BiComplex{BigInt})
end

@test begin
    z = random(BiComplex{BigInt})
    z == z + zero(z)
end

@test begin
    z = random(BiComplex{BigInt})
    z == z * one(z)
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    l = commutator(x, y)
    r = zero(BiComplex{BigInt})
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    l = alternatorL(x, y)
    r = zero(BiComplex{BigInt})
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    l = alternatorR(x, y)
    r = zero(BiComplex{BigInt})
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    l = flexator(x, y)
    r = zero(BiComplex{BigInt})
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    z = random(BiComplex{BigInt})
    l = associator(x, y, z)
    r = zero(BiComplex{BigInt})
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(BiComplex{BigInt})
    z = random(BiComplex{BigInt})
    l = jacobiator(x, y, z)
    r = 3 * x * y * z
    l == r
end

@test begin
    l = BiComplex(Cmplex(1,2))
    r = BiComplex(1,2,0)
    l == r
end

@test begin
    l = BiComplex(1, 2, 3.0)
    r = BiComplex(1.0, 2.0, 3.0)
    l == r
end

@test begin
    l = BiComplex(1, 2, 3, 4.0)
    r = BiComplex(1.0, 2.0, 3.0, 4.0)
    l == r
end


text:
# Script to generate Figures 18 and 19
using FileIO, Images, ImageFiltering, JLD2, MultiscaleGraphSignalTransforms, MultivariateStats, Plots, StatsBase
include("auxilaries.jl")

# Set up the resolution and display size
gr(dpi=200, size=(800,600))

# Load the original 512x512 composite texture image and its mastge
textures = FileIO.load("5block.png")
textures = Matrix{Float64}(textures) # convert it to a regular matrix

# Compute the mask image using the Gabor features + PCA
#
# Step 1: Parameter setups
# Orientations
ndir = 2
Θ = [ π/3; 5*π/6 ]

# Spatial frequencies
Ξ = [0.2; 0.3; 0.4; 0.5]
# Convert them to wavelengths
Λ = 1.0 ./ Ξ

# Gaussian bandwidth for Gabor filters
bw = 1.0
σ = 1.0/π * sqrt(log(2)/2) * (2^bw+1)/(2^bw-1)

# Spatial aspect ratio: σx = σ; σy = σ/γ
γ = 1.0

# Step 2: Apply Gabor filters; the same size as the input image.
m, n = size(textures)

K = length(Θ) *length(Λ) # K should be square number
F = zeros(m, n, K) # This is the space for the Gabor filtered images
Fg = zeros(m, n, K) # This is the space for the Gaussian smoothed version of F
kr = zeros(m+1, n+1, K)
ki = zeros(m+1, n+1, K)

k = 1
kl = (0, 0)
for θ in Θ
    for λ in Λ
        kernel = Kernel.gabor(m, n, λ*σ, θ, λ, γ, 0.0)
        kernel[1] ./= 2*π*σ*σ/γ
        kernel[2] ./= 2*π*σ*σ/γ
        global kl = size(kernel[1])
        kr[1:kl[1],1:kl[2],k] = kernel[1]
        ki[1:kl[1],1:kl[2],k] = kernel[2]
        cker = (centered(kernel[1]), centered(kernel[2]))
        F[:,:,k] = sqrt.(imfilter(textures, reflect(cker[1])).^2 +
            imfilter(textures, reflect(cker[2])).^2)
        # "reflect"ing the kernel is necessary for convolution;
        # otherwise it does the correlation.
        Fg[:,:,k] = imfilter(F[:,:,k], Kernel.gaussian(3.0*λ))
        global k += 1
    end
end

# Display those Gabor features
plot(K, layout=(length(Θ),length(Λ)), framestyle=:none, legend=false)
for k = 1:K
    heatmap!(kr[226:286,226:286,k], subplot=k, ratio=1, yaxis=:flip, showaxis=false,
        ticks=false, c=:grays, clims=(minimum(kr), maximum(kr)), colorbar=false)
end
display(current())

# Displaying the imaginary part of Gabor kernels may not be necessary...
#plot(K, layout=(Int(sqrt(K)),Int(sqrt(K))), framestyle=:none, legend=false)
#for k = 1:K
#    heatmap!(ki[:,:,k], subplot=k, ratio=1, yaxis=:flip, showaxis=false,
#      ticks=false, c=:grays, clims=(minimum(ki), maximum(ki)), colorbar=false)
#end
#display(current())

plot(K, layout=(length(Θ),length(Λ)), framestyle=:none, legend=false)
for k = 1:K
    heatmap!(Fg[:,:,k], subplot=k, ratio=1, yaxis=:flip, showaxis=false,
        ticks=false, c=:grays, clims=(0, maximum(F)), colorbar=false)
end
display(current())

# Step 3.  Compute the PCA and extract the first principal component
Xtr=reshape(Fg, (m*n, K))
Xtr=Xtr'
dt=StatsBase.fit(ZScoreTransform, Xtr) # Each variable is standardized to have
# mean 0, stddev 1
Xtr=StatsBase.transform(dt, Xtr)
M = MultivariateStats.fit(PCA, Xtr)
Ytr=MultivariateStats.transform(M, Xtr)
Ytr=Ytr'
Ytr=reshape(Ytr, (m, n, size(Ytr, 2)))
plot(K, layout=(length(Θ),length(Λ)), framestyle=:none, legend=false)
for k = 1:outdim(M)
    heatmap!(Ytr[:,:,k], subplot=k, ratio=1, yaxis=:flip, showaxis=false,
             ticks=false, c=:grays, clims=(minimum(Ytr), maximum(Ytr)), colorbar=false)
end
display(current())

mask = Ytr[:,:,1]
display(heatmap(mask, ratio=1, c=:grays, yaxis=:flip))
# End of 512 x 512 mask generation.

#
# Subsample both the original and mask images to 128 x 128
#
textures = deepcopy(textures[1:4:end,1:4:end])
mask = deepcopy(mask[1:4:end,1:4:end])
m, n = size(textures) # recapture the subsampled matrix size

#
# Normalize the mask
#
mask = (mask.-minimum(mask))./(maximum(mask)-minimum(mask));

#
# Generate Fig. 18a
#
display(heatmap(textures, ratio=1, yaxis =:flip, showaxis = false, ticks = false,
        color = :grays, colorbar = false))
savefig("textures_orig.pdf")
savefig("textures_orig.png")

#
# Generate Fig. 18b
#
display(heatmap(mask, ratio=1, yaxis =:flip, showaxis = false, ticks = false,
        color = :grays, colorbar = false))
savefig("textures_mask.pdf")
savefig("textures_mask.png")

#
# Now, let's compute the weight matrix based on the Gaussian affinity
# of the small windows (or radius r) around pixels and the pixel locations.
#
# Set up the key parameters
r = 3 # specify radius of neighbors; so far the best
σ = 0.0005 # the best so far with r=3 or r=5

# Do the weight matrix computation
W = image_Gaussian_affinity(mask, r, σ)

#
# Generate G (GraphSig object) and GP (GraphPart object) using the computed
# weight matrix W, and the subsampled texture image as a graph signal.
# Note that GraphSig requires a matrix data even if it is just a one vector, i.e., f = textures[:] does not work!
G = GraphSig(W, f = reshape(textures, (length(textures), 1)))
GP = partition_tree_fiedler(G)
dmatrix = ghwt_analysis!(G, GP=GP)

#
# Construct or search the specific basis
#
# Haar
BS_haar = bs_haar(GP)
dvec_haar = dmatrix2dvec(dmatrix, GP, BS_haar)

# Walsh
BS_walsh = bs_walsh(GP)
dvec_walsh = dmatrix2dvec(dmatrix, GP, BS_walsh)

# GHWT_c2f
dvec_c2f, BS_c2f = ghwt_c2f_bestbasis(dmatrix, GP)

# GHWT_f2c
dvec_f2c, BS_f2c = ghwt_f2c_bestbasis(dmatrix, GP)

# eGHWT
dvec_eghwt, BS_eghwt = eghwt_bestbasis(dmatrix, GP)

#
# Generate Figure 19a (Approximation error plot)
#
DVEC = [ dvec_haar[:], dvec_walsh[:], dvec_c2f[:], dvec_f2c[:], dvec_eghwt[:] ]
T = [ "Graph Haar", "Graph Walsh", "GHWT_c2f ", "GHWT_f2c", "eGHWT" ]
L = [ (:dashdot,:orange), (:dashdot,:blue), (:solid, :red), (:solid, :green),
      (:solid, :black) ]
approx_error2(DVEC, T, L, 0.5)
savefig("textures_approx_error_r3_sigma0005.pdf")
savefig("textures_approx_error_r3_sigma0005.png")

#
# Generate Figure 19b (Top 9 eGHWT basis vectors)
#
top_vectors_plot3(dvec_eghwt, m, n, BS_eghwt, GP)
savefig("textures_eghwt09_r3_sigma0005.pdf")
savefig("textures_eghwt09_r3_sigma0005.png")


text:
#Test code
#-------------------------------------------------------------------------------

using LibPSFC

#No real test code yet... just demonstrate use:

function testfileaccess(path::String)
	println("\n\nfile: $path")
	println("------------------")
	reader = LibPSFC._open(path)
	display(reader.properties)
	@show names(reader)
	close(reader)
end

sampledata(filename::String) = joinpath(LibPSFC.rootpath, "core/data", filename)

testfileaccess(sampledata("opBegin"))
testfileaccess(sampledata("pss0.fd.pss"))
testfileaccess(sampledata("timeSweep"))
testfileaccess(sampledata("srcSweep"))

signame = "INN"
println("\nRead in $signame vector:")
println("------------------")
reader = LibPSFC._open(sampledata("timeSweep"))
t = readsweep(reader)
y = read(reader, signame)
close(reader)

@show t[1], t[end]
flush(stdout)

:Test_Complete


text:
### Routines for constructing smooth functions that generate field data ###

import Base: *, -, +, show

abstract type AbstractSpatialField end
abstract type AbstractGeneratedField end

#include("motionprofiles.jl")
using SpaceTimeFields
import SpaceTimeFields: Abstract1DProfile, radius, strength

## Empty spatial field

"""
    EmptySpatialField()

Create a blank spatial field. This is primarily useful for initializing a
sum of spatial fields.

# Example
```jldoctest
julia> g = EmptySpatialField()
EmptySpatialField()

julia> g(2,3)
0.0
```
"""
struct EmptySpatialField <: AbstractSpatialField end
(g::EmptySpatialField)(a...) = Float64(0)

## Spatial Gaussian field ##

"""
    SpatialGaussian(σx,σy,x0,y0,A[,derivdir=0])

Set up a spatial field in the form of a Gaussian centered at `x0,y0` with
radii `σx` and `σy` in the respective directions and amplitude `A`. If the
optional parameter `deriv` is set to 1 or 2, then it returns the first
derivative of a Gaussian in that direction (`x` or `y`, respectively).

`SpatialGaussian(σx,σy,x0,y0,A,u,v[,derivdir=0])` generates a Gaussian
that convects at velocity `(u,v)`. It can be evaluated with an additional
argument for time.
"""
struct SpatialGaussian{CT,GX,GY} <: AbstractSpatialField
  gx :: GX
  gy :: GY
  A :: Float64
  u :: Float64
  v :: Float64
  SpatialGaussian(gx::Abstract1DProfile,gy::Abstract1DProfile,A,u,v) = new{true,typeof(gx),typeof(gy)}(gx,gy,A,u,v)
  SpatialGaussian(gx::Abstract1DProfile,gy::Abstract1DProfile,A) = new{false,typeof(gx),typeof(gy)}(gx,gy,A,0.0,0.0)
end

SpatialGaussian(σx::Real,σy::Real,x0::Real,y0::Real,A::Real;deriv::Int=0) =
                _spatialdgaussian(σx,σy,x0,y0,A,Val(deriv))
SpatialGaussian(σx::Real,σy::Real,x0::Real,y0::Real,A::Real,u::Real,v::Real;deriv::Int=0) =
                _spatialdgaussian(σx,σy,x0,y0,A,u,v,Val(deriv))

_spatialdgaussian(σx,σy,x0,y0,A,::Val{0}) = SpatialGaussian(Gaussian(σx,A) >> x0,
                                                            Gaussian(σy,1) >> y0,A)
_spatialdgaussian(σx,σy,x0,y0,A,::Val{1}) = SpatialGaussian(DGaussian(σx,A) >> x0,
                                                            Gaussian(σy,1) >> y0,A)
_spatialdgaussian(σx,σy,x0,y0,A,::Val{2}) = SpatialGaussian(Gaussian(σx,A) >> x0,
                                                            DGaussian(σy,1) >> y0,A)

_spatialdgaussian(σx,σy,x0,y0,A,u,v,::Val{0}) = SpatialGaussian(Gaussian(σx,A) >> x0,
                                                                Gaussian(σy,1) >> y0,A,u,v)
_spatialdgaussian(σx,σy,x0,y0,A,u,v,::Val{1}) = SpatialGaussian(DGaussian(σx,A) >> x0,
                                                                Gaussian(σy,1) >> y0,A,u,v)
_spatialdgaussian(σx,σy,x0,y0,A,u,v,::Val{2}) = SpatialGaussian(Gaussian(σx,A) >> x0,
                                                                DGaussian(σy,1) >> y0,A,u,v)

SpatialGaussian(σ,x0,y0,A;deriv::Int=0) = SpatialGaussian(σ,σ,x0,y0,A,deriv=deriv)

(g::SpatialGaussian{GX,GY})(x,y) where {CT,GX,GY} = g.gx(x)*g.gy(y)
# ignore the time argument if it is called with this...
(g::SpatialGaussian{false,GX,GY})(x,y,t) where {GX,GY} = g(x,y)
(g::SpatialGaussian{true,GX,GY})(x,y,t) where {GX,GY} = g.gx(x-g.u*t)*g.gy(y-g.v*t)

## Scaling spatial fields

struct ScaledField{N <: Real, P <: AbstractSpatialField} <: AbstractSpatialField
    s::N
    p::P
end
function show(io::IO, p::ScaledField)
    print(io, "$(p.s) × ($(p.p))")
end
s::Number * p::AbstractSpatialField = ScaledField(s, p)
-(p::AbstractSpatialField) = ScaledField(-1, p)

(p::ScaledField)(x,y) = p.s*p.p(x,y)
(p::ScaledField)(x,y,t) = p.s*p.p(x,y,t)

## Adding spatial fields together.

struct AddedFields{T <: Tuple} <: AbstractSpatialField
    ps::T
end
function show(io::IO, Σp::AddedFields)
    println(io, "AddedFields:")
    for p in Σp.ps
        println(io, "  $p")
    end
end

"""
    p₁::AbstractSpatialField + p₂::AbstractSpatialField

Add the fields so that `(p₁ + p₂)(x,y) = p₁(x,y) + p₂(x,y)`.

"""
+(p::AbstractSpatialField, Σp::AddedFields) = AddedFields((p, Σp.ps...))
+(Σp::AddedFields, p::AbstractSpatialField) = AddedFields((Σp.ps..., p))
function +(Σp₁::AddedFields, Σp₂::AddedFields)
    AddedFields((Σp₁..., Σp₂...))
end

-(p₁::AbstractSpatialField, p₂::AbstractSpatialField) = p₁ + (-p₂)

+(p::AbstractSpatialField...) = AddedFields(p)

# Evaluate at x,y , assuming that t = 0 for any time-varying member
function (Σp::AddedFields)(x,y)
    f = 0.0
    for p in Σp.ps
        f += p(x,y,0.0)
    end
    f
end

# Evaluate at x,y,t, ignoring t for any constant member
function (Σp::AddedFields)(x,y,t)
    f = 0.0
    for p in Σp.ps
        f += p(x,y,t)
    end
    f
end

## For generating an actual instance of a spatial field
"""
    GeneratedField(d::GridData,field::AbstractSpatialField...,grid::PhysicalGrid)

Create an instance of a spatial field function `field` on scalar grid data `d`,
based on a grid `grid`. After creating the instance `g = GeneratedField(d,field,grid)`,
then the resulting grid data can be accessed by typing `g()`. For vector grid data,
a separate `field` must be supplied for each component.

If the fields are time dependent, then you can also evaluate `g(t)` at the
desired time. The time argument is ignored if the fields are static.
"""
struct GeneratedField{T <: GridData}
    fielddata :: T
    fieldfcns :: Vector{AbstractSpatialField}
    grid :: PhysicalGrid
end

function GeneratedField(d::ScalarGridData,field::AbstractSpatialField,g::PhysicalGrid)
    #xg, yg = coordinates(d,g)
    #tmp = _generatedfield(xg,yg,field,d,0.0)
    tmp = similar(d)
    _generatedfield!(tmp,field,g,0.0)
    GeneratedField{typeof(tmp)}(tmp,AbstractSpatialField[field],g)
end

function GeneratedField(d::CollectedGridData,fields::Vector{AbstractSpatialField},g::PhysicalGrid)
  @assert length(fields) == _numberofcomponents(typeof(d))
  tmp = similar(d)
  _generatedfield!(tmp,fields,g,0.0)
  GeneratedField{typeof(tmp)}(tmp,fields,g)
end

function _generatedfield!(d::ScalarGridData,field::AbstractSpatialField,g::PhysicalGrid,t::Real)
    xg, yg = coordinates(d,g)
    d .= _generatedfield(xg,yg,field,d,t)
    return d
end

_generatedfield!(d::ScalarGridData,fields::Vector{AbstractSpatialField},g::PhysicalGrid,t::Real) =
    _generatedfield!(d,fields[1],g,t)

function _generatedfield!(d::CollectedGridData,fields::Vector{AbstractSpatialField},g::PhysicalGrid,t::Real)
  cnt = 0
  for (i,fname) in enumerate(propertynames(d))
    if typeof(getfield(d,fname)) <: GridData
        cnt += 1
        x, y = coordinates(getfield(d,fname),g)
        getfield(d,fname) .= _generatedfield(x,y,fields[cnt],getfield(d,fname),t)
    end
  end
  return d
end

GeneratedField(d::VectorGridData,
      fieldu::AbstractSpatialField,fieldv::AbstractSpatialField,
      g::PhysicalGrid) = GeneratedField(d,AbstractSpatialField[fieldu,fieldv],g)

_generatedfield(xg,yg,field::AbstractSpatialField,d::ScalarGridData,t) =
          typeof(d)(_evaluatefield(field,xg,yg,t))

_evaluatefield(field,xg,yg,t) = field.(xg*ones(1,length(yg)),ones(length(xg))*yg',t)

datatype(f::GeneratedField{T}) where {T} = T
grid(f::GeneratedField{T}) where {T} = f.grid

(f::GeneratedField{T})() where {T} = f.fielddata

function (f::GeneratedField)(t::Real) where {T<:ScalarGridData}
    _generatedfield!(f.fielddata,f.fieldfcns,f.grid,t)
end

## For generating a transient form of a spatial field
"""
    ModulatedField(g::GeneratedField,modfcn::Abstract1DProfile)

Create a time-modulated form of a generated spatial field, useful for
introducing a forcing field onto the grid. The supplied field `g`
is modulated by a function `modfcn` with a specified profle shape. The
resulting object can be evaluated with a single argument (time) and returns a
`GridData` type object of the same type contained in `g`.
"""
struct ModulatedField
  gfield :: GeneratedField
  modfcn :: Abstract1DProfile
end

(f::ModulatedField)(t) = f.modfcn(t)*f.gfield()

datatype(f::ModulatedField) = datatype(f.gfield)

"""
    PulseField(g::GeneratedField,t0,σ)

Create a transient form of a generated spatial field, useful for
introducing a short-lived forcing field onto the grid. The supplied field `g`
is modulated by a Gaussian pulse centered at time `t0` with width `σ`. This pulse's
maximum is unity, so that the amplitude of the overall field is set by `g`. The
resulting object can be evaluated with a single argument (time) and returns a
`ScalarGridData` type object of the same type contained in `g`.
"""
#=
struct PulseField
  gfield :: GeneratedField
  modfcn :: Abstract1DProfile
end
=#

PulseField(gfield::GeneratedField{T},t0::Real,σt::Real) where {T <: GridData} =
            ModulatedField(gfield,Gaussian(σt,sqrt(π*σt^2)) >> t0)

#(f::PulseField)(t) = f.modfcn(t)*f.gfield()

#datatype(f::PulseField) = datatype(f.gfield)

## Allow all spatial fields to accept a vector argument

#for T in InteractiveUtils.subtypes(AbstractSpatialField)
#    @eval (f::$T)(x::Vector{R}) where {R <: Float64} = f(x...)
#end

max_stars_repo_name: JuliaIBPM/CartesianGrids

max_issues_repo_name: JuliaIBPM/CartesianGrids

max_forks_repo_name: JuliaIBPM/CartesianGrids


text:
using Pkg, DelimitedFiles, MPIReco

# Install required packages
for P in [:HTTP,:TensorDecompositions, :PyPlot]
 try
   @eval using $P
 catch
   Pkg.add(String(P))
   @eval using $P
 end
end

# Download data
include("downloadData.jl")

include("./utils/subsampling.jl")
include("./utils/FRTruncation.jl")

# load system matrix
datadirSM = "./data/"
bSF = MPIFile(datadirSM*"SM.mdf")
# backup for local testing
# bSF = MPIFile("/opt/mpidata/20180228_180048_OpenMPIData_1_1/9")
nx,ny,nz = calibSize(bSF);
f = [2145,8629,60903]
S = ComplexF64.( getSystemMatrix(bSF,f,bgCorrection=true, tfCorrection=false) )

# sampling locations (linear index)
samplingIdx = vec( readdlm(joinpath(@__DIR__,"data/pdSamplingIdx.txt"), '\t', Int64) )

# get undersampled data
y = sfMeas(S, samplingIdx)

# CS-recovery
par_cs = Dict{Symbol,Any}()
par_cs[:shape] = (nx,ny,nz)   # size of the spatial measurement grid
par_cs[:λ_l1] = 0.02          # l1-regularization parameter
par_cs[:ρ_l1] = 2.0	      # parameter ρ associated with the l1-regularization
par_cs[:iterationsInner] = 50 # max. number of inner Split Bregman iterations
par_cs[:iterations] = 10      # number of outer Split Bregman iterations
par_cs[:relTol] = 1.e-2       # stopping criterion for the inner Split Bregman loop
par_cs[:absTol] = 1.e-4       # stopping criterion for the inner Split Bregman loop

S_cs = smRecovery(y, samplingIdx, par_cs)

# CSLR-recovery
par_cslr = Dict{Symbol,Any}()
par_cslr[:shape] = (nx,ny,nz)
par_cslr[:λ_l1] = 0.02
par_cslr[:ρ_l1] = 2.0
par_cslr[:λ_lr] = 0.02	        # TNN regularization parameter (for CSLR and CSFR)
par_cslr[:ρ_lr] = 0.2           # parameter ρ associated with the low rank regularization
par_cslr[:iterationsInner] = 50
par_cslr[:iterations] = 10
par_cslr[:relTol] = 1.e-2
par_cslr[:absTol] = 1.e-4

S_cslr = smRecovery(y, samplingIdx, par_cslr)

# reshape into tensor format
S = reshape(S,nx,ny,nz,:)
S_cs = reshape(S_cs,nx,ny,nz,:)
S_cslr = reshape(S_cslr,nx,ny,nz,:)

# CSFR
r = (6,6,6) # maximum HOSVD rank
S_csfr = FRTruncation3d(S_cslr,r)

# plot frequency components
figure()
subplot(4,3,1)
title("k=$(f[1])")
ylabel("measured")
imshow(abs.(S[:,:,17,1]),cmap="gray"); xticks([]); yticks([])
subplot(4,3,2)
title("k=$(f[2])")
imshow(abs.(S[:,:,17,2]),cmap="gray"); xticks([]); yticks([])
subplot(4,3,3)
title("k=$(f[3])")
imshow(abs.(S[:,:,17,3]),cmap="gray"); xticks([]); yticks([])
subplot(4,3,4)
ylabel("CS")
imshow(abs.(S_cs[:,:,17,1]),cmap="gray"); xticks([]); yticks([])
subplot(4,3,5)
imshow(abs.(S_cs[:,:,17,2]),cmap="gray"); xticks([]); yticks([])
subplot(4,3,6)
imshow(abs.(S_cs[:,:,17,3]),cmap="gray"); xticks([]); yticks([])
subplot(4,3,7)
ylabel("CSLR")
imshow(abs.(S_cslr[:,:,17,1]),cmap="gray"); xticks([]); yticks([])
subplot(4,3,8)
imshow(abs.(S_cslr[:,:,17,2]),cmap="gray"); xticks([]); yticks([])
subplot(4,3,9)
imshow(abs.(S_cslr[:,:,17,3]),cmap="gray"); xticks([]); yticks([])
subplot(4,3,10)
ylabel("CSFR")
imshow(abs.(S_csfr[:,:,17,1]),cmap="gray"); xticks([]); yticks([])
subplot(4,3,11)
imshow(abs.(S_csfr[:,:,17,2]),cmap="gray"); xticks([]); yticks([])
subplot(4,3,12)
imshow(abs.(S_csfr[:,:,17,3]),cmap="gray"); xticks([]); yticks([])


text:
using HDF5

HDF5File <: HDF5Object

HDF5Object


text:
abstract type AbstractStep end

# create a function that gets the resolved selection names (a step between resolve and select)

mutable struct Recipe{T}
    df::T
    roles::Dict{Symbol,Vector{Symbol}}
    steps::Vector{AbstractStep}
end

# @forward Recipe.df names getindex rename! resolve
Base.names(r::Recipe) = names(r.df)
Base.getindex(r::Recipe, i...) = getindex(r.steps, i...)

update!(roles::Dict, step::AbstractStep) = roles

function recipe(df; roles...)
    if isempty(roles)
        Recipe(df, Dict{Symbol,Vector{Symbol}}(), Vector{AbstractStep}())
    else
        Recipe(df, resolveroles(df, Dict(roles)), Vector{AbstractStep}())
    end
end

function fit!(r::Recipe, df)
    for step in r.steps
        fit!(step, df)
        update!(r.roles, step)
        transform!(step, df)
    end
    df
end

fit(r::Recipe, df) = fit!(r, copy(df))

function transform!(r::Recipe, df)
    for step in r.steps
        transform!(step, df)
    end
    df
end

transform(r::Recipe, df) = transform!(r, copy(df))


text:
using ActorInterfaces.Classic
using Test

struct Spawner end

struct SpawnTree
    childcount::UInt8
    depth::UInt8
end

@ctx function (me::Spawner)(msg::SpawnTree)
    if msg.depth > 0
        for i = 1:msg.childcount
            child = spawn(Spawner())
            send(child, SpawnTree(msg.childcount, msg.depth - 1))
        end
    end
    return nothing
end

const TREE_HEIGHT = 15
const TREE_SIZE = 2^(TREE_HEIGHT + 1) - 1

function run_spawnertest(lib::ActorLib)
    return @testset "ActorInterfaceTests/spawner: Spawning children and sending messages to them" begin
        @test ex_actorcount(lib) == 0
        root = Spawner()
        rootaddr = ex_spawn!(lib, root)
        @test ex_actorcount(lib) == 1
        ex_send!(lib, rootaddr, SpawnTree(2, TREE_HEIGHT))
        println("Building a tree of $TREE_SIZE actors and delivering the same amount of messages")
        ex_runtofinish(lib)
        @test ex_actorcount(lib) == TREE_SIZE
    end
end


text:
using Random
mutable struct Robot
    name
    Robot() = new(new_name())
end

const name_history = Set{String}()

function new_name()
    name = randstring('A':'Z', 2) * lpad(mod(rand(Int64), 1000), 3, "0")
    if !in(name, name_history)
        push!(name_history, name)
        return name
    end
    new_name()
end

function reset!(instance::Robot)
    instance.name = new_name()
end

function name(instance::Robot)
    instance.name
end

max_stars_repo_name: koladilip/julia-exercism

max_issues_repo_name: koladilip/julia-exercism

max_forks_repo_name: koladilip/julia-exercism


text:
const IUPACnumDict = Dict{Char,Int8}(Dict())
IUPACnumDict['A']=1
IUPACnumDict['C']=2
IUPACnumDict['G']=3
IUPACnumDict['T']=4
IUPACnumDict['U']=5
IUPACnumDict['R']=6
IUPACnumDict['Y']=7
IUPACnumDict['S']=8
IUPACnumDict['W']=9
IUPACnumDict['K']=10
IUPACnumDict['M']=11
IUPACnumDict['B']=12
IUPACnumDict['D']=13
IUPACnumDict['H']=14
IUPACnumDict['V']=15
IUPACnumDict['N']=16

IUPACnucLookup=[[true false false false false true false false true false true false true true true true]
[false true false false false false true true false false true true false true true true]
[false false true false false true false true false true false true true false true true]
[false false false true true false true false true true false true true true false true]
[false false false true true false true false true true false true true true false true]
[true false true false false true false true true true true true true true true true]
[false true false true true false true true true true true true true true true true]
[false true true false false true true true false true true true true true true true]
[true false false true true true true false true true true true true true true true]
[false false true true true true true true true true false true true true true true]
[true true false false false true true true true false true true true true true true]
[false true true true true true true true true true true true true true true true]
[true false true true true true true true true true true true true true true true]
[true true false true true true true true true true true true true true true true]
[true true true false false true true true true true true true true true true true]
[true true true true true true true true true true true true true true true true]];

function toIUPACnum(strin::String)
    return [IUPACnumDict[i] for i in strin]
end

function IUPAC_equals(a1,a2)
    return IUPACnucLookup[a1,a2]
end

"""Needleman-Wunch, with IUPAC codes, and end gaps penalized slightly less."""
function IUPAC_nuc_edit_dist(s1::String, s2::String; edge_reduction = 0.9999)
    # edge_reduction is a multiplicative score that gets multiplied to
    # the penalties along the edges, to prefer terminal gaps.
    ins_cost = -1.0  # -> : horizontal
    del_cost = -1.0  # V : vertical
    mismatch_cost = -1.0
    match_cost = 0.0

    s1numArr = toIUPACnum(s1)
    s2numArr = toIUPACnum(s2)
    s1arr = collect(s1)  # vertical
    s2arr = collect(s2)  # horizontal
    arr = zeros(length(s1arr)+1, length(s2arr)+1)
    traceArr = zeros(Int, length(s1arr), length(s2arr))

    # this will need to be generalized when we want to allow overhang
    arr[:, 1] = edge_reduction*del_cost*(0:length(s1arr))
    arr[1,:] = edge_reduction*ins_cost*(0:length(s2arr))

    for i in 2:length(s1arr)+1
        for j in 2:length(s2arr)+1
            if s1arr[i-1] == s2arr[j-1]
                diag = arr[i-1, j-1] + match_cost
            elseif IUPAC_equals(s1numArr[i-1],s2numArr[j-1])
                diag = arr[i-1, j-1] + match_cost-0.00001
            else
                diag = arr[i-1, j-1] + mismatch_cost
            end

            # to handle the lower edge penalties.
            delMult = 1
            if i == length(s1arr)+1
                delMult = edge_reduction
            end
            insMult=1
            if j == length(s2arr)+1
                insMult = edge_reduction
            end

            ins = arr[i-1, j]+(ins_cost*insMult)
            del = arr[i, j-1]+(del_cost*delMult)
            scores = [diag, del, ins]
            best = indmax(scores)
            if best == 1
                arr[i, j] = diag
            elseif best ==2
                arr[i, j] = del
            elseif best ==3
                arr[i, j] = ins
            end
            traceArr[i-1, j-1] = best
        end
    end
    alignedScore = arr[end, end]
    return alignedScore
end

"""Needleman-Wunch, with IUPAC codes, and end gaps penalized slightly less."""
function IUPAC_nuc_nw_align(s1::String, s2::String; edge_reduction = 0.9999)
    # edge_reduction is a multiplicative score that gets multiplied to
    # the penalties along the edges, to prefer terminal gaps.
    ins_cost = -1.0  # -> : horizontal
    del_cost = -1.0  # V : vertical
    mismatch_cost = -1.0
    match_cost = 0.0

    s1numArr = toIUPACnum(s1)
    s2numArr = toIUPACnum(s2)
    s1arr = collect(s1)  # vertical
    s2arr = collect(s2)  # horizontal
    arr = zeros(length(s1arr)+1, length(s2arr)+1)
    traceArr = zeros(Int, length(s1arr), length(s2arr))

    # this will need to be generalized when we want to allow overhang
    arr[:, 1] = edge_reduction*del_cost*(0:length(s1arr))
    arr[1,:] = edge_reduction*ins_cost*(0:length(s2arr))

    for i in 2:length(s1arr)+1
        for j in 2:length(s2arr)+1
            if s1arr[i-1] == s2arr[j-1]
                diag = arr[i-1, j-1] + match_cost
            elseif IUPAC_equals(s1numArr[i-1],s2numArr[j-1])
                diag = arr[i-1, j-1] + match_cost-0.00001
            else
                diag = arr[i-1, j-1] + mismatch_cost
            end

            # to handle the lower edge penalties.
            delMult = 1
            if i == length(s1arr)+1
                delMult = edge_reduction
            end
            insMult=1
            if j == length(s2arr)+1
                insMult = edge_reduction
            end

            ins = arr[i-1, j]+(ins_cost*insMult)
            del = arr[i, j-1]+(del_cost*delMult)
            scores = [diag, del, ins]
            best = indmax(scores)
            if best == 1
                arr[i, j] = diag
            elseif best ==2
                arr[i, j] = del
            elseif best ==3
                arr[i, j] = ins
            end
            traceArr[i-1, j-1] = best
        end
    end
    alignedScore = arr[end, end]

    # return arr

    # the trace endpoing will need to be generalized if we want to
    # allow overhang.
    trI, trJ = length(s1arr), length(s2arr)

    # First compute the trace running backwards. Initialized to the
    # maximum size. With unit penalties, is it possible to tell how
    # big this should be in advance?
    backtrace = Array{Int}(length(s1arr) + length(s2arr))
    # If you hit any boundary, you have to run all the way to the side.
    btInd = 1
    while (trI > 0) && (trJ > 0)
        backtrace[btInd] = traceArr[trI, trJ]
        if backtrace[btInd] == 1
            trI += -1
            trJ += -1
        elseif backtrace[btInd] == 2
            trJ += -1
        elseif backtrace[btInd] == 3
            trI += -1
        end
        btInd += 1
    end

    # If you hit the boundaries not at the top left corner.
    while trI > 0
        backtrace[btInd] = 3
        btInd += 1
        trI += -1
    end

    while trJ > 0
        backtrace[btInd] = 2
        btInd += 1
        trJ += -1
    end

    backtrace = backtrace[1:btInd-1]

    # This will need to be generalized to work on non-strings. Not
    # important.
    ali1arr = Array{Char}(length(backtrace))
    ali2arr = Array{Char}(length(backtrace))
    ind1 = 1
    ind2 = 1
    for i in 1:length(backtrace)
        tr = backtrace[length(backtrace)-(i-1)]
        if tr == 1
            ali1arr[i] = s1arr[ind1]
            ali2arr[i] = s2arr[ind2]
            ind1 += 1
            ind2 += 1
        elseif tr == 2
            ali2arr[i] = s2arr[ind2]
            ali1arr[i] = '-'
            ind2 += 1
        elseif tr == 3
            ali1arr[i] = s1arr[ind1]
            ali2arr[i] = '-'
            ind1 += 1
        end
    end
    return join(ali1arr), join(ali2arr), alignedScore
end

function banded_edit_dist(s1::String, s2::String; edge_reduction = 0.99, band_coeff = 1)
    # calculate band width
    avg_len = (length(s1) + length(s2)) / 2
    bandwidth = Int(ceil(band_coeff * sqrt(avg_len)))
    # edge_reduction is a multiplicative score that gets multiplied to
    # the penalties along the edges, to prefer terminal gaps.
    ins_cost = -1.0  # -> : horizontal
    del_cost = -1.0  # V : vertical
    mismatch_cost = -1.0
    match_cost = 1.0

    s1arr = collect(s1)  # vertical
    s2arr = collect(s2)  # horizontal

    # dimension difference added to band length on one side
    # to give padding for large sequence length differences
    dim_diff = length(s1arr) - length(s2arr)
    maxlen = max(length(s1), length(s2))
    arr = zeros(maxlen + 1, 2*bandwidth + 1 + abs(dim_diff))
    traceArr = zeros(Int, size(arr))
    for i in 1:bandwidth+1
        add_to_band!(arr, edge_reduction*del_cost*(i-1), i, 1, bandwidth, dim_diff)
        add_to_band!(arr, edge_reduction*ins_cost*(i-1), 1, i, bandwidth, dim_diff)
    end

    for i in 2:length(s1arr)+1
        # keep j within 2*bandwidth + 1 band of i, but keep bandwidth sized buffer on edges
        lo = 0
        # to keep a big enough padding around the edges when there's a big length mismatch
        if i <= bandwidth + 1 + max(0, dim_diff)
            lo = 2
        else
            lo = i - bandwidth - max(0, dim_diff)
        end
        for j in lo:length(s2arr)+1
            # keep j within band
            if !(in_band(i, j, bandwidth, dim_diff))
                break
            end

            if in_band(i-1, j-1, bandwidth, dim_diff) && s1arr[i-1] == s2arr[j-1]
                diag = get_band_val(arr, i-1, j-1, bandwidth, dim_diff) + match_cost
            else
                diag = get_band_val(arr, i-1, j-1, bandwidth, dim_diff) + mismatch_cost
            end

            # to handle the lower edge penalties.
            delMult = (i == length(s1arr)+1) ? edge_reduction : 1
            insMult = (j == length(s2arr)+1) ? edge_reduction : 1

            ins = in_band(i-1, j, bandwidth, dim_diff) ?
                    (get_band_val(arr, i-1, j, bandwidth, dim_diff)+(ins_cost*insMult)) : NaN
            del = in_band(i, j-1, bandwidth, dim_diff) ?
                    (get_band_val(arr, i, j-1, bandwidth, dim_diff)+(del_cost*delMult)) : NaN

            scores = [diag, del, ins]
            best = indmax(scores)
            add_to_band!(arr, scores[best], i, j, bandwidth, dim_diff)
            add_to_band!(traceArr, best, i-1, j-1, bandwidth, dim_diff)
        end
    end
    alignedScore = get_band_val(arr, length(s1arr)+1, length(s2arr)+1, bandwidth, dim_diff)

    return alignedScore
end

#Indmax in 2D, returning the coordinates
function maxcoord(mat)
    return ind2sub(mat,indmax(mat))
end

#THIS PRIMER MATCHING FUNCTION DOESN'T ASSUME WE KNOW PRIMER PAIRS
function matchPrimer(seq,forwardPrimers,revPrimers,forwardNames,revNames)
    fwdseq = seq;
    revseq = reverse_complement(fwdseq);
    fwdSeqFwdPrim = [IUPAC_nuc_edit_dist(fwdseq[1:length(prim)],prim) for prim in forwardPrimers];
    revSeqFwdPrim = [IUPAC_nuc_edit_dist(revseq[1:length(prim)],prim) for prim in forwardPrimers];
    fwdSeqRevPrim = [IUPAC_nuc_edit_dist(fwdseq[1:length(prim)],prim) for prim in revPrimers];
    revSeqRevPrim = [IUPAC_nuc_edit_dist(revseq[1:length(prim)],prim) for prim in revPrimers];

    fwdDirectionScores = zeros(length(forwardPrimers),length(revPrimers))
    revDirectionScores = zeros(length(forwardPrimers),length(revPrimers))

    for i in 1:length(forwardPrimers)
        for j in 1:length(revPrimers)
            fwdDirectionScores[i,j]=fwdSeqFwdPrim[i]+revSeqRevPrim[j];
            revDirectionScores[i,j]=revSeqFwdPrim[i]+fwdSeqRevPrim[j];
        end
    end

    fwd,back = maximum(fwdDirectionScores),maximum(revDirectionScores)

    #fwdDirectionScores = fwdSeqFwdPrim.+revSeqRevPrim;
    #revDirectionScores = revSeqFwdPrim.+fwdSeqRevPrim;
    fwd,back = maximum(fwdDirectionScores),maximum(revDirectionScores)
    fwdDir = fwd>back
    if fwdDir
        primInd = maxcoord(fwdDirectionScores)
        worstScore = minimum([fwdSeqFwdPrim[primInd[1]],revSeqRevPrim[primInd[2]]])
        retseq = fwdseq
    else
        primInd = maxcoord(revDirectionScores)
        worstScore = minimum([revSeqFwdPrim[primInd[1]],fwdSeqRevPrim[primInd[2]]])
        retseq = revseq
    end
    return primInd,forwardNames[primInd[1]]*"_"revNames[primInd[2]],retseq,round(worstScore),(!fwdDir)
end

function deMux(seqs,forwardPrimers,reversePrimers, forwardNames,reverseNames)
    matched = [matchPrimer(seq,forwardPrimers,reversePrimers,forwardNames,reverseNames) for seq in seqs];
    return [[i[2],i[3], i[5], i[4]] for i in matched]
end

function demux_fastx(fwdPrimers, revPrimers; minLen = 2300, errorRate = 0.01,
    path = "/Users/vrkumar/Downloads/",
    filename = "LP97.fastq",
    mismatchThresh = 2
)
    forwardNames = [i[1] for i in fwdPrimers];
    reverseNames = [i[1] for i in revPrimers];
    forwardPrimers = [i[2] for i in fwdPrimers];
    reversePrimers = [i[2] for i in revPrimers];

    newName = path*filename*"_"*string(minLen)*"_"*string(errorRate)*".fastq"
    usearch_filter(path*filename,newName, errorRate=errorRate,minLength=minLen,labelPrefix="seq",errorOut = true)

    if filename[end] == 'a'
        names, seqs = read_fasta_with_names(newName)
        phreds = nothing;
    else
        seqs, phreds, names = read_fastq(newName);
    end

    deMuxed = deMux(seqs,forwardPrimers,reversePrimers,forwardNames,reverseNames);

    primerPairs = [i[1] for i in deMuxed];
    primerMatchedSeqs = [i[2] for i in deMuxed];
    reversed = [i[3] for i in deMuxed];
    scores = [i[4] for i in deMuxed]

    if phreds != nothing
        out_phreds = [reversed[i] ? reverse(phreds[i]) : phreds[i] for i in 1:length(reversed)]
    else
        out_phreds = [0 for i in 1:length(reversed)]
    end

    out = Dict()
    for i in 1:length(primerPairs)
        if scores[i] > mismatchThresh
            continue
        end
        if !(haskey(out, primerPairs[i]))
            out[primerPairs[i]] = [[],[],[]]
        end
        push!(out[primerPairs[i]][1], primerMatchedSeqs[i])
        push!(out[primerPairs[i]][2], out_phreds[i])
        push!(out[primerPairs[i]][3], names[i])
    end

    for p in keys(out)
        if phreds != nothing
            write_fastq(path*p*".fastq", out[p][1], Array{Array{Int8,1},1}(out[p][2]), names = out[p][3])
        else
            write_fasta(path*p*".fastq", out[p][1], names = out[p][3])
        end
    end

end


text:
module HCModelKit

using OrderedCollections: OrderedDict
using StaticArrays: @SVector, @SMatrix

import LinearAlgebra: det, dot
import Latexify

export Expression, Constant, Variable, Operation

export @var,
       @unique_var,
       evaluate,
       subs,
       variables,
       differentiate,
       monomials,
       Compiled,
       CompiledSystem,
       CompiledHomotopy,
       compile,
       interpret,
       interpreted,
       System,
       Homotopy,
       evaluate,
       evaluate!,
       evaluate_gradient,
       evaluate_jacobian,
       evaluate_jacobian!,
       jacobian,
       jacobian!,
       dt,
       dt!,
       dt_jacobian,
       dt_jacobian!

include("expression.jl")
include("codegen.jl")

end # module


text:
insts = OrderedDict()
insts["minimal"] = [
    ((3, 2, true, 0),),
    ((3, 2, false, 0),),
    ((3, 2, true, 10),),
    ((3, 2, false, 10),),
    ]
insts["fast"] = [
    ((5, 3, true, 0),),
    ((5, 3, false, 0),),
    ((5, 3, true, 10),),
    ((5, 3, false, 10),),
    ((30, 10, true, 0),),
    ((30, 10, false, 0), (default_tol_relax = 1000,)),
    ((30, 10, true, 10),),
    ((30, 10, false, 10),),
    ]
insts["various"] = [
    ((10, 3, true, 10),),
    ((10, 3, false, 10),),
    ((20, 6, true, 10),),
    ((20, 6, false, 10),),
    ((40, 12, true, 10),),
    ((40, 12, false, 10),),
    ]
return (SparsePCANative, insts)


text:
# Autogenerated wrapper script for Mesa_jll for i686-w64-mingw32
export opengl32sw

using Zlib_jll
JLLWrappers.@generate_wrapper_header("Mesa")
JLLWrappers.@declare_library_product(opengl32sw, "opengl32sw.dll")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll)
    JLLWrappers.@init_library_product(
        opengl32sw,
        "bin\\opengl32sw.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    BondPricePs(y,c0,t,FaceValue=1)

Calculate a single coupon bond price from the spot rate curve and the payment stream
(arbitrary periods and coupons), m periods

# Input
- `y::Vector or Number`:    m vector or scalar, effective spot rates, e.g. [0.07;0.08], m periods
- `c0::Vector or Number`:   m vector or scalar, coupon rate, e.g. 0.06 or 0.09/2, m periods
- `t::Vector`:              m vector, time to the cash flows
- `FaceValue::Number`:      (optional) face value, [1]

# Output
- `Q::Number`:  scalar, bond price (eg. 1.01)

# Notice
The face value is paid at the same time as the last coupon.
Dates (t) should be expressed as fractions of the period.

The following is calculated:
```
Q =  c[1]/(1+y(1))^t(1) + c[2]/(1+y(2))^t(2) + ... + (1+c[m])/(1+y(m))^t(m)
```
The coupons and interest rates (optionally) depend on time to payment, so `y` and `c0`
can be either scalars or m-element vectors.

# Example
```
y = [0.0626;0.07]
c = [0.09; 0.09]
t = 1:2
gives
Q = 1.0367
```

Paul Soderlind (Paul.Soderlind@unisg.ch)

"""
function BondPricePs(y,c0,t,FaceValue=1)

  isa(c0,Number) ? c = fill(c0,length(t)) : c = copy(c0)
  c[end] = c[end] + FaceValue

  cdisc  = c./((1.0.+y).^t)        #c/(1+y(1))^t1 + c/(1+y(2))^t2 + ...+ c/(1+y(m))^tm
  Q      = sum(cdisc)              #price (sum(.,1)) to handle also m=1

  return Q

end
#-----------------------------------------------------------------------------

#-----------------------------------------------------------------------------
"""
    BondPriceQDyPs(y,c0,t,FaceValue=1)

Calculate coupon bond price and derivative.

# Input
- see BondPricePs

# Output
- `Q::Number`:      scalar, bond price (eg. 1.01)
- `dQ_dy::Number`:  scalar, derivative of bond price wrt y (valid only when y is a scalar)

"""
function BondPriceQDyPs(y,c0,t,FaceValue=1)

  isa(c0,Number) ? c = fill(c0,length(t)) : c = copy(c0)
  c[end] = c[end] + FaceValue

  cdisc  = c./((1.0.+y).^t)        #c/(1+y(1))^t1 + c/(1+y(2))^t2 + ...+ c/(1+y(m))^tm
  Q      = sum(cdisc)              #price (sum(.,1)) to handle also m=1

  dcdisc = cdisc .* (-t./(1.0.+y))        #derivative wrt ytm
  dQ_dy  = sum(dcdisc)                    #only strictly valid when y is a scalar

  return Q, dQ_dy

end
#-----------------------------------------------------------------------------

#-----------------------------------------------------------------------------
"""
    BondYieldToMatPs(Q,c0,t,FaceValue=1,method=1,yLH=[-0.1;0.5],tol=1e-7)

Calculate yield to maturity from bond price (several methods available).
Can also be used for general IRR calculations. Works with arbitrary coupon periods.

# Input
- `Q::Number`:              scalar
- `c::Vector or Number`:    see BondPricePs
- `t::Vector or Number`:    see BondPricePs
- `FaceValue::Number`:      optional, scalar or n vector, face value, [1]
- `method::Number`:         optional, 1: Newton-Raphson; 2: bisection [1]
- `yLH::Number`:            optional, if method==1: scalar, yLH[1] is initial guess of roots
                            if method==2: 2 vector, yLH[1] is lower boundary and yLH[2] upper boundary
- `tol::Number`:            optional, scalar, convergence criterion for y, [1e-7]

# Output
- `y::Number`           effective yield to maturity (per period)

# Notice
- method 1 (Newton-Raphson) is pretty fast.
- method 2 (bisection) is a bit slow, but very robust.
- To use for IRR calculations, set Q=0 and FaceValue=0
- For Newton-Raphson, we could use (c+(FaceValue-B)/n)/((FaceValue+B)/2) as a good starting value

# Requires
- BondPricePs, BondPriceQDyPs

"""
function BondYieldToMatPs(Q,c0,t,FaceValue=1,method=1,yLH=[-0.1;0.5],tol=1e-7)

  if method == 1                   #Newton-Raphson

    y = yLH[1]
    Dy = 1e+198
    while abs(Dy) > tol
      (F0,J0) = BondPriceQDyPs(y,c0,t,FaceValue)    #F(y0), derivative
      Dy = (Q - F0)/J0
      y  = y + Dy
    end

  elseif method == 2                   #bisection

    (yL,yH) = (yLH[1],yLH[2])                #lower and upper boundary for yield
    if yL > yH
      println("warning: lower bound greater than upper bound")
    end
    F0 = BondPricePs(yL,c0,t,FaceValue)     #create starting value, so [yL,yH] brackets the roots
    while F0 < Q                             #as long as theoretical < actual price
      yL = yL - 0.01
      F0 = BondPricePs(yL,c0,t,FaceValue)
    end
    F0 = BondPricePs(yH,c0,t,FaceValue)
    while F0 > Q
      yH = yH + 0.01
      F0 = BondPricePs(yH,c0,t,FaceValue)
    end

    y = (yL + yH)/2
    while (yH-yL) > tol                    #iteration loop
      y = (yL + yH)/2                      #mid point for yield
      F0 = BondPricePs(y,c0,t,FaceValue)  #price at guessed yield
      if F0 >= Q                           #logical where F0 >= Q  (decreasing function)
        yL = y                             # => root must be higher than y
      else
        yH = y                             # -> root must be lower than y
      end
    end

  end                                      #end different methods

  return y

end
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
"""
    BondDurationPs(Q,c0,t,FaceValue=1,method=1,yLH=[-0.1;0.5],tol=1e-7)

Calculate bond durations. See also BondYieldToMatPs

# Input
- see BondYieldToMatPs

# Output
- `D::Number`:        scalar, (dollar) duration
- `Da::Number`:       scalar, adjusted (or modified) duration
- `Dma::Numberc`:     scalar, Macaulays duration
- `ytm::Number`:      scalar

# Notice
The bond price satisfies (with yield expressed per period)

Q =  c/(1+y)^t(1) + c/(1+y)^t(2) + ... + (1+c)/(1+y)^t(m)

where y is the yield to maturity, c is per period coupon
and Macaulay's duration is

D =  [t(1)*c/(1+y)^t(1) + t(2)*c/(1+y)^t(2) + ... + t(m)*(1+c)/(1+y)^t(m)]/Q

# Requires
- BondPricePs, BondYieldToMatPs

"""
function BondDurationPs(Q,c0,t,FaceValue=1,method=1,yLH=[-0.1;0.5],tol=1e-7)

  y = BondYieldToMatPs(Q,c0,t,FaceValue,method,yLH,tol)  #ytm

  (_,Dy) = BondPriceQDyPs(y,c0,t,FaceValue)
  D      = -Dy                                    #duration
  Da     = D / Q                                  #adjusted duration
  Dmac   = D*(1+y)/ Q                             #Macaulays duration

  return D, Da, Dmac, y

end
#------------------------------------------------------------------------------

max_stars_repo_name: PaulSoderlind/PaulSoderlindCode

max_issues_repo_name: PaulSoderlind/PaulSoderlindCode

max_forks_repo_name: PaulSoderlind/PaulSoderlindCode


text:
# ------------------------------------------------------------------------------------------------------------------------------------------------
struct BackPath

    sides::NTuple{2,HorizonSide}
    path::Vector{Int}

end
function BackPath(robot::Robot, sides::NTuple{2,HorizonSide}=(Sud,West))
    local path=Int[]
    while ( not(isborder(robot,sides[1])) || not(isborder(robot,sides[2])) )
        for side in sides
            push!(path, movements!(robot,side))
        end
    end
    return new(reverse(inverse.(sides)),path)
end
# ------------------------------------------------------------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------------------------------------------------------------
function moving_back_to_start!(robot::Robot, vector::BackPath)
    i = 1
    for n in vector.path # n -- количество шагов робота в одном направлении
        movements!(robot, vector.sides[i])
        i = i % 2 + 1
    end
end
# ------------------------------------------------------------------------------------------------------------------------------------------------

back_bath = BackPath(r) # переносит робота в юго-западный угол. Возвращает массив с reverse путем, который робот проделал для того, чтобы оказаться в юго-заподном углу.
moving_back_to_start!(r, back_path) # возвращает робота из юго-заподного угла в начальное проложение

max_stars_repo_name: Igetsuga/AlgorithmicsHomework

max_issues_repo_name: Igetsuga/AlgorithmicsHomework

max_forks_repo_name: Igetsuga/AlgorithmicsHomework


text:
using DifferentialEquations
using ITensors
using ITensorTDVP
using KrylovKit
using LinearAlgebra
using Test

include(joinpath(pkgdir(ITensorTDVP), "examples", "03_models.jl"))
include(joinpath(pkgdir(ITensorTDVP), "examples", "03_solvers.jl"))

# Functions need to be defined in global scope (outside
# of the @testset macro)

ω₁ = 0.1
ω₂ = 0.2

ode_alg = Tsit5()
ode_kwargs = (; reltol=1e-8, abstol=1e-8)

ω⃗ = [ω₁, ω₂]
f⃗ = [t -> cos(ω * t) for ω in ω⃗]

function ode_solver(H⃗₀, time_step, ψ₀; kwargs...)
  return ode_solver(
    -im * TimeDependentSum(f⃗, H⃗₀),
    time_step,
    ψ₀;
    solver_alg=ode_alg,
    ode_kwargs...,
    kwargs...,
  )
end

krylov_kwargs = (; tol=1e-8, eager=true)

function krylov_solver(H⃗₀, time_step, ψ₀; kwargs...)
  return krylov_solver(
    -im * TimeDependentSum(f⃗, H⃗₀), time_step, ψ₀; krylov_kwargs..., kwargs...
  )
end

@testset "Time dependent Hamiltonian" begin
  n = 4
  J₁ = 1.0
  J₂ = 0.1

  time_step = 0.1
  time_stop = 1.0

  nsite = 2
  maxdim = 100
  cutoff = 1e-8

  s = siteinds("S=1/2", n)
  ℋ₁₀ = heisenberg(n; J=J₁, J2=0.0)
  ℋ₂₀ = heisenberg(n; J=0.0, J2=J₂)
  ℋ⃗₀ = [ℋ₁₀, ℋ₂₀]
  H⃗₀ = [MPO(ℋ₀, s) for ℋ₀ in ℋ⃗₀]

  ψ₀ = complex.(MPS(s, j -> isodd(j) ? "↑" : "↓"))

  ψₜ_ode = tdvp(ode_solver, H⃗₀, time_stop, ψ₀; time_step, maxdim, cutoff, nsite)

  ψₜ_krylov = tdvp(krylov_solver, H⃗₀, time_stop, ψ₀; time_step, cutoff, nsite)

  ψₜ_full, _ = ode_solver(prod.(H⃗₀), time_stop, prod(ψ₀))

  @test norm(ψ₀) ≈ 1
  @test norm(ψₜ_ode) ≈ 1
  @test norm(ψₜ_krylov) ≈ 1
  @test norm(ψₜ_full) ≈ 1

  ode_err = norm(prod(ψₜ_ode) - ψₜ_full)
  krylov_err = norm(prod(ψₜ_krylov) - ψₜ_full)

  @test krylov_err > ode_err
  @test ode_err < 1e-3
  @test krylov_err < 1e-3
end

nothing


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

using Plots;

"""
    democylind()

Demonstration of the construction of a cylinder.

# Examples:
```julia
julia> democylind()
```
"""
function democylind()
srf = nrbcylind(3.0,1.0,[0.0;0.0;0.0],deg2rad(270),deg2rad(180));
nrbplot(srf,[20;20],
        title = "Cylindrical section by extrusion of a circular arc.",
        legend = false);
end # democylind


text:
# This file was generated, do not modify it. # hide
#hideall
using Markdown

mdC_code = Markdown.htmlesc(raw"""#include <stdio.h>
int main(){
    printf("Hello Pikachu\n");
    return 0;
}""")
mdfile = joinpath(dirname(@OUTPUT), "how2embeddC.md")
open(mdfile,"w") do f
    print(f, mdC_code)
end

C_code=raw"""
#include <stdio.h>
int main(){
    printf("Hello Pikachu\n");
    return 0;
}
"""

exefile = tempname()

#=
This trick is taken from

https://discourse.julialang.org/t/how-to-make-a-c-function-compiled-by-myself-available-to-ccall/7972/26
=#

open(`gcc -Wall -O2 -march=native -xc -o $exefile -`, "w") do f
    print(f, C_code)
end

run(`$exefile`)


text:
using BenchmarkTools
using Bobby

b = Bobby.setBoard()
SUITE = BenchmarkGroup()
SUITE["perft"] = BenchmarkGroup()
SUITE["perft"]["starting position"] = @benchmarkable Bobby.perft(b,4)

fens = ["r3k2r/p1ppqpb1/bn2pnp1/3PN3/1p2P3/2N2Q1p/PPPBBPPP/R3K2R w KQkq - 0 0",
        "8/2p5/3p4/KP5r/1R3p1k/8/4P1P1/8 w - - 0 0",
        "r3k2r/Pppp1ppp/1b3nbN/nP6/BBP1P3/q4N2/Pp1P2PP/R2Q1RK1 w kq - 0 1",
        "r2q1rk1/pP1p2pp/Q4n2/bbp1p3/Np6/1B3NBn/pPPP1PPP/R3K2R b KQ - 0 1",
        "rnbq1k1r/pp1Pbppp/2p5/8/2B5/8/PPP1NnPP/RNBQK2R w KQ - 1 8",
        "r4rk1/1pp1qppp/p1np1n2/2b1p1B1/2B1P1b1/P1NP1N2/1PP1QPPP/R4RK1 w - - 0 10",
        "n1n5/PPPk4/8/8/8/8/4Kppp/5N1N b - - 0 1",
        "8/7p/p5pb/4k3/P1pPn3/8/P5PP/1rB2RK1 b - d3 0 28",
        "rnbqkb1r/ppppp1pp/7n/4Pp2/8/8/PPPP1PPP/RNBQKBNR w KQkq f6 0 3",
        "8/2p5/3p4/KP5r/1R3p1k/8/4P1P1/8 w - - 0 1"]

for (i, f) in enumerate(fens)
    b = Bobby.loadFen(f)
    SUITE["perft"]["perft"*string(i)] = @benchmarkable Bobby.perft(b, 4)
end


text:
const Prime = UInt64

struct Primes
    rangei :: Int64
    primebits :: BitArray{1}
    function Primes(n :: Int64)
        if n < 3
          if n < 2 return new(-1, falses(0)) # no elements
          else return new((0, trues(0))) end # n = 2: meaning is 1 element of 2
        end
        limi :: Int = (n - 1) ÷ 2 # calculate the required array size
        isprimes :: BitArray = trues(limi)
        @inbounds(
        for i in 1:limi
            p = i + i + 1
            start = (p * p - 1) >>> 1 # shift/divide if LLVM doesn't optimize
            if start > limi
                return new(limi, isprimes)
            end
            if isprimes[i]
                for j in start:p:limi
                  isprimes[j] = false
                end
            end
        end)
    end
end

Base.eltype(::Type{Primes}) = Prime

function Base.length(P::Primes)::Int64
    if P.rangei < 0 return 0 end
    return 1 + count(P.primebits)
end

function Base.iterate(P::Primes, state::Int = 0)::
                                        Union{Tuple{Prime, Int}, Nothing}
    lmt = P.rangei
    if state > lmt return nothing end
    if state <= 0 return (UInt64(2), 1) end
    let
        prmbts = P.primebits
        i = state
        @inbounds(
        while i <= lmt && !prmbts[i] i += 1 end)
        if i > lmt return nothing end
        return (i + i + 1, i + 1)
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode