text:
using Test
using LinearAlgebra
using PolynomialBasis
using ImplicitDomainQuadrature
using HDGElasticity

function curve_length(imap,quad)
    l = 0.0
    for (p,w) in quad
        detjac = HDGElasticity.determinant_jacobian(imap,p)
        l += detjac*w
    end
    return l
end

function curve_length(imap,quad,facescale)
    l = 0.0
    for (idx,(p,w)) in enumerate(quad)
        detjac = HDGElasticity.determinant_jacobian(imap,p)
        l += detjac*facescale[idx]*w
    end
    return l
end

function curve_normals(imap,points)
    npts = length(points)
    normals = zeros(2,npts)
    for (idx,x) in enumerate(points)
        t = gradient(imap,x)
        n = [-t[2],t[1]]
        normals[:,idx] .= n/norm(n)
    end
    return normals
end

imap = InterpolatingPolynomial(2,1,2)
icoeffs = [-1.,-1.,0.,0.,1.,-1.]
update!(imap,icoeffs)
quad = tensor_product_quadrature(1,5)
refl = curve_length(imap,quad)
testlength = sqrt(5) + 0.5*asinh(2.0)
@test isapprox(refl,testlength,atol=1e-2)

cellmap = HDGElasticity.CellMap([0.,0.],[2.,1.])
refnormals = curve_normals(imap,quad.points)
invjac = HDGElasticity.inverse_jacobian(cellmap)
transformed_normals = diagm(invjac)*refnormals
magnitude = mapslices(x->norm(x),transformed_normals,dims=1)
normals = hcat([transformed_normals[:,i]/magnitude[i] for i = 1:size(magnitude)[2]]...)
jac = HDGElasticity.jacobian(cellmap)
reftangents = hcat([gradient(imap,quad.points[:,i]) for i = 1:length(quad)]...)
tangents = diagm(jac)*reftangents

facescale = HDGElasticity.scale_area(cellmap,normals)
transformed_length = curve_length(imap,quad,facescale)
testlength = sqrt(2) + asinh(1)
@test isapprox(transformed_length,testlength,atol=1e-2)

ext: jl


text:
# ==========================================================================
# Unprocessed Solution for the primal problem

immutable PrimalSolution{T<:AbstractArray}
  minimizer::T
  minimum::Float64
  iterations::Int
  isconverged::Bool
end

function PrimalSolution{T<:AbstractArray}(
    minimizer::T,
    minimum::Float64,
    iterations::Int,
    isconverged::Bool)
  PrimalSolution{T}(minimizer, minimum, iterations, isconverged)
end

function Base.show(io::IO, r::PrimalSolution)
  println(io, typeof(r))
  _printvariable(io, 14, "minimizer()", minimizer(r))
  _printvariable(io, 14, "minimum()", minimum(r))
  _printvariable(io, 14, "iterations()", iterations(r))
  _printvariable(io, 14, "isconverged()", isconverged(r), newline = false)
end

coef(s::PrimalSolution) = s.minimizer
minimizer(s::PrimalSolution) = s.minimizer
Base.minimum(s::PrimalSolution) = s.minimum
iterations(s::PrimalSolution) = s.iterations
isconverged(s::PrimalSolution) = s.isconverged


text:
#Note: you have to ensure, when in an interactive session (e.g. atom)
# That BEDFunc isn't loaded from somewhere else and you are
#  not getting conflicts with this.
using Base.Test
DIR = dirname(dirname(@__FILE__()))
path=joinpath(DIR,"src/gtf.jl")
include(path)
filename = joinpath(DIR,"testdata/gencode.v25.annotation.100.gtf")
df=parse_gtf_gene_features_to_df(filename, [], ["gene"])
@test nrow(df) == 13

filename = joinpath(DIR,"testdata/gencode.v25.annotation.100.gtf.gz")
df=parse_gtf_gene_features_to_df(filename, [], ["gene"])
@test nrow(df) == 13


text:
using NumericalMethodsforEngineers
using Test

# write your own tests here

code_tests = [
  # Ch01 tests
  "ch01_nmex.jl",

  # Ch02 tests
  "ch02_lufac.jl",
  "ch02_ldlt.jl",
  #"ch02_SymmetricBandedMatrix.jl",
  #"ch02_SkylineMatrix.jl",
  "ch02_gaussian.jl",
  "ch02_penalized.jl",
  "ch02_jacobi.jl",
  "ch02_bicgstab.jl",                               # Initial version of bicgstab

  # Ch03 tests
  # Ch04 tests

  "ch04_symata.jl",

  # Ch05 tests

  # Ch06 tests
  "ch06_p6.2.jl",

  # Ch07 tests
  "ch07_ivp_ex7.3.jl",                           # Four direct methods
  "ch07_bvp_ex7.4a.jl",                        # Shooting method example 7.4a
  "ch07_bvp_ex7.4b.jl",                        # Shooting method example 7.4b

  # Symata Weighted Residual Method tests
  "wrm_with_using_symata.jl",
  "wrm_with_using_nmfe.jl",
  "../examples/ch07/WRM/Ex.7.13.jl",         # 1 point LagrangePolynomial
  "../examples/ch07/WRM/Ex.7.14.jl",     # 2 point LagrangePolynomial
  "../examples/ch07/WRM/Ex.7.15.jl",     # 1 point collocation
  "../examples/ch07/WRM/Ex.7.15b.jl",    # 2 point collocation
  "../examples/ch07/WRM/Ex.7.16.jl",     # 1 point subdomain
  "../examples/ch07/WRM/Ex.7.16b.jl",    # 2 point subdomain
  "../examples/ch07/WRM/Ex.7.17.jl",     # 1 point least-squares
  "../examples/ch07/WRM/Ex.7.17b.jl",    # 2 point least-squares
  "../examples/ch07/WRM/Ex.7.18.jl",     # 1 point galerkin
  "../examples/ch07/WRM/Ex.7.18b.jl",   # 2 point galerkin

  # Ch08 tests
  "ch08_ex01.jl",

  # Extra tests
  "vianello_symata.jl"                           # Vianello procedure

]

println("Running tests:")

@testset "NumericalMethodsforEngineers.jl" begin
  for my_test in code_tests
      println("\n  * $(my_test) *")
      include("$(my_test)")
  end
  println()
end


text:
using HPAT

@acc hpat function q05(category, education, gender, file_name)
    web_clickstreams = DataSource(DataTable{:wcs_user_sk=Int64,:wcs_item_sk=Int64}, HDF5, file_name)
    item = DataSource(DataTable{:item_sk=Int64,:i_category=ASCIIString,:i_category_id=Int64}, HDF5, file_name)
    customer = DataSource(DataTable{:c_customer_sk=Int64,:c_current_cdemo_sk=Int64}, HDF5, file_name)
    customer_demographics = DataSource(DataTable{:cd_demo_sk=Int64,:cd_education_status=ASCIIString,:cd_gender=ASCIIString}, HDF5, file_name)

    user_items = join(web_clickstreams, item, :wcs_item_sk==:i_item_sk, :item_sk)

    user_clicks_in_cat = aggregate(user_items, :item_sk, :clicks_in_category = sum(:i_category==category),
                                                         :clicks_in_1 = sum(:i_category_id==1),
                                                         :clicks_in_2 = sum(:i_category_id==2),
                                                         :clicks_in_3 = sum(:i_category_id==3),
                                                         :clicks_in_4 = sum(:i_category_id==4),
                                                         :clicks_in_5 = sum(:i_category_id==5),
                                                         :clicks_in_6 = sum(:i_category_id==6),
                                                         :clicks_in_7 = sum(:i_category_id==7))
    customer_clicks = join(user_clicks_in_cat, customer, :wcs_user_sk==:c_customer_sk)
    customer_demo_clicks = join(customer_clicks, customer_demographics, :c_current_cdemo_sk==:cd_demo_sk)
    data = hcat( convert(Vector{Int64}, customer_demo_clicks[:cd_education_status].==education),
            convert(Vector{Int64}, customer_demo_clicks[:cd_gender].==gender),
            customer_demo_clicks[:clicks_in_1],
            customer_demo_clicks[:clicks_in_2],
            customer_demo_clicks[:clicks_in_3],
            customer_demo_clicks[:clicks_in_4],
            customer_demo_clicks[:clicks_in_5],
            customer_demo_clicks[:clicks_in_6],
            customer_demo_clicks[:clicks_in_7])

    model = HPAT.API.LogisticRegression(data, customer_demo_clicks[:clicks_in_category])
    return model
end

println(q05("food", "colledge", "male", "data.hdf5"))


text:
include("core.jl")
include("arrayops.jl")
include("broadcast.jl")
include("fuzz.jl")


text:
using ReachabilityAnalysis
using ReachabilityAnalysis: center
using Plots
using LaTeXStrings

# w = 2pi/T
A = [0.0 1; -4π 0]
X0 = BallInf([0.0, 10.0], 0.01)
X = Universe(2)
b = [1.0, -1.0]
B = hcat(b)
u = -3.0
U = ConstantInput(linear_map(B, Singleton([u])))

P = @ivp(x' = A*x + u, x ∈ X, u ∈ U, x(0) ∈ X0)

# homogenize system
α = 1.0
A_homog = [A b*u/α; zeros(1, 3)]
X0_homog = cartesian_product(X0, Singleton([α]))
X_homog = Universe(3)
P_homog = @ivp(x' = A_homog*x, x ∈ X_homog, x(0) ∈ X0_homog)
P = P_homog

# flowpipe (tight approximation)
δ = 0.00001
T = 0.05
T_full = 0.15

alg_fw = GLGM06(δ=δ, approx_model=Forward())
sol_fw = solve(P; T=T, alg=alg_fw);

# discretizations
Ω0_chull = ConvexHullArray(map(set, sol_fw))
Ω0_ch = discretize(P, T, CorrectionHull(order=4)).x0

# random trajectories
x0s = vcat(sample(X0, 0, include_vertices=true),
           sample(X0, 50, sampler=LazySets.FaceSampler(1)))

P_orbit = [@ivp(x' = A*x + u, x ∈ X, u ∈ U, x(0) ∈ x0) for x0 in x0s]
sol_orbit = [solve(Pi, T=T_full, alg=ORBIT(δ=T/20.)) for Pi in P_orbit]

fig = plot(legend=:bottomleft, legendfontsize=15, xticks=[], yticks=[], size=(900, 400))

Φ = exp(A_homog * T)

sol_fw_ch = [Projection(Φ^i * Ω0_ch, [1, 2]) for i in 0:2]
sol_fw_chull = [ConvexHullArray([Projection(Φ^i * X, [1, 2]) for X in array(Ω0_chull)]) for i in 0:2]

[plot!(fig, s, vars=(1, 2), alpha=0.2, c=:magenta, seriestype=:path, marker=:none) for s in sol_orbit]
plot!(fig, sol_fw_ch[2:end], alpha=0, color=:blue, linecolor=:blue, lw=2, linealpha=1)
plot!(fig, sol_fw_chull[2:end], color=:purple, linecolor=:purple, lw=2, linealpha=1)
plot!(fig, sol_fw_ch[1], alpha=0, color=:blue, linecolor=:blue, lw=2, linealpha=1)
plot!(fig, sol_fw_chull[1], 1e-4, alpha=0.3, color=:purple, linecolor=:purple, lw=2, linealpha=1)

plot!(fig, X0, c=:yellow, alpha=1., lab="")
plot!(fig, sol_fw[1501:2000], vars=(1, 2), lw=0.1, c=:orange, alpha=0.8, lab="")
plot!(fig, sol_fw[end], vars=(1, 2), c=:green, alpha=0.8, lab="")

p = center(box_approximation(sol_fw_ch[1]))
plot!(fig, color=:black, annotations=(p[1], p[2] - .12, text(L"\Omega_0", 22)))
p = center(box_approximation(sol_fw_ch[2]))
plot!(fig, color=:black, annotations=(p[1] - .025, p[2] - .1, text(L"\Omega_1", 22)))
p = center(box_approximation(sol_fw_ch[3]))
plot!(fig, color=:black, annotations=(p[1] - .05, p[2] - .07, text(L"\Omega_2", 22)))

savefig(fig, "intro_recurrence.pdf")


text:
cd(@__DIR__); using Pkg; Pkg.activate(".");

using Revise, StableDQMC, GenericSVD, JacobiSVD, LinearAlgebra, Plots, Statistics, LaTeXStrings, DataStructures
using BenchmarkTools, StatsPlots, DataFrames, Infiltrator

pyplot();
# PyPlot.plt.style.use("default")
PyPlot.plt.style.use("publication_tex")
default(
    framestyle = :grid,
    grid = true,
    guidefont=font(15),
#     framestyle = :box,
#     grid = false,
#     guidefont=font(18),
    legend = true,
    size = (1.6*360, 360),
    xtickfont=font(13),
    ytickfont=font(13),
    legendfont=font(12),
    linewidth = 1.5
);

## Model

# Hubbard
L = 2
nsites = L^2
n = nsites * 2 # spin
t = 1
U = 6
μ = 0
Δτ = 0.1

# Kinetic part (PBC)
T = diagm(-1 => fill(-t, n-1), 1 => fill(-t, n-1), 0 => fill(μ, n))
T[1,end] = -t; T[end,1] = -t;

# on-site interaction
λ = acosh(exp(0.5*U*Δτ))
eV = Diagonal(vcat(fill(exp(λ), nsites), fill(exp(-λ), nsites)))

B = exp(-Δτ .* T) * eV

println()
@info "Hubbard model: U = $U, L = $L"
@info "Condition number of B: $(cond(B))"

## Bchain product

Bbig = BigFloat.(B)
N = 500

setprecision(100) do
    global svs_genericsvd = calc_Bchain_svd(Bbig, N; svdalg = genericsvd)[2];
end
svs_qr_udt = calc_Bchain_qr(B, N)[2];
svs_gesvd  = calc_Bchain_svd(B, N; svdalg = gesvd)[2];
svs_gesdd  = calc_Bchain_svd(B, N; svdalg = gesdd)[2];
svs_gesvj  = calc_Bchain_svd(B, N; svdalg = gesvj)[2];

betas = range(1, N * 0.1, length=N)
plot(betas, svs_genericsvd, color="#CA0020", labels = reshape(["exact", "", "","", "","", "", ""], (1,8)))
plot!(betas, svs_qr_udt, color="#E66101", labels = reshape(["QR", "", "","", "","", "", ""], (1,8)))
plot!(betas, svs_gesvj, color="#FDB863", labels = reshape(["SVD (Jacobi)", "", "","", "","", "", ""], (1,8)))
plot!(betas, svs_gesvd, color="#008837", labels = reshape(["SVD", "", "","", "","", "", ""], (1,8)))
plot!(betas, svs_gesdd, color="#5E3C99", labels = reshape(["SVD (D\\&C)", "", "","", "","", "", ""], (1,8)))

ylabel!("log singular values")
xlabel!(L"inverse temperature $\beta$")
# ylims!(-160, 100)
xlims!(0., 50.)
# savefig("../paper/figures/decomp_comparison_simple.pdf")

## Equal-time Green's function

function accuracy_greens(B; svdinversion = inv_one_plus)
    Bbig = BigFloat.(B)

    Δs_qr    = SortedDict{Int, Float64}()
    Δs_gesvd = SortedDict{Int, Float64}()
    Δs_gesdd = SortedDict{Int, Float64}()
    Δs_gesvj = SortedDict{Int, Float64}()

    for beta in range(5, 40, step=5)
        N = Int(beta * 10)

        local G_genericsvd
        setprecision(1000) do
            F_genericsvd = calc_Bchain_svd(Bbig, N; svdalg = genericsvd)[1];
            G_genericsvd = inv(I + Matrix(F_genericsvd))
        end

        # Bchain
        F_qr_udt = calc_Bchain_qr(B, N)[1];
        F_gesvd  = calc_Bchain_svd(B, N; svdalg = gesvd)[1];
        F_gesdd  = calc_Bchain_svd(B, N; svdalg = gesdd)[1];
        F_gesvj  = calc_Bchain_svd(B, N; svdalg = StableDQMC.gesvj)[1];

        # inversion
        G_qr_udt = inv_one_plus(F_qr_udt)
        G_gesvd  = svdinversion(F_gesvd; svdalg = gesvd!)
        G_gesdd  = svdinversion(F_gesdd; svdalg = gesdd!)
        G_gesvj  = svdinversion(F_gesvj; svdalg = StableDQMC.gesvj!)

        Δs_qr[beta]    = maximum(abs.(G_qr_udt - G_genericsvd))
        Δs_gesvd[beta] = maximum(abs.(G_gesvd - G_genericsvd))
        Δs_gesdd[beta] = maximum(abs.(G_gesdd - G_genericsvd))
        Δs_gesvj[beta] = maximum(abs.(G_gesvj - G_genericsvd))
    end

    return Δs_qr, Δs_gesvd, Δs_gesdd, Δs_gesvj
end

# Δs_qr, Δs_gesvd, Δs_gesdd, Δs_gesvj = Juno.@enter accuracy_greens(B);
Δs_qr, Δs_gesvd, Δs_gesdd, Δs_gesvj = accuracy_greens(B);

##

betas = collect(keys(Δs_qr))
p = plot(betas, log.(10, values(Δs_qr)), color="#E66101", marker=true, label="QR")
plot!(betas, log.(10, values(Δs_gesvd)), color="#008837", marker=true, label="SVD")
plot!(betas, log.(10, values(Δs_gesdd)), color="#5E3C99", marker=true, label="SVD (D\\&C)")
plot!(betas, log.(10, values(Δs_gesvj)), color="#FDB863", marker=true, label="SVD (Jacobi)")

xlabel!(L"inverse temperature $\beta$")
ylabel!(L"\log(\textrm{max}(\textrm{abs}(G - G_{\textrm{exact}})))")
# ylims!(0., 1e-3)
savefig("../paper/figures/accuracy_greens_svd_regularinv_U_$(U)_L_$(L)_t_$(t).pdf")

## Time-displaced Green's function

# G(slice, 1) = [B(slice, 1)^-1 + B(beta, slice)]^-1
# G(slice, 1) = [B(slice, 1)^-1 + B(beta, slice)]^-1
function accuracy_tdgf(B; loh = false)
    Bbig = BigFloat.(B)

#     betas = range(5, 40, step=5)
    beta = 40
    M = Int(beta * 10)
    slices = 1:10:M

    df = DataFrame(slice=Int[], qr=Float64[], gesdd=Float64[], gesvd=Float64[], gesvj=Float64[])

    inv_sum_method = loh ? inv_sum_loh : inv_sum

    for s in slices
        local TDGF_genericsvd
        setprecision(1000) do
            TDGF_genericsvd = calc_tdgf_svd(Bbig, s, M-s, inv_sum_method = inv_sum_method, svdalg_chain = genericsvd, svdalg_inv = genericsvd!)
        end

        TDGF_qr = calc_tdgf_qr(B, s, M-s, inv_sum_method = inv_sum_method)
        TDGF_gesdd = calc_tdgf_svd(B, s, M-s, inv_sum_method = inv_sum_method, svdalg_chain = gesdd, svdalg_inv = gesdd!)
        TDGF_gesvd = calc_tdgf_svd(B, s, M-s, inv_sum_method = inv_sum_method, svdalg_chain = gesvd, svdalg_inv = gesvd!)
        TDGF_gesvj = calc_tdgf_svd(B, s, M-s, inv_sum_method = inv_sum_method, svdalg_chain = StableDQMC.gesvj, svdalg_inv = StableDQMC.gesvj!)

        Δ_qr    = maximum(abs, TDGF_qr - TDGF_genericsvd)
        Δ_gesdd = maximum(abs, TDGF_gesdd - TDGF_genericsvd)
        Δ_gesvd = maximum(abs, TDGF_gesvd - TDGF_genericsvd)
        Δ_gesvj = maximum(abs, TDGF_gesvj - TDGF_genericsvd)

        push!(df, [s, Δ_qr, Δ_gesdd, Δ_gesvd, Δ_gesvj])
    end

    return df
end

df_tdgf = accuracy_tdgf(B; loh = true)

##

@df df_tdgf plot(:slice .* 0.1, (x -> log.(10, x)).([:qr, :gesvd, :gesdd, :gesvj]),
            marker=true,
            color=permutedims(["#E66101", "#008837", "#5E3C99", "#FDB863"]), # #E66101 for QR
            label=permutedims(["QR", "SVD", "SVD (D\\&C)", "SVD (Jacobi)"]),
#             xlims=(0,10)
           )
xlabel!(L"imaginary time $\tau$")
ylabel!(L"\log(\textrm{max}(\textrm{abs}(G(\tau, 0) - G_{\textrm{exact}}(\tau, 0))))")
savefig("../paper/figures/accuracy_tdgf_loh_U_$(U)_L_$(L)_t_$(t).pdf")


text:
# SAM
# ===

function checkkeytag(tag)
    if !isa(tag, AbstractString)
        throw(ArgumentError("tag must be a string"))
    elseif length(tag) != 2
        throw(ArgumentError("tag length must be 2"))
    end
    return
end

const auxtype = Dict{UInt8,DataType}(
    UInt8('A') => Char,
    UInt8('c') => Int8,
    UInt8('C') => UInt8,
    UInt8('s') => Int16,
    UInt8('S') => UInt16,
    UInt8('i') => Int32,
    UInt8('I') => UInt32,
    UInt8('f') => Float32,
    UInt8('d') => Float64,
    UInt8('Z') => String)

const auxtypechar = Dict{DataType,Char}(
    Char    => 'A',
    Int8    => 'c',
    UInt8   => 'C',
    Int16   => 's',
    UInt16  => 'S',
    Int32   => 'i',
    UInt32  => 'I',
    Float32 => 'f',
    Float64 => 'd',
    String  => 'Z')

# Parse hexadecimal byte format string.
function parse_hexbytearray(s)
    return [parse(UInt8, s[i:i+1], 16) for i in 1:2:endof(s)]
end

function parse_keyvals(s)
    ret = Dict{String,Any}()
    for pair in split(s, '\t')
        if pair[3] != ':'
            error("':' is expected")
        end
        ret[pair[1:2]] = pair[4:end]
    end
    return ret
end

include("flags.jl")
include("metainfo.jl")
include("header.jl")
include("record.jl")
include("reader.jl")
include("writer.jl")


text:
# TODO: this should go to ReSHOP!

"""
Add a node of type MP

# Arguments
- emp: the EMP object
- mdl: an opaque pointer to the ReSHOP.reshop_model
- node: the node to process
- reshop_mp_parent: an opaque pointer to the parent of the node
"""
function empnode(emp::EMPmaster, mdl::Ptr{ReSHOP.reshop_model}, node::MathPrgm, reshop_mp_parent::Ptr{ReSHOP.mathprgm}=Ptr{ReSHOP.mathprgm}(C_NULL))
    reshop_mp = ReSHOP.reshop_declare_mathprgm(node, mdl, emp.NLoffset)
    node.solverobj = reshop_mp

#    ReSHOP.reshop_add_mp(reshop_emp, reshop_node, reshop_mp)

    if reshop_mp_parent != C_NULL
        ReSHOP.emp_add_mp_mp(reshop_mp_parent, reshop_mp)
    else
        ReSHOP.emp_set_root(mdl, reshop_mp)
    end

    for mp_mp in node.mps
        empnode(emp, mdl, mp_mp, reshop_mp)
    end

    for mp_equil in node.equils
        empnode(emp, mdl, mp_equil, reshop_mp)
    end

end

"""
Add a node of type equilibrium

# Arguments
- emp: the EMP object
- mdl: an opaque pointer to the ReSHOP.reshop_model
- node: the node to process
- reshop_mp_parent: an opaque pointer to the parent of the node
"""
function empnode(emp::EMPmaster, mdl::Ptr{ReSHOP.reshop_model}, node::Vector{MathPrgm}, reshop_mp_parent::Ptr{ReSHOP.mathprgm}=Ptr{ReSHOP.mathprgm}(C_NULL))
    reshop_mpe = ReSHOP.emp_create_equil(length(node))

    reshop_mps = []
    for mp in node
        reshop_mp = ReSHOP.reshop_declare_mathprgm(mp, mdl, emp.NLoffset)
        mp.solverobj = reshop_mp
        ReSHOP.emp_equil_add(reshop_mpe, reshop_mp)
        push!(reshop_mps, reshop_mp)
    end

    if reshop_mp_parent != C_NULL
        ReSHOP.emp_add_mp_equil(mdl, reshop_node, reshop_mpe)
    else
        ReSHOP.emp_set_root(mdl, reshop_mpe)
    end

    for (mp, reshop_mp) in zip(node, reshop_mps)
        for mp_mp in mp.mps
            empnode(emp, mdl, mp_mp, reshop_mp)
        end

        for mp_equil in mp.equils
            empnode(emp, mdl, mp_equil, reshop_mp)
        end
    end
end

"""
    emptree(emp, mdl)

Pass the EMP graph structure to the EMP solver
"""
function emptree(emp::EMPmaster, mdl::Ptr{ReSHOP.reshop_model})
    if !isnothing(emp.root)
        ReSHOP.emp_mp_ensure(mdl, length(emp.mps))

        empnode(emp, mdl, emp.root)
    end
end


text:
const HTML_boolean_attributes = Set([
    :async
    :autofocus
    :autoplay
    :checked
    :contenteditable
    :controls
    :default
    :defer
    :disabled
    :formNoValidate
    :frameborder
    :hidden
    :ismap
    :itemscope
    :loop
    :multiple
    :muted
    :nomodule
    :novalidate
    :open
    :readonly
    :required
    :reversed
    :scoped
    :selected
    :typemustmatch
  ])

HTML_boolean = Dict{Symbol, Tuple{Vararg{Symbol,N}} where N}(
    :disabled => (:button , :fieldset, :input, :optgroup, :option, :select, :textarea),
    :checked => (:checkbox,)
  )


text:
# Silverman's rule of thumb for KDE bandwidth selection
function bandwidth(data::Vector, alpha::Float64 = 0.9)
    # Determine length of data
    ndata = length(data)

    # Calculate width using variance and IQR
    var_width = std(data)
    q25, q75 = quantile(data, [0.25, 0.75])
    quantile_width = (q75 - q25) / 1.34

    # Deal with edge cases with 0 IQR or variance
    width = min(var_width, quantile_width)
    if width == 0.0
        if var_width == 0.0
            width = 1.0
        else
            width = var_width
        end
    end

    # Set bandwidth using Silverman's rule of thumb
    return alpha * width * ndata^(-0.2)
end

# Store both grid and density for KDE over the real line
immutable UnivariateKDE
    x::Vector
    density::Vector
end

# Algorithm AS 176 for calculating univariate KDE
function kde(data::Vector, npoints::Integer = 512)
    # Determine length of data
    ndata = length(data)

    # Set bandwidth
    window = bandwidth(data)

    # Check that the window is a positive constant
    if window <= 0.0
        error("Window must be positive")
    end

    # Find interval that will contain almost all mass
    dlo = min(data) - 3 * window
    dhi = max(data) + 3 * window

    # Check that interval for estimation is valid
    if dlo >= dhi
        error("Interval must be specified as lower bound, upper bound")
    end

    # Set up a grid for discretized data
    grid = zeros(Float64, npoints)

    # Define some more constants
    step = (dhi - dlo) / npoints
    ainc = 1.0 / (ndata * step)
    npoints2 = fld(npoints, 2)
    hw = window / step
    fac1 = 32.0 * (atan(1.0) * hw / npoints)^2

    # Discretize the data using a histogram
    dlo1 = dlo - step
    for i in 1:ndata
        j = fld(data[i] - dlo1, step)
        if j >= 1 && j <= npoints
            grid[j] = grid[j] + ainc
        end
    end

    # Transform to Fourier basis
    ft = rfft(grid)

    # Find transform of KDE by convolving grid with the
    # Fourier transform of a Gaussian kernel
    for j = 2:length(ft)
        ft[j] *= exp(-fac1 * (j-1)^2)
    end

    # Invert the Fourier transform to get the KDE
    density = irfft(ft, npoints)

    # Fix any noise that crept in
    for j in 1:npoints
        if density[j] < 0.0
            density[j] = 0.0
        end
    end

    # Expand the grid over which KDE was calculated
    x = [dlo:step:(dhi - step)]

    return UnivariateKDE(x, density)
end


text:
function testselection(setup)
    N = setup.num_genes
    guides = Barcode[]

    for (gene, categ) in zip(1:N, [:decreasing, :inactive, :increasing])
        push!(guides, Barcode(gene, 1, gene-2, :linear, categ))
    end

    expand_to = setup.bottleneck_representation * length(guides)

    initial_cells = Array{Int}(expand_to)
    cell_phenotypes = Array{Float64}(size(initial_cells))
    for i in 1:length(initial_cells)
        initial_cells[i] = mod1(i, N)
        cell_phenotypes[i] = guides[initial_cells[i]].theo_phenotype
    end

    # ensure guide level phenotype isn't used so that CRISPRn works
    for guide in guides
        guide.theo_phenotype = -Inf
    end

    num_runs = 100
    results = Array{Float64}(N, num_runs)

    for i in 1:num_runs
        data = select(setup, initial_cells, cell_phenotypes, guides)
        data1 = StatsBase.counts(data[:bin1], 1:N)
        data2 = StatsBase.counts(data[:bin2], 1:N)
        results[:, i] = log2.(data2 ./ data1)
    end
    results
end

function testfacs(σ, quant, width)
    N = 3
    setup = FacsScreen()
    setup.num_genes = N
    setup.coverage = 1
    setup.bottleneck_representation = 20000
    setup.bin_info = Dict{Symbol, Tuple{Float64, Float64}}(:bin1 => (0, width), :bin2 => (1-width, 1))
    setup.σ = σ

    results = testselection(setup)

    ideal = log2(cdf.(Normal(-1, σ), quant)/cdf.(Normal(1,σ), quant))
    tocompare = collect(zip(mean(results, 2), [-ideal, 0.0, ideal]))
    all(Bool[isapprox(x[1], x[2], atol=0.1) for x in tocompare])
end

function testgrowth(ideal, num)
    N = 3
    setup = GrowthScreen()
    setup.num_genes = N
    setup.coverage = 1
    setup.num_bottlenecks = num
    setup.bottleneck_representation = 1000

    results = testselection(setup)
    tocompare = collect(zip(mean(results, 2), log2.(ideal)))
    all(Bool[isapprox(x[1], x[2], atol=0.1) for x in tocompare])
end

begin
    srand(1)
    @test testfacs(0.5, -0.50138209, 1/3)
    @test testfacs(0.5, 0, 1/2)
    @test testfacs(1, -0.58093999, 1/3)

    @test testgrowth([1/7, 2/7, 4/7]./(1/3), 1)
    @test testgrowth([1/21, 4/21,16/21]./(1/3), 2)
    srand()
end


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

using Markdown
using InteractiveUtils

# ╔═╡ de7a5e8a-1ee4-11eb-2bbf-0f0eb3e489ec
using Pkg, DrWatson

# ╔═╡ e4eeb13a-1ee4-11eb-0113-91f9a9c3b659
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, StatsBase, Combinatorics, DataFrames
	Random.seed!(1)
end

# ╔═╡ ac3ceabe-1ee4-11eb-3935-313e96aafd62
md"## Listing 2.8"

# ╔═╡ c4bef1a4-1ee9-11eb-0b58-8be033f6fbad
function bruteSetsProbabilityAllMiss(n)
    omega = collect(permutations(1:n))
    matchEvents = []
    for i in 1:n
        event = []
        for p in omega
            if p[i] == i
                push!(event,p)
            end
        end
        push!(matchEvents,event)
    end
    noMatch = setdiff(omega,union(matchEvents...))
    return length(noMatch)/length(omega)
end

# ╔═╡ e9f9616e-1eeb-11eb-177d-2fa9d0c18b13
formulaCalcAllMiss(n) = sum([(-1)^k/factorial(k) for k in 0:n])

# ╔═╡ e9f99800-1eeb-11eb-3333-59c96cd555d4
function mcAllMiss(n,N)
    function envelopeStuffer()
        envelopes = Random.shuffle!(collect(1:n))
        return sum([envelopes[i] == i for i in 1:n]) == 0
    end
    data = [envelopeStuffer() for _ in 1:N]
    return sum(data)/N
end

# ╔═╡ e9fa266e-1eeb-11eb-2630-59930c474ed3
begin
	N = 10^6
	df = DataFrame()

	for n in 1:6
		bruteForce = bruteSetsProbabilityAllMiss(n)
		fromFormula = formulaCalcAllMiss(n)
		fromMC = mcAllMiss(n,N)
		append!(df, DataFrame(
			:bruteForce => round(bruteForce,digits=4),
			:fromFormula => round(fromFormula,digits=4),
			:fromMC => round(fromMC,digits=4),
			:asymptotic => round(1/MathConstants.e,digits=4))
		)
	end
	df
end

# ╔═╡ ea086874-1eeb-11eb-2cb0-4186697d4294
md"## End of listing 2.8"

# ╔═╡ Cell order:
# ╟─ac3ceabe-1ee4-11eb-3935-313e96aafd62
# ╠═de7a5e8a-1ee4-11eb-2bbf-0f0eb3e489ec
# ╠═e4eeb13a-1ee4-11eb-0113-91f9a9c3b659
# ╠═c4bef1a4-1ee9-11eb-0b58-8be033f6fbad
# ╠═e9f9616e-1eeb-11eb-177d-2fa9d0c18b13
# ╠═e9f99800-1eeb-11eb-3333-59c96cd555d4
# ╠═e9fa266e-1eeb-11eb-2630-59930c474ed3
# ╟─ea086874-1eeb-11eb-2cb0-4186697d4294


text:
@testset "lx_lead" begin
    c = lxmock(raw"""\lead{Hello goodbye baz}""")
    @test lx_lead(c,0) == "@@lead Hello goodbye baz @@"
end

@testset "lx_alert" begin
    c = lxmock("""\\alert{foo **bar** baz}""")
    @test isapproxstr(lx_alert(c,0), """~~~

          ~~~
          foo **bar** baz
          ~~~

~~~""")
end

@testset "lx_figure" begin
    c = lxmock(raw"""\figure{path="/assets/img.png"}""")
    @test isapproxstr(lx_figure(c,0), """~~~

~~~""")
    c = lxmock(raw"""\figure{path="/assets/img.png",width="100%",caption="foo"}""")
    @test isapproxstr(lx_figure(c,0), """~~~

foo

~~~""")
end


text:
# This example is used to generate plots in the paper, "Reflection from a multi-species material and its transmitted effective wavenumber." Proc. R. Soc. (2018): 20170864.

# Everything related to ploting has been commented so that this example can run with requiring the Plots package. That way, this example can be checked everytime the package is modified.

using EffectiveWaves
# using LaTeXStrings

## choose material

  # concrete
  medium = LimeStone
  inclusion1 = AirDry
  inclusion2 = Brick

  ωfactor = 50.0;
  ωs = ωfactor.*LinRange(real(medium.c/10000),real(medium.c),100) # k from 0 to 1

  volfrac = 0.16
  r1 = 0.1/ωfactor; vol1 = 0.06
  r2 = 1.0/ωfactor; vol2 = volfrac - vol1

  sp1 = Specie(ρ=inclusion1.ρ ,r=r1, c=inclusion1.c, volfrac = vol1)
  sp2 = Specie(ρ=inclusion2.ρ ,r=r2, c=inclusion2.c, volfrac = volfrac-vol1)

  # True low partilce volume fraction wavenumber
  kTs  = wavenumber_low_volumefraction(ωs, medium, [sp1,sp2]);
  # Approximate wavenumber
  kTLSs = two_species_approx_wavenumber(ωs, medium, [sp1,sp2]);
  # Approximate challis
  kTCs = wavenumber_challis(ωs, medium, [sp1,sp2]);

  kTs_arr = [kTs,kTLSs,kTCs];
  speed_arr = [ ωs./real(ks) for ks in kTs_arr];
  push!(speed_arr, 0 .*ωs .+ real(medium.c))
  atten_arr = imag.(kTs_arr)

  # styles = [:solid :dashdot :dashdot :dot]
  # labs = [L"Effective $k_{*}$" L"Approximate $k_{*LS}$" L"Approximate $k_{*C}$" "Lime Stone"]
  ys_arr = speed_arr;
  xs = r1.*(ωs./real(medium.c));
  m =5;
  y1 = min(ys_arr[1][1:m]..., ys_arr[2][1:m]...);
  y2 = max(ys_arr[1][1:m]..., ys_arr[2][1:m]...);

  # using Plots
  # height=500
  # pyplot(linewidth=3, size=(2.6*height,height), border=false)

  # Plots.scalefontsizes(1.7)
  # plot(xs, ys_arr, xlabel=L"a_S k", ylabel="sound speed (m/s)", labels=labs, line = styles, xlims = (-0.002,maximum(xs))
        # , ylims = ( min(ys_arr[1]...,ys_arr[3]...,ys_arr[4]...)*0.995,  max(ys_arr[1]...,ys_arr[3]...,ys_arr[4]...)*1.005));
  # p1 = gray_square([xs[1],xs[m]],[y1,y2],plot!);

  ys_arr = atten_arr;
  # labs = [L"Effective $k_{*}$" L"Approximate $k_{*LS}$" L"Approximate $k_{*C}$"]
  # plot(xs, ys_arr, labels=labs, xlabel=L"a_S k", ylabel="attenuation (1/m)", line=styles, xlims = (-0.002,maximum(xs))
        # , ylims = ( min(ys_arr[1]...,ys_arr[3]...)*0.995,  max(ys_arr[1]...,ys_arr[3]...)*1.005));
  y1 = min(ys_arr[1][1:m]..., ys_arr[2][1:m]...);
  y2 = max(ys_arr[1][1:m]..., ys_arr[2][1:m]...);
  # p2 = gray_square([xs[1],xs[m]],[y1,y2],plot!);

  # plot(p1,p2)
  # try mkdir("media") end
  # savefig("media/compare_concrete.png")
  # savefig("media/compare_concrete.pdf")
  # gui()

## Zoomed in version
  ωs = LinRange(ωs[1],ωs[m],250) # k from 0 to 1
  m = length(ωs);
  kTs  = wavenumber_low_volumefraction(ωs, medium, [sp1,sp2]);
  kTLSs = two_species_approx_wavenumber(ωs, medium, [sp1,sp2]);
  kTCs  = wavenumber_challis(ωs, medium, [sp1,sp2]);

  kTs_arr = [kTs,kTLSs,kTCs];
  speed_arr = [ ωs./real(ks) for ks in kTs_arr];
  atten_arr = imag.(kTs_arr)

  # labs = [L"Effective $k_{*}$" L"Approximate $k_{*LS}$" L"Approximate $k_{*C}$"]

  ys_arr = speed_arr;
  xs = r1.*(ωs./real(medium.c));
  # plot(xs, ys_arr, xlabel=L"a_S k", ylabel="sound speed (m/s)", labels=labs
  #                , border = false, line = styles, xlims = (0,maximum(xs)));
  y1 = min(ys_arr[1][1:m]..., ys_arr[2][1:m]...);
  y2 = max(ys_arr[1][1:m]..., ys_arr[2][1:m]...);
  # p1 = gray_square([xs[1],xs[m]],[y1,y2],plot!);

  ys_arr = atten_arr;
  # labs = [L"Effective $k_{*}$" L"Approximate $k_{*LS}$" L"Approximate $k_{*C}$"]
  # plot(xs, ys_arr, xlabel=L"a_S k", ylabel="attenuation (1/m)", labels=labs
  #                , border=false, line = styles, xlims = (0,maximum(xs)));
  y1 = min(ys_arr[1][1:m]..., ys_arr[2][1:m]...);
  y2 = max(ys_arr[1][1:m]..., ys_arr[2][1:m]...);
  dy = abs(y2- y1)/120.0;
  # p2 = gray_square([xs[1],xs[m]],[y1,y2],plot!);

#   plot(p1,p2)
#   savefig("media/compare_concrete_zoom.png")
#   savefig("media/compare_concrete_zoom.pdf")
#
# Plots.scalefontsizes(1/1.7)

max_stars_repo_name: arturgower/EffectiveWaves


text:
"""
    const FactorTable = Dict{NamedTuple,Float64}

Type alias of a Dictionary mapping tuples of Variable name=>value assignments to probabilities.
"""
const FactorTable = Dict{NamedTuple,Float64}

"""
    struct Variable

A Type representing a single discrete variable.

Params:

    name::Symbol
name of Variable

    m::Int:
number of possible (discrete) values this Variable can take
"""
struct Variable
    name::Symbol
    m::Int # number of possible values
end

"""
    struct Factor

A Type representing a joint or conditional probability distribution over the variables in vars.

Fields:

    vars:
The variables represented by this distribution

    table:
Probability table mapping Variable assignments (NamedTuple) => probability (Float64)
"""
struct Factor
    vars::Vector{Variable}
    table::FactorTable
end

"""
    function variablenames(φ::Factor)

Get a vector of variable names (Symbols) for the Variables in a Factor
"""
variablenames(φ::Factor) = [var.name for var in φ.vars]

"""
    function assignments(vars::AbstractVector{Variable})

Get an array of all possible variable=>value assignment tuples for Variables in vars.
"""
function assignments(vars::AbstractVector{Variable})

    # get a vector of variable names
    n = [var.name for var in vars]

    # product(1:v1.m, 1:v2.m, 1:v3.m) -> cartesian product of integer indicies
    #                                 -> iterator over tuples of indicies
    return [namedtuple(n, v) for v in product((1:v.m for v in vars)...)]
end

"""
    function normalize!(φ::Factor)

Update φ::Factor so all values sum to 1.0 while maintaining relative value (so values represent probabilities).
"""
function normalize!(ϕ::Factor)

    # sum all probabilities in Factor table
    z = sum(p for (a,p) in ϕ.table)

    # divide each value by total sum
    for (a,p) in ϕ.table
            ϕ.table[a] = p/z
    end

    return ϕ
end

"""
    struct BayesianNetwork

Represents a Bayesian Network.

Fields:

    vars::Vector{Variable}:
Vector of Variables contained in network

    factors::Vector{Factor}:
Vector of Factors, one for each Variable in vars, representing conditional probability distribution over that variable given parents of var in network.

    graph::SimpleDiGraph{Int64}:
Represents the structure of the network variables and implicitly defines independence relationships.

"""
struct BayesianNetwork

    """Vector of every Variable in the network"""
    vars::Vector{Variable}

    """
        Network conditional independence specified by a `Vector` of `Factor`s.
        Each Factor speifies the conditional joint PMF over its subset of variables, given the
    """
    factors::Vector{Factor}

    """
        DAG that specifies the independence relationships among variables.
        Nodes
    """
    graph::SimpleDiGraph{Int64}
end

"""
    function probability(bn::BayesianNetwork, assignment)

Calculate the probability of a given variable assignment given the bayesian network bn.

# Arguments

* bn: The BayesianNetwork describing the joint distribution
* assignment: a NamedTuple of :Symbol=>Integer pairs giving the
            categorical index value of each variable.

# Examples
```julia
julia> X = Variable(:x, 2);

julia> Y = Variable(:y, 2);

julia> Z = Variable(:z, 2);

julia> ϕx = Factor([X, Y, Z], FactorTable(
           (x=1, y=1, z=1) => 0.08,
           (x=1, y=1, z=2) => 0.31,
           (x=1, y=2, z=1) => 0.09,
           (x=1, y=2, z=2) => 0.37,
           (x=2, y=1, z=1) => 0.01,
           (x=2, y=1, z=2) => 0.05,
           (x=2, y=2, z=1) => 0.02,
           (x=2, y=2, z=2) => 0.07,
       ));

julia> ϕy = Factor([Y], FactorTable(
           (y=1,) => 0.20,
           (y=2,) => 0.80,
       ));

julia> ϕz = Factor([Z], FactorTable(
           (z=1,) => 0.40,
           (z=2,) => 0.60,
       ));

julia> g = SimpleDiGraph(3);

julia> add_edge!(g,1,3);

julia> add_edge!(g,2,3);

julia> bn = BayesianNetwork([Y, Z, X], [ϕy, ϕz, ϕ], g);

julia> assignment = (x=1, y=2, z=1);

julia> probability(bn, assignment)
0.028800000000000006
```
"""
function probability(bn::BayesianNetwork, assignment)

        # Function to map Factor to a NamedTuple containing only the name=>value pairs for names specified in the (conditionally independent) Factor φ and values for them (and all other variables in bn) specified in `assignment`
        __select(φ) = select(assignment, variablenames(φ))

        # Function mapping Factor to probability of factor assignment
        probability(φ) = φ.table[__select(φ)]

        # Final probability is product of individual factor's probabilities (since they are implicitly independent)
        return prod(probability(φ) for φ in bn.factors)

end


text:
module mkcell
using MkCell

function parse_args(args)
    # 11 elements, whitespace separated
    # 9 for 3x3 matrix, [a b c; d e f; h i j]
    mat = reshape(parse.(Float64, args[1:9]), (3,3))'
    # scaling size
    scale = parse(Int, args[10])
    # search limits (can set automatically)
    lims = parse(Int, args[11])

    return mat, scale, lims

end

Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
    mat, scale, lims = parse_args(ARGS)
    @show MkCell.cellopt(mat, scale, lims)
    return 0
end

end


text:
using Distributions

"""
Converts a year value into an integer corresponding to fund's time index.
"""
function getindexfromyear(year::Int)
    baseyear = 1950
    return year - baseyear + 1
end

"""
Reads parameter csvs from data directory into a dictionary (parameter_name => default_value).
For parameters defined as distributions, this sets the value to their mode.
"""
function load_default_parameters(datadir = joinpath(dirname(@__FILE__), "..", "data"))
    files = readdir(datadir)
    filter!(i -> i != "desktop.ini", files)
    parameters = Dict{Any, Any}(splitext(file)[1] => readdlm(joinpath(datadir,file), ',', comments = true) for file in files)

    prepparameters!(parameters)

    return parameters
end

# For Truncated Gamma distributions, fund uses the mode of the untrucated distribution as it's default value.
import StatsBase.mode
function mode(d::Truncated{Gamma{Float64},Continuous})
    return mode(d.untruncated)
end

"""
Returns the mode for a distributional parameter; returns the value if it's not a distribution.
"""
getbestguess(p) = isa(p, ContinuousUnivariateDistribution) ? mode(p) : p

"""
Converts the original parameter dictionary loaded from the data files into a dictionary of default parameter values.

Original dictionary: parameter_name => string of distributions or values from csv file
Final dictionary: parameter_name => default value
"""
function prepparameters!(parameters)
    for (param_name, p) in parameters
        column_count = size(p,2)
        if column_count == 1
            parameters[param_name] = getbestguess(convertparametervalue(p[1,1]))
        elseif column_count == 2
            parameters[param_name] = Float64[getbestguess(convertparametervalue(p[j,2])) for j in 1:size(p,1)]
        elseif column_count == 3
            length_index1 = length(unique(p[:,1]))
            length_index2 = length(unique(p[:,2]))
            new_p = Array{Float64}(undef, length_index1, length_index2)
            cur_1 = 1
            cur_2 = 1
            for j in 1:size(p,1)
                new_p[cur_1,cur_2] = getbestguess(convertparametervalue(p[j,3]))
                cur_2 += 1
                if cur_2 > length_index2
                    cur_2 = 1
                    cur_1 += 1
                end
            end
            parameters[param_name] = new_p
        end
    end
end

"""
Takes as input a single parameter value.
If the parameter value is a string containing a distribution definition, it returns the distribtion.
If the parameter value is a number, it returns the number.
"""
function convertparametervalue(pv)
    if isa(pv,AbstractString)
        if startswith(pv,"~") && endswith(pv,")")
            args_start_index = something(findfirst(isequal('('), pv), 0)
            dist_name = pv[2:args_start_index-1]
            args = split(pv[args_start_index+1:end-1], ';')
            fixedargs = filter(i->!occursin("=", i),args)
            optargs = Dict(split(i,'=')[1]=>split(i,'=')[2] for i in filter(i->occursin("=", i),args))

            if dist_name == "N"
                if length(fixedargs)!=2 error() end
                if length(optargs)>2 error() end

                basenormal = Normal(parse(Float64, fixedargs[1]),parse(Float64, fixedargs[2]))

                if length(optargs)==0
                    return basenormal
                else
                    return Truncated(basenormal,
                        haskey(optargs,"min") ? parse(Float64, optargs["min"]) : -Inf,
                        haskey(optargs,"max") ? parse(Float64, optargs["max"]) : Inf)
                end
            elseif startswith(pv, "~Gamma(")
                if length(fixedargs)!=2 error() end
                if length(optargs)>2 error() end

                basegamma = Gamma(parse(Float64, fixedargs[1]),parse(Float64, fixedargs[2]))

                if length(optargs)==0
                    return basegamma
                else
                    return Truncated(basegamma,
                        haskey(optargs,"min") ? parse(Float64, optargs["min"]) : -Inf,
                        haskey(optargs,"max") ? parse(Float64, optargs["max"]) : Inf)
                end
            elseif startswith(pv, "~Triangular(")
                triang = TriangularDist(parse(Float64, fixedargs[1]), parse(Float64, fixedargs[2]), parse(Float64, fixedargs[3]))
                return triang
            else
                error("Unknown distribution")
            end
        elseif pv=="true"
            return true
        elseif pv=="false"
            return false
        elseif endswith(pv, "y")
            return parse(Int, strip(pv,'y'))
        else
            try
                return parse(Float64, pv)
            catch e
                error(pv)
            end
        end
        return pv
    else
        return pv
    end
end


text:
######### Stan program example  ###########

using StanSample, MCMCChains

ProjDir = @__DIR__

bernoullimodel = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
"

datatheta = joinpath(ProjDir, "bernoulli.data.json")
inittheta = joinpath(ProjDir , "bernoulli.init.json")

sm = SampleModel("bernoulli", bernoullimodel);

rc = stan_sample(sm, data=datatheta, init=inittheta)
#rc = stan_sample(sm, data=datatheta)

if success(rc)
  chn = read_samples(sm, :mcmcchains)
  chn |> display
end


text:
module Platforms

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict()

for texture in Maple.wood_platform_textures
    placements["Platform (Moving, $(uppercasefirst(texture)))"] = Ahorn.EntityPlacement(
        Maple.MovingPlatform,
        "rectangle",
        Dict{String, Any}(
            "texture" => texture
        ),
        function(entity)
            x, y = Int(entity.data["x"]), Int(entity.data["y"])
            width = Int(get(entity.data, "width", 8))
            entity.data["x"], entity.data["y"] = x + width, y
            entity.data["nodes"] = [(x, y)]
        end
    )

    placements["Platform (Sinking, $(uppercasefirst(texture)))"] = Ahorn.EntityPlacement(
        Maple.SinkingPlatform,
        "rectangle",
        Dict{String, Any}(
            "texture" => texture
        )
    )
end

Ahorn.editingOptions(entity::Maple.MovingPlatform) = Dict{String, Any}(
    "texture" => Maple.wood_platform_textures
)

Ahorn.editingOptions(entity::Maple.SinkingPlatform) = Dict{String, Any}(
    "texture" => Maple.wood_platform_textures
)

Ahorn.nodeLimits(entity::Maple.MovingPlatform) = 1, 1

Ahorn.resizable(entity::Maple.MovingPlatform) = true, false
Ahorn.resizable(entity::Maple.SinkingPlatform) = true, false

Ahorn.minimumSize(entity::Maple.MovingPlatform) = 8, 0
Ahorn.minimumSize(entity::Maple.SinkingPlatform) = 8, 0

function Ahorn.selection(entity::Maple.SinkingPlatform)
    x, y = Ahorn.position(entity)
    width = Int(get(entity.data, "width", 8))

    return Ahorn.Rectangle(x, y, width, 8)
end

function Ahorn.selection(entity::Maple.MovingPlatform)
    width = Int(get(entity.data, "width", 8))
    startX, startY = Int(entity.data["x"]), Int(entity.data["y"])
    stopX, stopY = Int.(entity.data["nodes"][1])

    return [Ahorn.Rectangle(startX, startY, width, 8), Ahorn.Rectangle(stopX, stopY, width, 8)]
end

outerColor = (30, 14, 25) ./ 255
innerColor = (10, 0, 6) ./ 255

function renderConnection(ctx::Ahorn.Cairo.CairoContext, x::Number, y::Number, nx::Number, ny::Number, width::Number)
    cx, cy = x + floor(Int, width / 2), y + 4
    cnx, cny = nx + floor(Int, width / 2), ny + 4

    length = sqrt((x - nx)^2 + (y - ny)^2)
    theta = atan(cny - cy, cnx - cx)

    Ahorn.Cairo.save(ctx)

    Ahorn.translate(ctx, cx, cy)
    Ahorn.rotate(ctx, theta)

    Ahorn.setSourceColor(ctx, outerColor)
    Ahorn.set_antialias(ctx, 1)
    Ahorn.set_line_width(ctx, 3);

    Ahorn.move_to(ctx, 0, 0)
    Ahorn.line_to(ctx, length, 0)

    Ahorn.stroke(ctx)

    Ahorn.setSourceColor(ctx, innerColor)
    Ahorn.set_antialias(ctx, 1)
    Ahorn.set_line_width(ctx, 1);

    Ahorn.move_to(ctx, 0, 0)
    Ahorn.line_to(ctx, length, 0)

    Ahorn.stroke(ctx)

    Ahorn.Cairo.restore(ctx)
end

function renderPlatform(ctx::Ahorn.Cairo.CairoContext, texture::String, x::Number, y::Number, width::Number)
    tilesWidth = div(width, 8)

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, "objects/woodPlatform/$texture", x + 8 * (i - 1), y, 8, 0, 8, 8)
    end

    Ahorn.drawImage(ctx, "objects/woodPlatform/$texture", x, y, 0, 0, 8, 8)
    Ahorn.drawImage(ctx, "objects/woodPlatform/$texture", x + tilesWidth * 8 - 8, y, 24, 0, 8, 8)
    Ahorn.drawImage(ctx, "objects/woodPlatform/$texture", x + floor(Int, width / 2) - 4, y, 16, 0, 8, 8)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.MovingPlatform, room::Maple.Room)
    width = Int(get(entity.data, "width", 8))

    x, y = Int(entity.data["x"]), Int(entity.data["y"])
    nx, ny = Int.(entity.data["nodes"][1])

    texture = get(entity.data, "texture", "default")

    renderConnection(ctx, x, y, nx, ny, width)
    renderPlatform(ctx, texture, x, y, width)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.SinkingPlatform, room::Maple.Room)
    width = Int(get(entity.data, "width", 8))

    x, y = Int(entity.data["x"]), Int(entity.data["y"])

    texture = get(entity.data, "texture", "default")

    renderConnection(ctx, x, y, x, Int(Ahorn.loadedState.room.size[2]), width)
    renderPlatform(ctx, texture, x, y, width)
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.MovingPlatform, room::Maple.Room)
    width = Int(get(entity.data, "width", 8))

    startX, startY = Int(entity.data["x"]), Int(entity.data["y"])
    stopX, stopY = Int.(entity.data["nodes"][1])

    texture = get(entity.data, "texture", "default")

    renderPlatform(ctx, texture, startX, startY, width)
    renderPlatform(ctx, texture, stopX, stopY, width)

    Ahorn.drawArrow(ctx, startX + width / 2, startY, stopX + width / 2, stopY, Ahorn.colors.selection_selected_fc, headLength=6)
end

end


text:
using QSM: Poisson2
using QSM: wrapped_laplacian
using QSM: wrapped_laplacian_boundary_neumann!
using QSM: wrapped_laplacian_boundary_periodic!
using QSM: solve_poisson_dct, solve_poisson_fft, solve_poisson_mgpcg

L2norm(x) = norm(x) / sqrt(length(x))

function test_problem(::Type{T}, sz::NTuple{3, Integer}) where {T<:AbstractFloat}
    f = sinpi

    R = CartesianIndices(map(n -> 2:n-1, sz))
    m = zeros(Bool, sz)
    u = zeros(T, sz)
    d2u = zeros(T, sz)
    d2uh = zeros(T, sz)

    X = range(0, 1, length=sz[1])
    Y = range(0, 1, length=sz[2])
    Z = range(0, 1, length=sz[3])

    dx = (X.step.hi, Y.step.hi, Z.step.hi)
    _u = [f(x).*f(y).*f(z) for x in X, y in Y, z in Z]

    m[R] .= true
    u[R] .= _u[R]
    d2u[R] .= -3*π^2*_u[R]

    A = Poisson2(T, m, dx)
    d2uh = mul!(d2uh, A, u)

    return A, u, d2u, d2uh
end

@testset "Multigrid" begin
    include("multigrid/multigrid.jl")
end

@testset "DCT with wrapped Laplacian" for T in (Float32, Float64)
    Random.seed!(501)
    for sz in [(11, 12, 13), (12, 13, 14, 3)]
        x, dx = rand(T, sz), Tuple(rand(T, 3))
        x .-= mean(x, dims=(1,2,3))
        d2x = wrapped_laplacian(x, dx)
        d2x = wrapped_laplacian_boundary_neumann!(d2x, x, dx)
        @test solve_poisson_dct(d2x, dx) ≈ x
    end
end

@testset "FFT with wrapped Laplacian" for T in (Float32, Float64)
    Random.seed!(501)
    for sz in [(11, 12, 13), (12, 13, 14, 3)]
        x, dx = rand(T, sz), Tuple(rand(T, 3))
        x .-= mean(x, dims=(1,2,3))
        d2x = wrapped_laplacian(x, dx)
        d2x = wrapped_laplacian_boundary_periodic!(d2x, x, dx)
        @test solve_poisson_fft(d2x, dx) ≈ x
    end

end

@testset "MGPCG with wrapped_laplacian" for T in (Float32, Float64)
    Random.seed!(501)
    for sz in [(127, 128, 129), (130, 131, 132)]
        A, u, d2u, d2uh = test_problem(T, sz)
        @test wrapped_laplacian(u, A.dx) ≈ d2uh
        @test solve_poisson_mgpcg(d2uh, A.interior, A.dx; rtol=eps(T)) ≈ u
    end
end


text: stringtonumber(str) = parse( Int, str )

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
#=
  Binary IO
=#

using Base.Threads: @threads, nthreads, threadid
using LinearAlgebra.BLAS: set_num_threads
using Base.Filesystem: mkdir, readdir

# Convert a matrix to a list of matrices
function matrixToBlock(mat::Array{T, 2}, nBlocks::Int64) where {T <: AbstractFloat}
  n = size(mat)[1]
  ret = Array{Array{T, 2}, 1}(undef, nBlocks)
  for i in 0:(nBlocks - 1)
    start = div((n * i), nBlocks) + 1
    finish = div((n*(i + 1)), nBlocks)
    ret[i + 1] = mat[start:finish, :]
  end

  return ret
end

"""
  Function to write 2D array to disk
"""
function write2DArray(arr::Array{T}, fileName::String) where {T <: AbstractFloat}

  # Opens Binary File
  io = open(fileName, "w")
  _size = size(arr)

  for i in _size
    write(io, Int64(i))
  end

  for el in arr
    write(io, el)
  end

  close(io)
  return 0

end

function writeBlockMatrix(blockMatrix::Array{Array{T, 2}, 1}, path::String) where {T <: AbstractFloat}
  mkdir(path)
  nchar::Int64 = length(path)
  # Unix dependent
  if path[nchar] != '/'
    path *= "/"
  end
  nBlocks::Int64 = length(blockMatrix)
  for i in 1:nBlocks
    write2DArray(blockMatrix[i], path * "block_" * string(i) * ".bin")
  end
  return
end

"""
  Function to read 2D array from disk
"""
function read2DArray(::Type{T}, fileName::String) where {T <: AbstractFloat}

  io = open(fileName, "r")

  _size = zeros(Int64, 2)

  for i in 1:2
    _size[i] = read(io, Int64)
  end

  arr = Array{T}(undef, tuple(_size...))

  for i in eachindex(arr)
    arr[i] = read(io, T)
  end
  close(io)
  return arr
end

function readBlockMatrix(::Type{T}, path::String) where {T <: AbstractFloat}
  nchar::Int64 = length(path)
  # Unix dependent
  if path[nchar] != '/'
    path *= "/"
  end
  files::Array{String, 1} = readdir(path)
  nBlocks::Int64 = length(files)
  blockMatrix = Array{Array{T, 2}, 1}(undef, nBlocks)
  for i in 1:nBlocks
    blockMatrix[i] = read2DArray(T, path * files[i])
  end
  return blockMatrix
end

"""
  Function to write an array to disk
"""
function writeNDArray(arr::Array{T}, fileName::String) where {T <: AbstractFloat}

  # Opens Binary File
  io = open(fileName, "w")
  _size = size(arr)

  write(io, Int64(length(_size)))

  for i in _size
    write(io, Int64(i))
  end

  for el in arr
    write(io, el)
  end

  close(io)
  return 0

end

"""
  Function to read an array from disk
"""
function readNDArray(::Type{T}, fileName::String) where {T <: AbstractFloat}
  io = open(fileName, "r")
  _len = read(io, Int64)
  _size = tuple(zeros(Int64, _len)...)

  for i in 1:_len
    _size[i] = read(io, Int64)
  end

  arr = Array{T}(undef, _size)

  for i in eachindex(arr)
    arr[i] = read(io, T)
  end

  return arr
end

max_stars_repo_name: dataPulverizer/glmSolver

max_issues_repo_name: dataPulverizer/glmSolver

max_forks_repo_name: dataPulverizer/glmSolver


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

mean(d, v::Symbol, w::AbstractWeighter) = mean(d[v], weight(d, w))
mean(d, v::Symbol, s::Number) = mean(d, v, BlockWeighter(ntuple(i->s,ncoords(d))))
mean(d, v::Symbol) = mean(d, v, median_heuristic(d))
mean(d, w::AbstractWeighter) = Dict(v => mean(d, v, w) for v in name.(variables(d)))
mean(d, s::Number) = mean(d, BlockWeighter(ntuple(i->s,ncoords(d))))

"""
    mean(sdata)
    mean(sdata, v)
    mean(sdata, v, s)

Spatial mean of spatial data `sdata`. Optionally,
specify the variable `v` and the block side `s`.
"""
mean(d::AbstractData) = mean(d, median_heuristic(d))

var(d, v::Symbol, w::AbstractWeighter) = var(d[v], weight(d, w), mean=mean(d, v, w), corrected=false)
var(d, v::Symbol, s::Number) = var(d, v, BlockWeighter(ntuple(i->s,ncoords(d))))
var(d, v::Symbol) = var(d, v, median_heuristic(d))
var(d, w::AbstractWeighter) = Dict(v => var(d, v, w) for v in name.(variables(d)))
var(d, s::Number) = var(d, BlockWeighter(ntuple(i->s,ncoords(d))))

"""
    var(sdata)
    var(sdata, v)
    var(sdata, v, s)

Spatial variance of spatial data `sdata`. Optionally,
specify the variable `v` and the block side `s`.
"""
var(d::AbstractData) = var(d, median_heuristic(d))

quantile(d, v::Symbol, p, w::AbstractWeighter) = quantile(d[v], weight(d, w), p)
quantile(d, v::Symbol, p, s::Number) = quantile(d, v, p, BlockWeighter(ntuple(i->s,ncoords(d))))
quantile(d, v::Symbol, p) = quantile(d, v, p, median_heuristic(d))
quantile(d, p, w::AbstractWeighter) = Dict(v => quantile(d, v, p, w) for v in name.(variables(d)))
quantile(d, p::T, s::Number) where {T<:Union{Number,AbstractVector}} = quantile(d, p, BlockWeighter(ntuple(i->s,ncoords(d))))

"""
    quantile(sdata, p)
    quantile(sdata, v, p)
    quantile(sdata, v, p, s)

Spatial quantile of spatial data `sdata` at probability `p`.
Optionally, specify the variable `v` and the block side `s`.
"""
quantile(d::AbstractData, p) = quantile(d, p, median_heuristic(d))

EmpiricalHistogram(d, v::Symbol, w::AbstractWeighter) = fit(Histogram, d[v], weight(d, w))
EmpiricalHistogram(d, v::Symbol, s::Number) = EmpiricalHistogram(d, v, BlockWeighter(ntuple(i->s,ncoords(d))))
EmpiricalHistogram(d, v::Symbol) = EmpiricalHistogram(d, v, median_heuristic(d))
EmpiricalHistogram(d, w::AbstractWeighter) = Dict(v => EmpiricalHistogram(d, v, w) for v in name.(variables(d)))
EmpiricalHistogram(d, s::Number) = EmpiricalHistogram(d, BlockWeighter(ntuple(i->s,ncoords(d))))

"""
    EmpiricalHistogram(sdata)
    EmpiricalHistogram(sdata, v)
    EmpiricalHistogram(sdata, v, s)

Spatial histogram of spatial data `sdata`. Optionally,
specify the variable `v` and the block side `s`.
"""
EmpiricalHistogram(d) = EmpiricalHistogram(d, median_heuristic(d))

function median_heuristic(d)
  # select at most 1000 points at random
  nel = nelms(d)
  inds = sample(1:nel, min(nel, 1000), replace=false)
  X = coordinates(d, inds)
  D = pairwise(Euclidean(), X, dims=2)

  # median heuristic
  n = size(D, 1)
  m = median(D[i,j] for i in 1:n for j in 1:n if i > j)

  # bounding box constraint
  l = minimum(sides(boundbox(d)))

  min(m, l)
end


text:
# test the parts in optTopLevel
# Claudia February 2015

# ------------------5taxon network 1 hybridization-----------------
# starting topology: Case G
include("../examples/case_g_example.jl");
currT = deepcopy(net);
printEdges(currT)

# real network: Case H
df = readtable("CaseH_output.csv")
d = readDataCF(df)
#df2 = readtable("CaseG_output.csv")

optBL!(currT,d)
newT = deepcopy(currT);
count = 0
N = 100
move = whichMove(currT)
move = :CHdir
move = :MVorigin
move = :MVtarget
move = :nni

flag = proposedTop!(move,newT,true,count,N, zeros(Int,18), zeros(Int,6))
printEdges(newT)
printNodes(newT)
sum([e.hybrid?1:0 for e in newT.edge]) == 2 || error("there are not 2 hybrid edges")
newT.hybrid[1].k

optBL!(newT,d)
newloglik - currloglik
changeDirectionUpdate!(newT.node[5],newT,false)

currT = deepcopy(newT);
currloglik = newloglik
currxmin = newxmin

# ------------------5taxon network 1 hybridization-----------------
# starting topology: Case F
include("../case_f_example.jl");
currT = deepcopy(net);
printEdges(currT)

# real network: Case H
df = readtable("CaseH_output.csv")
d = readDataCF(df)

currloglik,currxmin = optBL!(currT,d)
updateParameters!(currT)
updateLik!(currT,currloglik)
newT = deepcopy(currT);
count = 0
N = 100
move = whichMove(currT)
move = :CHdir
move = :MVorigin
move = :MVtarget
move = :nni

flag = proposedTop!(move,newT,true,count,N)
printEdges(newT)
printNodes(newT)
sum([e.hybrid?1:0 for e in newT.edge]) == 2 || error("there are not 2 hybrid edges")
newT.hybrid[1].k

newloglik, newxmin = optBL!(newT,d)
newloglik - currloglik

currT = deepcopy(newT);
currloglik = newloglik
currxmin = newxmin

# -------------------5taxon tree------------------

include("../src/types.jl")
include("../src/functions.jl")

df = readtable("Tree_output.csv")
d = readDataCF(df)

# starting tree:
tree = "((6,4),(7,8),10);"
f = open("prueba_tree.txt","w")
write(f,tree)
close(f)
currT = readTopologyUpdate("prueba_tree.txt");
printEdges(currT)

optBL!(currT,d)
newT = deepcopy(currT);
count = 0
N = 10
move = :nni

flag = proposedTop!(move,newT,true,count,N, zeros(Int,18), zeros(Int,6))
flag
printEdges(newT)
printNodes(newT)

newloglik, newxmin = optBL!(newT,d)
newloglik - currloglik
`
currT = deepcopy(newT);
currloglik = newloglik
currxmin = newxmin


text:
"""
    function test_2021_03_26()
---
To test:
- if the right loci is found
- if the editing is right
"""
function test_2021_03_26()
    @load "dat/run/base.jld2" base
    qtl = sim_QTL(base, 10)[1]
    println(qtl.pos)
    snp = copy(base[:hap])
    ep = gedit(base[:hap], qtl, .8)
    println(sum(snp[ep, :]), ' ', sum(base[:hap][ep, :]))
end

"""
    function edit_2021_03_26()
---
1. simulate a base with `macs`, of 29 chromosomes
2. simulate a QTL setup
3. compare
   - random mating
   - selection with GS
   - selection with GEAS
"""
function edit_2021_03_26()
    Parameters = Dict(:nSire=> 125,
                      :nDam => 250, # male:female = 1:2
                      :nSib => 100,
                      :nC7e => 40, # number of sib to be challenged
                      :nG8n => 5,  # number of generations
                      :nQTL => [1000, 1000],
                      :h²   => [.5, .5],
                      :p8e  => .5, # percentage of affected in the binary trait
                      :e19e => .8, # edit_successsful_rate
                      :w4t  => 2   # weight on the binary trait EBV
                      )
    par = (; Parameters...)     # named tuple.  contents as above

    base = sim_base(1200, 29, 1e8, 1725)
    qtl = sim_QTL(base, par.nQTL...)

    p1, _ = breeding_program(base, par, qtl, edit=false)
    p2, _ = breeding_program(base, par, qtl, edit=true)
end


text:
using ModelingToolkit, Test
@parameters α β δ
expr = (((1 / β - 1) + δ) / α) ^ (1 / (α - 1))
ref = [β, δ, α]
sol = Num.(ModelingToolkit.get_variables(expr))
@test all(simplify, sol[i] == ref[i] for i in 1:3)

@parameters γ
s = α => γ
expr = (((1 / β - 1) + δ) / α) ^ (1 / (α - 1))
sol = ModelingToolkit.substitute(expr, s)
new = (((1 / β - 1) + δ) / γ) ^ (1 / (γ - 1))
@test iszero(sol - new)

# test namespace_expr
@parameters t a p(t)
pterm = p.val
pnsp = ModelingToolkit.namespace_expr(pterm, :namespace, :t)
@test typeof(pterm) == typeof(pnsp)
@test ModelingToolkit.getname(pnsp) == Symbol("namespace₊p")
asym = a.val
ansp = ModelingToolkit.namespace_expr(asym, :namespace, :t)
@test typeof(asym) == typeof(ansp)
@test ModelingToolkit.getname(ansp) == Symbol("namespace₊a")


text:
using JuliaInterp
using BenchmarkTools

const SUITE = BenchmarkGroup()

options = JuliaInterp.options(false, false, [@__MODULE__], [], UInt(0))

# Recursively call itself
f(i, j) = i == 0 ? j : f(i - 1, j + 1)
SUITE["recursive self 1_000"] = @benchmarkable JuliaInterp.interpret_ast(@__MODULE__, Meta.parse("f(1_000, 0)"), options)

# Long stack trace calling other functions
f0(i) = i
for i in 1:1_000
    @eval $(Symbol("f", i))(i) = $(Symbol("f", i-1))(i)
end
SUITE["recursive other 1_000"] = @benchmarkable Base.@invokelatest JuliaInterp.interpret_ast(@__MODULE__, Meta.parse("f1000(1)"), options)

# Tight loop
function f(X)
    s = 0
    for x in X
        s += x
    end
    return s
end
const X = rand(1:10, 10_000)
SUITE["tight loop 10_000"] = @benchmarkable JuliaInterp.interpret_ast(@__MODULE__, Meta.parse("f(X)"), options)

# Throwing and catching an error over a large stacktrace
function g0(i)
    try
        g1(i)
    catch e
        e
    end
end
for i in 1:1_000
    @eval $(Symbol("g", i))(i) = $(Symbol("g", i+1))(i)
end
g1001(i) = error()
SUITE["throw long 1_000"] = @benchmarkable JuliaInterp.interpret_ast(@__MODULE__, Meta.parse("g0(1)"), options)

# Function with many statements
macro do_thing(expr, N)
    e = Expr(:block)
    for i in 1:N
        push!(e.args, esc(expr))
    end
    return e
end

function counter()
    a = 0
    @do_thing(a = a + 1, 5_000)
    return a
end
SUITE["long function 5_000"] = @benchmarkable JuliaInterp.interpret_ast(@__MODULE__, Meta.parse("counter()"), options)

# Ccall
function ccall_ptr(ptr, x, y)
    ccall(ptr, Int, (Int, Int), x, y)
end
const ptr = @cfunction(+, Int, (Int, Int))
SUITE["ccall ptr"] = @benchmarkable JuliaInterp.interpret_ast(@__MODULE__, Meta.parse("ccall_ptr(ptr, 1, 5)"), options)

function powf(a, b)
    ccall(("powf", Base.Math.libm), Float32, (Float32,Float32), a, b)
end
SUITE["ccall library"] = @benchmarkable JuliaInterp.interpret_ast(@__MODULE__, Meta.parse("powf(2, 3)"), options)


text:
export iauGd2gc
"""
Transform geodetic coordinates to geocentric using the specified
reference ellipsoid.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  canonical transformation.

Given:
   n       int        ellipsoid identifier (Note 1)
   elong   double     longitude (radians, east +ve)
   phi     double     latitude (geodetic, radians, Note 3)
   height  double     height above ellipsoid (geodetic, Notes 2,3)

Returned:
   xyz     double[3]  geocentric vector (Note 2)

Returned (function value):
            int        status:  0 = OK
                              -1 = illegal identifier (Note 3)
                              -2 = illegal case (Note 3)

Notes:

   1. The identifier n is a number that specifies the choice of
   reference ellipsoid.  The following are supported:

      n    ellipsoid

      1     WGS84
      2     GRS80
      3     WGS72

   The n value has no significance outside the SOFA software.  For
   convenience, symbols WGS84 etc. are defined in sofam.h.

   2. The height (height, given) and the geocentric vector (xyz,
   returned) are in meters.

   3. No validation is performed on the arguments elong, phi and
   height.  An error status -1 means that the identifier n is
   illegal.  An error status -2 protects against cases that would
   lead to arithmetic exceptions.  In all error cases, xyz is set
   to zeros.

   4. The inverse transformation is performed in the function iauGc2gd.

Called:
   iauEform     Earth reference ellipsoids
   iauGd2gce    geodetic to geocentric transformation, general
   iauZp        zero p-vector

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# int iauGd2gc ( int n, double elong, double phi, double height,
#                double xyz[3] )

function iauGd2gc(n::Int, elong::Real, phi::Real, height::Real)
   xyz = zeros(Float64, 3)

   status = ccall((:iauGd2gc, libsofa_c), Cint,
         (Cint, Cdouble, Cdouble, Cdouble,
         Ptr{Cdouble}),
         convert(Int32, n),
         convert(Float64, elong),
         convert(Float64, phi),
         convert(Float64, height),
         xyz)

   return status, SVector{3}(xyz)
 end


text:
using SIIPExamples
pkgpath = dirname(dirname(pathof(SIIPExamples)))
include(joinpath(pkgpath, "test", "2_PowerSystems_examples", "parse_matpower.jl"))

folder = mktempdir()
path = joinpath(folder, "system.json")
@info "Serializing to $path"
to_json(sys, path)

filesize(path) / 1000000 #MB

sys2 = System(path)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
arguement = int(ARGS[1])
fib(n) = n < 2 ? n : fib(n - 1) + fib(n - 2)
result = fib(arguement)

println(result)

max_stars_repo_name: csenkbeil/RecursiveFibonacci

max_issues_repo_name: csenkbeil/RecursiveFibonacci

max_forks_repo_name: csenkbeil/RecursiveFibonacci


text:
using DataFrames, CSV, Statistics, StatsBase

df = CSV.read("./preprocessing/input/land-change-impacts.csv")

# remove last (empty) column
df = df[.!ismissing.(df.Scenario),:]

select!(df, Not(:Model))
select!(df, Not(:Region))

df = stack(df, 4:size(df)[2])

rename!(df, [:year, :value, :scenario, :variable, :unit])

println(first(df, 8))

df[df.variable .== "Primary Energy|Biomass",:variable] .= "Bioenergy"
df[df.variable .== "Carbon Sequestration|CCS|Biomass",:variable] .= "Carbon Sequestration from Biomass CCS"
df[df.variable .== "Price|Agriculture|Non-Energy Crops and Livestock|Index",:variable] .= "Food Price"

df = unstack(df, :year, :value)

CSV.write("./src/assets/data/land-change-impact.csv", df)

max_stars_repo_name: SensesProject/land-use

max_issues_repo_name: SensesProject/land-use

max_forks_repo_name: SensesProject/land-use


text:
"""
    SumAudioGenerator

An audio generator used to combine (sum) multiple other sources.
"""
mutable struct SumAudioGenerator <: AbstractAudioGenerator
    generators::Vector{AbstractAudioGenerator}
end

function reset!(gen::SumAudioGenerator)
    for g in gen.generators
        reset!(g)
    end
    nothing
end
function has_samples(gen::SumAudioGenerator, sample_rate::Float64)
    any(
        has_samples(g, sample_rate)
        for g in gen.generators
    )
end
function next_samples!(gen::SumAudioGenerator, max_n::Int,
                       sample_rate::Float64)
    samples = zeros(Float64, max_n, 2)
    max_len = 0
    for g in gen.generators
        s = next_samples!(g, max_n, sample_rate)
        len = size(s)[1]
        samples[1:len, :] .+= s
        max_len = max(max_len, len)
    end
    samples[1:max_len]
end


text:
type AnalyticEuropeanEngine{P <: AbstractBlackScholesProcess} <: PricingEngine
  process::P
end

function _calculate!(pe::AnalyticEuropeanEngine, opt::EuropeanOption)
  #TODO check non striked payoff
  payoff = opt.payoff
  variance = black_variance(pe.process.blackVolatility, opt.exercise.dates[end], payoff.strike)

  dividendDiscount = discount(pe.process.dividendYield, opt.exercise.dates[end])

  riskFreeDiscount = discount(pe.process.riskFreeRate, opt.exercise.dates[end])

  spot = state_variable(pe.process).value
  forwardPrice = spot *  dividendDiscount / riskFreeDiscount

  black = BlackCalculator(payoff, forwardPrice, sqrt(variance), riskFreeDiscount)

  opt.results.value = value(black)
  opt.results.delta = delta(black, spot)
  opt.results.deltaForward = delta_forward(black)
  opt.results.elasticity = elasticity(black, spot)
  opt.results.gamma = gamma(black, spot)

  rfdc = pe.process.riskFreeRate.dc
  divdc = pe.process.dividendYield.dc
  voldc = pe.process.blackVolatility.dc

  t = year_fraction(rfdc, reference_date(pe.process.riskFreeRate), opt.exercise.dates[end])
  opt.results.rho = rho(black, t)

  t = year_fraction(divdc, reference_date(pe.process.dividendYield), opt.exercise.dates[end])
  opt.results.dividendRho = dividend_rho(black, t)

  t = year_fraction(voldc, reference_date(pe.process.blackVolatility), opt.exercise.dates[end])
  opt.results.vega = vega(black, t)

  opt.results.theta = theta(black, spot, t)
  opt.results.thetaPerDay = theta_per_day(black, spot, t)

  opt.results.strikeSensitivity = strike_sensitivity(black)
  opt.results.itmCashProbability = itm_cash_probability(black)

  return opt
end


text:
Z = Strang(1)
@test full(Z) == reshape([2.0],(1,1))

n = rand(1:10)
Z = Strang(n)

#Matvec product
b = randn(n)
@test_approx_eq Z*b full(Z)*b

m = rand(1:10)
A = randn(m, n)
@test_approx_eq A*Z A*full(Z)


text:
#################### Double Metropolis-Hastings Sampler ####################

#################### Types and Constructors ####################

mutable struct DMHTune <: SamplerTune
  seed_data::AbstractArray
  m::Int64

  DMHTune() = new()
  function DMHTune(seed_data::A, m::Int64) where A<:AbstractArray
	@assert m >= length(seed_data)
    new(seed_data, m)
  end
end

const DMHVariate = SamplerVariate{DMHTune}

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

function DMH(params::ElementOrVector{Symbol}, outer::Symbol,
             inner::Symbol)
  samplerfx = function(model::Model, block::Integer)
    block = SamplingBlock(model, block, true)
	targets = keys(model, :target, block)
    f = x -> logpdf!(block, x)
	fp = (x, y) -> pseudologpdf!(block, x, y)
	cl = x -> conditional_likelihood(block, x)
    v = SamplerVariate(block, f, NullFunction(); args...)

    DMH_sample!(v::DMHVariate, f, fp, cl)

    relist(block, v)
  end
  Sampler(params, samplerfx, DMHTune())
end

function DMH_sample!(v::DMHVariate, lf::Function, pslf::Function, cl::Function)
	# 1. propose theta prime
	# 2. Generate the auxillary variable using theta prime

	# this way of generating theta_prime from the current values of theta
	# takes care of the transition probability from theta_prime to theta and vice versa
	# the values equal and will cancel out.
	θ_prime = v + v.tune.scale .* rand(v.tune.proposal(0.0, 1.0), length(v))

	# calculate logpdf values
	lf_xt = lf(v.value)
	lf_xtp = lf(θ_prime)

	# store old value for possible future reference
	θ = v.value

	# prematurely assign, to allow computations to go through properly
	v[:] = θ_prime

	# Sample new pseudo observations based on the θ_prime
	y = inner_sampler(v, cl)

	# calculate logpdfs based on pseudo observations
	lf_ytp = pslf(v.value, y)
	lf_yt = pslf(θ, y)

	# calculate acceptance probability
	r = exp((lf_yt + lf_xtp) - (lf_xt + lf_ytp))

	# RWM acceptance logic is a bit reversed here. Above the new value is
	# prematurely assigned to allow computations in the inner sampler to go through.
	# If the sample is accepted nothing needs to be done anymore, otherwise the
	# old value will be reassigend.
	if rand() > r
		# sample is rejected, so use the old value.
		v[:] = θ
	end
	v
end

function inner_sampler(v::DMHTune, cond_prob::Function)
	#X::Array{Int64,1}, v_params, h_params,
	#u_params, spatial_sums::Array{Float64,2}, ling_sums::Array{Float64,2}, m::Int64)
	m = v.m
	X = v.seed_data

	feature_vals = sort!(unique(X))
	N = length(X)
	samples = deepcopy(X)
	idx_list = Vector(1:N)
	counter = 0
	while true
		random_index_order = shuffle(1:N)
		for i in random_index_order
			if ismissing(X[i]) # no missing values yet, and this part needs improvement;
				# To do: add the neighbour and global majority methods
				throw("We should not end up here")
				missing_probs = cond_prob.(Ref(samples), i, feature_vals, Ref(spatial_sums),
					Ref(ling_sums), Ref(v_params), Ref(h_params), Ref(u_params))
				missing_probs ./= sum(missing_probs)
				new_x = sample(feature_vals, StatsBase.weights(missing_probs))
				samples[i] = new_x
			end
				probs = cond_prob.(Ref(samples), i, feature_vals, Ref(spatial_sums), Ref(ling_sums),
					Ref(v_params), Ref(h_params), Ref(u_params))
				probs ./= sum(probs)
				new_x = sample(feature_vals, StatsBase.weights(probs))
				samples[i] = new_x
			counter += 1
			if counter > m
				return samples
			end
		end
	end
end


text:
module Purple

using Mixtape
import Mixtape: CompilationContext, transform, allow, optimize!
using Core.Compiler: compact!, getfield_elim_pass!, adce_pass!
using CodeInfoTools
using CodeInfoTools: resolve
using SymbolicUtils
using SymbolicUtils: Symbolic, term, Sym, similarterm, promote_symtype

#####
##### Exports
#####

export lift, deactivate, activate

#####
##### Purple
#####

# Utility functions for lifting / de-lifting to and from Sym types.
unwrap(s::Type{Sym{T, K}}) where {T, K} = T
unwrap(s::Type) = s
wrap(t::Type{T}) where T = Sym{T, Nothing}

# `activate` is a "semantic stub" function which gets filled in by
# the src of Tuple{typeof(f), syms...}.
# In short, `activate` activates the transform at that call site.
function activate(f::Function, syms::S) where S <: Tuple
    return f(syms...)
end

# `deactivate` is a fallback stub which just constructs a new Term.
function deactivate(f::Function, syms::S) where S <: Tuple
    return term(f, syms...)
end

# Compilation context to control the Mixtape pipeline.
struct LiftContext <: CompilationContext
    opt::Bool
end

# Turn on the transform for `activate`.
allow(ctx::LiftContext, mod::Module, fn::typeof(activate), inner, args...) = true

# Wrap each call with `deactivate`.
# Otherwise, if the Expr has .head == :(=), try and wrap the RHS.
function swap!(b, v, st::Expr, slotmap)
    if st.head == :call
        args = walk(l -> resolve(get(slotmap, l, l)), st).args[2 : end]
        atup = push!(b, Expr(:call, tuple, args...))
        return Expr(:call, deactivate, st.args[1], atup)
    elseif st.head == :(=) &&
        st.args[2] isa Expr &&
        st.args[2].head == :call
        callexpr = st.args[2]
        args = walk(l -> resolve(get(slotmap, l, l)), callexpr).args[2 : end]
        atup = push!(b, Expr(:call, tuple, args...))
        return Expr(:(=), st.args[1],
                    Expr(:call, deactivate, callexpr.args[1], atup))
    else
        return st
    end
end

swap!(b, v, st, slotmap) = walk(l -> get(slotmap, l, l), st)
swap!(b, v, st::Core.NewvarNode, slotmap) = st

function transform(mix::LiftContext, src, sig)

    if !(sig[2] <: Function) ||
        sig[2] === Core.IntrinsicFunction
        return src
    end # If target is not a function, just return src.

    # Create Builder and try to get the src from the inner
    # call using the unwrapped Sym arg types.
    b = CodeInfoTools.Builder(src)
    inner = sig[2].instance
    symtypes = sig[3]
    argtypes = Tuple{map(unwrap, sig[3].parameters)...}
    largtypes = length(argtypes.parameters)
    try
        inner = CodeInfoTools.code_info(inner, argtypes)
    catch
        return src
    end

    # If we can't get the inner src, just return the original src.
    # Because we are inside of `activate` --
    # this is just applying the function to the splatted syms arg.
    inner === nothing && return src

    # Slots as args.
    slotmap = Dict()
    for (ind, a) in enumerate(inner.slotnames[2 : 1 + largtypes])
        o = push!(b, Expr(:call, Base.getindex, slot(3), ind))
        setindex!(slotmap, o, get_slot(inner, a))
    end

    # Extra slots.
    for (ind, a) in enumerate(inner.slotnames[2 + largtypes : end])
        s = slot!(b, a)
    end

    # Now "inline" the src code from inner into the builder.
    # Make the right corrections to slots, etc.
    for (v, st) in enumerate(inner.code)
        ex = swap!(b, v, st, slotmap)
        setindex!(slotmap, push!(b, ex), var(v))
    end
    return CodeInfoTools.finish(b)
end

# If `opt` is turned on (= true) inside the `ctx`,
# run `ssa_inlining_pass!` with the rest of the passes.
function optimize!(ctx::LiftContext, b)
    if ctx.opt
        ir = julia_passes!(b)
    else
        ir = get_ir(b)
        ir = compact!(ir)
        ir = getfield_elim_pass!(ir)
        ir = adce_pass!(ir)
        ir = compact!(ir)
    end
    if any(==(Union{}), ir.stmts.type)
        @debug "Inference failed to infer -- found Union{}" ir
    end
    return ir
end

# Interface to emit transformed CodeInfo or a callable.
function lift(fn, tt::Type{T};
        clean = false, jit = false, opt = false) where T
    clean && return Mixtape.emit(fn, tt; opt = opt)
    symtypes = Tuple{map(tt.parameters) do p
                         Sym{p, Nothing}
                     end...}
    if !jit
        return Mixtape.emit(activate,
                            Tuple{typeof(fn), symtypes};
                            ctx = LiftContext(opt),
                            opt = opt)
    else
        entry = Mixtape.jit(activate,
                            Tuple{typeof(fn), symtypes};
                            ctx = LiftContext(opt))
        return function (symargs...)
            entry(fn, symargs)
        end
    end
end

end # module


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(csc(a+b*x), x) == :(-1*b^-1*arctanh(cos(a+b*x)))
@test integrate(csc(a+b*x)^2, x) == :(-1*b^-1*cot(a+b*x))
@test integrate(csc(a+b*x)^3, x) == :(-1//2*b^-1*arctanh(cos(a+b*x))+-1//2*b^-1*cot(a+b*x)*csc(a+b*x))
@test integrate(csc(a+b*x)^4, x) == :(-1*b^-1*cot(a+b*x)+-1//3*b^-1*cot(a+b*x)^3)
@test integrate(csc(a+b*x)^5, x) == :(-3//8*b^-1*arctanh(cos(a+b*x))+-3//8*b^-1*cot(a+b*x)*csc(a+b*x)+-1//4*b^-1*csc(a+b*x)^3*cot(a+b*x))
@test integrate(csc(a+b*x)^6, x) == :(-1*b^-1*cot(a+b*x)+-2//3*b^-1*cot(a+b*x)^3+-1//5*b^-1*cot(a+b*x)^5)
@test integrate(csc(a+b*x)^7, x) == :(-5//16*b^-1*arctanh(cos(a+b*x))+-5//16*b^-1*cot(a+b*x)*csc(a+b*x)+-5//24*b^-1*csc(a+b*x)^3*cot(a+b*x)+-1//6*b^-1*csc(a+b*x)^5*cot(a+b*x))
@test integrate(csc(a+b*x)^8, x) == :(-1*b^-1*cot(a+b*x)^3+-1*b^-1*cot(a+b*x)+-3//5*b^-1*cot(a+b*x)^5+-1//7*b^-1*cot(a+b*x)^7)
@test integrate(csc(a+b*x)^7//2, x) == :(-6//5*b^-1*csc(a+b*x)^(1/2)*cos(a+b*x)+-2//5*b^-1*csc(a+b*x)^5//2*cos(a+b*x)+-6//5*b^-1*csc(a+b*x)^(1/2)*sin(a+b*x)^(1/2)*Elliptic.E((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate(csc(a+b*x)^5//2, x) == :(-2//3*b^-1*csc(a+b*x)^3//2*cos(a+b*x)+2//3*b^-1*csc(a+b*x)^(1/2)*sin(a+b*x)^(1/2)*Elliptic.F((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate(csc(a+b*x)^3//2, x) == :(-2*b^-1*csc(a+b*x)^(1/2)*cos(a+b*x)+-2*b^-1*csc(a+b*x)^(1/2)*sin(a+b*x)^(1/2)*Elliptic.E((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate(csc(a+b*x)^(1/2), x) == :(2*b^-1*csc(a+b*x)^(1/2)*sin(a+b*x)^(1/2)*Elliptic.F((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate(csc(a+b*x)^-1//2, x) == :(2*b^-1*csc(a+b*x)^(1/2)*sin(a+b*x)^(1/2)*Elliptic.E((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate(csc(a+b*x)^-3//2, x) == :(-2//3*b^-1*csc(a+b*x)^-1//2*cos(a+b*x)+2//3*b^-1*csc(a+b*x)^(1/2)*sin(a+b*x)^(1/2)*Elliptic.F((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate(csc(a+b*x)^-5//2, x) == :(-2//5*b^-1*csc(a+b*x)^-3//2*cos(a+b*x)+6//5*b^-1*csc(a+b*x)^(1/2)*sin(a+b*x)^(1/2)*Elliptic.E((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate(csc(a+b*x)^-7//2, x) == :(-10//21*b^-1*csc(a+b*x)^-1//2*cos(a+b*x)+-2//7*b^-1*csc(a+b*x)^-5//2*cos(a+b*x)+10//21*b^-1*csc(a+b*x)^(1/2)*sin(a+b*x)^(1/2)*Elliptic.F((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate((c*csc(a+b*x))^7//2, x) == :(-6//5*b^-1*c^3*(c*csc(a+b*x))^(1/2)*cos(a+b*x)+-2//5*c*b^-1*(c*csc(a+b*x))^5//2*cos(a+b*x)+-6//5*b^-1*c^4*(c*csc(a+b*x))^-1//2*sin(a+b*x)^-1//2*Elliptic.E((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate((c*csc(a+b*x))^5//2, x) == :(-2//3*c*b^-1*(c*csc(a+b*x))^3//2*cos(a+b*x)+2//3*b^-1*c^2*(c*csc(a+b*x))^(1/2)*sin(a+b*x)^(1/2)*Elliptic.F((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate((c*csc(a+b*x))^3//2, x) == :(-2*c*b^-1*(c*csc(a+b*x))^(1/2)*cos(a+b*x)+-2*b^-1*c^2*(c*csc(a+b*x))^-1//2*sin(a+b*x)^-1//2*Elliptic.E((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate((c*csc(a+b*x))^(1/2), x) == :(2*b^-1*(c*csc(a+b*x))^(1/2)*sin(a+b*x)^(1/2)*Elliptic.F((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate((c*csc(a+b*x))^-1//2, x) == :(2*b^-1*(c*csc(a+b*x))^-1//2*sin(a+b*x)^-1//2*Elliptic.E((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate((c*csc(a+b*x))^-3//2, x) == :(-2//3*b^-1*c^-1*(c*csc(a+b*x))^-1//2*cos(a+b*x)+2//3*b^-1*c^-2*(c*csc(a+b*x))^(1/2)*sin(a+b*x)^(1/2)*Elliptic.F((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate((c*csc(a+b*x))^-5//2, x) == :(-2//5*b^-1*c^-1*(c*csc(a+b*x))^-3//2*cos(a+b*x)+6//5*b^-1*c^-2*(c*csc(a+b*x))^-1//2*sin(a+b*x)^-1//2*Elliptic.E((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate((c*csc(a+b*x))^-7//2, x) == :(-10//21*b^-1*c^-3*(c*csc(a+b*x))^-1//2*cos(a+b*x)+-2//7*b^-1*c^-1*(c*csc(a+b*x))^-5//2*cos(a+b*x)+10//21*b^-1*c^-4*(c*csc(a+b*x))^(1/2)*sin(a+b*x)^(1/2)*Elliptic.F((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate((csc(x)^2)^7//2, x) == :(-5//16*arcsinh(cot(x))+-5//16*(csc(x)^2)^(1/2)*cot(x)+-5//24*(csc(x)^2)^3//2*cot(x)+-1//6*(csc(x)^2)^5//2*cot(x))
@test integrate((csc(x)^2)^5//2, x) == :(-3//8*arcsinh(cot(x))+-3//8*(csc(x)^2)^(1/2)*cot(x)+-1//4*(csc(x)^2)^3//2*cot(x))
@test integrate((csc(x)^2)^3//2, x) == :(-1//2*arcsinh(cot(x))+-1//2*(csc(x)^2)^(1/2)*cot(x))
@test integrate((csc(x)^2)^(1/2), x) == :(-1*arcsinh(cot(x)))
@test integrate((csc(x)^2)^-1//2, x) == :(-1*(csc(x)^2)^-1//2*cot(x))
@test integrate((csc(x)^2)^-3//2, x) == :(-2//3*(csc(x)^2)^-1//2*cot(x)+-1//3*(csc(x)^2)^-3//2*cot(x))
@test integrate((csc(x)^2)^-5//2, x) == :(-8//15*(csc(x)^2)^-1//2*cot(x)+-4//15*(csc(x)^2)^-3//2*cot(x)+-1//5*(csc(x)^2)^-5//2*cot(x))
@test integrate((csc(x)^2)^-7//2, x) == :(-16//35*(csc(x)^2)^-1//2*cot(x)+-8//35*(csc(x)^2)^-3//2*cot(x)+-6//35*(csc(x)^2)^-5//2*cot(x)+-1//7*(csc(x)^2)^-7//2*cot(x))
@test integrate((a*csc(x)^2)^7//2, x) == :(-5//16*a^7//2*arctanh(a^(1/2)*(a*csc(x)^2)^-1//2*cot(x))+-5//16*a^3*(a*csc(x)^2)^(1/2)*cot(x)+-5//24*a^2*(a*csc(x)^2)^3//2*cot(x)+-1//6*a*(a*csc(x)^2)^5//2*cot(x))
@test integrate((a*csc(x)^2)^5//2, x) == :(-3//8*a^5//2*arctanh(a^(1/2)*(a*csc(x)^2)^-1//2*cot(x))+-3//8*a^2*(a*csc(x)^2)^(1/2)*cot(x)+-1//4*a*(a*csc(x)^2)^3//2*cot(x))
@test integrate((a*csc(x)^2)^3//2, x) == :(-1//2*a^3//2*arctanh(a^(1/2)*(a*csc(x)^2)^-1//2*cot(x))+-1//2*a*(a*csc(x)^2)^(1/2)*cot(x))
@test integrate((a*csc(x)^2)^(1/2), x) == :(-1*a^(1/2)*arctanh(a^(1/2)*(a*csc(x)^2)^-1//2*cot(x)))
@test integrate((a*csc(x)^2)^-1//2, x) == :(-1*(a*csc(x)^2)^-1//2*cot(x))
@test integrate((a*csc(x)^2)^-3//2, x) == :(-1//3*(a*csc(x)^2)^-3//2*cot(x)+-2//3*a^-1*(a*csc(x)^2)^-1//2*cot(x))
@test integrate((a*csc(x)^2)^-5//2, x) == :(-1//5*(a*csc(x)^2)^-5//2*cot(x)+-8//15*a^-2*(a*csc(x)^2)^-1//2*cot(x)+-4//15*a^-1*(a*csc(x)^2)^-3//2*cot(x))
@test integrate((a*csc(x)^2)^-7//2, x) == :(-1//7*(a*csc(x)^2)^-7//2*cot(x)+-16//35*a^-3*(a*csc(x)^2)^-1//2*cot(x)+-8//35*a^-2*(a*csc(x)^2)^-3//2*cot(x)+-6//35*a^-1*(a*csc(x)^2)^-5//2*cot(x))
@test integrate((a*csc(x)^3)^5//2, x) == :(-154//585*a^2*(a*csc(x)^3)^(1/2)*cot(x)+-154//195*a^2*(a*csc(x)^3)^(1/2)*cos(x)*sin(x)+-22//117*a^2*(a*csc(x)^3)^(1/2)*csc(x)^2*cot(x)+-2//13*a^2*(a*csc(x)^3)^(1/2)*csc(x)^4*cot(x)+154//195*a^2*(a*csc(x)^3)^(1/2)*sin(x)^3//2*Elliptic.E(-1//2*x+1//4*pi,2))
@test integrate((a*csc(x)^3)^3//2, x) == :(-10//21*a*(a*csc(x)^3)^(1/2)*cos(x)+-10//21*a*(a*csc(x)^3)^(1/2)*sin(x)^3//2*Elliptic.F(-1//2*x+1//4*pi,2)+-2//7*a*(a*csc(x)^3)^(1/2)*cot(x)*csc(x))
@test integrate((a*csc(x)^3)^(1/2), x) == :(-2*(a*csc(x)^3)^(1/2)*cos(x)*sin(x)+2*(a*csc(x)^3)^(1/2)*sin(x)^3//2*Elliptic.E(-1//2*x+1//4*pi,2))
@test integrate((a*csc(x)^3)^-1//2, x) == :(-2//3*(a*csc(x)^3)^-1//2*cot(x)+-2//3*(a*csc(x)^3)^-1//2*sin(x)^-3//2*Elliptic.F(-1//2*x+1//4*pi,2))
@test integrate((a*csc(x)^3)^-3//2, x) == :(-14//45*a^-1*(a*csc(x)^3)^-1//2*cos(x)+-14//15*a^-1*(a*csc(x)^3)^-1//2*sin(x)^-3//2*Elliptic.E(-1//2*x+1//4*pi,2)+-2//9*a^-1*(a*csc(x)^3)^-1//2*sin(x)^2*cos(x))
@test integrate((a*csc(x)^3)^-5//2, x) == :(-26//77*a^-2*(a*csc(x)^3)^-1//2*cot(x)+-78//385*a^-2*(a*csc(x)^3)^-1//2*cos(x)*sin(x)+-26//77*a^-2*(a*csc(x)^3)^-1//2*sin(x)^-3//2*Elliptic.F(-1//2*x+1//4*pi,2)+-26//165*a^-2*(a*csc(x)^3)^-1//2*sin(x)^3*cos(x)+-2//15*a^-2*(a*csc(x)^3)^-1//2*sin(x)^5*cos(x))
@test integrate((a*csc(x)^4)^7//2, x) == :(-1*a^3*(a*csc(x)^4)^(1/2)*cos(x)*sin(x)+-3*a^3*(a*csc(x)^4)^(1/2)*cos(x)^2*cot(x)^3+-2*a^3*(a*csc(x)^4)^(1/2)*cos(x)^2*cot(x)+-20//7*a^3*(a*csc(x)^4)^(1/2)*cos(x)^2*cot(x)^5+-6//11*a^3*(a*csc(x)^4)^(1/2)*cos(x)^2*cot(x)^9+-5//3*a^3*(a*csc(x)^4)^(1/2)*cos(x)^2*cot(x)^7+-1//13*a^3*(a*csc(x)^4)^(1/2)*cos(x)^2*cot(x)^11)
@test integrate((a*csc(x)^4)^5//2, x) == :(-1*a^2*(a*csc(x)^4)^(1/2)*cos(x)*sin(x)+-6//5*a^2*(a*csc(x)^4)^(1/2)*cos(x)^2*cot(x)^3+-4//3*a^2*(a*csc(x)^4)^(1/2)*cos(x)^2*cot(x)+-4//7*a^2*(a*csc(x)^4)^(1/2)*cos(x)^2*cot(x)^5+-1//9*a^2*(a*csc(x)^4)^(1/2)*cos(x)^2*cot(x)^7)
@test integrate((a*csc(x)^4)^3//2, x) == :(-1*a*(a*csc(x)^4)^(1/2)*cos(x)*sin(x)+-2//3*a*(a*csc(x)^4)^(1/2)*cos(x)^2*cot(x)+-1//5*a*(a*csc(x)^4)^(1/2)*cos(x)^2*cot(x)^3)
@test integrate((a*csc(x)^4)^(1/2), x) == :(-1*(a*csc(x)^4)^(1/2)*cos(x)*sin(x))
@test integrate((a*csc(x)^4)^-1//2, x) == :(-1//2*(a*csc(x)^4)^-1//2*cot(x)+(1/2)*x*(a*csc(x)^4)^-1//2*csc(x)^2)
@test integrate((a*csc(x)^4)^-3//2, x) == :(-5//16*a^-1*(a*csc(x)^4)^-1//2*cot(x)+-5//24*a^-1*(a*csc(x)^4)^-1//2*cos(x)*sin(x)+-1//6*a^-1*(a*csc(x)^4)^-1//2*sin(x)^3*cos(x)+5//16*x*a^-1*(a*csc(x)^4)^-1//2*csc(x)^2)
@test integrate((a*csc(x)^4)^-5//2, x) == :(-63//256*a^-2*(a*csc(x)^4)^-1//2*cot(x)+-21//128*a^-2*(a*csc(x)^4)^-1//2*cos(x)*sin(x)+-21//160*a^-2*(a*csc(x)^4)^-1//2*sin(x)^3*cos(x)+-9//80*a^-2*(a*csc(x)^4)^-1//2*sin(x)^5*cos(x)+-1//10*a^-2*(a*csc(x)^4)^-1//2*sin(x)^7*cos(x)+63//256*x*a^-2*(a*csc(x)^4)^-1//2*csc(x)^2)

max_stars_repo_path: test/4 Trig functions/4.6 Cosecant/4.6.0 (a csc)^m (b trg)^n.jl


text:
const RED     = "\x1b[31m"
const GREEN   = "\x1b[32m"
const BOLD    = "\x1b[1m"
const DEFAULT = "\x1b[0m"

colored(s::AbstractString, color) = string(color, s, DEFAULT)

red(s::AbstractString)   = colored(s, RED)
green(s::AbstractString) = colored(s, GREEN)
bold(s::AbstractString)  = colored(s, BOLD)


text:
using Test

# Identify files in test/ that are testing matching files in src/
#  - src/Source.jl will be matched by test/test_Source.jl

filebase = String[]
for (root, dirs, files) in walkdir("../src")
    append!(filebase,
            map(file -> replace(file, r"(.*).jl" => s"\1"),
                filter(file -> occursin(r".*\.jl", file), files)))
end

testbase = map(file -> replace(file, r"test_(.*).jl" => s"\1"),
                filter(str -> occursin(r"^test_.*\.jl$", str), readdir()))

@testset "EcoBase.jl" begin
    println()
    @info "Running tests for files:"
    for t in testbase
        println("    = $t.jl")
    end
    println()

    @info "Running tests..."
    @testset for t in testbase
        fn = "test_$t.jl"
        println("    * Testing $t.jl ...")
        include(fn)
    end
end

# Identify tests with no matching file
superfluous = filter(f -> f ∉ filebase, testbase)
if length(superfluous) > 0
    println()
    @info "Potentially superfluous tests:"
    for f in superfluous
        println("    + $f.jl")
    end
end

# Identify files with no matching test
notest = filter(f -> f ∉ testbase, filebase)
if length(notest) > 0
    println()
    @info "Potentially missing tests:"
    for f in notest
        println("    - $f.jl")
    end
end


text:
#SkyTEM 10 ohm-m response for verify_skytem_td.jl
looparea = 499.34
rRx = 17.
zRxLM = -31.9
zRxHM = -30.19
zTx = -30.0
rTx = 12.607
lowpassfcs = [300000, 450000.0]
## model
zfixed   = [-1e6,   0]
rho      = [1e12,   10]
## response provided
dLM=[ 5.738414E-06
 4.236351E-06
 3.111256E-06
 2.243019E-06
 1.608307E-06
 1.161711E-06
 8.251093E-07
 5.775482E-07
 3.972306E-07
 2.682531E-07
 1.798862E-07
 1.183093E-07
 7.655095E-08
 4.898255E-08
 3.071389E-08
 1.892388E-08
 1.144417E-08
 6.789652E-09
 3.946301E-09
]/looparea

dHM = [
8.457780E-07
4.794991E-07
2.861063E-07
1.729966E-07
1.053307E-07
6.465972E-08
3.950176E-08
2.406247E-08
1.459279E-08
8.806743E-09
5.283412E-09
3.147092E-09
1.863403E-09
1.094271E-09
6.370838E-10
3.673768E-10
2.096507E-10
1.181879E-10
6.198442E-11
]/looparea

## LM times and ramp
LM_times = [ 1.300E-05
 1.650E-05
 2.100E-05
 2.700E-05
 3.450E-05
 4.350E-05
 5.500E-05
 6.950E-05
 8.800E-05
 1.115E-04
 1.405E-04
 1.775E-04
 2.240E-04
 2.820E-04
 3.555E-04
 4.480E-04
 5.645E-04
 7.110E-04
 8.955E-04]

LM_ramp =[-1.000E-03   0.00
-9.800E-04   0.46
-9.550E-04   0.79
-9.200E-04   0.94
-9.100E-04   0.95
-8.000E-04   1.00
0.000E+00    1.00
4.400E-07    0.98
8.000E-07    0.93
1.100E-06    0.87
1.580E-06    0.69
2.160E-06    0.46
2.940E-06    0.22
3.900E-06    0.07
5.340E-06    0.01
6.480E-06    0.00]
## HM times and ramp
HM_times = [8.570E-05
 1.092E-04
 1.382E-04
 1.752E-04
 2.217E-04
 2.797E-04
 3.532E-04
 4.457E-04
 5.622E-04
 7.087E-04
 8.932E-04
 1.126E-03
 1.418E-03
 1.786E-03
 2.250E-03
 2.833E-03
 3.568E-03
 4.493E-03
 5.808E-03]

HM_ramp = [-1.000E-02   0.00
-8.700E-03   0.36
-7.100E-03   0.64
-5.500E-03   0.81
-3.300E-03   0.93
-1.200E-03   0.98
-5.000E-05   1.00
 0.000E+00   1.00
 3.700E-06   0.94
 1.540E-05   0.72
 3.080E-05   0.43
 4.220E-05   0.21
 5.090E-05   0.01
 5.500E-05   0.00 ]

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
@everywhere using BioSymbols
@everywhere using BioSequences
using ArgParse
using Glob

@everywhere using Rifraf

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        # script-only arguments
        "--phred-cap"
        help = "maximum PHRED score"
        arg_type = Int8
        default = Int8(0)

        "--prefix"
        help = "prepended to each filename to make name"
        arg_type = String
        default = ""

        "--keep-unique-name"
        help = "keep only unique middle part of filename"
        action = :store_true

        # RIFRAF arguments
        "--reference"
        help = string("reference fasta file.",
                      " uses first sequence unless --reference-map is given.")
        arg_type = String
        default = ""

        "--reference-map"
        help = "file mapping input filename to reference id"
        arg_type = String
        default = ""

        "--ref-errors"
        help = "comma-seperated reference error ratios - mm, codon ins, codon del"
        arg_type = String
        default = "10,0.1,0.1,1,1"

        "--max-iters"
        help = "maximum iterations before giving up"
        arg_type = Int
        default = 100

        "--verbose", "-v"
        help = "print progress"
        arg_type = Int
        default = 0

        # required arguments
        "seq-errors"
        help = "comma-seperated sequence error ratios - mismatch, insertion, deletion"
        arg_type = String
        required = true

        "input"
        help = "a single file or a glob. each filename should be unique."
        required = true

        "output"
        help = "output fasta file"
        required = true

    end
    return parse_args(s)
end

@everywhere function dofile(file, reffile, refid, args)
    if args["verbose"] >= 1
        println(STDERR, "reading sequences from '$(file)'")
        if length(reffile) > 0
            println(STDERR, "reading reference from '$(reffile)'")
        end
    end
    reference = DNASeq()
    ref_records = []
    if length(reffile) > 0
        ref_records = Rifraf.read_fasta_records(reffile)
    end
    if length(refid) > 0
        filt_records = collect(filter(r -> FASTA.identifier(r) == refid, ref_records))
        if length(filt_records) == 0
            error("reference '$refid' not found in `$reffile`")
        end
        if length(filt_records) > 1
            error("multiple references with id '$refid' found in `$reffile`")
        end
        ref_record = filt_records[1]
        reference = DNASeq(sequence(ref_record))
    elseif length(ref_records) > 0
        ref_record = ref_records[1]
        reference = DNASeq(sequence(ref_record))
    end

    score_args = map(x -> parse(Float64, x), split(args["seq-errors"], ","))
    scores = Scores(ErrorModel(score_args...))

    ref_score_args = map(x -> parse(Float64, x), split(args["ref-errors"], ","))
    ref_scores = Scores(ErrorModel(ref_score_args...))

    sequences, phreds, _ = Rifraf.read_fastq(file, seqtype=DNASeq)
    if args["verbose"] >= 1
        println(STDERR, "starting run")
    end
    phred_cap = args["phred-cap"]
    if phred_cap > 0
        phreds = Vector{Phred}[Rifraf.cap_phreds(p, phred_cap)
                               for p in phreds]
    end
    params = RifrafParams(scores=scores,
                          ref_scores=ref_scores,
                          max_iters=args["max-iters"],
                          verbose=args["verbose"])
    return rifraf(sequences, phreds;
                  reference=reference,
                  params=params)
end

function common_prefix(strings)
    minlen = minimum([length(s) for s in strings])
    cand = strings[1][1:minlen]
    x = 0
    for i = 1:minlen
        if all([s[i] == cand[i] for s in strings])
            x = i
        else
            break
        end
    end
    return cand[1:x]
end

function common_suffix(strings)
    return reverse(common_prefix([reverse(s) for s in strings]))
end

function main()
    args = parse_commandline()
    input = args["input"]

    dir, pattern = splitdir(input)
    infiles = glob(pattern, dir)
    if length(infiles) == 0
        if args["verbose"] >= 1
            println(STDERR, "warning: no input files found.")
        end
        return
    end
    infiles = sort(infiles)
    basenames = [basename(f) for f in infiles]
    if length(Set(basenames)) != length(basenames)
        error("Files do not have unique names")
    end

    if length(args["reference"]) > 0
        if !isfile(args["reference"])
            error("reference file not found")
        end
        if length(args["reference-map"]) > 0
            if !isfile(args["reference-map"])
                error("reference map file not found")
            end
        end
    else
        if length(args["reference-map"]) > 0
            error("--reference-map is invalid without --reference")
        end
    end

    reffile = args["reference"]
    refmapfile = args["reference-map"]
    refids = ["" for i in 1:length(infiles)]
    if length(refmapfile) > 0
        refmap_handle = open(refmapfile)
        name_to_ref = Dict()
        for line in readlines(refmap_handle)
            name, refid = split(line)
            name_to_ref[name] = refid
        end
        close(refmap_handle)
        # only process seqs with a reference
        present_names = Set(keys(name_to_ref))
        infiles = sort(collect(filter(f -> basename(f) in present_names, infiles)))
        basenames = [basename(f) for f in infiles]
        refids = [name_to_ref[name] for name in basenames]
    end
    results = pmap((f, rid) -> dofile(f, reffile, rid, args),
                   infiles, refids)

    plen = 0
    slen = 0
    if args["keep-unique-name"]
        plen = length(common_prefix(basenames))
        snames = [n[plen+1:end] for n in basenames]
        slen = length(common_suffix(snames))
    end

    n_converged = 0
    prefix = args["prefix"]
    stream = open(FASTA.Writer, args["output"])
    for (i, result) in enumerate(results)
        if typeof(result) == RemoteException
            throw(result)
        end
        if result.state.converged
            n_converged += 1
            name = basenames[i]
            if args["keep-unique-name"]
                name = name[plen + 1:end - slen]
            end
            seqname = string(prefix, name)
            record = FASTA.Record(seqname, result.consensus)
            write(stream, record)
        end
    end
    close(stream)
    if args["verbose"] >= 1
        println(STDERR, "done. $n_converged / $(length(results)) converged.")
    end
end

main()


text:
using AtomGraphs
using Test

@testset "AtomGraphs.jl" begin
    include("AtomGraph_tests.jl")
    include("graph_building_tests.jl")
end