text:
function rotation_matrix(theta, ax)
    x,y,z = LinearAlgebra.normalize(ax)
    c = cos(theta)
    s = sin(theta)
    R = [(c+x^2*(1-c)) (x*y*(1-c)-s*z) (z*x*(1-c)+s*y);
         (x*y*(1-c)+s*z) (c+y^2*(1-c)) (z*y*(1-c)-s*x);
         (x*z*(1-c)-s*y) (y*z*(1-c)+s*x) (c+z^2*(1-c))]
    return R
end
function parallel_transport(coordinate::Matrix{T_Number}) where {T_Number<:Number}
    r = copy(coordinate)
    n_points = size(r, 1)
    if size(r, 2) == 2
        r = cat(r, zeros(eltype(r), size(r, 1)), dims = 2)
    end
    t = 1.0:size(r,1)
    dt = (maximum(t) - minimum(t)) / length(r)
    dr = gradient(r)
    ddr = gradient(dr)
    dddr = gradient(ddr)

    T = similar(dr)
    N = similar(dr)
    B = similar(dr)
    T[1, :] = LinearAlgebra.normalize(dr[1, :])
    B[1, :] = LinearAlgebra.normalize([T[1,2], -T[1,1], 0 ])
    N[1, :] = LinearAlgebra.normalize(cross(B[1, :], T[1, :]))
    for i = 2:n_points
        T[i, :] = LinearAlgebra.normalize(dr[i, :])
        b_ = cross(T[i-1,:],T[i,:])
        if isapprox(norm(b_),0)
            B[i, :] = B[i-1,:]
        else
            b_ = LinearAlgebra.normalize(b_)
            phi = acos(dot(T[i-1,:],T[i,:]))
            R = rotation_matrix(phi,b_)
            B[i, :] = R * B[i-1,:]
        end
        N[i, :] = LinearAlgebra.normalize(cross(B[i, :], T[i, :]))
    end

    ds = Array{eltype(dr),1}(undef, n_points)
    curvature = similar(ds)
    torsion = similar(ds)

    for i in eachindex(ds)
        selected_dr = dr[i, :]
        selected_ddr = ddr[i, :]
        selected_dddr = dddr[i, :]
        ds[i] = norm(selected_dr)

        dr_cross_ddr = cross(selected_dr, selected_ddr)
        curvature_num = norm(dr_cross_ddr) # || dr X ddr ||
        curvature_denom = norm(selected_dr)^3 # || dr ||^3
        curvature[i] = curvature_num / curvature_denom

        torsion_num = dot(dr_cross_ddr, selected_dddr) # (dr X ddr) . dddr
        torsion_denom = norm(dr_cross_ddr)^2 # || dr X ddr ||^2
        torsion[i] = torsion_num / torsion_denom
    end
    speed = ds ./ dt
    return T,N,B,curvature,torsion,speed
end


text:
using ApproxOperator

P = 1000.0
Ẽ = 3e6
# ν̃ = 0.3
ν̃ = 0.49999
E = Ẽ/(1-ν̃^2)
ν = ν̃/(1-ν̃)
L = 48.0
D = 12.0
I = D^3/12
EI = E*I

u = (x,y,z) -> (
    -P*y/6/EI*((6*L-3x)*x + (2+ν)*(y^2-D^2/4)),
     P/6/EI*(3*ν*y^2*(L-x) + (4+5*ν)*D^2*x/4 + (3*L-x)*x^2),
    -P/EI*(L-x)*y,
    -P/6/EI*((6*L-3*x)*x + (2+ν)*(3*y^2-D^2/4)),
     P/6/EI*((6*L-3*x)*x - 3*ν*y^2 + (4+5*ν)*D^2/4),
     P/EI*(L-x)*y*ν
)
t = (x,y,z) -> (0.0,P/2/I*(D^2/4-y^2))
b = (x,y,z) -> (0.0,0.0)
n = (x,y,z) -> (1.0,0.0,1.0)

aps = import_msh("./msh/cantilever.msh")
nₚ = length(aps["Domain"][1].nodes)

ops = [
    # 1
    PlaneStress_Ω(b,E,ν),
    # 2
    PlaneStrain_Ωᵛ(Ẽ,ν̃),
    # 3
    PlaneStrain_Ωᵈ(b,Ẽ,ν̃),
    # 4
    PlaneStress_Γᵗ(t),
    # 5
    PlaneStress_Γᵍ_penalty(u,n,E*1e7),
    # 6
    HₑError_PlaneStress(u,E,ν),
    # 7
    VTKExport(
                filename="vtk/cantilever_locking.vtk",
                topic="incompressible test for cantilever beam problem",
                pointdata=Dict{String,Symbol}("disp"=>:u_2D),
                celldata=Dict{String,Symbol}("stress"=>:σ_PlaneStrain),
                parameters=Dict{String,Float64}("E"=>Ẽ,"ν"=>ν̃)
            )
]

## locking
k = zeros(2nₚ,2nₚ)
f = zeros(2nₚ)

ops[1](aps["Domain"],k,f)
ops[4](aps["Traction"],f)
ops[5](aps["EssentialBC"],k,f)

d = k\f

Hₑ_Locking,L₂_Locking = ops[6](aps["Domain"],d)

ops[7](aps["Domain"],d)

## locking-free
k = zeros(2nₚ,2nₚ)
f = zeros(2nₚ)

ops[3](aps["Domain"],k,f)
set_integration_rule!(aps["Domain"],:QuadGI1)
ops[2](aps["Domain"],k,f)
ops[4](aps["Traction"],f)
ops[5](aps["EssentialBC"],k,f)

d = k\f

set_integration_rule!(aps["Domain"],:QuadGI2)
Hₑ_Free,L₂_Free = ops[6](aps["Domain"],d)

ops[7].filename = "vtk/cantilever_free.vtk"
ops[7](aps["Domain"],d)

ext: jl

max_stars_repo_name: jmwjc/ApproxOperator

max_issues_repo_name: jmwjc/ApproxOperator

max_forks_repo_name: jmwjc/ApproxOperator


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
A = 10(1:9) .+ (1:9)'

# %%
i = 4
@views A00 = A[begin:i-1, begin:i-1]

# %%
@views a10 = A[begin:i-1, i:i]

# %%
@views A02 = A[begin:i-1, i+1:end]

# %%
@views A10 = A[i:i, begin:i-1]

# %%
@views a11 = A[i:i, i:i]

# %%
@views a12 = A[i:i, i+1:end]

# %%
@views A20 = A[i+1:end, begin:i-1]

# %%
@views a21 = A[i+1:end, i:i]

# %%
@views A22 = A[i+1:end, i+1:end]

# %%
?"subarray"

# %%
?view

# %%
?@view

# %%
?@views

# %%

max_stars_repo_path: 0016/How to split a matrix.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0016/How to split a matrix.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0016/How to split a matrix.jl

max_forks_repo_name: genkuroki/public


text:
import Pkg
Pkg.activate("..")

using Distributions: Beta, Bernoulli
using CrimsonSkyline
using Random: seed!
using Logging
using StatsBase: mean

seed!(2021)

function coin_model(t::Trace, data::Vector{Bool})
    bias = sample(t, "bias", Beta(3.0, 3.0))
    coin = Bernoulli(bias)
    for (i, d) in enumerate(data)
        observe(t, "flip $i", coin, d)
    end
end

function main()
    @warn "~~~Coin flip problem~~~"
    data = [true, true, false, false, true, true, true, true]
    @info "Demonstrating trace structure"
    t = trace()
    coin_model(t, data[1:2])
    @info "Trace with data subset: $t"
    @info "Running default mh inference (no user tuning required)"
    @time inference_results = mh(coin_model; params = (data,))
    # posterior parameters:
    # alpha' = alpha + sum x_i = 3.0 + 6 = 9.0
    # beta' = beta + n - sum x_i = 3.0 + 8 - 6 = 5.0
    # posterior = Beta(alpha', beta') = Beta(9.0, 5.0)
    # posterior mean = alpha' / (alpha' + beta') = 0.6429
    true_posterior_mean = 0.6429
    @info "True posterior bias: $true_posterior_mean"
    @info "Mean posterior bias estimate: $(mean(inference_results, "bias"))"
    plot_marginal(inference_results, "bias", "plots/coin_flip", "bias.png")
end

main()


text:
##############################################################################
#
# TD-IDF
#
##############################################################################

tf_idf{T <: Real}(dtm::Matrix{T}) = tf_idf!(dtm, Array(Float64, size(dtm)...))

tf_idf{T <: Real}(dtm::SparseMatrixCSC{T}) =  tf_idf!(dtm, sparse(Int[], Int[], 0.0, size(dtm)...))

tf_idf!{T <: Real}(dtm::AbstractMatrix{T}) = tf_idf!(dtm, dtm)

tf_idf!{T <: Real}(dtm::SparseMatrixCSC{T}) = tf_idf!(dtm, dtm)

tf_idf(dtm::DocumentTermMatrix) = tf_idf(dtm.dtm)

# This does not make sense, since DocumentTermMatrix is based on an array of integers
#tf_idf!(dtm::DocumentTermMatrix) = tf_idf!(dtm.dtm)

# The second Matrix will be overwritten with the result
# Will work correctly if dtm and tfidf are the same matrix
function tf_idf!{T1 <: Real, T2 <: AbstractFloat}(dtm::AbstractMatrix{T1}, tfidf::AbstractMatrix{T2})
    n, p = size(dtm)

    # TF tells us what proportion of a document is defined by a term
    for i in 1:n
        words_in_document = 0
        for j in 1:p
            words_in_document += dtm[i, j]
        end
        tfidf[i, :] = dtm[i, :] ./ words_in_document
    end

    # IDF tells us how rare a term is in the corpus
    documents_containing_term = vec(sum(dtm .> 0, 1))
    idf = log(n ./ documents_containing_term)

    # TF-IDF is the product of TF and IDF
    for i in 1:n
        for j in 1:p
           tfidf[i, j] = tfidf[i, j] * idf[j]
        end
    end

    return tfidf
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ConnectivityIssue <: SwaggerModel
    origin::Any # spec type: Union{ Nothing, String } # spec name: origin
    severity::Any # spec type: Union{ Nothing, String } # spec name: severity
    type::Any # spec type: Union{ Nothing, String } # spec name: type
    context::Any # spec type: Union{ Nothing, Vector{IssueContext} } # spec name: context

    function ConnectivityIssue(;origin=nothing, severity=nothing, type=nothing, context=nothing)
        o = new()
        validate_property(ConnectivityIssue, Symbol("origin"), origin)
        setfield!(o, Symbol("origin"), origin)
        validate_property(ConnectivityIssue, Symbol("severity"), severity)
        setfield!(o, Symbol("severity"), severity)
        validate_property(ConnectivityIssue, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        validate_property(ConnectivityIssue, Symbol("context"), context)
        setfield!(o, Symbol("context"), context)
        o
    end
end # type ConnectivityIssue

const _property_map_ConnectivityIssue = Dict{Symbol,Symbol}(Symbol("origin")=>Symbol("origin"), Symbol("severity")=>Symbol("severity"), Symbol("type")=>Symbol("type"), Symbol("context")=>Symbol("context"))
const _property_types_ConnectivityIssue = Dict{Symbol,String}(Symbol("origin")=>"String", Symbol("severity")=>"String", Symbol("type")=>"String", Symbol("context")=>"Vector{IssueContext}")
Base.propertynames(::Type{ ConnectivityIssue }) = collect(keys(_property_map_ConnectivityIssue))
Swagger.property_type(::Type{ ConnectivityIssue }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ConnectivityIssue[name]))}
Swagger.field_name(::Type{ ConnectivityIssue }, property_name::Symbol) =  _property_map_ConnectivityIssue[property_name]

const _allowed_ConnectivityIssue_origin = ["Local", "Inbound", "Outbound"]

const _allowed_ConnectivityIssue_severity = ["Error", "Warning"]

const _allowed_ConnectivityIssue_type = ["Unknown", "AgentStopped", "GuestFirewall", "DnsResolution", "SocketBind", "NetworkSecurityRule", "UserDefinedRoute", "PortThrottled", "Platform"]

function check_required(o::ConnectivityIssue)
    true
end

function validate_property(::Type{ ConnectivityIssue }, name::Symbol, val)
    if name === Symbol("origin")
        Swagger.validate_param(name, "ConnectivityIssue", :enum, val, _allowed_ConnectivityIssue_origin)
    end
    if name === Symbol("severity")
        Swagger.validate_param(name, "ConnectivityIssue", :enum, val, _allowed_ConnectivityIssue_severity)
    end
    if name === Symbol("type")
        Swagger.validate_param(name, "ConnectivityIssue", :enum, val, _allowed_ConnectivityIssue_type)
    end
end


text:
############################################
"""
    printtopoESRI(topo::VisClaw.Topo, filename::AbstractString; nodatavalue::Integer=-9999, center::Bool=true)

print topo data as a file (ESRI ASCII format)
"""
function printtopoESRI(topo::VisClaw.Topo, filename::AbstractString="topo.asc"; nodatavalue::Integer=-9999, center::Bool=true)
    nrows, ncols = size(topo.elevation)
    ## print
    open(filename, "w") do fileIO
        ## header
        @printf(fileIO, "%d    ncols\n", ncols)
        @printf(fileIO, "%d    nrows\n", nrows)
        if center
            @printf(fileIO, "%e    xllcenter\n", topo.x[1])
            @printf(fileIO, "%e    yllcenter\n", topo.y[1])
        else
            @printf(fileIO, "%e    xllcorner\n", topo.x[1])
            @printf(fileIO, "%e    yllcorner\n", topo.y[1])
        end
        @printf(fileIO, "%e    cellsize\n", topo.dx)
        @printf(fileIO, "%d    nodatavalue\n", nodatavalue)
        ## elevation
        if isa(topo.elevation, BitArray) || isa(topo.elevation, Array{Int,2}) || isa(topo.elevation, Array{UInt,2})
            [(if j != ncols
                @printf(fileIO, "%d ", topo.elevation[i,j])
            else
                @printf(fileIO, "%d\n", topo.elevation[i,j])
            end) for j=1:ncols, i=nrows:-1:1]
        else
            [(if j != ncols
                @printf(fileIO, "%14.6e ", topo.elevation[i,j])
            else
                @printf(fileIO, "%14.6e\n", topo.elevation[i,j])
            end) for j=1:ncols, i=nrows:-1:1]
        end
    end ## close
    return nothing
end
############################################
const printtopo = printtopoESRI

############################################
"""
    printdtopo(dtopo::VisClaw.DTopo, filename::AbstractString; center::Bool=true)

print dtopo data as a topotype-3-file
"""
function printdtopo(dtopo::VisClaw.DTopo, filename::AbstractString="dtopo.asc"; center::Bool=true)
    ## print
    open(filename, "w") do fileIO
        ## header
        @printf(fileIO, "%d    mx\n", dtopo.mx)
        @printf(fileIO, "%d    my\n", dtopo.my)
        @printf(fileIO, "%d    mt\n", dtopo.mt)
        if center
            @printf(fileIO, "%e    xllcenter\n", dtopo.x[1])
            @printf(fileIO, "%e    yllcenter\n", dtopo.y[1])
        else
            @printf(fileIO, "%e    xllcorner\n", dtopo.x[1])
            @printf(fileIO, "%e    yllcorner\n", dtopo.y[1])
        end
        @printf(fileIO, "%e    t0\n", dtopo.t0)
        @printf(fileIO, "%e    dx\n", dtopo.dx)
        @printf(fileIO, "%e    dy\n", dtopo.dy)
        @printf(fileIO, "%e    dt\n", dtopo.dt)
        ## dtopo
        if dtopo.mt == 1
            [(if j != dtopo.mx
                  @printf(fileIO, "%14.6e ", dtopo.deform[i,j])
              else
                  @printf(fileIO, "%14.6e\n", dtopo.deform[i,j])
              end) for j=1:dtopo.mx, i=dtopo.my:-1:1]
        else
            [(if j != dtopo.mx
                  @printf(fileIO, "%14.6e ", dtopo.deform[i,j,k])
              else
                  @printf(fileIO, "%14.6e\n", dtopo.deform[i,j,k])
              end) for j=1:dtopo.mx, i=dtopo.my:-1:1, k=1:dtopo.mt]
        end
    end ## close
    return nothing
end
############################################


text:
const known_directions = [:x,   :y,   :z,
                          :yz,  :xz,  :xy,
                          :zy,  :zx,  :yx,
                          :xyz, :yxz, :xzy, :zyx]

"""
    relerr(x, truex)

Calculate relative error between calculated value `x` and theoretical
value `truex`.
"""
relerr(x, truex) = abs((x - truex) / truex)

"""
    scan(x)

Scan an array `x`, that is calculate running sum of its elements.
"""
function scan(x :: AbstractVector)
    result = copy(x)

    for i in 2:length(result)
        result[i] += result[i-1]
    end

    return result
end

"""
    integrate(f, range)

Integrate a function of one argument `f` across the range `range`
using Simpson's 1/3 rule
"""
function integrate(f :: Function, r :: AbstractRange)
    iter = zip(r, drop(r, 1))
    sum = 0.0

    for (a, b) in iter
        middle = (a + b) / 2
        sum += step(r)/6 * (f(a) + 4f(middle) + f(b))
    end

    return sum
end

# Value noise for tests

function lolrng(x    :: Integer,
                y    :: Integer,
                z    :: Integer,
                seed :: Integer)
    r1 = UInt32(x) * UInt32(0x1B873593)
    r2 = UInt32(y) * UInt32(0x19088711)
    r3 = UInt32(z) * UInt32(0xB2D05E13)

    r  = UInt32(seed) + r1 + r2 + r3
    r ⊻= r >> UInt32(5)
    r *= UInt32(0xCC9E2D51)
    return r / 0xffffffff
end

function octave(x    :: Number,
                y    :: Number,
                z    :: Number,
                oct  :: Integer,
                seed :: Integer)
    divisor = 2.0^(-oct)

    xi = (x / divisor) |> floor |> Int
    yi = (y / divisor) |> floor |> Int
    zi = (z / divisor) |> floor |> Int

    δx = rem(x, divisor) / divisor
    δy = rem(y, divisor) / divisor
    δz = rem(z, divisor) / divisor

    v000 = lolrng(xi,     yi,     zi, seed)
    v001 = lolrng(xi + 1, yi,     zi, seed)
    v010 = lolrng(xi,     yi + 1, zi, seed)
    v011 = lolrng(xi + 1, yi + 1, zi, seed)

    v100 = lolrng(xi,     yi,     zi + 1, seed)
    v101 = lolrng(xi + 1, yi,     zi + 1, seed)
    v110 = lolrng(xi,     yi + 1, zi + 1, seed)
    v111 = lolrng(xi + 1, yi + 1, zi + 1, seed)

    inter(v1, v2, x) = v1 + (v2 - v1)*x
    v00 = inter(v000, v001, δx)
    v01 = inter(v010, v011, δx)
    v10 = inter(v100, v101, δx)
    v11 = inter(v110, v111, δx)

    v0 = inter(v00, v01, δy)
    v1 = inter(v10, v11, δy)

    v = inter(v0, v1, δz)
    return v
end

value_noise(x        :: Number,
            y        :: Number,
            z        :: Number,
            octaves  :: Integer,
            seed     :: Integer) =
                mapreduce((x, octave) -> x ./ 2.0^octave, +,
                          (octave(x, y, z, o, seed) for o in 0:octaves-1),
                          countfrom(0)) ./ 2*(1 - 2.0^(-octaves))

two_phase_noise_3d() =
    let noise = [value_noise(x/10, y/10, z/10, 6, rand(UInt32))
                 for x in 1:50, y in 1:50, z in 1:50]
        noise .< 0.5
    end

two_phase_noise_2d() =
    let noise = [value_noise(x/10, y/10, 0, 6, rand(UInt32))
                 for x in 1:50, y in 1:50]
        noise .< 0.5
    end


text:
abstract type UFDivisionPotential <: EoSModel end

is_splittable(::UFDivisionPotential) = false

struct WCADivision <: UFDivisionPotential end
struct BHDivision <: UFDivisionPotential end

abstract type UFTheoryModel{T<:UFDivisionPotential} <: EoSModel end

function ϵₓ(model::UFTheoryModel,V,T,z)

end

function σₓ(model::UFTheoryModel,V,T,z)

end
function I_f(model::UFTheoryModel,V,T,z,ϵₓ,σₓ,C)

end

function I_i(model::UFTheoryModel,V,T,z,ϵₓ,σₓ,C)

end
const UFTheoryConsts = (
    Iu_LJ_WCA =[
-0.37447216 -0.00081593 -0.53140336 0.06813345
 -0.01631296 0.01735490 -0.48321813 -0.10811097
 -0.11968918 0.06900995 -0.09297811 -0.05479239
 0.42113581 -0.11298284 0.06698083 0.13611539
 -0.13660913 0.02263133 0.00933081 0.21970435],

    Iu_LJ_BH =[
-0.88956910 0.00410326 -0.00176185 0.12604517
-0.45301169 0.09609255 -0.00939335 0.13116768
-0.02908313 0.14360795 -0.04132516 0.14642522
0.12098751 0.02496398 -0.00397194 0.12907731
0.16797647 -0.15518686 0.03169724 0.14408187],

Iu_Mie_WCA =[
-1.78056093 -3.16686373 7.68360475 -93.60419002
-0.01668085 0.36495829 -2.35505504 -14.95211664
1.43501852 -3.17823271 -2.65975346 -1.68622036
0.00771078 -0.71067776 3.37028434 1.73371125
-0.75639404 -4.01519074 6.52257588 -34.61784393
0.04103887 -0.56727066 2.18602934 14.15314345
0.48329816 0.19213347 7.61859331 -9.92811362
-0.05943504 1.70951783 5.12258334 -49.19660269
0.29800680 -0.04245160 11.33527039 -5.89983990
-0.03863233 0.63090719 -0.33588296 -24.96073802
-0.40370239 -0.63858832 9.01192466 -41.33876761
-0.12403302 0.18998693 -14.73978776 151.80982268
0.07958646 1.79956210 1.22984038 -32.09818586
-0.03238907 -0.02204735 -0.27406920 -4.38812980
-0.00158420 -0.92021684 15.17385932 2.83803911
0.77885655 -22.86944530 110.24039194 205.53391110
0.10810701 1.32984049 -6.75001640 -76.40272812
0.15279280 -1.40020406 1.84877442 -2.70192234
-0.26259091 0.76348918 0.25624963 24.26218002
0.25097389 -0.09398157 -2.42094525 9.05317197
],

If_Mie_WCA = [
    -0.00676215 0.35802345 -4.68894827 21.71301236
0.00576239 -0.19869030 3.52670932 -21.52329250
-0.01585189 -0.02353289 -0.38046094 -2.22940520
0.01688067 -0.14880487 1.27641823 2.81305965
0.46266889 4.14762330 16.29605943 40.66463742
0.21097529 0.29711213 -10.88510499 108.50376227
0.16099600 0.07601413 -14.07795812 68.10396204
-0.17351007 0.40133563 -3.10110242 -13.06800562
-0.00259994 1.38844039 -5.72096750 29.61608733
-0.06312554 6.15404500 -2.00412098 -55.95330391
0.28768461 -2.07560237 -34.73499062 -16.16371594
-0.04182885 -1.73566128 -13.79223815 304.28754909
0.12463286 -1.06919900 0.00000075 51.67574708
-0.18247949 1.07222692 -4.09891917 -0.00001846
0.09499683 0.38741975 12.47877072 100.23721592
0.00802157 0.99105527 -14.01113384 -124.75984338
-0.04486163 1.30492578 -3.81392000 -6.69041253
0.08329549 0.61403648 -4.27136752 -10.04750594
],
If_LJ_WCA = [
0.00326607 -0.00631528 -0.00134063 0.00803155 -0.00412890
0.98177233 0.37428833 -0.02613570 -0.23353006 0.11235306
0.09812966 0.44943620 0.18812166 -0.04076140 -0.19168200
0.10917110 -0.10567363 -0.01328055 0.04651761 0.01508000
0.00246771 0.06797539 -0.00376730 0.08330469 -0.06683759
],

If_LJ_BH = [
0.00048683 -0.00971268 0.00731929 -0.00108245 0.00299320
0.87965112 0.29643192 -0.15355144 -0.18000878 0.04834741
0.22925638 0.24842618 0.10330378 -0.02903055 -0.09036307
0.01554676 -0.01605316 0.02748238 -0.00863994 -0.09759044
0.02629016 0.02591591 -0.01051095 0.01659820 0.0116228
],
)

#mie fluid,BH division. LJ fluid is just a particular case
function dHS_BH(T̄ ,m)
    a = (72564575e-01,1.44851940e-02,-4.47700974e-04)
    b = (2.41366577e-01,9.94732182e-03,-2.83801421e-04,-2.48237845e-06)
    c = (1.44734716e-03,-4.97470110e-05,2.32705654e-07)
    A = evalpoly(a,m)
    B = evalpoly(b,m)
    C = evalpoly(c,m)
    return (one(T̄) + A*T̄)/evalpoly((one(B),B,C),T̄)
end

#mie fluid, WCA (BFC),
function dHS_WCA(T̄,m)
    m == 12 && return 2/(1+sqrt(T̄))^(1/6)
    a = (0.0,7.08059224e-01, 5.52970256e-01, 5.33848228e-01)
    b = (0.0,2.37463060e-01, 1.91956811e-01, -3.99938968e-02)
    c = (0.0,5.75476615e-01, 8.20345330e-01, 5.03543986e-01)
    d = (0.0,4.90928991e-01, 1.49870438e-01, -3.71426331e-02)
    e = (0.0,8.15769030e-03, -3.43788879e-03, 4.81941979e-04)
    minv = 1/m
    A = evalpoly(a,minv)
    B = evalpoly(b,minv)
    C = evalpoly(c,minv)
    D = evalpoly(d,minv)
    E = evalpoly(e,minv)
    y = T̄-1
    num = (one(m),A,B)
    den = (one(m),C,D,E)
    return evalpoly(num,y)/evalpoly(den,y)
end


text:
module mtest2
using StaticArrays
using MeshCore: VecAttrib
import MeshSteward: vselect, initbox, updatebox!
using Test
function test()
    xyz = [0.0 0.0; 633.3333333333334 0.0; 1266.6666666666667 0.0; 1900.0 0.0; 0.0 400.0; 633.3333333333334 400.0; 1266.6666666666667 400.0; 1900.0 400.0; 0.0 800.0; 633.3333333333334 800.0; 1266.6666666666667 800.0; 1900.0 800.0]
    N, T = size(xyz, 2), eltype(xyz)
    locs =  VecAttrib([SVector{N, T}(xyz[i, :]) for i in 1:size(xyz, 1)])
    @test length(locs) == 12
    x = locs[SVector{2}([2, 4])]
    @test x[1] == SVector{2}([633.3333333333334 0.0])
    @test length(locs[1]) == 2
    @test eltype(locs[1]) == Float64

    box = [0.0 0.0 0.0 0.0]
    inflate = 0.01
    outputlist = vselect(locs; box = box, inflate = inflate)
    @test isapprox(outputlist, [1])

    box = initbox(xyz)
    outputlist = vselect(locs; box = box, inflate = inflate)
    @test isapprox(outputlist, collect(1:size(xyz, 1)))

    box = initbox(xyz[2, :])
    updatebox!(box, xyz[6, :])
    outputlist = vselect(locs; box = box, inflate = inflate)
    @test isapprox(outputlist, [2, 6])

    box = initbox(xyz[5, :])
    updatebox!(box, xyz[6, :])
    outputlist = vselect(locs; box = box, inflate = inflate)
    @test isapprox(outputlist, [5, 6])
    true
end
end
using .mtest2
mtest2.test()

include("samplet4.jl")

module mt4topo1
using StaticArrays
using MeshCore: P1, T4, ShapeColl,  manifdim, nvertices, nshapes, indextype
using MeshCore: VecAttrib
using MeshCore: IncRel, nrelations, nentities, ir_boundary
using MeshSteward: connectedv
using ..samplet4: samplet4mesh
using Test
function test()
    xyz, cc = samplet4mesh()
    # Construct the initial incidence relation
    N, T = size(xyz, 2), eltype(xyz)
    locs =  VecAttrib([SVector{N, T}(xyz[i, :]) for i in 1:size(xyz, 1)])
    vrts = ShapeColl(P1, length(locs))
    tets = ShapeColl(T4, size(cc, 1))
    ir30 = IncRel(tets, vrts, cc)
    vl = connectedv(ir30)
    @test length(vl) == size(xyz, 1)
    ir20 = ir_boundary(ir30)
    vl = connectedv(ir20)
    @test length(vl) == 54
    true
end
end
using .mt4topo1
mt4topo1.test()

module mt4topv1
using StaticArrays
using LinearAlgebra
using MeshCore: P1, T4, ShapeColl,  manifdim, nvertices, nshapes, indextype
using MeshCore: VecAttrib
using MeshCore: IncRel, nrelations, nentities, ir_boundary
using MeshSteward: connectedv, vselect, vtkwrite
using ..samplet4: samplet4mesh
using Test
function test()
    xyz, cc = samplet4mesh()
    # Construct the initial incidence relation
    N, T = size(xyz, 2), eltype(xyz)
    locs =  VecAttrib([SVector{N, T}(xyz[i, :]) for i in 1:size(xyz, 1)])
    vrts = ShapeColl(P1, length(locs))
    tets = ShapeColl(T4, size(cc, 1))
    ir30 = IncRel(tets, vrts, cc)
    ir30.right.attributes["geom"] = locs
    outputlist = vselect(locs; distance = 0.3, from = locs[1], inflate = 0.001)
    nmatched = 0
    for i in 1:length(locs)
        if norm(locs[i]) < 0.3 + 0.001
            nmatched = nmatched + 1
        end
    end
    @test nmatched == length(outputlist)
    vtkwrite("delete_me", ir30)
    try rm("delete_me" * ".vtu"); catch end
    true
end
end
using .mt4topv1
mt4topv1.test()

module mt4topv2
using StaticArrays
using LinearAlgebra
using MeshCore: P1, T4, ShapeColl,  manifdim, nvertices, nshapes, indextype
using MeshCore: VecAttrib
using MeshCore: IncRel, nrelations, nentities, ir_boundary
using MeshSteward: connectedv, vselect, vtkwrite
using ..samplet4: samplet4mesh
using Test
function test()
    xyz, cc = samplet4mesh()
    # Construct the initial incidence relation
    N, T = size(xyz, 2), eltype(xyz)
    locs =  VecAttrib([SVector{N, T}(xyz[i, :]) for i in 1:size(xyz, 1)])
    vrts = ShapeColl(P1, length(locs))
    tets = ShapeColl(T4, size(cc, 1))
    ir30 = IncRel(tets, vrts, cc)
    ir30.right.attributes["geom"] = locs
    outputlist = vselect(locs; plane = [1.0, 0.0, 0.0, 5.0], thickness = 0.001)
    nmatched = 0
    for i in 1:length(locs)
        if abs(locs[i][1] - 5.0) < 0.001
            nmatched = nmatched + 1
        end
    end
    @test nmatched == length(outputlist)
    # vtkwrite("delete_me", ir30)
    # try rm("delete_me" * ".vtu"); catch end
    true
end
end
using .mt4topv2
mt4topv2.test()

module mt4topv3
using StaticArrays
using LinearAlgebra
using MeshCore: P1, T4, ShapeColl,  manifdim, nvertices, nshapes, indextype
using MeshCore: VecAttrib
using MeshCore: IncRel, nrelations, nentities, ir_boundary
using MeshSteward: connectedv, vselect, vtkwrite
using ..samplet4: samplet4mesh
using Test
function test()
    xyz, cc = samplet4mesh()
    # Construct the initial incidence relation
    N, T = size(xyz, 2), eltype(xyz)
    locs =  VecAttrib([SVector{N, T}(xyz[i, :]) for i in 1:size(xyz, 1)])
    vrts = ShapeColl(P1, length(locs))
    tets = ShapeColl(T4, size(cc, 1))
    ir30 = IncRel(tets, vrts, cc)
    ir30.right.attributes["geom"] = locs
    outputlist = vselect(locs; nearestto = locs[13])
    @test length(outputlist) == 1
    @test outputlist[1] == 13
    true
end
end
using .mt4topv3
mt4topv3.test()


text:
# be sure they are installed Pkg.add
using Random, Distributions, StatsBase, Plots;

Random.seed!(123)

function cost_calculation_rds_postgres_per_user(v_invocations, v_gb)
    out_costs = 0.09
    storage_cost = 0.23
    out_costs_per_user = out_costs * v_gb * v_invocations + storage_cost * v_gb
end

function cost_calculation_rds_postgres(v_invocations, v_gb)
    user_costs = cost_calculation_rds_postgres_per_user.(v_invocations, v_gb)
    backup_costs = sum(v_gb) * 0.010 +  sum(v_gb) * 0.022 # s3 1 backup
    backup_costs + sum(user_costs)
end

function cost_calculation_s3(v_invocations, v_gb)
    storage_costs = sum(v_gb) * 0.023
    post_cost = length(v_gb) * 0.005 / 1000 * 31 * 1 # once per day store day transaction file
    select_costs = length(v_gb) * v_invocations * 0.0004 / 1000
    storage_costs + post_cost + select_costs
end

function gb_dist(users)
    data_distribution = SkewNormal(2.0, 3.0, 6.0)
    #plot(x->pdf(data_distribution,x+4))
    sample_values = rand(data_distribution, users)
    dt = fit(UnitRangeTransform, sample_values, dims = 1)
    # normalise 0:1
    StatsBase.transform(dt, sample_values)
end

max_users= 100000
v_gb = 0.0001 .+ gb_dist(max_users) .* 1.5;    # Distribution 100kb -> 1.5 Gb
v_invocations = 1 # download

u=1000
p,s = cost_calculation_rds_postgres(v_invocations, v_gb[1:u]), cost_calculation_s3(v_invocations, v_gb[1:u])

max_stars_repo_name: rutjes-dev/storage-calculations

max_issues_repo_name: rutjes-dev/storage-calculations

max_forks_repo_name: rutjes-dev/storage-calculations


text:
export readgoterms

# read in go ontology data
# 2 columns
# 1: node name, 2: ontology assocated to node
# ignore ontology pairs that are not in nodesidx
function readgoterms(fd::IO,nodes::Vector)
    gobases = ["GO:0008372","GO:0008150","GO:0003674"]
    nodesidx = indexmap(nodes)
    goknown = [Vector{String}() for i in 1:length(nodesidx)]
    goterms = Set{String}()
    for line = eachline(fd)
        nodename,goterm = split(strip(line))
        #goterm = goterm[4:end]
        (!haskey(nodesidx,nodename) || findfirst(gobases,goterm)>0) && continue
        node = nodesidx[nodename]
        push!(goknown[node],goterm)
        push!(goterms,goterm)
    end
    goterms = collect(goterms)
    goknown,goterms
end
readgoterms(file::AbstractString,args...;vargs...) =
    open(fd->readgoterms(fd,args...;vargs...),file,"r")


text:
"""
metaformattedclassresult(result::Array,Xtest::Array)

LightGBM実行形式で出力される行列フォーマットに変換

Converts to matrix format output in LightGBM executable format
# Arguments
`result::Array`:prediction
`Xtest::Array`:the features data.
"""
function metaformattedclassresult(result::Array,X::Array)
    rowsize=size(X,1)
    colsize=size(result,1)/rowsize
    colsize=Int(colsize)
    metaformattedresult=zeros(rowsize,colsize)

    for i in 0:rowsize-1
        for j in 1:colsize
            metaformattedresult[i+1,j]=result[i*colsize+j,1]
        end
    end

    return metaformattedresult
end

"""
metaformattedclassresult(metaformattedresult::Array)

分類予測結果から、予測精度の一番高い結果だけをLightGBM実行形式で出力される形式で出力

From the classification prediction result, only the result with the highest prediction accuracy is output in the output format in the LightGBM execution format
# Arguments
`metaformattedresult::Array`:formatted prediction result.
"""
function metaformattedclassresult(metaformattedresult::Matrix)
    rowsize,colsize=size(metaformattedresult)
    metaformattedclassresult=zeros(Int,rowsize,1)
    metaformattedclassresult=Array(metaformattedclassresult)

    for i in 1:rowsize
        work=-Inf
        for j in 1:colsize
            if(metaformattedresult[i,j] >= work)
                work =metaformattedresult[i,j]
                metaformattedclassresult[i,1]=j-1
            end
        end
    end

    return metaformattedclassresult
end

"""
formattedclassfit(result::Array,Xtest::Array)

予測精度の一番高い結果だけをLightGBM実行形式で出力される形式で出力

Only the result with the highest prediction accuracy is output in the output format in the LightGBM execution format
# Arguments
`result::Array`:prediction result.
`Xtest::Array`:the features data.
"""
function formattedclassfit(result::Array,X::Array)
    return metaformattedclassresult(metaformattedclassresult(result,X))
end

"""
predict2(estimator::LGBMEstimator, Xtest::Array)

分類予測の場合、予測精度の一番高い結果だけを出力する。その他の予測は結果をそのまま出力

 In the case of classification prediction, only the result with the highest prediction accuracy is output. Other predictions output the results as they are
# Arguments
* `estimator::LGBMEstimator`: the estimator to use in the prediction.
* `X::Matrix{T<:Real}`: the features data.
* `predict_type::Integer`: keyword argument that controls the prediction type. `0` for normal
    scores with transform (if needed), `1` for raw scores, `2` for leaf indices.
* `num_iterations::Integer`: keyword argument that sets the number of iterations of the model to
    use in the prediction. `< 0` for all iterations.
* `verbosity::Integer`: keyword argument that controls LightGBM's verbosity. `< 0` for fatal logs
    only, `0` includes warning logs, `1` includes info logs, and `> 1` includes debug logs.
* `is_row_major::Bool`: keyword argument that indicates whether or not `X` is row-major. `true`
    indicates that it is row-major, `false` indicates that it is column-major (Julia's default).
"""
function predict2(estimator::LGBMEstimator, X::Matrix{TX}; predict_type::Integer = 0,
                           num_iterations::Integer = -1, verbosity::Integer = 1,
                           is_row_major::Bool = false) where TX<:Real
    result=LightGBM.predict(estimator, X, predict_type = predict_type,
                           num_iterations = num_iterations, verbosity = verbosity,
                           is_row_major = is_row_major)

    if(typeof(estimator) == LGBMMulticlass )
        result=formattedclassfit(result,X)
    elseif(typeof(estimator) == LGBMRegression )
        result=result
    elseif(typeof(estimator) == LGBMBinary )
        result=result
    else
        println("Error")
    end

    return result
end


text:
"""
    Code for convenience functions for calculating CCFs for an AbstractChunkList (see RvSpectMLBase)
Author: Eric Ford
Created: August 2020
"""

"""  `calc_ccf_chunklist ( chunklist, ccf_plans )`
Convenience function to compute CCF based on a spectrum's chunklist.
# Inputs:
- chunklist
- vector of ccf plans (one for each chunk)
# Optional Arguments:
- `assume_sorted`:  if true, skips checking the line_list is sorted by wavelength
# Return:
CCF summed over all chunks in a spectrum's chunklist, evaluated using the
line list and mask_shape from the ccf plan for each chunk.
"""
function calc_ccf_chunklist(chunk_list::AbstractChunkList,
                                plan_for_chunk::AbstractVector{PlanT};
                                Δfwhm::Real = 0,
                                assume_sorted::Bool = false ) where {
                                            PlanT<:AbstractCCFPlan }
   num_vels = calc_length_ccf_v_grid(first(plan_for_chunk))
   ccf_out = Array{Float64,1}(undef, num_vels)
   if Δfwhm > 0
      this_plan_for_chunk = copy(plan_for_chunk)
      #increase_mask_fwhm!(this_plan_for_chunk,Δfwhm)
      map(chid->increase_mask_fwhm!(this_plan_for_chunk[chid],Δfwhm), 1:num_chunks )
   else
      this_plan_for_chunk = plan_for_chunk
   end
   calc_ccf_chunklist!(ccf_out, chunk_list, this_plan_for_chunk, assume_sorted=assume_sorted )
   return ccf_out
end

"""  `calc_ccf_chunklist! ( ccfs_out, chunklist, ccf_plans )`
Convenience function to compute CCF based on a spectrum's chunklist.
# Inputs:
- chunklist
- vector of ccf plans (one for each chunk)
# Optional Arguments:
- `assume_sorted`:  if true, skips checking the line_list is sorted by wavelength
# Return:
CCF summed over all chunks in a spectrum's chunklist, evaluated using the
line list and mask_shape from the ccf plan for each chunk.
"""
function calc_ccf_chunklist!(ccf_out::AbstractArray{T1,1}, chunk_list::AbstractChunkList,
                                plan_for_chunk::AbstractVector{PlanT};
                                assume_sorted::Bool = false ) where {
                                            T1<:Real, PlanT<:AbstractCCFPlan }
   num_chunks = length(chunk_list)
   @assert length(plan_for_chunk) == num_chunks
   num_vels = calc_length_ccf_v_grid(first(plan_for_chunk))
   @assert size(ccf_out,1) == num_vels
   @assert all(map(chid->calc_length_ccf_v_grid(plan_for_chunk[chid]), 1:num_chunks ) .== num_vels )
   @assert assume_sorted || all(map(chid->issorted( plan_for_chunk[chid].line_list.λ ), 1:num_chunks))

   num_lines = mapreduce(chid->length(plan_for_chunk[chid].line_list), +, 1:length(chunk_list.data) )
   ccf_out .= mapreduce(chid->calc_ccf_chunk(chunk_list.data[chid], plan_for_chunk[chid],
                    assume_sorted=true ), +, 1:length(chunk_list.data) )
   ccf_out .*= num_lines
end

"""  `calc_ccf_and_var_chunklist ( chunklist, ccf_plans )`
Convenience function to compute CCF based on a spectrum's chunklist.
# Inputs:
- chunklist
- vector of ccf plans (one for each chunk)
# Optional Arguments:
- `assume_sorted`:  if true, skips checking the line_list is sorted by wavelength
# Return:
CCF summed over all chunks in a spectrum's chunklist, evaluated using the
line list and mask_shape from the ccf plan for each chunk.
"""
function calc_ccf_and_var_chunklist(chunk_list::AbstractChunkList,
                                plan_for_chunk::AbstractVector{PlanT};
                                ccf_var_scale::Real = 1.0, Δfwhm::Real = 0,
                                assume_sorted::Bool = false ) where {
                                            PlanT<:AbstractCCFPlan }
  @assert length(chunk_list) == length(plan_for_chunk)
  num_chunks = length(chunk_list)
  num_vels = maximum(map(chid->calc_length_ccf_v_grid(plan_for_chunk[chid]), 1:num_chunks ))
  ccf_out = Array{Float64,1}(undef, num_vels)
  ccf_var_out = Array{Float64,1}(undef, num_vels)

  if Δfwhm > 0
     this_plan_for_chunk = copy(plan_for_chunk)
     map(chid->increase_mask_fwhm!(this_plan_for_chunk[chid],Δfwhm), 1:num_chunks )
  else
     this_plan_for_chunk = plan_for_chunk
  end

  #(ccf_out, ccf_var_out ) = mapreduce(chid->calc_ccf_and_var_chunk(chunk_list.data[chid], plan_for_chunk[chid],
  #                 assume_sorted=assume_sorted  ), add_tuple_sum, 1:length(chunk_list.data) )
  calc_ccf_and_var_chunklist!(ccf_out, ccf_var_out, chunk_list,this_plan_for_chunk, ccf_var_scale=ccf_var_scale, assume_sorted=assume_sorted)

  return (ccf=ccf_out, ccf_var=ccf_var_out)
end

"""  `calc_ccf_and_var_chunklist! ( ccf_out, ccf_var_out, chunklist, ccf_plans )`
Convenience function to compute CCF based on a spectrum's chunklist.
# Inputs:
- chunklist
- vector of ccf plans (one for each chunk)
# Optional Arguments:
- `assume_sorted`:  if true, skips checking the line_list is sorted by wavelength
# Return:
CCF summed over all chunks in a spectrum's chunklist, evaluated using the
line list and mask_shape from the ccf plan for each chunk.
"""
function calc_ccf_and_var_chunklist!(ccfs_out::AbstractArray{T1,1}, ccf_vars_out::AbstractArray{T2,1},
                                chunk_list::AbstractChunkList, plan_for_chunk::AbstractVector{PlanT};
                                ccf_var_scale::Real = 1.0, assume_sorted::Bool = false ) where {
                                            T1<:Real, T2<:Real, PlanT<:AbstractCCFPlan }
  num_chunks = length(chunk_list)
  @assert length(plan_for_chunk) == num_chunks
  @assert assume_sorted || all(map(chid->issorted( plan_for_chunk[chid].line_list.λ ), 1:num_chunks))
  num_vels = maximum(map(chid->calc_length_ccf_v_grid(plan_for_chunk[chid]), 1:num_chunks ))
  @assert size(ccfs_out,1) == num_vels
  @assert size(ccf_vars_out,1) == num_vels

  #(ccf_out, ccf_var_out ) = mapreduce(chid->calc_ccf_and_var_chunk(chunk_list.data[chid], plan_for_chunk[chid],
  #                  assume_sorted=assume_sorted  ), add_tuple_sum, 1:length(chunk_list.data) )

  ccfs_out .= 0
  ccf_vars_out .= 0
  ccfs_tmp = zeros(num_vels)
  ccf_vars_tmp = zeros(num_vels)
  for chid in 1:num_chunks
     calc_ccf_and_var_chunk!(ccfs_tmp, ccf_vars_tmp,
        chunk_list.data[chid], plan_for_chunk[chid], ccf_var_scale=ccf_var_scale, assume_sorted=true )
     ccfs_out .+= ccfs_tmp
     ccf_vars_out .+= ccf_vars_tmp
  end

  return (ccfs=ccfs_out, ccf_vars=ccf_vars_out)
end

"""  `calc_ccf_and_covar_chunklist! ( ccf_out, ccf_cocovar_out, chunklist, ccf_plans )`
Convenience function to compute CCF based on a spectrum's chunklist.
# Inputs:
- chunklist
- vector of ccf plans (one for each chunk)
# Optional Arguments:
- `assume_sorted`:  if true, skips checking the line_list is sorted by wavelength
# Return:
CCF summed over all chunks in a spectrum's chunklist, evaluated using the
line list and mask_shape from the ccf plan for each chunk.
"""
function calc_ccf_and_covar_chunklist!(ccfs_out::AbstractArray{T1,1}, ccf_covars_out::AbstractArray{T2,2},
                                chunk_list::AbstractChunkList, plan_for_chunk::AbstractVector{PlanT};
                                ccf_var_scale::Real = 1.0, assume_sorted::Bool = false ) where {
                                            T1<:Real, T2<:Real, PlanT<:AbstractCCFPlan }
  num_chunks = length(chunk_list)
  @assert length(plan_for_chunk) == num_chunks
  @assert assume_sorted || all(map(chid->issorted( plan_for_chunk[chid].line_list.λ ), 1:num_chunks))
  num_vels = maximum(map(chid->calc_length_ccf_v_grid(plan_for_chunk[chid]), 1:num_chunks ))
  @assert size(ccfs_out,1) == num_vels
  @assert size(ccf_covars_out,1) == size(ccf_covars_out,2) == num_vels

  #(ccf_out, ccf_covar_out ) = mapreduce(chid->calc_ccf_and_covar_chunk(chunk_list.data[chid], plan_for_chunk[chid],
  #                  assume_sorted=assume_sorted  ), add_tuple_sum, 1:length(chunk_list.data) )

  ccfs_out .= 0
  ccf_covars_out .= 0
  ccfs_tmp = zeros(num_vels)
  ccf_covars_tmp = zeros(num_vels,num_vels)
  for chid in 1:num_chunks
     calc_ccf_and_covar_chunk!(ccfs_tmp, ccf_covars_tmp,
        chunk_list.data[chid], plan_for_chunk[chid], ccf_var_scale=ccf_var_scale, assume_sorted=true )
     ccfs_out .+= ccfs_tmp
     ccf_covars_out .+= ccf_covars_tmp
  end

  return (ccfs=ccfs_out, ccf_covars=ccf_covars_out)
end

"""  `calc_ccf_and_covar_chunklist ( chunklist, ccf_plans )`
Convenience function to compute CCF based on a spectrum's chunklist.
# Inputs:
- chunklist
- vector of ccf plans (one for each chunk)
# Optional Arguments:
- `assume_sorted`:  if true, skips checking the line_list is sorted by wavelength
# Return:
CCF summed over all chunks in a spectrum's chunklist, evaluated using the
line list and mask_shape from the ccf plan for each chunk.
"""
function calc_ccf_and_covar_chunklist(chunk_list::AbstractChunkList,
                                plan_for_chunk::AbstractVector{PlanT};
                                ccf_var_scale::Real = 1.0,
                                Δfwhm::Real = 0,
                                assume_sorted::Bool = false ) where {
                                            PlanT<:AbstractCCFPlan }
  @assert length(chunk_list) == length(plan_for_chunk)
  num_chunks = length(chunk_list)
  num_vels = maximum(map(chid->calc_length_ccf_v_grid(plan_for_chunk[chid]), 1:num_chunks ))
  ccf_out = Array{Float64,1}(undef, num_vels)
  ccf_covar_out = Array{Float64,2}(undef, num_vels, num_vels)

  if Δfwhm > 0
     this_plan_for_chunk = copy(plan_for_chunk)
     #increase_mask_fwhm!(this_plan_for_chunk,Δfwhm)
     map(chid->increase_mask_fwhm!(this_plan_for_chunk[chid],Δfwhm), 1:num_chunks )
  else
     this_plan_for_chunk = plan_for_chunk
  end

  #(ccf_out, ccf_covar_out ) = mapreduce(chid->calc_ccf_and_covar_chunk(chunk_list.data[chid], plan_for_chunk[chid],
  #                 assume_sorted=assume_sorted  ), add_tuple_sum, 1:length(chunk_list.data) )
  calc_ccf_and_covar_chunklist!(ccf_out, ccf_covar_out, chunk_list,this_plan_for_chunk, ccf_var_scale=ccf_var_scale, assume_sorted=assume_sorted)

  return (ccf=ccf_out, ccf_covar=ccf_covar_out)
end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
function main(args)
    return args
end
# Load inline arguments for worker unit
num_workers, this = main(ARGS)
num_workers = parse(Int64, num_workers)
this        = parse(Int64, this)

#
# Important directories
simdir          = "/remote/ceph2/group/gerda/data/mpik/gerda-simulations/gerda-gems-db/pss-dev"
simdir_old      = "/remote/ceph/group/gerda/data/simulation/gerda-mage-sim"
mapping_file    = simdir_old*"/UTILS/det-data/ged-mapping.json"
parameters_file = simdir_old*"/UTILS/det-data/ged-parameters.json"

#
# Load packages and functions
include("../src/init.jl")
include("../src/fct.jl")
include("../src/worker_fct.jl")

#
# Output paths and filter settings
plots_base_path = "../../waveforms/sim/plots/raw/"
base_path_raw   = "../../waveforms/sim/raw/"
base_path       = "../../waveforms/sim/processed/"
n_sim_events    = 5000
isTP = 0;
isBL = 0;
mult = 1;
sampling_time   = 1u"ns"
hits_threshold  = 0.005; # MeV
E = 400;

#
# Load workers log file
log_file = base_path_raw * "log-" * string(this) * ".json"
if !isdir(base_path_raw)
    mkpath(base_path_raw)
end
if !isfile(log_file)
    open(log_file, "w") do f
        JSON.print(f, Dict(), 4)
    end
    global log = Dict()
else
    global log = JSON.parsefile(log_file)
end;

#
# Collect simulation filepaths
current_dir = pwd()
cd(simdir)
filenames = glob("calib/*/*/*/*.root")
cd(current_dir)

filenames = get_share_for_worker(filenames, num_workers, this)

for filename in filenames
    E_lim = E
    output_path = base_path_raw*basename(filename)*"/";
    if !haskey(log, basename(filename))
        Base.run(`clear`)
        @info(string(findfirst(x->x == filename, filenames)) * " of " * string(length(filenames)) * " in progress!")
        file = TFile(joinpath(simdir, filename))
        tree = file["fTree"];

        tt = Table(eventnumber = tree.eventnumber[:],
            hits_iddet   = tree.hits_iddet[:],
            hits_edep    = tree.hits_edep[:],
            hits_xpos    = tree.hits_xpos[:] .* 10, # mm
            hits_ypos    = tree.hits_ypos[:] .* 10, # mm
            hits_zpos    = tree.hits_zpos[:] .* 10) # mm

        tt = tt |> @filter(length(_.hits_edep) != 0 && sum(_.hits_edep) >= E_lim/1000 && length(unique(_.hits_iddet)) == 1) |> Table;

        data = Table(
            evtno        = [],
            multiplicity = [],
            detno        = [],
            hits_totnum  = [],
            edep         = [],
            pos          = []
        )
        prog = Progress(size(tt,1), dt=0.5,
            barglyphs=BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
            barlen=10)
        for i in eachindex(tt)
            multi = 1
            det = tt[i].hits_iddet[1]
            parameters[sim_to_channel[det][2]]["upside_down"] == true ? upside_down = -1 : upside_down = 1
            append!(data, Table(
            evtno        = [tt[i].eventnumber],
            multiplicity = [multi],
            detno        = [Array{Int64,1}(zeros(length(tt[i].hits_iddet)) .+ sim_to_channel[det][1])],
            hits_totnum  = [length(tt[i].hits_edep)],
            edep         = [tt[i].hits_edep],
            pos          = [[ SVector{3}(([ tt[i].hits_xpos[k] .- parameters[sim_to_channel[det][2]]["detcenter_x"],
                                        tt[i].hits_ypos[k] .- parameters[sim_to_channel[det][2]]["detcenter_y"],
                                        upside_down .* (tt[i].hits_zpos[k] .- parameters[sim_to_channel[det][2]]["detcenter_z"] .+ upside_down * parameters[sim_to_channel[det][2]]["height"]/2)
                            ] * u"mm")...) for k in eachindex(tt[i].hits_xpos) ]]
                )
            )
            next!(prog)
        end
        tt = nothing
        dets = []
        for detno in data.detno
            push!(dets, unique(detno)[1])
        end
        dets = unique(dets)
        prog = Progress(length(dets), dt=0.5,
            barglyphs=BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
            barlen=10)
        @info("Creating files for each detector")
        for detno in dets
            det = channel_to_name[detno]
            t = data |> @filter(unique(_.detno)[1] == detno) |> Table;
            !isdir(output_path) ? mkpath(output_path) : ""
            HDF5.h5open(output_path * lpad(detno, 2, "0") * "-" * det * ".h5", "w") do h5f
                LegendHDF5IO.writedata( h5f, "data", Table(
                    evtno       = t.evtno,
                    event_mult  = t.multiplicity,
                    detno       = VectorOfArrays(t.detno),
                    hits_totnum = t.hits_totnum,
                    edep        = VectorOfArrays(t.edep),
                    pos         = VectorOfArrays(t.pos)))
            end
            next!(prog)
        end
        #
        # Log the progress
        log[basename(filename)] = "Done!"
        open(log_file, "w") do f
            JSON.print(f, log, 4)
        end
    end
end


text:
__precompile__()

module UnitFloat

import Base:
    +,
    -,
    *,
    bits,
    convert,
    exponent,
    isfinite,
    isinf,
    isless,
    one,
    show,
    signbit,
    significand,
    typemax,
    typemin,
    zero

export UFloat

"""
    UFloat

The 64-bit type storing the *fraction* and the *exponent* of the unit-interval
floating point number.
```
         |-----------------------------------------|----------------------|
         | Exponent (41 bits)                      | Fraction (23 bits)   |
         |-----------------------------------------|----------------------|
bit no   |63                                     23|22                   0|
```
"""
primitive type UFloat 64 end # TODO can we extend AbstractFloat, like BigFloat?

"""
    _UFloat

Construct a `UFloat` given the exponent and the fraction. The exponent must be
strictly negative and >= `-2^41`. Only the last 23 bits of the given fraction
are used.
"""
@inline function _UFloat(exponent::Int64, fraction::Float32)
    bits = UInt64(reinterpret(UInt32, fraction)) & 0x00000000007fffff # ~UInt(0) >> 41
    bits |= ~reinterpret(UInt64, -exponent) << 23
    reinterpret(UFloat, bits)
end

@inline _is_zero(x)   = x == zero(typeof(x))
@inline _is_one(x)    = x == one(typeof(x))

@inline zero(::Type{UFloat}) = reinterpret(UFloat, 0x0000000000000000)
@inline one(::Type{UFloat})  = reinterpret(UFloat, 0xffffffffff800000)
zero(::UFloat) = zero(UFloat)
one(::UFloat) = one(UFloat)
typemin(::Type{UFloat}) = zero(UFloat)
typemax(::Type{UFloat}) = one(UFloat)
typemin(::UFloat) = zero(UFloat)
typemax(::UFloat) = one(UFloat)
isinf(::UFloat) = false
isfinite(::UFloat) = true

bits(u::UFloat) = bits(reinterpret(UInt64, u))
signbit(u::UFloat) = false

@inline function _exponent_unsafe(u::UFloat)
    # gives incorrect results for zero(UFloat) and one(UFloat)
    -reinterpret(Int64, ~reinterpret(UInt64, u) >> 23)
end
@inline function _exponent_unsafe(f::Float32)
    Int64((reinterpret(UInt32, f) & 0x7f800000) >> 23) - 127
end
@inline function _exponent_unsafe(f::Float64)
    Int64((reinterpret(UInt64, f) & 0x7ff0000000000000) >> 52) - 1023
end

function exponent(u::UFloat)
    _is_zero(u) && throw(DomainError())
    _exponent_unsafe(u)
end

@inline function _significand_unsafe(u::UFloat)
    # incorrect for zero(UFloat) and one(UFloat)
    bits = UInt32((reinterpret(UInt64, u) & 0x00000000007fffff)) # zero first 41 bits
    bits |= 0x3f800000 # add zero exponent to Float32
    reinterpret(Float32, bits)
end
@inline function _significand_unsafe(f::Float32)
    reinterpret(Float32, (reinterpret(UInt32, f) & 0x007fffff) | 0x3f800000)
end
@inline function _significand_unsafe(f::Float64)
    _significand_unsafe(Float32(f))
end

function significand(u::UFloat)
    _is_zero(u) && return 0.0f0
    _is_one(u) && return 1.0f0
    _significand_unsafe(u)
end

function convert(::Type{UFloat}, f::F) where {F<:AbstractFloat}
    f < zero(F) && throw(InexactError())
    f == zero(F) && return zero(UFloat)
    f == one(F) && return one(UFloat) # one special case because exponent==0
    exp, frac = exponent(f), significand(f)
    exp >= 0 && throw(InexactError())
    _UFloat(exp, Float32(frac))
end

function convert(::Type{F}, u::UFloat) where {F<:AbstractFloat}
    _is_zero(u) && return zero(F)
    _is_one(u) && return one(F)
    frac = _significand_unsafe(u)
    exp = _exponent_unsafe(u)
    ldexp(F(frac), exp)
end

convert(::Type{UFloat}, b::Bool) = b ? one(UFloat) : zero(UFloat)
convert(::Type{Bool}, u::UFloat) = u == one(UFloat)

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

function _multiply(a, b)::UFloat
    sa = _significand_unsafe(a)
    sb = _significand_unsafe(b)
    ea = _exponent_unsafe(a)
    eb = _exponent_unsafe(b)

    sn = sa * sb
    en = ea + eb + _exponent_unsafe(sn) # cannot be <= 1<<41

    en > -2199023255552 ? _UFloat(en, sn) : zero(UFloat)
end

@inline function _add(a, b)::UFloat
    _is_zero(a) && return UFloat(b)
    _is_zero(b) && return UFloat(a)

    sa = _significand_unsafe(a)
    sb = _significand_unsafe(b)
    ea = _exponent_unsafe(a)
    eb = _exponent_unsafe(b)

    ed = ea - eb

    sn = sb + ldexp(sa, ed)
    en = eb + _exponent_unsafe(sn)

    en >= 0 ? one(UFloat) : _UFloat(en, sn)
end

+(a::UFloat, b::UFloat) = _add(a, b)
*(a::UFloat, b::UFloat) = _multiply(a, b)
for T in [Float32, Float64], (op, fn) in [(:+, :_add), (:*, :_multiply)]
    @eval begin
        ($op)(a::UFloat, b::$T) = ($fn)(a, b)
        ($op)(a::$T, b::UFloat) = ($fn)(a, b)
    end
end

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

@inline function _isless(a, b)::Bool
    a == one(typeof(a)) && return false
    b == one(typeof(b)) && return true # assumed a ≠ 1.0
    b == zero(typeof(b)) && return false
    a == zero(typeof(a)) && return true # assumed b ≠ 0.0

    ea = _exponent_unsafe(a)
    eb = _exponent_unsafe(b)

    ea < eb && return true
    ea > eb && return false

    # ea == eb
    sa = _significand_unsafe(a)
    sb = _significand_unsafe(b)

    sa < sb
end

isless(a::UFloat, b::UFloat) = _isless(a, b)

begin
    local T = [UFloat, Float32, Float64]
    for T1 in T, T2 in T
        if T1 != T2
            @eval begin
                isless(a::$T1, b::$T2) = _isless(a, b)
            end
        end
    end
end

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

function show(io::IO, u::UFloat)
    if _is_zero(u)
        print(io, "0.0uf")
    elseif _is_one(u)
        print(io, "1.0uf")
    else
        # represent the number as a×10^b, with b int
        e2 = exponent(u)
        e10 = e2 * log10(2.0)
        e10i = round(e10)
        e10r = e10 - e10i
        s = Float32(significand(u) * exp10(e10r))
        while round(s, 5) >= 10.0f0
            s /= 10.0f0; e10i += 1
        end
        while round(s, 5) < 1.0f0
            s *= 10.0f0; e10i -= 1
        end
        @printf(io, "%.5fuf%d", s, e10i)
    end
end

end # module


text:
using SimplePartitions: find_part
import SimpleGraphs

export parents
parents(g::SimpleDigraph, v) = g.NN[v] |> collect

export children
children(g::SimpleDigraph, v) = g.N[v] |> collect

export partners
function partners(g::SimpleDigraph, v::Symbol)
    s = map(collect(children(g,v))) do x
        parents(g,x)
    end

    isempty(s) && return []

    setdiff(union(s...),[v]) |> collect
end

# function stochParents(m::Model, g::SimpleDigraph, v::Symbol, acc=Symbol[])
#     pars = parents(g,v)

#     result = union(pars, acc)
#     for p in pars
#         union!(result, _stochParents(m, g, findStatement(m,p)))
#     end
#     result
# end

# _stochParents(m::Model, g::SimpleDigraph, st::Sample, acc=Symbol[]) = union([st.x], acc)
# _stochParents(m::Model, g::SimpleDigraph, st::Assign, acc=Symbol[]) = stochParents(m, g, st.x, union([st.x],acc))
# _stochParents(m::Model, g::SimpleDigraph, st::Arg, acc=Symbol[]) = [st.x]
# _stochParents(m::Model, g::SimpleDigraph, st::Bool, acc=Symbol[]) = []

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

function stochChildren(m::Model, g::SimpleDigraph, v::Symbol, acc=Symbol[])
    pars = children(g,v)

    result = union(pars, acc)
    for p in pars
        union!(result, _stochChildren(m, g, findStatement(m,p)))
    end
    result
end

_stochChildren(m::Model, g::SimpleDigraph, st::Sample, acc=Symbol[]) = union([st.x], acc)
_stochChildren(m::Model, g::SimpleDigraph, st::Assign, acc=Symbol[]) = stochChildren(m, g, st.x, union([st.x],acc))
_stochChildren(m::Model, g::SimpleDigraph, st::Bool, acc=Symbol[]) = []

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

function stochPartners(m::Model, g::SimpleDigraph, v::Symbol)
    s = map(collect(stochChildren(m,g,v))) do x
        parents(g,x)
    end

    isempty(s) && return []

    setdiff(union(s...),[v]) |> collect
end

function markovBlanketVars(m::Model, g::SimpleDigraph, v::Symbol)
    markovBlanketVars(m,g,findStatement(m,v))
end

function markovBlanketVars(m::Model, g::SimpleDigraph,st::Sample)
    p = poset(m)
    g = digraph(m)

    body = Symbol[st.x]

    args = union(parents(g,st.x), stochPartners(m,g,st.x))
    for arg in args
        union!(body, interval(p, arg, st.x))
    end

    ys = stochChildren(m,g,st.x)
    union!(body, ys)
    for y in ys
        union!(body, interval(p, st.x, y))
    end

    setdiff!(body, args)
    (args, body)
end

function markovBlanketVars(m::Model, g::SimpleDigraph,st::Assign)
    (parents(g,st.x), Symbol[st.x])
end

# markovBlanket(g,v) = [v] ∪ stochParents(g,v) ∪ stochChildren(g,v) ∪ stochPartners(g,v)

# function Base.convert(::Type{SimpleGraph}, d::SimpleDigraph)
#     g = SimpleGraph{Symbol}()
#     add_edges!(g, elist(d))
#     g
# end

# # Type piracy, this should really go in SimpleGraphs
# # OTOH there's not anything else this could really mean
# SimpleGraphs.components(g::SimpleDigraph) = SimpleGraphs.components(convert(SimpleGraph, g))

export markovBlanket
function markovBlanket(m::Model, x :: Symbol)
    g = digraph(m) #Creates a new SimpleDigraph, so no need to copy before mutating

    # Find the strongly connected component containing x
    # part = find_part(SimpleGraphs.components(g), x) |> union

    # newargs = (arguments(m) ∪ [x]) ∩ part
    # setdiff!(part, newargs)

    (args, vars) = markovBlanketVars(m,g,x)

    M = getmodule(m)
    m_init = Model(M, args, NamedTuple(), NamedTuple(), nothing)
    m_init = merge(m_init, Model(M,findStatement(m,x)))
    m = foldl(vars; init= m_init) do m0,v
        merge(m0, Model(M,findStatement(m, v)))
    end
    m = merge(m, Model(M,findStatement(m,x)))

end


text:
using LinearAlgebra
using SparseArrays
using Plots
using ForwardDiff

"This code performs a convergence test for two finite difference methods applied to
Laplace's equation with mixed Dirichlet-Neumann boundary conditions."

# method of manufactured solutions: define the exact solution, then figure out
# what f,α,β should be to make it the solution to your PDE
uexact(x) = log(2+sin(pi*x))
dudx_exact(x) = ForwardDiff.derivative(uexact,x)
f(x) = -ForwardDiff.derivative(dudx_exact,x) # -d^2u/dx^2 = f

α = uexact(-1)
β = dudx_exact(1)

# The O(h) accurate method
function solve1(m,f,α,β)
    x = LinRange(-1,1,m+2) # add 2 endpoints
    xint = x[2:end]
    h = x[2]-x[1]
    A = (1/h^2) * spdiagm(0=>2*ones(m+1),-1=>-ones(m),1=>-ones(m))
    A[end,end] = 1/h^2
    b = f.(xint)
    b[1] += α/h^2

    b[end] += β/h # approach 1 - O(h) accuracy only

    u = A\b
    return u,x,xint,h
end

# The O(h^2) accurate method
function solve2(m,f,α,β)
    x = LinRange(-1,1,m+2) # add 2 endpoints
    xint = x[2:end] # leave x_{m+1} in the list of points!
    h = x[2]-x[1]
    A = (1/h^2) * spdiagm(0=>2*ones(m+1),-1=>-ones(m),1=>-ones(m))
    A[end,end] = 1/h^2
    b = f.(xint)
    b[1] += α/h^2

    b[end] = .5*f(x[end]) + β/h

    u = A\b
    return u,x,xint,h
end

# # this code plots the solution
# m = 40
# u1,x,xint1,h = solve1(m,f,α,β)
# u2,x,xint2,h = solve2(m,f,α,β)
# plot(xint1,u1,mark=:dot,markersize=2,label="O(h) solution")
# plot!(xint2,u2,mark=:dot,markersize=2,label="O(h^2) solution")
# plot!(x,uexact.(x),label="Exact solution")
# plot!(leg=:bottomright)

mvec = 2 .^ (2:9)
hvec = zeros(length(mvec))
err1 = zeros(length(mvec))
err2 = zeros(length(mvec))
for (i,m) in enumerate(mvec)
    u1,x,xint,h = solve1(m,f,α,β)
    u2,x,xint,h = solve2(m,f,α,β)
    hvec[i] = h
    err1[i] = maximum(@. abs(uexact(xint) - u1))
    err2[i] = maximum(@. abs(uexact(xint) - u2))
end
plot(hvec,err1,marker=:dot,label="Max error for method 1")
plot!(hvec,err2,marker=:dot,label="Max error for method 2")
plot!(hvec,hvec.^2,linestyle=:dash,label="O(h^2)")
plot!(xaxis=:log,yaxis=:log,legend=:bottomright)


text:
# Linear Systems STO
function simulate(m::linSTO)

  # Define Time
  t = collect(linspace(m.STOev.t0, m.STOev.tf, 10000))

  # Perform Simulation
  x, xpts, J = simulatelinsto(m.tau, m.STOev.x0, m.STOev.Q, m.STOev.E, m.STOev.A, t)

  return x, xpts, J, t

end

function simulate(m::linSTO, tau::Array{Float64,1})

  # Define Time
  t = collect(linspace(m.STOev.t0, m.STOev.tf, 10000))

  # Perform Simulation
  x, xpts, J = simulatelinsto(tau, m.STOev.x0, m.STOev.Q, m.STOev.E, m.STOev.A, t)

  return x, xpts, J, t

end

function simulate(m::linSTO, tau::Array{Float64,1}, t::Array{Float64,1})

  # Perform Simulation
  x, xpts, J = simulatelinsto(tau, m.STOev.x0, m.STOev.Q, m.STOev.E, m.STOev.A, t)

  return x, xpts, J, t

end

# Nonlinear Systems STO
function simulate(m::nlinSTO)

  # Define Time
  t = collect(linspace(m.STOev.t0, m.STOev.tf, 10000))

  # Get original Q, x0
  Q = m.STOev.Q[1:end-1, 1:end-1]
  E = m.STOev.E[1:end-1, 1:end-1]
  x0 = m.STOev.x0[1:end-1]

  # Perform Simulation
  x, xpts, J = simulatenlinsto(m.STOev.nonlin_dyn, m.tau, x0,  Q, E,  m.STOev.uvec, t)

  return x, xpts, J, t

end

function simulate(m::nlinSTO, tau::Array{Float64, 1})

  # Define Time
  t = collect(linspace(m.STOev.t0, m.STOev.tf, 10000))

  # Get original, Q, x0
  Q = m.STOev.Q[1:end-1, 1:end-1]
  E = m.STOev.E[1:end-1, 1:end-1]
  x0 = m.STOev.x0[1:end-1]

  # Perform Simulation
  x, xpts, J = simulatenlinsto(m.STOev.nonlin_dyn, tau, x0,  Q, E,  m.STOev.uvec, t)

  return x, xpts, J, t

end

function simulate(m::nlinSTO, tau::Array{Float64, 1}, t::Array{Float64, 1})

  # Get original, Q, x0
  Q = m.STOev.Q[1:end-1, 1:end-1]
  E = m.STOev.E[1:end-1, 1:end-1]
  x0 = m.STOev.x0[1:end-1]

  # Perform Simulation
  x, xpts, J = simulatenlinsto(m.STOev.nonlin_dyn, tau, x0,  Q, E,  m.STOev.uvec, t)

  return x, xpts, J, t

end

# Linearized Nonlinear System STO
function simulatelinearized(m::nlinSTO)

  # Define Time
  t = collect(linspace(m.STOev.t0, m.STOev.tf, 10000))

  # Get original Q, x0
  Q = m.STOev.Q[1:end-1, 1:end-1]
  E = m.STOev.E[1:end-1, 1:end-1]
  x0 = m.STOev.x0[1:end-1]

  # Perform Simulation
  x, xpts, J = simulatelinearizedsto(m.STOev.nonlin_dyn, m.STOev.nonlin_dyn_deriv, m.tau, m.STOev.tgrid, m.STOev.tfdelta, m.STOev.uvec,  x0, Q, E,  t)

  return x, xpts, J, t

end

# Linearized Nonlinear System STO
function simulatelinearized(m::nlinSTO, tau::Array{Float64,1})

  # Define Time
  t = collect(linspace(m.STOev.t0, m.STOev.tf, 10000))

  # Get original Q, x0
  Q = m.STOev.Q[1:end-1, 1:end-1]
  E = m.STOev.E[1:end-1, 1:end-1]
  x0 = m.STOev.x0[1:end-1]

  # Perform Simulation
  x, xpts, J = simulatelinearizedsto(m.STOev.nonlin_dyn, m.STOev.nonlin_dyn_deriv, tau, m.STOev.tgrid, m.STOev.uvec,  x0, Q, E,  t)

  return x, xpts, J, t

end

# Linearized Nonlinear System STO
function simulatelinearized(m::nlinSTO, tau::Array{Float64,1}, t::Array{Float64, 1})

  # Get original Q, x0
  Q = m.STOev.Q[1:end-1, 1:end-1]
  E = m.STOev.E[1:end-1, 1:end-1]
  x0 = m.STOev.x0[1:end-1]

  # Perform Simulation
  x, xpts, J = simulatelinearizedsto(m.STOev.nonlin_dyn, m.STOev.nonlin_dyn_deriv, tau, m.STOev.tgrid, m.STOev.uvec, x0, Q, E, t)

  return x, xpts, J, t

end

### Lower Level Functions

# Linear Systems
function simulatelinsto(tau::Array{Float64,1}, x0::Array{Float64, 1}, Q::Array{Float64,2}, E::Array{Float64,2}, A::Array{Float64,3}, t::Array{Float64,1})
  # Get dimensions
  nx = size(A, 1)  # Number of States
  N = length(tau)  # Number of switches

  tau = [t[1]; tau; t[end]]  # Extend tau vector to simplify numbering

  # Compute Initial States
  xpts = Array{Float64}(nx, N+1)
  xpts[:,1] = x0
  for i = 2:N+1
    xpts[:,i] = expm(A[:, :, i-1]*(tau[i] - tau[i-1]))*xpts[:, i-1]
  end

  # Compute State Trajectory
  x = Array{Float64}(nx, length(t))
  x[:, 1] = x0
  tauind = 1  # Index to keep track of the current mode

  for i = 2:length(t)
    # Check if we are still in the current switching interval. Otherwise Change
    if tauind < N+1
      if t[i] > tau[tauind + 1]
        tauind += 1
      end
    end

    # Compute State
    x[:, i] = expm(A[:, :, tauind]*(t[i] - tau[tauind]))*xpts[:, tauind]

  end

  # Numerically Integrate Cost Function
  Jtoint = diag(x'*Q*x)
  J = trapz(t, Jtoint) + (x[:, end]'*E*x[:, end])[1]

  return x, xpts, J

end

# Linearized Nonlinear System
  function simulatelinearizedsto(nonlin_dyn::Function, nonlin_dyn_deriv::Function, tau::Array{Float64,1}, tgrid::Array{Float64,1}, uvec::Array{Float64, 2}, x0::Array{Float64, 1}, Q::Array{Float64,2},  E::Array{Float64,2}, t::Array{Float64,1})

    # Get dimensions
    nx = length(x0)  # Number of States
    N = length(tau)  # Number of switches
    ngrid = length(tgrid)  # Number of elements in time grid

    # Create merged and sorted time vector with grid and switching times
    tvec, tauIdx = mergeSortFindIndex(tgrid, tau)

    # Create matrix of Linearized Dynamics
    A = Array{Float64}(nx+1, nx+1, N+ngrid - 1)

    # Compute Initial States
    xpts = Array{Float64}(nx+1, N+ngrid)  # Augmented State for Linearization
    xpts[:,1] = [x0; 1]

      uIdx = 1  # Initialize index for current u

      # Compute Matrix Exponentials
      for i =1:N+ngrid-1  # Iterate over all grid (incl sw times)

        # Verify which U input applies
        if uIdx <= N
          if i>= tauIdx[uIdx + 1]
            uIdx += 1
          end
        end

        # Linearize Dynamics
        A[:,:,i] = linearizeDyn(nonlin_dyn, nonlin_dyn_deriv, xpts[1:end-1,i], uvec[:,uIdx])

        # Compute Next Point in Simulation
        xpts[:,i+1] = expm(A[:, :, i]*(tvec[i+1] - tvec[i]))*xpts[:, i]

      end

      # Compute State Trajectory
      x = Array{Float64}(nx+1, length(t))
      x[:, 1] = [x0; 1]
      tauind = 1  # Index to keep track of the current mode

      for i = 2:length(t)

        # Check if we are still in the current switching interval. Otherwise Change
        if tauind < N+ngrid-1
          if t[i] > tvec[tauind + 1]
            tauind += 1
          end
        end

        # Compute State
        x[:, i] = expm(A[:, :, tauind]*(t[i] - tvec[tauind]))*xpts[:, tauind]

      end

      # Remove Extra State for Cost Function Evaluation
      x = x[1:end-1, :]

      # Numerically Integrate Cost Function
      Jtoint = diag(x'*Q*x)
      J = trapz(t, Jtoint) + (x[:, end]'*E*x[:, end])[1]

      return x, xpts, J

  end

# Noninear Systems
function simulatenlinsto(nonlin_dyn::Function, tau::Array{Float64,1}, x0::Array{Float64, 1}, Q::Array{Float64,2}, E::Array{Float64,2}, uvec::Array{Float64,2}, t::Array{Float64,1})
  # Get dimensions
  nx = length(x0)  # Number of States
  N = length(tau)  # Number of switches

  tau = [t[1]; tau; t[end]]  # Extend tau vector to simplify numbering

  # Define empty state trajectory
  x = zeros(nx, length(t))
  x[:,1] = x0

  # Define indeces to determine current switching mode
  tempInd1 = 1
  tempInd2 = 1
  xprevSwitch = x0

  # Create Vector of Points at the switching instants
  xpts = zeros(nx, N+1)
  xpts[:, 1] = x0

  for i = 1:N+1 # Integrate over all the intervals

    # redefine Dynamic function
    nldyn(x, p, t) = nonlin_dyn(x, uvec[:, i])

    while t[tempInd2] < tau[i+1]
      tempInd2 = tempInd2 + 1  # Increase time index
    end

    if tempInd2>tempInd1  # There has been a progress and the switching times are not collapsed. So we integrate.
      prob = DiffEqBase.ODEProblem(nldyn,xprevSwitch,(t[tempInd1],t[tempInd2]))
      if tempInd2 == tempInd1 + 1  # Only one step progress. Handle as special case for the integrator
        sol = DiffEqBase.solve(prob,OrdinaryDiffEq.Tsit5();save_everystep=false)
        xtemp = [sol[1] sol[end]]
        #_, xmap = ode45(nldyn, xprevSwitch, [t[tempInd1]; (t[tempInd1]+t[tempInd2])/2; t[tempInd2]], points=:specified)  # Integrate
        #xmap = hcat(xmap...)
        #xtemp = [xmap[:,1] xmap[:,3]]          # Take only two points
      else
        sol = DiffEqBase.solve(prob,OrdinaryDiffEq.Tsit5();saveat=t[tempInd1:tempInd2])
        xtemp = hcat((sol[i] for i in eachindex(sol))...)
      end

      x[:, tempInd1:tempInd2] = xtemp

      # Update indeces for next iteration
      xprevSwitch = x[:, tempInd2]
      tempInd1 = tempInd2

      # Update vector of switching states
      if i < N+1
        xpts[:, i+1] = x[:, tempInd2]
      end
    end

  end

  # Numerically Integrate Cost Function
  Jtoint = diag(x'*Q*x)
  J = trapz(t, Jtoint) + (x[:,end]'*E*x[:,end])[1]

  return x, xpts, J

end

#### Auxiliary functions
function simulateinput(m::nlinSTO)

  # Define Time
  t = collect(linspace(m.STOev.t0, m.STOev.tf, 10000))

  # u = computenlswinput(m.taucomplete, m.STOev.uvec, t)
  u = computenlswinput(m.tau, m.STOev.uvec, t)  # Fix all TAUCOMPLETE

  return u, t

end

function simulateinput(m::nlinSTO, t::Array{Float64, 1})

  # u = computenlswinput(m.taucomplete, m.STOev.uvec, t)
  u = computenlswinput(m.tau, m.STOev.uvec, t)  # Fix all TAUCOMPLETE

  return u, t

end

# Compute Actual Inputs from Artificial Ones
function computenlswinput(tauopt, uvec, t)
  N = length(tauopt)
  nu = size(uvec, 1)
  tau = [t[1]; tauopt; t[end]]  # Add initial and final switching to simplify numbering

  # Define empty input vector
  usim = zeros(nu, length(t))

  # Define indeces to span time vector
  tempInd1 = 1
  tempInd2 = 1

  for i = 1:N+1  # Iterate over all intervals
    while t[tempInd2] < tau[i+1]
      tempInd2 += 1
    end

    if tempInd2 > tempInd1
      usim[:, tempInd1:tempInd2] = repmat(uvec[:, i], 1, tempInd2 - tempInd1 + 1)
      tempInd1 = tempInd2
    end

  end
  return usim

end

# Trapezoidal integration rule
function trapz{Tx<:Number, Ty<:Number}(x::Vector{Tx}, y::Vector{Ty})
    local n = length(x)
    if (length(y) != n)
        error("Vectors 'x', 'y' must be of same length")
    end
    r = zero(zero(Tx) + zero(Ty))
    if n == 1; return r; end
    for i in 2:n
        r += (x[i] - x[i-1]) * (y[i] + y[i-1])
    end
    return r/2
end


text:
module sequenceCompare
    export DTWDistance

    # http://en.wikipedia.org/wiki/Dynamic_time_warping
    function DTWDistance{T}(s::Array{T}, t::Array{T})
        n, m = size(s,2), size(t, 2)
        DTW = zeros(n+1,m+1)
        #println("sizeof DTW: ", size(DTW))
        DTW[:,1], DTW[1,:] = Inf, Inf
        DTW[1,1] = 0

        for i=2:n+1
            for j=2:m+1
                #println("i=", i, "/", n, " j=", j, "/", m, " size(s)=", size(s), " size(t)=", size(t))
                cost = sum((s[:,i-1] - t[:,j-1]).^2)

                DTW[i,j] = cost + min(DTW[i-1, j  ], #insertion
                                      DTW[i  , j-1], #deletion
                                      DTW[i-1, j-1]) #match
                #println("cost is ", cost, " i=", i, "/", n, " j=", j, "/", m, " DTW=", DTW[i,j])
            end
        end

        return DTW[n+1,m+1]
    end

    # note. there is also a windowed version for better performance!  See the wikipedia article
    function DTWDistance{T}(s::Array{T}, t::Array{T}, w::Int64)
        n, m = size(s,2), size(t, 2)
        w = max(w, abs(n-m))
        DTW = zeros(n+1,m+1)
        #println("sizeof DTW: ", size(DTW))
        DTW[:,:] = Inf
        DTW[1,1] = 0

        for i=2:n+1
            for j=max(2, i-w):min(m+1, i+w)
                #println("i=", i, "/", n, " j=", j, "/", m, " size(s)=", size(s), " size(t)=", size(t))
                cost = sum((s[:,i-1] - t[:,j-1]).^2)

                DTW[i,j] = cost + min(DTW[i-1, j  ], #insertion
                                      DTW[i  , j-1], #deletion
                                      DTW[i-1, j-1]) #match
                #println("cost is ", cost, " i=", i, "/", n, " j=", j, "/", m, " DTW=", DTW[i,j])
            end
        end

        return DTW[n+1,m+1]
    end

    function HaversineDistance{T}(lat1::T, lon1::T, lat2::T, lon2::T)
        REarth = 6371.
        lat, lon = (lat1-lat2)*π/180., (lon1-lon2)*π/180.
        lat1, lat2 = lat1*π/180., lat2*π/180.
        a = sin(lat/2)*sin(lat/2)+cos(lat1)*cos(lat2)*sin(lon/2)*sin(lon/2)
        d = 2*atan2(sqrt(a), sqrt(1-a))
        d = REarth*d

        return d
    end

    function RMSE{T}(pred::T, real::T)
        return sqrt(mean((pred-real).*(pred-real)))
    end

    function RMSLE{T}(pred::T, real::T)
    return sqrt(mean((log(pred+1)-log(real+1)).^2))
    end

    function meanHaversideDistance{T}(lat1::T, lon1::T, lat2::T, lon2::T)
        return mean(HaversideDistance{T}(lat1, lon1, lat2, lon2))
    end

    function TravelTimePredict(submission, answers)
        return RMSLE(submission, answers)
    end

    function DestinationPredict(submission, answers)
        lat_sub, lon_sub = submission[:,1], submission[:,2]
        lat_real, lon_real = answers[:,1], answers[:,2]
        return meanHaversineDistance(lat_sub, lon_sub, lat_real, lon_real)
    end

end


text:
#=
Given two singly linked lists that intersect at some point, find the intersecting node. The lists are non-cyclical.

For example, given A = 3 -> 7 -> 8 -> 10 and B = 99 -> 1 -> 8 -> 10, return the node with value 8.

In this example, assume nodes with the same value are the exact same node objects.

Do this in O(M + N) time (where M and N are the lengths of the lists) and constant space.
=#

using Test
include("Solutions/problem20_find_intersecting_node.jl")

intersecting_node = Node(8, Node(10, nothing))
A = Node(3, Node(7, intersecting_node))
B = Node(99, Node(1, intersecting_node))

@test findintersect(A, B) == intersecting_node

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
module ParameterEstimateScatterPlots

using RecipesBase

@recipe function scatterestimatesrecipe(
    ::Type{Val{:scatterestimates}},
    truth,
    estimates,
    z;
    name = "",
    units = ""
)

    legend := :topleft

    # Scatter plot
    @series begin
        seriestype := :scatter
        markercolor := :yellow
        markeralpha := 0.7
        markersize := 5
        markerstrokewidth := 1
        label := ""
        xguide := "True $name" * (units == "" ? "" : " ($units)")
        yguide := "Estimated $name" * (units == "" ? "" : " ($units)")
        x := truth
        y := estimates
    end

    # Line of identity
    @series begin
        seriestype := :straightline
        linecolor := :black
        linewidth := 2
        label := "Ideal"
        x := [0, 1]
        y := [0, 1]
    end

    # Trend line
    (m, b) = [truth ones(eltype(truth), length(truth))] \ estimates
    @series begin
        seriestype := :straightline
        linecolor := :blue
        linewidth := 4
        label := "Trend"
        x := [0, 1]
        y := [b, m + b]
    end

end

export scatterestimates
scatterestimates(args...; kwargs...) = RecipesBase.plot(args...; kwargs..., seriestype = :scatterestimates)

end


text:
# Julia translation of http://nbviewer.jupyter.org/github/barbagroup/AeroPython/blob/master/lessons/03_Lesson03_doublet.ipynb
# Lesson 3 doublet

using PyPlot
using Distributions

close("all")
meshgrid(x,y) = (repmat(x',length(y),1),repmat(y,1,length(x)))

N = 50                                # number of points in each direction
x_start, x_end = -2.0, 2.0            # boundaries in the x-direction
y_start, y_end = -1.0, 1.0            # boundaries in the y-direction
x = linspace(x_start, x_end, N)    # creates a 1D-array with the x-coordinates
y = linspace(y_start, y_end, N)    # creates a 1D-array with the y-coordinates

X,Y=meshgrid(x,y)

kappa = 1.0

x_doublet, y_doublet = 0.0, 0.0   # location of the source

function get_velocity_doublet(strength, xd, yd, X, Y):
    """Returns the velocity field generated by a doublet.

    Arguments
    ---------
    strength -- strength of the doublet.
    xd, yd -- coordinates of the doublet.
    X, Y -- mesh grid.
    """
    u = - strength/(2*pi).*((X-xd).^2-(Y-yd).^2)./((X-xd).^2+(Y-yd).^2).^2
    v = - strength/(2*pi)*2.0.*(X-xd).*(Y-yd)./((X-xd).^2+(Y-yd).^2).^2

    return u, v
  end

function get_stream_function_doublet(strength, xd, yd, X, Y):
    """Returns the stream-function generated by a doublet.

    Arguments
    ---------
    strength -- strength of the doublet.
    xd, yd -- coordinates of the doublet.
    X, Y -- mesh grid.
    """
    psi = - strength/(2*pi).*(Y-yd)./((X-xd).^2+(Y-yd).^2)

    return psi
  end

  # computes the velocity field on the mesh grid
u_doublet, v_doublet = get_velocity_doublet(kappa, x_doublet, y_doublet, X, Y)

# computes the stream-function on the mesh grid
psi_doublet = get_stream_function_doublet(kappa, x_doublet, y_doublet, X, Y)

# Freestream and Doublet

  u_inf = 1.0

  # computes the freestream velocity field
  u_freestream = u_inf * ones((N, N))
  v_freestream = zeros((N, N))

  # computes the stream-function
  psi_freestream = u_inf * Y

  u = u_freestream + u_doublet
  v = v_freestream + v_doublet
  psi = psi_freestream + psi_doublet

  # plots the streamlines
  size = 10
  PyPlot.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))
  PyPlot.grid(true)
  PyPlot.xlabel("x", fontsize=16)
  PyPlot.ylabel("y", fontsize=16)
  PyPlot.xlim(x_start, x_end)
  PyPlot.ylim(y_start, y_end)
  PyPlot.streamplot(X, Y, u, v, density=2, linewidth=1, arrowsize=1, arrowstyle="->")
  PyPlot.contour(X, Y, psi, levels=[0.], colors="#CD2305", linewidths=2, linestyles="solid");
  PyPlot.scatter(x_doublet, y_doublet, color="#CD2305", s=80, marker="o")

  # calculates the stagnation points
x_stagn1, y_stagn1 = +sqrt(kappa/(2*pi*u_inf)), 0
x_stagn2, y_stagn2 = -sqrt(kappa/(2*pi*u_inf)), 0

# adds the stagnation points to the figure
PyPlot.scatter([x_stagn1, x_stagn2], [y_stagn1, y_stagn2], color="g", s=80, marker="o");

# computes the pressure coefficient field
cp = 1.0 - (u.^2+v.^2)./u_inf^2

# plots the pressure coefficient field
size = 10
PyPlot.figure(figsize=(1.1*size, (y_end-y_start)/(x_end-x_start)*size))
PyPlot.xlabel("x", fontsize=16)
PyPlot.ylabel("y", fontsize=16)
PyPlot.xlim(x_start, x_end)
PyPlot.ylim(y_start, y_end)
contf = PyPlot.contourf(X, Y, cp, levels=linspace(-2.0, 1.0, 100), extend="both")
cbar = PyPlot.colorbar(contf)
cbar[:set_label](L"$C_p$", fontsize=16)
cbar[:set_ticks]([-2.0, -1.0, 0.0, 1.0])
PyPlot.scatter(x_doublet, y_doublet, color="#CD2305", s=80, marker="o")
PyPlot.contour(X, Y, psi, levels=[0.], colors="#CD2305", linewidths=2, linestyles="solid");
PyPlot.scatter([x_stagn1, x_stagn2], [y_stagn1, y_stagn2], color="g", s=80, marker="o");

gamma = 4.0                      # strength of the vortex
x_vortex, y_vortex = 0.0, 0.0    # location of the vortex

function get_velocity_sink(strength, xs, ys, X, Y):
    """Returns the velocity field generated by a sink.

    Arguments
    ---------
    strength -- strength of the sink.
    xs, ys -- coordinates of the sink.
    X, Y -- mesh grid.
    """
    u = strength/(2*pi).*(X-xs)./((X-xs).^2+(Y-ys).^2)
    v = strength/(2*pi).*(Y-ys)./((X-xs).^2+(Y-ys).^2)

    return u, v
  end

function get_stream_function_sink(strength, xs, ys, X, Y):
      """Returns the stream-function generated by a sink.

      Arguments
      ---------
      strength -- strength of the sink.
      xs, ys -- coordinates of the sink.
      X, Y -- mesh grid.
      """
      psi = strength/(2*pi)*atan2((Y-ys), (X-xs))

      return psi
    end

    function get_velocity_vortex(strength, xv, yv, X, Y):
          """Returns the velocity field generated by a vortex.

          Arguments
          ---------
          strength -- strength of the vortex.
          xv, yv -- coordinates of the vortex.
          X, Y -- mesh grid.
          """
          u = + strength/(2*pi).*(Y-yv)./((X-xv).^2+(Y-yv).^2)
          v = - strength/(2*pi).*(X-xv)./((X-xv).^2+(Y-yv).^2)

          return u, v
        end

    function get_stream_function_vortex(strength, xv, yv, X, Y):
        """Returns the stream-function generated by a vortex.

        Arguments
        ---------
        strength -- strength of the vortex.
        xv, yv -- coordinates of the vortex.
        X, Y -- mesh grid.
        """
        psi = strength/(4*pi).*log((X-xv).^2+(Y-yv).^2)

        return psi
      end

    # computes the velocity field on the mesh grid
u_vortex, v_vortex = get_velocity_vortex(gamma, x_vortex, y_vortex, X, Y)

# computes the stream-function on the mesh grid
psi_vortex = get_stream_function_vortex(gamma, x_vortex, y_vortex, X, Y)

# superposition of the doublet and the vortex on the freestream flow
u = u_freestream + u_doublet + u_vortex
v = v_freestream + v_doublet + v_vortex
psi = psi_freestream + psi_doublet + psi_vortex

# calculates cylinder radius
R = sqrt(kappa/(2*pi*u_inf))

# calculates the stagnation points
x_stagn1, y_stagn1 = +sqrt(R^2-(gamma/(4*pi*u_inf))^2), -gamma/(4*pi*u_inf)
x_stagn2, y_stagn2 = -sqrt(R^2-(gamma/(4*pi*u_inf))^2), -gamma/(4*pi*u_inf)

# plots the streamlines
size = 10
PyPlot.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))
PyPlot.xlabel("x", fontsize=16)
PyPlot.ylabel("y", fontsize=16)
PyPlot.xlim(x_start, x_end)
PyPlot.ylim(y_start, y_end)
PyPlot.streamplot(X, Y, u, v, density=2, linewidth=1, arrowsize=1.5, arrowstyle="->")
#circle = PyPlot.Circle((0, 0), radius=R, color="#CD2305", alpha=0.5)
#PyPlot.gca().add_patch(circle)
PyPlot.scatter(x_vortex, y_vortex, color="#CD2305", s=13000, marker="o")
PyPlot.scatter([x_stagn1, x_stagn2], [y_stagn1, y_stagn2], color="g", s=80, marker="o");

# calculates the surface tangential velocity on the cylinder
theta = linspace(0, 2*pi, 100)
u_theta = -2*u_inf*sin(theta) - gamma/(2*pi*R)

# computes the surface pressure coefficient
cp = 1.0 - (u_theta./u_inf).^2

# if there was no vortex
u_theta_no_vortex = -2*u_inf*sin(theta)
cp_no_vortex = 1.0 - (u_theta_no_vortex./u_inf).^2

q=0.5*1.225*u_inf^2
drag1 = - cp.*q .* 2*pi .* sin(2*pi) .* R
DRAG = sum(drag1)

lift1  = cp.*q .* 2*pi .* cos(2*pi) .* R
LIFT = sum(lift1)

println(DRAG)
println(LIFT)

# plots the surface pressure coefficient
size = 6
PyPlot.figure(figsize=(size, size))
PyPlot.grid(true)
PyPlot.xlabel(L"$theta$", fontsize=18)
PyPlot.ylabel(L"$C_p$", fontsize=18)
PyPlot.xlim(minimum(theta), maximum(theta))
PyPlot.plot(theta, cp, color="#CD2305", linewidth=2, linestyle="-", label = "With Vortex")
PyPlot.plot(theta, cp_no_vortex, color="g", linewidth=2, linestyle="-",label = "Without Vortex")
PyPlot.legend(loc = "best");


text:
using RecyclePermanent
using RNGPool
import NonUniformRandomVariateGeneration.sampleCategorical
using BenchmarkTools
using Combinatorics

include("compare.jl")

n = 100
N = 1000

A = exp.(randn(n, N)) ./ N

trials = 1000

setRNGs(12345)

function permanent(A::Matrix)
  n, N = size(A)
  perms = permutations(1:N, n)
  v::Float64 = 0.0
  for σ in perms
    u::Float64 = 1.0
    for i in 1:n
      u *= A[i, σ[i]]
    end
    v += u
  end
  return v
end

function kuznetsovLogPermanent(A::Matrix, ks::Vector{Int64},
  probs::Vector{Float64})
  n, N = size(A)
  v::Float64 = 0.0
  for p in 1:n
    for i in 1:N
      @inbounds probs[i] = A[p,i]
    end
    for j in 1:p-1
      @inbounds probs[ks[j]] = 0.0
    end
    sp::Float64 = sum(probs)
    if sp == 0 return -Inf end
    k::Int64 = sampleCategorical(probs, sp)
    v += log(sp)
    p < n && (ks[p] = k)
  end
  return v
end

function kuznetsovLogPermanent(A::Matrix)
  n, N = size(A)
  ks = Vector{Int64}(undef, n-1)
  probs = Vector{Float64}(undef, N)
  return kuznetsovLogPermanent(A, ks, probs)
end

# too slow for large n, N
# lp = log(permanent(A))

kuznetsovOutput = testEstimator(() -> kuznetsovLogPermanent(A), trials)
kmOutput = testEstimator(() -> RecyclePermanent.kmLogPermanent(A), trials)

@btime kuznetsovLogPermanent($A)
@btime RecyclePermanent.kmLogPermanent($A)


text:
using RandomStuff
using Test
using Random
using TestSetExtensions
# using SafeTestsets

@testset ExtendedTestSet "RandomStuff.jl" begin

  # @testset "Random Strings" begin
  #   @test isa(rand(String), String)
  #   @test length(rand(String)) == 8
  #   @test length(rand(String, 10)) == 10

  #   @test length(rand(String)) == 10
  #   @test rand(String) == "abcdefgh"
  # end

  # @testset "Random Symbols" begin
  #   @test isa(rand(Symbol), Symbol)
  #   @test length(rand(Symbol) |> string) == 8
  #   @test length(rand(Symbol, 10) |> string) == 10
  # end

  # @testset "Random Pairs" begin
  #   @test isa(rand(Pair{String,Int}), Pair{String,Int})
  # end

  @includetests ARGS
end

# @safetestset "Value Tests" begin include("tests_values.jl") end
# @safetestset StringTests = "String Tests" begin include("tests_strings.jl") end


text:
using GPUifyLoops: @launch, @loop
using Oceananigans.Operators: Δy, ΔzF
using Oceananigans.Utils: @loop_xy, @loop_xz, launch_config

#####
##### Algorithm for adding fluxes associated with non-trivial flux boundary conditions.
##### Inhomogeneous Value and Gradient boundary conditions are handled by filling halos.
#####

"""
    apply_z_bcs!(Gc, arch, grid, args...)

Apply flux boundary conditions to a field `c` by adding the associated flux divergence to
the source term `Gc` at the top and bottom.
"""
function apply_z_bcs!(Gc, c, arch, args...)
    @launch(device(arch), config=launch_config(Gc.grid, :xy),
            _apply_z_bcs!(Gc.data, Gc.grid, c.boundary_conditions.z.bottom,
                          c.boundary_conditions.z.top, args...))
    return nothing
end

"""
    apply_y_bcs!(Gc, arch, grid, args...)

Apply flux boundary conditions to a field `c` by adding the associated flux divergence to
the source term `Gc` at the left and right.
"""
function apply_y_bcs!(Gc, c, arch, args...)
    @launch(device(arch), config=launch_config(Gc.grid, :xz),
            _apply_y_bcs!(Gc.data, Gc.grid, c.boundary_conditions.y.left,
                          c.boundary_conditions.y.right, args...))
    return nothing
end

"""
    _apply_z_bcs!(Gc, grid, bottom_bc, top_bc, args...)

Apply a top and/or bottom boundary condition to variable `c`.
"""
function _apply_z_bcs!(Gc, grid, bottom_bc, top_bc, args...)
    @loop_xy i j grid begin
        apply_z_bottom_bc!(Gc, bottom_bc, i, j, grid, args...)
           apply_z_top_bc!(Gc, top_bc,    i, j, grid, args...)
    end
end

"""
    _apply_z_bcs!(Gc, grid, left_bc, right_bc, args...)

Apply a left and/or right boundary condition to variable `c`.
"""
function _apply_y_bcs!(Gc, grid, left_bc, right_bc, args...)
    @loop_xz i k grid begin
         apply_y_left_bc!(Gc, left_bc,  i, k, grid, args...)
        apply_y_right_bc!(Gc, right_bc, i, k, grid, args...)
    end
end

# Fall back functions for boundary conditions that are not of type Flux.
@inline apply_z_top_bc!(args...) = nothing
@inline apply_z_bottom_bc!(args...) = nothing

@inline apply_y_right_bc!(args...) = nothing
@inline apply_y_left_bc!(args...) = nothing

# Shortcuts for 'no-flux' boundary conditions.
@inline apply_z_top_bc!(Gc, ::NFBC, args...) = nothing
@inline apply_z_bottom_bc!(Gc, ::NFBC, args...) = nothing

@inline apply_y_right_bc!(Gc, ::NFBC, args...) = nothing
@inline apply_y_left_bc!(Gc, ::NFBC, args...) = nothing

"""
    apply_z_top_bc!(Gc, top_flux::BC{<:Flux}, i, j, grid, args...)

Add the part of flux divergence associated with a top boundary condition on `c`.
Note that because

    `tendency = ∂c/∂t = Gc = - ∇ ⋅ flux`

a positive top flux is associated with a *decrease* in `Gc` near the top boundary.
If `top_bc.condition` is a function, the function must have the signature

    `top_bc.condition(i, j, grid, boundary_condition_args...)`
"""
@inline apply_z_top_bc!(Gc, top_flux::BC{<:Flux}, i, j, grid, args...) =
    @inbounds Gc[i, j, grid.Nz] -= getbc(top_flux, i, j, grid, args...) / ΔzF(i, j, grid.Nz, grid)

@inline apply_y_right_bc!(Gc, right_flux::BC{<:Flux}, i, k, grid, args...) =
    @inbounds Gc[i, grid.Ny, k] -= getbc(right_flux, i, k, grid, args...) / Δy(i, grid.Ny, k, grid)

"""
    apply_z_bottom_bc!(Gc, bottom_flux::BC{<:Flux}, i, j, grid, args...)

Add the flux divergence associated with a bottom flux boundary condition on `c`.
Note that because

    `tendency = ∂c/∂t = Gc = - ∇ ⋅ flux`

a positive bottom flux is associated with an *increase* in `Gc` near the bottom boundary.
If `bottom_bc.condition` is a function, the function must have the signature

    `bottom_bc.condition(i, j, grid, boundary_condition_args...)`
"""
@inline apply_z_bottom_bc!(Gc, bottom_flux::BC{<:Flux}, i, j, grid, args...) =
    @inbounds Gc[i, j, 1] += getbc(bottom_flux, i, j, grid, args...) / ΔzF(i, j, 1, grid)

@inline apply_y_left_bc!(Gc, left_flux::BC{<:Flux}, i, k, grid, args...) =
    @inbounds Gc[i, 1, k] += getbc(left_flux, i, k, grid, args...) / Δy(i, 1, k, grid)


text:
#
# children() and related functions
#

"""
    children(node)

Get the immediate children of `node`.

This is the primary function that needs to be implemented for custom tree types. It should return an
iterable object for which an appropriate implementation of `Base.pairs` is available.

The default behavior is to assume that if an object is iterable, iterating over
it gives its children. Non-iterable types are treated as leaf nodes.
"""
children(node) = Base.isiterable(typeof(node)) ? node : ()

"""
    ischild(node1, node2)

Check if `node1` is a child of `node2`.

By default this iterates through ``children(node2)``, so performance may be improved by adding a
specialized method for given node type.
"""
ischild(node1, node2) = any(node -> node === node1, children(node2))

#
# Subtree-related functions
#

"""
    intree(node, root)

Check if `node` is a member of the tree rooted at `root`.

By default this traverses through the entire tree in search of `node`, and so may be slow if a
more specialized method has not been implemented for the given tree type.

See also: [`isdescendant`](@ref)
"""
intree(node, root) = any(n -> n === node, PreOrderDFS(root))

"""
    isdescendant(node1, node2)

Check if `node1` is a descendant of `node2`. This isequivalent to checking whether `node1` is a
member of the subtree rooted at `node2` (see [`intree`](@ref)) except that a node cannot be a
descendant of itself.

Internally this calls `intree(node1, node2)` and so may be slow if a specialized method of that
function is not available.
"""
isdescendant(node1, node2) = node1 !== node2 && intree(node1, node2)

"""
    treesize(node)

Get the size of the tree rooted at `node`.

By default this recurses through all nodes in the tree and so may be slow if a more specialized
method has not been implemented for the given type.
"""
treesize(node) = 1 + mapreduce(treesize, +, children(node), init=0)

"""
    treebreadth(node)

Get the number of leaves in the tree rooted at `node`. Leaf nodes have a breadth of one.

By default this recurses through all nodes in the tree and so may be slow if a more specialized
method has not been implemented for the given type.
"""
treebreadth(node) = isempty(children(node)) ? 1 : mapreduce(treebreadth, +, children(node))

"""
    treeheight(node)

Get the maximum depth from `node` to any of its descendants. Leaf nodes have a height of zero.

By default this recurses through all nodes in the tree and so may be slow if a more specialized
method has not been implemented for the given type.
"""
treeheight(node) = isempty(children(node)) ? 0 : 1 + mapreduce(treeheight, max, children(node))


text:
# reference: https://www.learndatasci.com/tutorials/reinforcement-q-learning-scratch-python-openai-gym/

using PyCall
gym = pyimport("gym")

env = gym.make("Taxi-v3").env

# Training

Q = zeros(env.observation_space.n, env.action_space.n)

alpha = 0.1
gamma = 0.6
epsilon = 0.1

all_epochs = []
all_penalties = []

for i = 1:10000
    state = env.reset()

    epochs, penalties, reward, = 0, 0, 0
    done = false

    while !done
        if rand() < epsilon
            action = env.action_space.sample() # Explore action space
        else
            action = argmax(Q[state+1,:])-1 # Exploit learned values
        end

        next_state, reward, done, info = env.step(action)

        next_max = maximum(Q[next_state+1, :])

        Q[state+1, action+1] = (1-alpha) * Q[state+1, action+1] + alpha * (reward + gamma * next_max)

        if reward == -10
            penalties += 1
        end

        state = next_state
        epochs += 1

    end

    if mod(i, 100)==0
        println("Episode = $i")
    end
end

# Testing the learned Q function on a new environment

state = env.reset()
epochs, penalties, reward = 0, 0, 0

done = false

frames = []
while !done
    action = argmax(Q[state+1,:]) - 1
    state, reward, done, info = env.step(action)

    if reward == -10
        penalties += 1
    end

    push!(frames, Dict(
            "frame"=> env.render(mode="ansi"),
            "state"=> state,
            "action"=> action,
            "reward"=> reward
        )
    )
    epochs += 1
end

# Visualize the result
for (k,f) in enumerate(frames)
    for _ = 1:9
        print("\033[F\033[K")
    end
    print(f["frame"])
    print("Time = $k\n")
    sleep(0.1)
end