text:
using CGPrism3D
using Test
using LinearAlgebra
#using StatsBase

#const K = 3
#const x = K+1
#const y = K/2

@time dataT = dataTet();
@time dataF = dataFsymb();
@time dataA = dataPrsmA();

# Testing Penstagon Identity -- related to 3-2 Pachner move

#@time f2 = tensorGlue(dataF,dataF,[2,4,6],[2,3,1])
#@time f3 = tensorGlue(f2,dataF,[5,8,7,3,4],[1,2,4,5,6])
#@time fs1 = tensorSum(f3,[4])
#@show sum(f2[2]-fs1[2])

#@time t2 = tensorGlueTet3D(dataT,dataT,[1,2,3],[1,2,3])
#@time t3 = tensorGlueTet3D(t2,dataT,[1,5,6,9,8],[1,2,3,5,6])
#@time ts1 = tensorSumTet3D(t3,[1])
#@show sum(t2[2]-ts1[2])

# Tests on Prisms -- Full Algorithm

# Create a prism from three tetrahedra -- equivalent to dataA
@time tet2 = tensorGlue(dataT,dataT,[1,2,3],[1,2,3])
@time prAn = tensorGlue(tet2,dataT,[2,7,9],[1,2,3])
sp = sortperm(prAn[1])
prA = prAn[1][sp], prAn[2][sp]
@show sum(prA[2]-dataA[2])

#@time fullSplitTet3D(prA,[4,5,6],[7,8,9,10,11,12],[1,2,3]);
#@time fullSplit3D(prA,[4,5,6],[7,8,9,10,11,12],[1,2,3]);

## 2-2 move on faces [5,6,9,8,1] and [4,6,11,12,2] to switch diagonals -- gets prism B

#@time fa1 = tensor22moveB(prA,[5,6,8,9,1]) # face [i,l,j,k,n]
#@time prB = tensor22moveB(fa1,[4,6,12,11,2])
#@show sum(prB[2]-prA[2])

# Split prisms into tetrahedra and pyramid

#@time spltA = fullSplitTet3D(prA,[4,5,6],[7,8,9,10,11,12],[1,2,3]) # split into U, V
#ua, va = spltA[1], spltA[2]

#@time spltB = fullSplitTet3D(prB,[6,9,11],[3,4,5,7,8,12],[10,2,1]) # split into U, V
#ub, vb = spltB[1], spltB[2]

#@show (sum(ub[2]-dataT[2]),sum(vb[2]-tet2[2]))
#@show (sum(ua[2]-dataT[2]),sum(va[2]-tet2[2]))

# Glue tetrahedra and pyramids back into effective prism

#@time uavb = tensorGlueTet3D(ua,vb,[1,3,5],[2,6,8])
#@time ubva = tensorGlueTet3D(ub,va,[1,3,5],[6,5,8])
#@show (sum(uavb[2]-prA[2]),sum(ubva[2]-prA[2]))

#Using only Prism A to coarse grain

@time pruvA = tensorGlue(prA,prA,[1,5,6,8,9],[2,4,6,12,11])
@time pruvAA = tensorSum(pruvA,[1])

# 2-2 move /F-move on faces 0,1,2,(3,4)- not necessary ?
#@time fA0FN = tensor22moveA(pruvAA,[18,15,6,2,7])
@time fA0F = tensor22moveA(pruvAA,[18,15,6,2,7])
@time fA0F = tensor22moveA(pruvAA,[18,15,6,2,7])
#@time fA1F = tensor22move(pruvAA,[17,18,10,9,8])
#@time fA2F = tensor22move(pruvAA,[14,13,3,2,4])
#@time fA3F = tensor22move(pruvAA,[17,16,14,5,12])
#@time fA4F = tensor22move(pruvAA,[10,11,5,3,1])
#@show (length(fA0F[2]),length(fA1F[2]),length(fA2F[2]),length(fA3F[2]),length(fA4F[2]))

@time fA01F = tensor22moveB(fA0F,[17,18,10,9,8])
@time fA012F = tensor22moveB(fA01F,[14,13,3,2,4])
#@time fA0123F = tensor22move(fA012F,[17,16,14,5,12])
#@time fA01234F = tensor22move(fA0123F,[10,11,5,3,1])
#@show (length(fA012F[2]),length(fA0123F[2]),length(fA01234F[2]))

@time f133A = fullSplit3D(fA012F,[9,6,18],[1,2,3,4,5,10,12,13,14,15,16,17],[7,8,11])
#@time f233A = fullSplitTet3D(fA012F,[13,15,2],[1,3,5,6,8,9,10,11,12,14,17,18],[7,4,16])
#@show (length(f133A[1][2]),length(f133A[2][2]),length(f233A[1][2]),length(f233A[2][2]))

@time CGPrsmA1 = fullSplit3D(f133A[2],[8,10,2],[1,3,5,6,7,9,12,14,15],[13,4,11])
#@time CGPrsmA2 = fullSplitTet3D(f233A[2],[6,4,12],[1,2,3,7,9,10,11,14,15],[13,5,8])
#@show (length(CGPrsmA1[1][2]),length(CGPrsmA1[2][2]),length(CGPrsmA2[1][2]),length(CGPrsmA2[2][2]))
@show length(CGPrsmA1[1][2]),length(CGPrsmA1[2][2])
@show ( sum( CGPrsmA1[2][2] - prA[2]) ) #, sum( CGPrsmA2[2][2] - prA[2]) )

#@show sort(permuteInd.(CGPrsmA1[2][1])) == (prA[1])
@time cgprsm = tensorPermute(CGPrsmA1[2])
@show length(cgprsm[1])
@show sort(cgprsm[1]) == sort(prA[1])

# Now Using prsm A and B to glue and coarsegrain

#@time pruv = tensorGlueTet3D(uavb,ubva,[1,2,6,7,8],[3,2,4,10,9]) # get effective prism
#@time pruvAB = tensorSumTet3D(pruv,[1])# sum over bulk index

# Coarse graining procedures
# 2-2 move /F-move on faces 0,1,2

#@time Face0 = tensor22move(pruvAB,[12,15,8,9,6])
#@time Face1 = tensor22move(pruvAB,[14,12,3,2,1])
#@time Face2 = tensor22move(pruvAB,[17,16,11,9,7])
#@show (length(Face0[2]),length(Face1[2]),length(Face2[2]))

# It doesn't matter the order of moves
#@time Face01 = tensor22move(Face0,[14,12,3,2,1])
#@time Face10 = tensor22move(Face1,[12,15,8,9,6])
#@time Face02 = tensor22move(Face0,[17,16,11,9,7])
#@time Face20 = tensor22move(Face2,[12,15,8,9,6])
#@time Face12 = tensor22move(Face1,[17,16,11,9,7])
#@time Face21 = tensor22move(Face2,[14,12,3,2,1])
#@show (length(Face01[2]),length(Face10[2]),length(Face02[2]),length(Face20[2]),length(Face12[2]),length(Face21[2]))

#@time Face012 = tensor22move(Face01,[17,16,11,9,7])
#@time Face102 = tensor22move(Face10,[17,16,11,9,7])
#@show (length(Face012[2]),length(Face102[2]))

# Final Pachner moves-- 3-1 pachner moves using SVD

#@time CGtet1 = fullSplitTet3D(Face012,[8,12,2],[3,4,5,7,9,11,13,14,15,16,17,18],[1,6,10])
#@time CGtet2 = fullSplitTet3D(Face012,[15,16,9],[1,2,3,4,5,8,10,11,12,13,14,17],[7,6,18])
#@show (length(CGtet1[1][1]),length(CGtet1[2][1]),length(CGtet2[1][1]),length(CGtet2[2][1]))

#@time CGtetfull1 = fullSplitTet3D(CGtet2[2],[6,2,9],[3,4,5,8,10,11,12,13,15],[1,14,7])
#@time CGtetfull2 = fullSplitTet3D(CGtet1[2],[9,10,5],[1,2,3,6,7,8,11,13,15],[4,14,12])
#@show (length(CGtetfull1[2][1]),length(CGtetfull1[2][2]),length(CGtetfull2[2][1]),length(CGtetfull2[2][2]) )
#@show fans = sum( CGtetfull2[2][2] - dataA[2])
#if numchop(fans) == 0
#    println("SUCCESS -- :)")
#end


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

using CSV
using DataFrames
using Plots
using StatsPlots
using Pipe
using Statistics
using RCall
using ProgressMeter
using LinearAlgebra
using StatsFuns
using MCPhylo
using Random
Random.seed!(8074252001070266643)

##

wals = CSV.read("../data/charMtx.csv", DataFrame)

famFreqs = combine(groupby(wals, :glot_fam), nrow)
sort!(famFreqs, :nrow, rev=true)

families = filter(x -> x.nrow>1, famFreqs).glot_fam
isolates = filter(x -> x.nrow==1, famFreqs).glot_fam

lineages = vcat(families, isolates)

isoDict = @pipe wals |>
      filter(x -> x.glot_fam ∈ isolates, _) |>
      zip(_.glot_fam, _.longname) |>
      Dict

##

d = CSV.read("../data/fpairMtx.csv", DataFrame)
d = innerjoin(d, wals[:,[:longname, :glot_fam]], on=:taxon => :longname)
##

function myvar(x)
    if length(x) == 1
        return 0.
    end
    var(x, corrected=false)
end

##

function getEmpiricalStatistics(charNum)
    fp = names(d)[charNum+1]
    f1, f2 = Symbol.(split(fp, "-"))
    binaryRepresentation = select(
        innerjoin(
            d[:, [1, charNum + 1]],
            wals[:, [:glot_fam, :longname, f1, f2]],
            on = :taxon => :longname,
        ) |> dropmissing,
        ["taxon", "glot_fam", fp],
    )
    binaryRepresentation[!, :a] = Int.(binaryRepresentation[:, 3] .== "a")
    binaryRepresentation[!, :b] = Int.(binaryRepresentation[:, 3] .== "b")
    binaryRepresentation[!, :c] = Int.(binaryRepresentation[:, 3] .== "c")
    binaryRepresentation[!, :d] = Int.(binaryRepresentation[:, 3] .== "d")
    totalVariance = sum(mapslices(myvar, Array(binaryRepresentation[:, 4:end]), dims = 1))
    lineagewiseVariance = combine(
        x -> sum(mapslices(myvar, Array(x[:, 4:end]), dims = 1)),
        groupby(binaryRepresentation, :glot_fam),
    )
    rename!(lineagewiseVariance, [:glot_fam, :variance])
    lvar = mean(lineagewiseVariance.variance)
    varianceBetweenLineages = @pipe binaryRepresentation |>
          groupby(_, :glot_fam) |>
          combine(x -> mapslices(mean, Array(x[:, 4:end]), dims = 1), _) |>
          Array(_[:, 2:end]) |>
          mapslices(myvar, _, dims = 1) |>
          sum
    totalVariance, lvar, varianceBetweenLineages
end

empiricalStats = DataFrame(getEmpiricalStatistics.(1:28))
rename!(empiricalStats, [:totalVariance, :lineagewiseVariance, :crossLineageVariance])

##

treefiles = readdir("../data/posteriorTrees", join=true);

##

fm2trees = Dict()
@showprogress for fm in families
      fm2trees[fm] = MCPhylo.ParseNewick("../data/posteriorTrees/$fm.posterior.tree")
end

for fm in isolates
      fm2trees[fm] = repeat([Node(isoDict[fm])], 1000)
end

##

offDiag = [CartesianIndex(x,y) for y in 1:4 for x in 1:4 if x!=y];
taxa = d.taxon;

##

nSims = 1000

stats = []

##
@showprogress for g in 1:nSims
    treeI = sample(1:1000)

    simChar = Dict{String,Int}()

    global fullSrates = rand(LogNormal(0, 1), 12)
    q = zeros(4, 4)
    q[offDiag] = fullSrates
    for i = 1:4
        q[i, i] -= sum(q[i, :])
    end

    D, U = eigen(q)
    Uinv = inv(U)

    eq = real.(Uinv[4, :])
    eq /= sum(eq)

    for fmI = 1:length(lineages)
        fmI
        tree = fm2trees[lineages[fmI]][treeI]
        preo = pre_order(tree)
        tStates = Vector{Int}(undef, length(preo))
        tStates[tree.num] = rand(Categorical(eq))
        for nd in preo[2:end]
            mState = tStates[nd.mother.num]
            p = real.(U * diagm(exp.(nd.inc_length * D)) * Uinv)
            dState = rand(Categorical(p[mState, :]))
            tStates[nd.num] = dState
        end
        for nd in preo
            if nd.nchild == 0
                simChar[nd.name] = tStates[nd.num]
            end
        end
    end

    x = [simChar[l] for l in taxa]
    x_bin = DataFrame(hcat([x .== i for i = 1:4]...))
    x_bin[!, :glot_fam] = d.glot_fam

    tv = mapslices(myvar, Array(x_bin[:, 1:4]), dims = 1) |> sum
    lmv = @pipe x_bin |>
          groupby(_, :glot_fam) |>
          combine(x -> sum(mapslices(myvar, Array(x[:, 1:4]), dims = 1)), _) |>
          mean(_[:, 2])

    vc = @pipe x_bin |>
          groupby(_, :glot_fam) |>
          combine(x -> mapslices(mean, Array(x[:, 1:4]), dims = 1), _) |>
          mapslices(myvar, Array(_[:, 2:end]), dims = 1) |>
          sum

    push!(stats, [tv, lmv, vc])
end

##
uSim = DataFrame(Tuple.(stats))
rename!(uSim, names(empiricalStats))

##

df = copy(empiricalStats)
df[:,:method] .= "empirical"
df[:,:model] .= "nature"

uSim[:,:method] .= "simulated"
uSim[:,:model] .= "universal"

df = vcat(df, uSim)

##

nSims = 1000

lstats = zeros(nSims, 3)

@showprogress for g in 1:nSims
    global q, fullSrates, p
    treeI = sample(1:1000)
    simChar = Dict{String,Int}()
    for fmI in 1:length(lineages)
        fullSrates = rand(LogNormal(0,1), 12)
        q = zeros(4,4)
        q[offDiag] = fullSrates
        for i in 1:4
            q[i,i] -= sum(q[i,:])
        end

        D, U = eigen(q)
        Uinv = inv(U)

        eq = abs.(Uinv[4,:])
        eq /= sum(eq)
        tree = fm2trees[lineages[fmI]][treeI]
        preo = pre_order(tree)
        tStates = Vector{Int}(undef, length(preo))
        tStates[tree.num] = rand(Categorical(eq/sum(eq)))
        for nd in preo[2:end]
              mState = tStates[nd.mother.num]
              p = real.(U * diagm(exp.(nd.inc_length * D)) * Uinv)
              ps = abs.(p[mState,:]/sum(p[mState,:]))
              ps /= sum(ps)
              dState = rand(Categorical(ps))
              tStates[nd.num] = dState
        end
        for nd in preo
            if nd.nchild == 0
                simChar[nd.name] = tStates[nd.num]
            end
        end
    end

    x = [simChar[l] for l in taxa]
    x_bin = DataFrame(hcat([x .== i for i in 1:4]...))
    x_bin[!,:glot_fam] = d.glot_fam
    x_bin

    tv = mapslices(myvar, Array(x_bin[:,1:4]), dims=1) |> sum
    lmv = @pipe x_bin |> groupby(_, :glot_fam) |>
        combine(x -> sum(mapslices(myvar, Array(x[:,1:4]), dims=1)), _) |>
        mean(_[:,2])

    vc = @pipe x_bin |> groupby(_, :glot_fam) |>
        combine(x -> mapslices(mean, Array(x[:,1:4]), dims=1), _) |>
    mapslices(myvar, Array(_[:,2:end]), dims=1) |> sum

    lstats[g,:] = [tv, lmv, vc]
end
##
lSim = DataFrame(lstats)
rename!(lSim, names(empiricalStats))

lSim[:,:method] .= "simulated"
lSim[:,:model] .= "lineageSpecific"

df = vcat(df, lSim)

##

theme(:solarized_light)

upscale = 1 #8x upscaling in resolution
fntsm = Plots.font("sans-serif", pointsize=round(20.0*upscale))
fntlg = Plots.font("sans-serif", pointsize=round(25.0*upscale))
default(titlefont=fntlg, guidefont=fntlg, tickfont=fntsm, legendfont=fntsm)
default(size=(800*upscale,600*upscale)) #Plot canvas size

##

u_tv_df = @pipe df |>
      filter(x -> x.model != "lineageSpecific", _) |>
      select(_, [:totalVariance, :method]) |>
      stack

p1 = @df u_tv_df boxplot(
    :method,
    :value,
    legend = false,
    grouping = :method,
    title = "total variance",
    linewidth=2*upscale
)

@df u_tv_df dotplot!(p1, :method, :value, markersize=2*upscale, color=:black, alpha=0.4)

##

u_lv_df = @pipe df |>
      filter(x -> x.model != "lineageSpecific", _) |>
      select(_, [:lineagewiseVariance, :method]) |>
      stack

p2 = @df u_lv_df boxplot(
    :method,
    :value,
    legend = false,
    grouping = :method,
    title = "lineage-wise variance",
    linewidth=2*upscale,
)
@df u_lv_df dotplot!(p2, :method, :value, markersize=2*upscale, color=:black, alpha=0.4)

##

u_clv_df = @pipe df |>
      filter(x -> x.model != "lineageSpecific", _) |>
      select(_, [:crossLineageVariance, :method]) |>
      stack

p3 = @df u_clv_df boxplot(
    :method,
    :value,
    legend = false,
    grouping = :method,
    title = "cross-lineage variance",
    linewidth=2*upscale
)
@df u_clv_df dotplot!(p3, :method, :value, markersize=2*upscale, color=:black, alpha=0.4)

##

try
    mkdir("../img")
catch e
end

savefig(p1, "../img/priorSimulationUniversal_1.png")
savefig(p2, "../img/priorSimulationUniversal_2.png")
savefig(p3, "../img/priorSimulationUniversal_3.png")

##

l_tv_df = @pipe df |>
      filter(x -> x.model != "universal", _) |>
      select(_, [:totalVariance, :method]) |>
      stack

p4 = @df l_tv_df boxplot(
    :method,
    :value,
    legend = false,
    grouping = :method,
    title = "total variance",
    linewidth=2*upscale
)

@df l_tv_df dotplot!(p4, :method, :value, markersize=2*upscale, color=:black, alpha=0.4)

##

l_lv_df = @pipe df |>
      filter(x -> x.model != "universal", _) |>
      select(_, [:lineagewiseVariance, :method]) |>
      stack

p5 = @df l_lv_df boxplot(
    :method,
    :value,
    legend = false,
    grouping = :method,
    title = "lineage-wise variance",
    linewidth=2*upscale,
)
@df l_lv_df dotplot!(p5, :method, :value, markersize=2*upscale, color=:black, alpha=0.4)

##

l_clv_df = @pipe df |>
      filter(x -> x.model != "universal", _) |>
      select(_, [:crossLineageVariance, :method]) |>
      stack

p6 = @df l_clv_df boxplot(
    :method,
    :value,
    legend = false,
    grouping = :method,
    title = "cross-lineage variance",
    linewidth=2*upscale
)
@df l_clv_df dotplot!(p6, :method, :value, markersize=2*upscale, color=:black, alpha=0.4)

##

savefig(p4, "../img/priorSimulationLineage_1.png")
savefig(p5, "../img/priorSimulationLineage_2.png")
savefig(p6, "../img/priorSimulationLineage_3.png")

ext: jl

max_stars_repo_name: gerhardJaeger/phylogeneticTypology

max_issues_repo_name: gerhardJaeger/phylogeneticTypology

max_forks_repo_name: gerhardJaeger/phylogeneticTypology


text:
using Documenter

makedocs(
    sitename = "AppliGate",
    format = Documenter.HTML(),
    pages = [
        "AppliGate" => "index.md",
        "Blog" => "blog.md",
        "Projects" => "projects.md",
        "Contact" => "contact.md",
        "Employees" => "employees.md"
    ]
)


text:
# find midpoint of a range without causing an overflow to Inf
function robustmidpoint(a::Float64,b::Float64)
	l = min(a,b)
	u = max(a,b)
	if sign(l) == sign(u)
		m = l + (u-l)/2
	else
		m = (u+l)/2
	end
	m
end


text:
using IntervalWavelets
using Test

# Side, index, resolution
expected_boundary_function_values = Dict{Tuple{Sides, Int, Int}, Vector{Float64}}(
    (LEFT, 0, 0) => [0.0 ; -1.1266 ; 0.0],
    (LEFT, 1, 0) => [0.0 ; 1.40586 ; -0.365497 ; 0.0],
    (LEFT, 0, 1) => [0.0 ; 0.412364 ; -1.1266 ; 0.206182 ; 0.0],
    (LEFT, 1, 1) => [0.0 ; 0.849475 ; 1.40586 ; -0.00765026 ; -0.365497 ; 0.0668906 ; 0.0],
    (RIGHT, 0, 0) => [0.0 ; 0.651639 ; 0.0],
    (RIGHT, 1, 0) => [0.0 ; 1.25341 ; 0.142971 ; 0.0],
    (RIGHT, 0, 1) => [0.0 ; 0.445077 ; 0.651639 ; 0.890155 ; 0.0],
    (RIGHT, 1, 1) => [0.0 ; 0.856098 ; 1.25341 ; 0.327042 ; 0.142971 ; -0.14055 ; 0.0]
)

@testset "$side boundary scaling functions" for side in [IntervalWavelets.LEFT, IntervalWavelets.RIGHT]
    @testset "Specific boundary scaling functions at resolution $R" for R in 0:1, index in 0:1
        phi = boundary_scaling_functions(side, 2, R)

        boundary_function_values = phi[index].(support(phi[index]))
        @test boundary_function_values ≈ expected_boundary_function_values[(side, index, R)] atol = 10.0^-4
    end

    @testset "Function values at 0" begin
        phi = boundary_scaling_functions(side, 2, 2)

        x = DyadicRational(0,0)
        if side == LEFT
            @test phi[0](x) ≈ 1.95133 atol = 10.0^-4
            @test phi[1](x) ≈ 0.293094 atol = 10.0^-4
        else
            @test phi[0](x) ≈ 1.128671 atol = 10.0^-4
            @test phi[1](x) ≈ -0.42407 atol = 10.0^-4
        end
    end

    @testset "Number of vanishing moments must be sufficiently high" begin
        @test_throws DomainError boundary_scaling_functions(side, 0, 0)
        @test_throws DomainError boundary_scaling_functions(side, 1, 0)
    end

    @testset "Increase resolution of boundary scaling function" for p in 2:3, index in 0:p - 1
        phi = boundary_scaling_functions(side, p, 0)

        phi0_support = support(phi[index])
        expected_support = IntervalWavelets.all_dyadic_rationals(support_boundaries(phi[index])..., 0)
        @test phi0_support == expected_support

        phi1 = IntervalWavelets.increase_resolution(phi)
        expected_support = IntervalWavelets.all_dyadic_rationals(support_boundaries(phi1[index])..., 1)
        @test support(phi1[index]) == expected_support

        @test phi1[index].(phi0_support) == phi[index].(phi0_support)
    end
end


text:
struct Visual
    plottype::PlotFunc
    attributes::NamedArguments
end
Visual(plottype=Any; kwargs...) = Visual(plottype, NamedArguments(kwargs))

function (v::Visual)(input::ProcessedLayer)
    plottype = Makie.plottype(v.plottype, input.plottype)
    attributes = merge(input.attributes, v.attributes)
    return ProcessedLayer(input; plottype, attributes)
end

visual(plottype=Any; kwargs...) = transformation(Visual(plottype; kwargs...))


text:
export rand_lines, rand_network, rand_network_circuit, ensemble, ensemble_circuit

using Random
using Distributions: Normal

"""
```julia
rand_lines(dims,n,l)
```
Generates a random collection of lines with initial point
laying within the rectangle between the origin and `dims`.
Determines second point of line to be distance 'l' and random
direction (as sampled unifromly from the N-1 sphere) from
the initial point.
"""
function rand_lines(dims::SVector{N,T},n::Int,l::T) where {N,T}
    x₀s = map(x->x.*dims,rand(SVector{N,T},n))
    A = rand(Normal(),N,n)
    # normalize columns of A to leave columns
    # as uniformly sampled on the N-1 sphere
    A .= A./sqrt.(sum(A.^2,dims=1)).*l
    x₁s = [SVector{N}(c) for c ∈ eachcol(A)]
    Line.(x₀s,x₀s+x₁s)
end

"""
```julia
rand_network(dims,n,l)
```
Generates a random `StickNetwork` using `rand_lines(dims,n,l)`
and adds the properties of resistivity and wire diameter.
"""
function rand_network(dims::SVector{N,T},n,l) where {N,T}
    lines = rand_lines(dims::SVector{N,T},n::Int,l::T)
    StickNetwork(lines, dims)
end

"""
```julia
rand_network_circuit(dims,n,l,ρ,D,V₊,Rⱼ,Rₑ,M)
```
Generates a random nanowire network using model `M` with dimensions `dims`,
number of wires `n`, wire resistivity `ρ`, wire diameter `D`, voltage `V₊`,
junction resistance `Rⱼ`, and electrode resistance `Rₑ`. Electrode is along
the leftmost extent of the region and ground along the rightmost extent.
Electrode and ground are positioned `eps(T)` in from the edges to minimize
rejected intersections.
"""
function rand_network_circuit(dims::SVector{2,T},n,l,ρ,D,V₊,Rⱼ,Rₑ,M) where {T}
    lines = rand_lines(dims::SVector{2,T},n::Int,l::T)
    wprops = repeat([WireProp(ρ,D)],n)
    elec = Line(SA[eps(T),zero(T)],SA[eps(T),dims[2]])
    grnd = Line(SA[dims[1]-eps(T),zero(T)],SA[dims[1]-eps(T),dims[2]])
    NWN_circuit{M}(lines, wprops, [elec], [V₊], [grnd], dims, Rⱼ, Rₑ)
end

"""
```julia
rand_network_circuit(dims,n,l,ρ,D,V₊,Rⱼ,Rₑ,M,Val(:conducting))
```
Generates a random nanowire network using model `M` which is gaurenteed to conduct. Otherwise the same as `rand_network_circuit`.
"""
function rand_network_circuit(dims::SVector{2,T},n,l,ρ,D,V₊,Rⱼ,Rₑ,M,::Val{:conducting}) where {T}
    lines = rand_lines(dims::SVector{2,T},n::Int,l::T)
    wprops = repeat([WireProp(ρ,D)],n)
    elec = Line(SA[eps(T),zero(T)],SA[eps(T),dims[2]])
    grnd = Line(SA[dims[1]-eps(T),zero(T)],SA[dims[1]-eps(T),dims[2]])
    nwn = NWN_circuit{M}(lines, wprops, [elec], [V₊], [grnd], dims, Rⱼ, Rₑ)
    while !is_conducting(nwn)
        lines = rand_lines(dims::SVector{2,T},n::Int,l::T)
        nwn = NWN_circuit{M}(lines, wprops, [elec], [V₊], [grnd], dims, Rⱼ, Rₑ)
    end
    nwn
end

"""
```julia
ensemble(n,dims,ρₙ,l)
```
Creates an ensemble of `StickNetwork`s with dimensions `dims`, stick densites `ρₙ`, and stick lengths `l`.
"""
function ensemble(n,dims::SVector{N,T},ρₙ,l) where {N,T}
    [rand_network(dims, round(Int, ρₙ*prod(dims)), l) for i∈1:n]
end

"""
```julia
ensemble_circuit(n,dims,ρₙ,l,ρ,D,V₊,Rⱼ,Rₑ,M,args...)
```
Generates an ensemble of nanowire networks by calling `rand_network_circuit(dims,round(Int,ρₙ*prod(dims)),l,ρ,D,V₊,Rⱼ,Rₑ,M,args...)` `n` times.
"""
function ensemble_circuit(n,dims::SVector{2,T},ρₙ,l,ρ,D,V₊,Rⱼ,Rₑ,M,args...) where {T}
    [rand_network_circuit(dims,round(Int,ρₙ*prod(dims)),l,ρ,D,V₊,Rⱼ,Rₑ,M,args...) for i∈1:n]
end


text:
len = length
dc  = deepcopy

macro clear()
  run(`clear`)
end

macro timeit( expression )
  maxtime = 2. # seconds
  tic()
  expression
  dt = toq()
  maxiters = int( maxtime / dt )
  times = zeros(maxiters)
  println("Running for maximum of $maxtime seconds...")
  for k = 1:maxiters
    tic()
    expression
    times[k] = toq();
  end
  println("Number of iterations: $maxiters")
  println("Units:                seconds          ")
  println("Best time:            $(minimum(times))")
  println("Worst time:           $(maximum(times))")
  println("Average time:         $(sum(times)/len(times))")
end

macro LStrings()
  include("/Users/terra/.julia/mymods/latex.jl")
end

macro MyPyPlotHelp()
  println("plt.title(titlestr, y = 1.05 )")
end

function Id( n )
  #= Makes a square Identity Matrix (n by n) =#
  return diagm(ones(n))
end

function symplecticMatrix( n )
  #= Makes a square Symplectic Block Matrix Omega (2n by 2n)
  #  Ω = [ 0 I_n; -I_n 0 ]
  =#
  I = Id( n )
  Z = zeros(n,n);
  Ω = [ Z I; -I Z ]
  return Ω
end

function symplecticTest( A, tol=1e-10 )
  #=
    symplecticTest( A[,tolerance=1e-10])

    Tests if matrix is Symplectic with the SBM Omega.
      A must be a square matrix of size 2n by 2n.
      Also, A'*Ω*A - Ω < tol = 1e-10 by default.
    Returns Boolean.
  =#
  M,N = size(A)
  if M != N
    println("ERROR: MATRIX NOT SQUARE")
    return false
  elseif mod(M,2) != 0
    println("ERROR: MATRIX NOT 2n BY 2n")
    return false
  end
  tol = 1e-10
  n   = int(M / 2)
  Ω   = symplecticMatrix( n )

  if maximum(maximum(abs(A' * Ω * A - Ω))) < tol
    return true
  end
  return false
end

function tayCoefMatrix( q, j=q )
  A = ones(q+1,j+1)
  J = collect(0:j)
  for k = 1:q+1
    p = k - 1
    f = factorial(p)
    A[k,:] = J.^p ./ f                             # Note 0^0 := 1
  end
  return A
end

function minor( A, j, k )
  s = size(A);
  if s[1] != s[2]
    print("Error -- Matrix is nonsquare")
    return None
  end
  small_size = (2,2);
  if s == small_size
    return A
  end
  t = (s[1]-1,s[2]-1)
  B = zeros(t)
  vc = 1
  for v = 1:s[2]
    if v != k
      uc = 1
      for u = 1:s[1]
        if u != j
          B[uc,vc] = A[u,v]
          uc += 1
        end
      end
    vc += 1
    end
  end
  return B
end

function diagDet( A )
  s = size(A)
  if s[1] != s[2]
    print("Error -- Matrix is nonsquare")
    return None
  end
  d = 0
  for k = 1:s[2]
    d += A[k,k] * det(minor(A,k,k))
  end
  return d
end

function welfertCheck( A )
  s = size(A)
  if s[1] != s[2]
    print("Error -- Matrix is nonsquare")
    return None
  end
  dd = diagDet(A)
  d  = det(A)
  return d - dd
end

function avg( a, slice=1 )
  return sum(avg,slice) / size(avg,slice)
end

function average( a, slice=1 )
  return avg(a, slice)
end

max_stars_repo_name: jyalim/myconf

max_issues_repo_name: jyalim/myconf

max_forks_repo_name: jyalim/myconfig


text:
@inline getbody(mechanism::Mechanism, id::Integer) = mechanism.bodies[id]
@inline getbody(mechanism::Mechanism, id::Nothing) = mechanism.origin
function getbody(mechanism::Mechanism, name::String)
    if mechanism.origin.name == name
        return mechanism.origin
    else
        for body in mechanism.bodies
            if body.name == name
                return body
            end
        end
    end
    return
end
@inline geteqconstraint(mechanism::Mechanism, id::Integer) = mechanism.eqconstraints[id]
function geteqconstraint(mechanism::Mechanism, name::String)
    for eqc in mechanism.eqconstraints
        if eqc.name == name
            return eqc
        end
    end
    return
end
@inline getineqconstraint(mechanism::Mechanism, id::Integer) = mechanism.ineqconstraints[id]
function getineqconstraint(mechanism::Mechanism, name::String)
    for ineqc in mechanism.eqconstraints
        if ineqc.name == name
            return ineqc
        end
    end
    return
end

getcomponent(mechanism::Mechanism, id::Nothing) = mechanism.origin
function getcomponent(mechanism::Mechanism, id::Integer)
    if haskey(mechanism.bodies, id)
        return getbody(mechanism, id)
    elseif haskey(mechanism.eqconstraints, id)
        return geteqconstraint(mechanism, id)
    elseif haskey(mechanism.ineqconstraints, id)
        return getineqconstraint(mechanism, id)
    else
        return
    end
end

function getcomponent(mechanism::Mechanism, name::String)
    component = getbody(mechanism,name)
    if component === nothing
        component = geteqconstraint(mechanism,name)
    end
    if component === nothing
        component = getineqconstraint(mechanism,name)
    end
    return component
end

@inline function normf(mechanism::Mechanism, body::Body)
    f = dynamics(mechanism, body)
    return dot(f, f)
end

@inline function normf(mechanism::Mechanism, eqc::EqualityConstraint)
    f = g(mechanism, eqc)
    return dot(f, f)
end

@inline function normf(mechanism::Mechanism, ineqc::InequalityConstraint)
    f = gs(mechanism, ineqc)
    d = h(ineqc)
    return dot(f, f) + dot(d, d)
end

@inline function normfμ(mechanism::Mechanism, ineqc::InequalityConstraint)
    f = gs(mechanism, ineqc)
    d = hμ(ineqc, mechanism.μ)
    return dot(f, f) + dot(d, d)
end

@inline function normf(mechanism::Mechanism)
    mechanism.normf = 0

    foreachactive(addNormf!, mechanism.bodies, mechanism)
    foreachactive(addNormf!, mechanism.eqconstraints, mechanism)
    foreachactive(addNormf!, mechanism.ineqconstraints, mechanism)

    return sqrt(mechanism.normf)
end

@inline function normf(mechanism::LinearMechanism)
    mechanism.normf = norm([fdynamics(mechanism); fconstraints(mechanism)])

    return sqrt(mechanism.normf)
end

@inline function meritf(mechanism::Mechanism)
    mechanism.normf = 0

    foreachactive(addNormf!, mechanism.bodies, mechanism)
    foreachactive(addNormf!, mechanism.eqconstraints, mechanism)
    foreachactive(addNormfμ!, mechanism.ineqconstraints, mechanism)

    return sqrt(mechanism.normf)
end

@inline function normΔs(mechanism::Mechanism)
    mechanism.normΔs = 0

    foreachactive(addNormΔs!, mechanism.bodies, mechanism)
    foreachactive(addNormΔs!, mechanism.eqconstraints, mechanism)
    foreachactive(addNormΔs!, mechanism.ineqconstraints, mechanism)

    return sqrt(mechanism.normΔs)
end

@inline function addNormf!(component::Component, mechanism::Mechanism)
    mechanism.normf += normf(mechanism, component)
    return
end

@inline function addNormfμ!(ineqc::InequalityConstraint, mechanism::Mechanism)
    mechanism.normf += normfμ(mechanism, ineqc)
    return
end

@inline function addNormΔs!(component::Component, mechanism::Mechanism)
    mechanism.normΔs += normΔs(component)
    return
end

@inline function discretizestate!(mechanism::Mechanism)
    foreach(discretizestate!, mechanism.bodies, mechanism.Δt)
    return
end

@inline function activate!(mechanism::Mechanism, id::Integer)
    component = getcomponent(mechanism, id)
    activate!(component)
    activate!(mechanism.graph,id)
    return
end

@inline function deactivate!(mechanism::Mechanism, id::Integer)
    component = getcomponent(mechanism, id)
    deactivate!(component)
    deactivate!(mechanism.graph,id)
    return
end

@inline function fdynamics(mechanism::LinearMechanism)
    A = mechanism.A
    Bu = mechanism.Bu
    Bλ = mechanism.Bλ
    z0 = mechanism.z
    z1 = mechanism.zsol[2]
    λ1 = mechanism.λsol[2]
    u = mechanism.u

    return A*z0 + Bu*u + Bλ*λ1 - z1
end

@inline function fconstraints(mechanism::LinearMechanism)
    G = mechanism.G
    z1 = mechanism.zsol[2]

    return G*z1
end

function disassemble(mechanism::Mechanism{T}) where T
    origin = mechanism.origin
    bodies = mechanism.bodies.values
    eqconstraints = mechanism.eqconstraints.values
    ineqconstraints = mechanism.ineqconstraints.values

    # Flip component ids
    for body in bodies
        body.id *= -1
    end
    for eqc in eqconstraints
        eqc.id *= -1
        if eqc.parentid !== nothing
            eqc.parentid *= -1
        end
        eqc.childids *= -1
    end
    for ineqc in ineqconstraints
        ineqc.id *= -1
        if ineqc.parentid !== nothing
            ineqc.parentid *= -1
        end
        ineqc.childids *= -1
    end

    # Set CURRENTID
    global CURRENTID = -1
    for body in bodies
        if body.id <= CURRENTID
            CURRENTID = body.id-1
        end
    end
    for eqc in eqconstraints
        if eqc.id <= CURRENTID
            CURRENTID = eqc.id-1
        end
    end
    for ineqc in ineqconstraints
        if ineqc.id <= CURRENTID
            CURRENTID = ineqc.id-1
        end
    end

    # Set origin to next id
    oldoid = origin.id
    origin.id = getGlobalID()
    for eqc in eqconstraints
        if eqc.parentid === nothing
            eqc.parentid = origin.id
        end
    end
    for ineqc in ineqconstraints
        if ineqc.parentid === nothing
            ineqc.parentid = origin.id
        end
    end

    return origin, bodies, eqconstraints, ineqconstraints
end


text:
using Econometrics, Statistics, Random

# version which generates shock internally
function SVmodel(θ, n, burnin)
    ϕ = θ[1]
    ρ = θ[2]
    σ = θ[3]
    hlag = 0.0
    ys = zeros(n)
    for t = 1:burnin+n
        # bound the log variance at a high value,
        # too bullet-proof things which use this
        h = min(ρ*hlag + σ*randn(),20.0)
        y = ϕ*exp(h/2.0)*randn()
        if t > burnin
            ys[t-burnin] = y
        end
        hlag = h
    end
    ys, ϕ*exp(hlag/2.0) # return the sample of returns, plus the final period volatility
end

# auxiliary model: HAR-RV
# Corsi, Fulvio. "A simple approximate long-memory model
# of realized volatility." Journal of Financial Econometrics 7,
# no. 2 (2009): 174-196.
function HAR(y)
    ylags = lags(y,10)
    X = [ones(size(y,1)) ylags[:,1]  mean(ylags[:,1:4],dims=2) mean(ylags[:,1:10],dims=2)]
    # drop missings
    y = y[11:end]
    X = X[11:end,:]
    βhat = X\y
    σhat = std(y-X*βhat)
    vcat(βhat,σhat)
end

function aux_stat(y)
    s = std(y)
    y = abs.(y)
    m = mean(y)
    s2 = std(y)
    k = std(y.^2.0)
    c = cor(y[1:end-1],y[2:end])
    # ratios of quantiles of moving averages to detect clustering
    q = try
        q = quantile(ma(y/s2,4), [0.1, 0.25, 0.75, 0.9])
    catch
        q = [1.0, 1.0, 0.0, 0.0]
    end
        c1 = q[4]/q[1]
        c2 = q[3]/q[2]
    vcat(m, s, s2, k, c, c1, c2, HAR(y))
    # four noise stats
    #vcat(m, s, s2, k, c, c1, c2, HAR(y), randn(4)./sqrt(size(y,1)))
end

# asymptotic Gaussian likelihood function of statistic
function logL(θ, m, n, η, ϵ, withdet=true)
    S = size(η,2)
    k = size(m,1)
    ms = zeros(S, k)
    Threads.@threads for s = 1:S
        y, junk = SVmodel(θ, n, η[:,s], ϵ[:,s])
        y = min.(y, 100.0)
        y = max.(y,-100.0)
        ms[s,:] = sqrt(n)*aux_stat(y)
    end
    mbar = mean(ms,dims=1)[:]
    if ~any(isnan.(mbar))
        Σ = cov(ms)
        x = (m .- mbar)
        lnL = try
            if withdet
                lnL = -0.5*log(det(Σ)) - 0.5*x'*inv(Σ)*x # for Bayesian
            else
                lnL = 0.5*x'*inv(Σ)*x # for classic indirect inference (note sign change)
            end
        catch
            lnL = -Inf
        end
     else
         lnL = -Inf
     end
     return lnL
end

# prior just checks that we're in the bounds
function prior(theta, lb, ub)
    a = 0.0
    if(all((theta .>= lb) .& (theta .<= ub)))
        a = 1.0
    end
    return a
end

# uniform random walk in one dimension
function proposal1(current, tuning, lb, ub)
    trial = copy(current)
    if rand() > 0.1
        i = rand(1:size(trial,1))
        trial[i] += tuning[i].*randn()
    else
        trial = lb + (ub - lb).*rand(size(lb,1))
    end
    return trial
end

# MVN random walk, or occasional draw from prior
function proposal2(current, cholV, lb, ub)
    trial = copy(current)
    if rand() > 0.1
        trial += cholV'*randn(size(trial))
    else
        trial = lb + (ub - lb).*rand(size(lb,1))
    end
    return trial
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

module GlassCat

using Polynomials
using Plots
using StringEncodings
using Unitful
using StaticArrays
using Base: @.
import Unitful: Length, Temperature, Quantity, Units
using Unitful.DefaultSymbols
using Pkg
using ForwardDiff

include("constants.jl")

include("GlassTypes.jl")
export GlassID, info, glassid, glassname, glassforid
include("Air.jl")
export Air, isair

# include built glass cat source files
@assert AGFGLASSCAT_PATH === joinpath(@__DIR__, "data", "jl", "AGFGlassCat.jl")
if !isfile(AGFGLASSCAT_PATH)
    @warn "$(basename(AGFGLASSCAT_PATH)) not found! Running build steps."
    Pkg.build("OpticSim"; verbose=true)
end
include("data/jl/AGFGlassCat.jl") # this needs to be literal for intellisense to work
include("data/jl/CARGILLE.jl")
export CARGILLE

# include functionality for managing runtime (dynamic) glass cats: MIL_GLASSES and MODEL_GLASSES
include("runtime.jl")
export glassfromMIL, modelglass

# include functions for searching the glass cats
include("search.jl")
export glasscatalogs, glassnames, findglass

include("utilities.jl")
export plot_indices, index, polyfit_indices, absairindex, absorption, drawglassmap

# include utility functions for maintaining the AGF source list
include("sources.jl")
export add_agf

# include build utility scripts to make testing them a bit easier
include("generate.jl")

end # module
export GlassCat


text:
## Testing
# These routines are for the unit tests

using Test

## Supports @test_approx_eq

Test.approx_full(f::Fun) = f

## Spaces Tests

function testtransforms(S::Space;minpoints=1,invertibletransform=true)
    # transform tests
    v = rand(max(minpoints,min(100,ApproxFun.dimension(S))))
    plan = plan_transform(S,v)
    @test transform(S,v)  == plan*v

    iplan = plan_itransform(S,v)
    @test itransform(S,v)  == iplan*v

    if invertibletransform
        for k=max(1,minpoints):min(5,dimension(S))
            v = [zeros(k-1);1.0]
            @test transform(S,itransform(S,v)) ≈ v
        end

        @test transform(S,itransform(S,v)) ≈ v
        @test itransform(S,transform(S,v)) ≈ v
    end
end

function testcalculus(S::Space;haslineintegral=true,hasintegral=true)
    for k=1:min(5,dimension(S))
        v = [zeros(k-1);1.0]
        f = Fun(S,v)
        @test abs(DefiniteIntegral()*f-sum(f)) < 100eps()
        if haslineintegral
            @test DefiniteLineIntegral()*f ≈ linesum(f)
        end
        @test norm(Derivative()*f-f') < 100eps()
        if hasintegral
            @test norm(differentiate(integrate(f))-f) < 100eps()
            @test norm(differentiate(cumsum(f))-f) < 200eps()
            @test norm(first(cumsum(f))) < 100eps()
        end
    end
end

function testmultiplication(spa,spb)
    for k=1:10
        a = Fun(spa,[zeros(k-1);1.])
        M = Multiplication(a,spb)
        pts = ApproxFun.checkpoints(rangespace(M))
        for j=1:10
            b = Fun(spb,[zeros(j-1);1.])
            @test (M*b).(pts) ≈ a.(pts).*b.(pts)
        end
    end
end

function testspace(S::Space;
                    minpoints=1,invertibletransform=true,haslineintegral=true,hasintegral=true,
                    dualspace=S)
    testtransforms(S;minpoints=minpoints,invertibletransform=invertibletransform)
    testcalculus(S;haslineintegral=haslineintegral,hasintegral=hasintegral)
    if dualspace ≠ nothing
        testmultiplication(dualspace,S)
    end
end

## Operator Tests

function backend_testfunctional(A)
    @test rowstart(A,1) ≥ 1
    @test colstop(A,1) ≤ 1
    @test bandwidth(A,1) ≤ 0
    @test blockbandwidth(A,1) ≤ 0

    B=A[1:10]
    @test eltype(B) == eltype(A)
    for k=1:5
        @test B[k] ≈ A[k]
        @test isa(A[k],eltype(A))
    end
    @test isa(A[1,1:10],Vector)
    @test isa(A[1:1,1:10],AbstractMatrix)
    @test B ≈ A[1,1:10]
    @test transpose(B) ≈ A[1:1,1:10]
    @test B[3:10] ≈ A[3:10]
    @test B ≈ [A[k] for k=1:10]

    co=cache(A)
    @test co[1:10] ≈ A[1:10]
    @test co[1:10] ≈ A[1:10]
    @test co[20:30] ≈ A[1:30][20:30] ≈ A[20:30]
end

# Check that the tests pass after conversion as well
function testfunctional(A::Operator{T}) where T<:Real
    backend_testfunctional(A)
    backend_testfunctional(Operator{Float64}(A))
    backend_testfunctional(Operator{Float32}(A))
    backend_testfunctional(Operator{ComplexF64}(A))
end

function testfunctional(A::Operator{T}) where T<:Complex
    backend_testfunctional(A)
    backend_testfunctional(Operator{ComplexF32}(A))
    backend_testfunctional(Operator{ComplexF64}(A))
end

function backend_testinfoperator(A)
    @test isinf(size(A,1))
    @test isinf(size(A,2))
    B=A[1:5,1:5]
    @test eltype(B) == eltype(A)

    for k=1:5,j=1:5
        @test B[k,j] ≈ A[k,j]
        @test isa(A[k,j],eltype(A))
    end

    @test A[1:5,1:5][2:5,1:5] ≈ A[2:5,1:5]
    @test A[1:5,2:5] ≈ A[1:5,1:5][:,2:end]
    @test A[1:10,1:10][5:10,5:10] ≈ [A[k,j] for k=5:10,j=5:10]
    @test A[1:10,1:10][5:10,5:10] ≈ A[5:10,5:10]
    @test A[1:30,1:30][20:30,20:30] ≈ A[20:30,20:30]

    @test Matrix(A[Block(1):Block(3),Block(1):Block(3)]) ≈ Matrix(A[blockstart(rangespace(A),1):blockstop(rangespace(A),3),blockstart(domainspace(A),1):blockstop(domainspace(A),3)])
    @test Matrix(A[Block(3):Block(4),Block(2):Block(4)]) ≈ Matrix(A[blockstart(rangespace(A),3):blockstop(rangespace(A),4),blockstart(domainspace(A),2):blockstop(domainspace(A),4)])

    for k=1:10
        @test isfinite(colstart(A,k)) && colstart(A,k) > 0
        @test isfinite(rowstart(A,k)) && colstart(A,k) > 0
    end

    co=cache(A)
    @test co[1:10,1:10] ≈ A[1:10,1:10]
    @test co[1:10,1:10] ≈ A[1:10,1:10]
    @test co[20:30,20:30] ≈ A[1:30,1:30][20:30,20:30]

    let C=cache(A)
        resizedata!(C,5,35)
        resizedata!(C,10,35)
        @test C.data[1:10,1:C.datasize[2]] ≈ A[1:10,1:C.datasize[2]]
    end
end

# Check that the tests pass after conversion as well
function testinfoperator(A::Operator{T}) where T<:Real
    backend_testinfoperator(A)
    backend_testinfoperator(convert(Operator{Float64}, A))
    backend_testinfoperator(convert(Operator{Float32}, A))
    backend_testinfoperator(convert(Operator{ComplexF64}, A))
end

function testinfoperator(A::Operator{T}) where T<:Complex
    backend_testinfoperator(A)
    backend_testinfoperator(convert(Operator{ComplexF32}, A))
    backend_testinfoperator(convert(Operator{ComplexF64}, A))
end

function testraggedbelowoperator(A)
    @test israggedbelow(A)
    for k=1:20
        @test isfinite(colstop(A,k))
    end

    R = RaggedMatrix(view(A, 1:10, 1:min(10,size(A,2))))
    for j=1:size(R,2)
        @test colstop(R,j) == min(colstop(A,j),10)
    end

    testinfoperator(A)
end

function testbandedbelowoperator(A)
    @test isbandedbelow(A)
    @test isfinite(bandwidth(A,1))
    testraggedbelowoperator(A)

    for k=1:10
        @test colstop(A,k) ≤ max(0,k + bandwidth(A,1))
    end
end

function testalmostbandedoperator(A)
    testbandedbelowoperator(A)
end

function testbandedoperator(A)
    @test isbanded(A)
    @test isfinite(bandwidth(A,2))
    testalmostbandedoperator(A)
    for k=1:10
        @test rowstop(A,k) ≤ k + bandwidth(A,2)
    end

    @test isa(A[1:10,1:10],BandedMatrix)
end

function testblockbandedoperator(A)
    @test isblockbanded(A)
    testraggedbelowoperator(A)
    @test isfinite(blockbandwidth(A,2))
    @test isfinite(blockbandwidth(A,1))

    if -blockbandwidth(A,1) ≤ blockbandwidth(A,2)
        for K=1:10
            @test K - blockbandwidth(A,2) ≤ blockcolstop(A,K).n[1] ≤ K + blockbandwidth(A,1) < ∞
            @test K - blockbandwidth(A,1) ≤ blockrowstop(A,K).n[1] ≤ K + blockbandwidth(A,2) < ∞
        end
    end
end

function testbandedblockbandedoperator(A)
    @test isbandedblockbanded(A)
    testblockbandedoperator(A)
    @test isfinite(subblockbandwidth(A,1))
    @test isfinite(subblockbandwidth(A,2))

    @test isa(A[Block.(1:4),Block.(1:4)], BandedBlockBandedMatrix)
end


text:
struct NeuralNetRelaxation <: AbstractRelaxation
    num_layers::Int
end

function reformulate_unary_quadratic_term!(
    model::MOI.ModelLike,
    reformulation::Reformulation{NeuralNetRelaxation},
    _x::VI,
    _y::VI,
    l::Float64,
    u::Float64,
)
    x, y = _map_to_interval_if_needed(model, _x, _y, l, u)
    # Secant cut: y <= x
    MOI.add_constraint(model, 1.0SV(x) - SV(y), MOI.GreaterThan(0.0))
    return _formulate_zero_one_quadratic!(
        model,
        x,
        y,
        reformulation.relaxation.num_layers,
        !reformulation.soc_lower_bound,
    )
end

function _formulate_zero_one_quadratic!(
    model::MOI.ModelLike,
    x::VI,
    y::VI,
    num_layers::Int,
    impose_lower_bound::Bool,
)
    z_v, z_c = MOI.add_constrained_variables(
        model,
        [MOI.ZeroOne() for i in 1:num_layers],
    )
    g_v, g_c = MOI.add_constrained_variables(
        model,
        [MOI.Interval(0.0, 1.0) for i in 1:num_layers],
    )

    gv = Dict(i => SV(g_v[i]) for i in 1:num_layers)
    gv[0] = SV(x)

    agg_aff = 1.0 * SV(y) - SV(x) + sum(gv[s] / 2.0^(2s) for s in 1:num_layers)

    if impose_lower_bound
        agg_c = MOI.add_constraint(model, agg_aff, MOI.EqualTo(0.0))
    else
        agg_c = MOI.add_constraint(model, agg_aff, MOI.LessThan(0.0))
    end

    for i in 1:num_layers
        α_i = SV(z_v[i])
        MOI.add_constraint(model, gv[i] - 2.0gv[i-1], MOI.LessThan(0.0))
        MOIU.normalize_and_add_constraint(
            model,
            gv[i] + 2.0gv[i-1] - 2.0,
            MOI.LessThan(0.0),
        )
        MOI.add_constraint(
            model,
            gv[i] + 2.0gv[i-1] - 2.0 * α_i,
            MOI.GreaterThan(0.0),
        )
        MOIU.normalize_and_add_constraint(
            model,
            gv[i] - 2.0gv[i-1] + 2.0 * α_i,
            MOI.GreaterThan(0.0),
        )
    end
    return z_v
end


text:
# Defining the fitness function
# It calls fopa.jl which returns the gain profile and some performance parameters

# The input of fitness is vector with first element corrsponding to L, second element corresponding to λ_p1, and third element corresponding to λ_p2
# Output of the fitness function is defined as weighted sum of bandwidth (BW), peak gain (Gpeak) and gain ripply (G_var)
function fitness(A,fnc)
    if fnc == "ackley"
        x = A[1];
        y = A[2];
        f = -20*exp(-0.2* ( 0.5*(x^2 + y^2) )^0.5 ) - exp( 0.5*( cos(2*pi*x) + cos(2*pi*y) )) + exp(1) + 20;
        f = -f;
    end

    if fnc == "sphere"
        f = sum(A.^2);
        f = -f;
    end

    if fnc == "rosenbrock"
        D = length(A);
        sum = 0;
        for m = 1:D-1
            sum = sum + 100*(A[m+1]-A[m]^2)^2 + (1-A[m])^2;
        end
        f = -sum;
    end

    if fnc == "rastrigin"
        D = length(A);
        sum = 0;
        for m = 1:D
            sum = sum + A[m]^2 - 10*cos(2*pi*A[m])
        end
        f = -1*(sum + 10*D);
    end

    if fnc == "double_gaussian"
        x0 = 0.1;
        y0 = -0.2;
        z0 = 0.3;

        f1 = exp(-(A[1]-x0).^2 -(A[2]-y0).^2 -(A[3]-z0).^2);

        x0 = -0.2;
        y0 = 0.2;
        z0 = -0.1;

        f2 = exp(-(A[1]-x0).^2 -(A[2]-y0).^2 -(A[3]-z0).^2);
        f =f1 + 0.6*f2;
    end
    return f
end

# A = [.1 -.2 .3];     # test input

# println("Testing the fitness function:")
# println("Test input = ", A)
# println("Fitness value = ", fitness(A))

max_stars_repo_path: Simulated Annealing code/fitness.jl

max_issues_repo_path: Simulated Annealing code/fitness.jl

max_forks_repo_path: Simulated Annealing code/fitness.jl


text:
##########################################################
# Rafał Nowak
#
# Instytut Informatyki
# Uniwersytet Wrocławski
#
#
# Pomocniczne materiały do wykładu z analizy numerycznej
#
##########################################################
# Wyganae pakiety:
# sudo apt-get install python python-matplotlib python3-matplotlib
##############################

Pkg.add("PyPlot")
using PyPlot

f(x) = ( 1.0 - cos(x) ) ./ x.^2

x = linspace(0,2*pi,1000);
y = f(x);
PyPlot.plot(x, y, color="red", linewidth=2.0, linestyle="--")
PyPlot.show()

for x = [ 4.0 3.0 2.0 1.50 1.30 1.15 1.10 1.05 1.01 1.00 ]*1e-8
  y = f(x);
  @printf("f(%.2e) = %.2e\n", x, y);
end

max_stars_repo_name: Magikis/Uniwersity

max_issues_repo_name: Magikis/Uniwersity

max_forks_repo_name: Magikis/Uniwersity


text:
import HttpServer.Client
import HttpServer.HttpHandler
import HttpServer.Server
import HttpServer.ClientParser
import HttpServer.HttpParser
import HttpServer.Request  # Not sure if necessary, but otherwise dispatching on HttpCommon.Request can fail.
import HttpServer.Response
import HttpCommon.Cookie
import HttpCommon.STATUS_CODES
import URIParser.URI
export printstartinfo
show(io::IO, client::Client) = directto_abstractdevice(io, client)
function _show(d::AbstractDevice, client::Client)
    _log(d, typeof(client),  " id ", :bold, client.id, :normal)
    _log(d, " ", client.sock, "\n")
    _log(d, "\t\t\t", client.parser)
    nothing
end

# TODO is this good enough?
function show(io::IO, p::ClientParser)
    print(io, "HttpServer.ClientParser.HttpParser.libhttp-parser: v", HttpParser.version())
    if p.parser.http_major > 0
        print(io, " HTTP/", p.parser.http_major, ".", p.parser.http_minor)
    end
    nothing
end

"Response already has a decent show method, we're not overwriting that.
This metod is called only when logging to an Abstractdevice"
function _show(d::AbstractDevice, response::HttpServer.Response)
    _log(d, :green, "Response status: ", :bold, response.status," ")
    _log(d, get(STATUS_CODES, response.status, "--"), " ")
    if !isempty(response.headers)
        _log(d, :green, " Headers: ", :bold, response.headers.count)
        _log(d, :green, "\n", response.headers)
    end
    if !isempty(response.cookies)
        _log(d, " Cookies: ", :bold, length(response.cookies))
        _log(d, "\n", response.cookies)
    end
    if !isempty(response.data)
        _log(d, "\t", DataDispatch(response.data, get(response.headers, "Content-Type","---")))
    end
    nothing
end

"Request already has a decent show method, we're not overwriting that.
This metod is called only when logging to an Abstractdevice"
function _show(d::AbstractDevice, request::Request)
    _log(d,  :normal, :light_yellow, "Request ", :normal)
    _log(d, :bold, request.method, " ", :cyan, request.resource, "\n", :normal)
    if !isempty(request.data)
        _log(d, "\t", DataDispatch(request.data, get(request.headers, "Content-Type", "text/html; charset=utf-8")))
    end
    if request.uri != URI("")
        _log(d, :cyan, "\tUri:", :bold, _string(request.uri), :normal, "\n\t")
    end
    if !isempty(request.headers)
        _log(d, "\t", :cyan, " .headers: ", request.headers.count)
        _log(d, :cyan, "\n", request.headers)
    end
    nothing
end

"HttpServer does not define a show method for its server type. Defining this is not piracy."
show(io::IO, server::Server) = directto_abstractdevice(io, server)
function _show(d::AbstractDevice, server::Server)
    _log(d, :bold , :green,  typeof(server), "(\n", :normal)
    server.http != nothing  && _log(d, "\t", server.http)
    server.websock != nothing  && _log(d, "\t", server.websock)
    _log(d, :bold, :green, ")")
    nothing
end

"HttpServer does not define a show method for its HttpHandler type. Defining this is not piracy."
show(io::IO, httphandler::HttpHandler) = directto_abstractdevice(io, httphandler)
function _show(d::AbstractDevice, httphandler::HttpHandler)
    _log(d,  :bold, Base.info_color(), typeof(httphandler), "( " , :normal)
    _log(d,  ".handle: ", :blue, :bold, httphandler.handle, :normal, "\n")
    _log(d, "\t\t\t.events:\n")
    _log(d,  httphandler.events)
    _log(d,  "\t\t\t", ".socket:\t", httphandler.sock, :bold, Base.info_color(), ")\n")
    nothing
end

nothing


text:
include("structs.jl")
include("utils.jl")

using .Structs: Query
using .SimpleDist: dist
using .SimpleDist.MidDist: mid_dist

a = Query("a", 1)
b = Query("b", -1)

println(dist(a, b))
println(mid_dist(a, b))


text:
function hcat(x::Vector{T}) where T<:AbstractVecOrMat
    l    = length(x)
    if l == 0
        throw(ArgumentError("cannot flatten empty vector"))
    else
        x1   = x[1]
        m, n = size(x1, 1), size(x1, 2)
        B    = similar(x1, eltype(x1), (m, l*n))
        for i = 1:l
            B[:, (i - 1)*n + 1:i*n] = x[i]
        end
        return B
    end
end

if VERSION < v"1.3.0-alpha.115"
mul!(y::StridedVector{T},
     A::StridedMatrix{T},
     x::StridedVector{T},
     α::Number,
     β::Number) where {T<:BlasFloat} =
    BLAS.gemv!('N', convert(T, α), A, x, convert(T, β), y)

mul!(y::StridedVector{T},
     A::Adjoint{T,<:StridedMatrix{T}},
     x::StridedVector{T},
     α::Number,
     β::Number) where {T<:BlasReal} =
    BLAS.gemv!('T', convert(T, α), parent(A), x, convert(T, β), y)

mul!(y::StridedVector{T},
     A::Adjoint{T,<:StridedMatrix{T}},
     x::StridedVector{T},
     α::Number,
     β::Number) where {T<:BlasComplex} =
    BLAS.gemv!('C', convert(T, α), parent(A), x, convert(T, β), y)

function mul!(y::StridedVector, A::StridedMatrix, x::StridedVector, α::Number, β::Number)
    n = length(y)
    for i = 1:n
        y[i] *= β
    end
    for i = 1:n
        for l = 1:size(A,2)
            y[i] += α*A[i,l]*x[l]
        end
    end
    return y
end

function mul!(y::StridedVector, A::Adjoint{<:StridedMatrix}, x::StridedVector, α::Number, β::Number)
    n = length(y)
    for i = 1:n
        y[i] *= β
    end
    for i = 1:n
        for l = 1:size(A,2)
            y[i] += α*A[l,i]'*x[l]
        end
    end
    return y
end
end

function qr!(x::AbstractArray)
    nm = norm(x)
    rmul!(x, inv(nm))
    return x, nm
end


text:
mutable struct ListItem
    next::ListItem
    prev::ListItem
    function ListItem()
        x = new()
        x.next = x
        x.prev = x
        x
    end
    function ListItem(n::ListItem, p::ListItem)
        x=new()
        x.next = x
        x.prev = p
        x
    end

end

a = ListItem()
b = ListItem(b, a)
println(b.next == b)
println(b.prev == a)

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
export RectangularBinning

import DelayEmbeddings: Dataset, minima, maxima
import StaticArrays: SVector, MVector

"""
    BinningScheme

The supertype of all binning schemes in the CausalityTools ecosystem.
"""
abstract type BinningScheme end

"""
    RectangularBinningScheme <: BinningScheme

The supertype of all rectangular binning schemes in the CausalityTools ecosystem.
"""
abstract type RectangularBinningScheme <: BinningScheme end

"""
    RectangularBinning(ϵ) <: RectangularBinningScheme

Instructions for creating a rectangular box partition using the binning scheme `ϵ`.
Binning instructions are deduced from the type of `ϵ`.

Rectangular binnings may be automatically adjusted to the data in which the `RectangularBinning`
is applied, as follows:

1. `ϵ::Int` divides each coordinate axis into `ϵ` equal-length intervals,
    extending the upper bound 1/100th of a bin size to ensure all points are covered.

2. `ϵ::Float64` divides each coordinate axis into intervals of fixed size `ϵ`, starting
    from the axis minima until the data is completely covered by boxes.

3. `ϵ::Vector{Int}` divides the i-th coordinate axis into `ϵ[i]` equal-length
    intervals, extending the upper bound 1/100th of a bin size to ensure all points are
    covered.

4. `ϵ::Vector{Float64}` divides the i-th coordinate axis into intervals of fixed size `ϵ[i]`, starting
    from the axis minima until the data is completely covered by boxes.

Rectangular binnings may also be specified on arbitrary min-max ranges.

5. `ϵ::Tuple{Vector{Tuple{Float64,Float64}},Int64}` creates intervals
    along each coordinate axis from ranges indicated by a vector of `(min, max)` tuples, then divides
    each coordinate axis into an integer number of equal-length intervals. *Note: this does not ensure
    that all points are covered by the data (points outside the binning are ignored)*.

# Example 1: Grid deduced automatically from data (partition guaranteed to cover data points)

## Flexible box sizes

The following binning specification finds the minima/maxima along each coordinate axis, then
split each of those data ranges (with some tiny padding on the edges) into `10` equal-length
intervals. This gives (hyper-)rectangular boxes, and works for data of any dimension.

```julia
using Entropies
RectangularBinning(10)
```

Now, assume the data consists of 2-dimensional points, and that we want a finer grid along
one of the dimensions than over the other dimension.

The following binning specification finds the minima/maxima along each coordinate axis, then
splits the range along the first coordinate axis (with some tiny padding on the edges)
into `10` equal-length intervals, and the range along the second coordinate axis (with some
tiny padding on the edges) into `5` equal-length intervals.
This gives (hyper-)rectangular boxes.

```julia
using Entropies
RectangularBinning([10, 5])
```

## Fixed box sizes

The following binning specification finds the minima/maxima along each coordinate axis,
then split the axis ranges into equal-length intervals of fixed size `0.5` until the all data
points are covered by boxes. This approach yields (hyper-)cubic boxes, and works for
data of any dimension.

```julia
using Entropies
RectangularBinning(0.5)
```

Again, assume the data consists of 2-dimensional points, and that we want a finer grid along
one of the dimensions than over the other dimension.

The following binning specification finds the minima/maxima along each coordinate axis, then splits
the range along the first coordinate axis into equal-length intervals of size `0.3`,
and the range along the second axis into equal-length intervals of size `0.1` (in both cases,
making sure the data are completely covered by the boxes).
This approach gives a (hyper-)rectangular boxes.

```julia
using Entropies
RectangularBinning([0.3, 0.1])
```

# Example 2: Custom grids (partition not guaranteed to cover data points):

Assume the data consists of 3-dimensional points `(x, y, z)`, and that we want a grid
that is fixed over the intervals `[x₁, x₂]` for the first dimension, over `[y₁, y₂]`
for the second dimension, and over `[z₁, z₂]` for the third dimension. We when want to
split each of those ranges into 4 equal-length pieces. *Beware: some points may fall
outside the partition if the intervals are not chosen properly (these points are
simply discarded)*.

The following binning specification produces the desired (hyper-)rectangular boxes.

```julia
using Entropies, DelayEmbeddings

D = Dataset(rand(100, 3));

x₁, x₂ = 0.5, 1 # not completely covering the data, which are on [0, 1]
y₁, y₂ = -2, 1.5 # covering the data, which are on [0, 1]
z₁, z₂ = 0, 0.5 # not completely covering the data, which are on [0, 1]

ϵ = [(x₁, x₂), (y₁, y₂), (z₁, z₂)], 4 # [interval 1, interval 2, ...], n_subdivisions

RectangularBinning(ϵ)
```
"""
struct RectangularBinning{E} <: RectangularBinningScheme
    ϵ::E
end

"""
    get_minima(pts) → SVector

Return the minima along each axis of the dataset `pts`.
"""
function get_minima end

"""
    get_maxima(pts) → SVector

Return the maxima along each axis of the dataset `pts`.
"""
function get_maxima end

"""
    get_minmaxes(pts) → Tuple{Vector{Float}, Vector{Float}}

Return a vector of tuples containing axis-wise (minimum, maximum) values.
"""
function get_minmaxes end

function get_minima(pts::AbstractDataset)
    minima(pts)
end

function get_minima(pts::Vector{T}) where {T <: Union{SVector, MVector, Vector}}
    minima(Dataset(pts))
end

function get_maxima(pts::AbstractDataset)
    maxima(pts)
end

function get_maxima(pts::Vector{T}) where {T <: Union{SVector, MVector, Vector}}
    maxima(Dataset(pts))
end

function get_minmaxes(pts::AbstractDataset)
    mini, maxi = minima(pts), maxima(pts)
    minmaxes = [(mini[i], maxi[i]) for i = 1:length(mini)]
end

function get_minmaxes(pts::Vector{T}) where {T <: Union{SVector, MVector, Vector}}
    get_minmaxes(Dataset(pts))
end

"""
    get_minima_and_edgelengths(points,
        binning_scheme::RectangularBinning) → (Vector{Float}, Vector{Float})

Find the minima along each axis of the embedding, and computes appropriate
edge lengths given a rectangular `binning_scheme`, which provide instructions on how to
grid the space. Assumes the input is a vector of points.

See documentation for [`RectangularBinning`](@ref) for details on the
binning scheme.

# Example

## Example 1: Grid deduced automatically from data (partition guaranteed to cover data points)

```julia
using Entropies, DelayEmbeddings
pts = Dataset([rand(4) for i = 1:1000])

get_minima_and_edgelengths(pts, RectangularBinning(0.6))
get_minima_and_edgelengths(pts, RectangularBinning([0.5, 0.3, 0.4, 0.4]))
get_minima_and_edgelengths(pts, RectangularBinning(10))
get_minima_and_edgelengths(pts, RectangularBinning([10, 8, 5, 4]))
```

## Example 2: Custom grids (partition not guaranteed to cover data points):

```julia
using Entropies, DelayEmbeddings

# Here, we split the first, second and third coordinate axes in
# five equal-length pieces, defined over the intervals [a, b],
# [c, d], and [e, f], respectively (beware: some points may fall
# outside the partition if the intervals are not chosen properly).
D = Dataset(rand(100, 3));
a, b = 0.5, 1 # not completely covering the data, which are on [0, 1]
c, d = -2, 1.5 # covering the data, which are on [0, 1]
e, f = 0, 0.5 # not completely covering the data, which are on [0, 1]
ϵ = [(a, b), (c, d), (e, f)], 3
bin = RectangularBinning(ϵ)
get_minima_and_edgelengths(D, bin)
```
"""
function get_minima_and_edgelengths(points, binning_scheme::RectangularBinning)
    ϵ = binning_scheme.ϵ

    D = length(points[1])
    n_pts = length(points)

    axisminima = minimum.([minimum.([pt[i] for pt in points]) for i = 1:D])
    axismaxima = maximum.([maximum.([pt[i] for pt in points]) for i = 1:D])

    edgelengths = Vector{Float64}(undef, D)

    # Dictated by data ranges
    if ϵ isa Float64
        edgelengths = [ϵ for i in 1:D]
    elseif ϵ isa Vector{<:AbstractFloat}
        edgelengths .= ϵ
    elseif ϵ isa Int
        edgeslengths_nonadjusted = (axismaxima  - axisminima) / ϵ
        edgelengths = ((axismaxima + (edgeslengths_nonadjusted ./ 100)) - axisminima) ./ ϵ
    elseif ϵ isa Vector{Int}
        edgeslengths_nonadjusted = (axismaxima  .- axisminima) ./ ϵ
        edgelengths = ((axismaxima .+ (edgeslengths_nonadjusted ./ 100)) .- axisminima) ./ ϵ

    # Custom data ranges
    elseif ϵ isa Tuple{Vector{Tuple{T1, T2}},Int64} where {T1 <: Real, T2 <: Real}
        ranges = ϵ[1]
        n_bins = ϵ[2]
        length(ranges) == D || error("Tried to apply a $(length(ranges))-dimensional binning scheme to $(D)-dimensional data. Please provide ranges for all $(D) dimensions.")

        # We have predefined axis minima and axis maxima.
        stepsizes = zeros(Float64, D)
        edgelengths = zeros(Float64, D)
        for i = 1:D
            ranges[i][1] < ranges[i][2] || error("Range along dimension $i is not given as (minimum, maximum), got $(ranges[i])")
            edgelengths[i] = (maximum(ranges[i]) - minimum(ranges[i]))/n_bins
            axisminima[i] = minimum(ranges[i])
        end
    else
       error("get_minima_and_edgelengths not implemented for RectangularBinning(ϵ) with ϵ of type $(typeof(ϵ))")
    end

    axisminima, edgelengths
end

"""
    get_edgelengths(pts, binning_scheme::RectangularBinning) → Vector{Float}

Return the box edge length along each axis resulting from discretizing `pts` on a
rectangular grid specified by `binning_scheme`.

# Example

```julia
using Entropies, DelayEmbeddings
pts = Dataset([rand(5) for i = 1:1000])

get_edgelengths(pts, RectangularBinning(0.6))
get_edgelengths(pts, RectangularBinning([0.5, 0.3, 0.3, 0.4, 0.4]))
get_edgelengths(pts, RectangularBinning(8))
get_edgelengths(pts, RectangularBinning([10, 8, 5, 4, 22]))
```
"""
function get_edgelengths end

get_edgelengths(points, ϵ::RectangularBinning) = get_minima_and_edgelengths(points, ϵ)[2]

# need this for type stability when getting minima and edgelength for a Dataset for _non0hist
# otherwise, we're subtracting and dividing static vectors with regular vectors, which allocates
# all over the place.
import StaticArrays: SVector

function minima_edgelengths(points::AbstractDataset{D, T}, binning_scheme::RectangularBinning) where {D, T<:Real}
    mini, edgelengths = get_minima_and_edgelengths(points, binning_scheme)
    return SVector{D, T}(mini), SVector{D, Float64}(float.(edgelengths))
end
get_edgelengths(points::AbstractDataset, ϵ::RectangularBinning) = minima_edgelengths(points, ϵ)[2]


text:
@testset "AbsAlgAss" begin

  @testset "Decomposition" begin
    Fp = GF(3)
    G = small_group(8, 4)
    FpG = group_algebra(Fp, G)

    dec = decompose(FpG)
    @test length(dec) == 5
    dim1 = 0
    dim4 = 0
    for (A, toFpG) in dec
      if dim(A) == 1
        dim1 += 1
      elseif dim(A) == 4
        dim4 += 1
      end
    end
    @test dim1 == 4
    @test dim4 == 1

    # Check whether the maps "work"
    A, toFpG = dec[1]
    e = toFpG(one(A))
    @test e^2 == e
    ee = toFpG\one(FpG)
    @test ee^2 == ee

    # And now the same for AlgAss
    FpG = AlgAss(FpG)[1]
    dec = decompose(FpG)
    @test length(dec) == 5
    dim1 = 0
    dim4 = 0
    for (A, toFpG) in dec
      if dim(A) == 1
        dim1 += 1
      elseif dim(A) == 4
        dim4 += 1
      end
    end
    @test dim1 == 4
    @test dim4 == 1

    # Check whether the maps "work"
    A, toFpG = dec[1]
    e = toFpG(one(A))
    @test e^2 == e
    ee = toFpG\one(FpG)
    @test ee^2 == ee

    # And now for AlgMat
    A = matrix_algebra(Fp, 2)

    dec = decompose(A)
    @test length(dec) == 1

    B, BtoA = dec[1]
    @test dim(B) == 4

    e = BtoA(one(B))
    @test e^2 == e
    ee = BtoA\one(A)
    @test ee^2 == ee

    A = matrix_algebra(Fp, [ matrix(Fp, [ 1 1; 0 1 ]) ]) # not semisimple!
    @test_throws AssertionError decompose(A)

    Qx, x = FlintQQ["x"]
    A = AlgAss((x^2 + 1)*(x^2 + 3))
    dec = Hecke.as_number_fields(A)

    @test length(dec) == 2
    @test degree(dec[1][1]) == 2
    @test degree(dec[2][1]) == 2

    K, AtoK = dec[1]
    e = AtoK(one(A))
    @test e^2 == e
    ee = AtoK\one(K)
    @test ee^2 == ee
  end

  @testset "Generators" begin
    Qx, x = FlintQQ["x"]
    A = AlgAss((x^2 + 1)*(x^2 + 3))
    g, full_basis, v = gens(A, Val{true})

    @test length(full_basis) == dim(A)

    M = zero_matrix(FlintQQ, dim(A), dim(A))
    for i = 1:dim(A)
      Hecke.elem_to_mat_row!(M, i, full_basis[i])
    end
    @test rank(M) == dim(A)

    for i = 1:dim(A)
      b = full_basis[i]
      a = one(A)
      for (j, k) in v[i]
        a *= g[j]^k
      end
      @test b == a
    end
  end

end


text:
using GLVisualize, GeometryTypes, Reactive

if !isdefined(:runtests)
    window = glscreen()
    timesignal = bounce(linspace(0,1,360))
end
N = 20
# generate some rotations
function rotation_func(t)
    t = (t == 0f0 ? 0.01f0 : t)
    Vec2f0[(sin(x/t), cos(y/(t/2f0))) for x=1:N, y=1:N]
end

# us Reactive.map to transform the timesignal signal into the arrow flow
flow = map(rotation_func, timesignal)

# create a visualisation
vis = visualize(flow, ranges=(50:800,50:500))
_view(vis, window, camera=:orthographic_pixel)

if !isdefined(:runtests)
    renderloop(window)
end


text:
using MLStyle

flatten_row(expr) =
    @match expr begin
        Expr(:row, args...) => args
        a => [a]
    end

lisp(expr) =
    @match expr begin
        :[] => nothing

        :[f] => Expr(:call, lisp(f))

        Expr(:vcat, args...) =>
            lisp(Expr(:hcat, vcat(map(flatten_row, args)...)...))

        :[($(hd::QuoteNode)) ($(args...))] =>
            Expr(hd.value, map(lisp, args)...)

        :[($hd) ($(args...))] =>
            Expr(:call, hd, map(lisp, args)...)

        Expr(head, args...) =>
            Expr(head, map(lisp, args)...)

        a => a
    end

macro lisp(expr)
    lisp(expr) |> esc
end

max_stars_repo_name: thautwarm/MLStyle-Playground

max_issues_repo_name: thautwarm/MLStyle-Playground

max_forks_repo_name: thautwarm/MLStyle-Playground


text:
module AbstractDomains

import Base: convert, promote_rule
import Base: string, print, show, rand
import Base: abs, zero, one, in, inv, ndims, issubset, union, intersect, isequal

# FIXME: ifelse is built-in. New methods are not allowed
# import Base: ifelse

import Base: !
import Base: ==
import Base: ==
import Base: ==
import Base: ==
import Base: |
import Base: &
import Base: ==
import Base: !=
import Base: ==
import Base: ==
import Base: >
import Base: <
import Base: <
import Base: <
import Base: <=
import Base: >=
import Base: +
import Base: -
import Base: *
import Base: /
import Base: //
import Base: getindex
import Base: setindex!

export Interval,
       Domain,
       AbstractBool,
       Boxes, HyperBox,LazyBox,
       LazyRandomVector,
       t,f,tf,
       ⊔, ⊓, subsumes, isintersect, intersect, isrelational, domaineq,
       ndims,
       getindex,
       mid, mid_split, partial_split_box, mid_partial_split,
       unit,
       dims,
       hasdim,
       isrelational

include("common.jl")
include("domains.jl")

end


text:
type Distribution{FS<:FeatureSpace}
  space::FS
  features::Number
  total::Number
  smooth::Function
  smooth_data::Array
  mdata::Any
  Distribution() = new(FS(),0,0,_no_smoothing,[])
  Distribution(fv::FeatureVector) = new(fv,length(fv),get_total(fv),_no_smoothing,[])
  Distribution(c::Cluster) = new(c,length(c.vector_sum),get_total(c.vector_sum) ,_no_smoothing,[])
  Distribution(ds::DataSet) = new(ds,length(ds.vector_sum),get_total(ds.vector_sum) ,_no_smoothing,[])
  function get_total(fv::FeatureVector)
    total = 0
    for value in values(fv)
      total += value
    end
    return total
  end
end
Distribution(fv::FeatureVector) = Distribution{FeatureVector}(fv)
Distribution(c::Cluster) = Distribution{Cluster}(c)
Distribution(ds::DataSet) = Distribution{DataSet}(ds)

function getindex(d::Distribution, key)
  return d.space[key]
end

function probability(d::Distribution, feature)
  return d.smooth(d, feature, d.smooth_data)
end

function keys(d::Distribution)
  return keys(d.space)
end

function features(d::Distribution{FeatureVector})
  return keys(d.space)
end

function features(d::Distribution)
  return keys(d.space.vector_sum)
end

function isempty(d::Distribution)
  return isempty(d.space)
end

function entropy(d::Distribution)
  ent = 0
  for feature in features(d)
    ent -= probability(d,feature)*log2(probability(d,feature))
  end
  return ent
end

function info_gain(d1::Distribution, d2::Distribution)
  return entropy(d1)-entropy(d2)
end

function perplexity(d::Distribution)
  return 2^entropy(d)
end

function set_smooth!(d::Distribution{FeatureVector}, f::Function, sd::Array)
  d.smooth = f
  d.smooth_data = sd
end

#no smoothing default
function remove_smoothing!(d::Distribution)
  set_smooth!(d,_no_smoothing,[])
end

function _no_smoothing(d::Distribution{FeatureVector}, key, data::Array)
  return d.space[key] / d.total
end

function _no_smoothing(d::Distribution, feature, data::Array)
  return d.space.vector_sum[feature] / d.total
end

#add-delta smoothing, default to add-one smoothing
function delta_smoothing!(d::Distribution, δ::Number=1)
  if δ <= 0
    Base.warn("δ must be greater than 0")
  end
  set_smooth!(d,_δ_smoothing,[δ,d.features,d.total])
end

function _δ_smoothing(d::Distribution, key, data::Array)
  if !haskey(d.space, key)
    return (data[1])/(data[1]*(data[2]+1)+data[3])
  end
  return (d.space[key]+data[1])/(data[1]*(data[2]+1)+data[3])
end

function display(dist::Distribution)
  display(dist.space)
end

max_stars_repo_name: drawkward/SLU-TextMining

max_issues_repo_name: drawkward/SLU-TextMining

max_forks_repo_name: drawkward/SLU-TextMining


text:
@doc """
       QuasiNewton method for local optimization (single variable)

       Supported by:
       Ph D Kelvyn B. Sánchez
       e-mail: kelvyn.baruc@gmail.com
 """ ->
function QuasiNewton(f, var; x0 = [1 2], error = 0.0001, maxiter = 100, iter = true)
      try
            # Symbolic derivatives using Module SymPy
            dfk1 = diff(f, var)

            err = Inf
            i = 1
            xk = 0
            x_num = 0

            # Numeric iterations
            iter == true ? ti = time_ns() : nothing # To measure algorothm performance

            while err > error
                  x_num = x0

                  xk = x_num[2]-(N(dfk1 |> subs(var,x_num[2]))/(N(dfk1 |> subs(var,x_num[2])) -
                        N(dfk1 |> subs(var,x_num[1]))))*(x_num[2] - x_num[1])

                  if xk > x_num[2]
                        err = abs(x_num[1] - xk)
                        x0 = [x_num[2] xk]

                  elseif xk < x_num[1]
                        err = abs(x_num[2] - xk)
                        x0 = [xk x_num[1]]

                  elseif xk < x_num[2] && xk > x_num[1]
                         err = abs(x_num[1] - xk)
                         x0 = [xk x_num[2]]
                  end

                        if iter == true
                              if i == 1
                                    println("===========================================================")
                                    println(" QuasiNewton Solver for Unconstrained Optimization (0.0.2)")
                              end

                              if i == 1 || ceil((i/1)/10) == (i/1)/10
                                    println("===========================================================")
                                    @printf("Iter\t X\t\t f(x)\t\t |err|\n")
                                    println("===========================================================")
                              end

                              @printf("%i\t %4.6e\t %4.6e\t %4.6e\n", i, xk, N(f |> subs(var,xk)), err)
                        end

                  if i == maxiter
                        throw(ErrorException("Not Solved!. Iteration limit reached."))
                  end

                  i += 1
            end

            println("===========================================================")
            println(" ")

            iter == true ? tf = time_ns() : nothing # To measure algorithm performance
            iter == true ? println("elapse time: ", (tf - ti)/1.0e9, " seconds") : nothing

            return xk,  N(f |> subs(var,xk))

      catch Msg
            println("===========================================================")
            println(Msg)
            return NaN, NaN
      end
end


text:
using Mixins
using Documenter

DocMeta.setdocmeta!(Mixins, :DocTestSetup, :(using Mixins); recursive=true)

makedocs(;
    modules=[Mixins],
    authors="Owen Lynch <root@owenlynch.org> and contributors",
    repo="https://github.com/olynch/Mixins.jl/blob/{commit}{path}#{line}",
    sitename="Mixins.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)


text:
using ModelingToolkit, StaticArrays, LinearAlgebra, LabelledArrays
using DiffEqBase, ForwardDiff
using Test

# Define some variables
@parameters t σ ρ β
@variables x(t) y(t) z(t)
D = Differential(t)

# Define a differential equation
eqs = [D(x) ~ σ*(y-x),
       D(y) ~ t*x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

de = ODESystem(eqs)
ff = ODEFunction(de, [x,y,z], [σ,ρ,β], jac=true)

a = @SVector [1.0,2.0,3.0]
b = SLVector(x=1.0,y=2.0,z=3.0)
c = [1.0,2.0,3.0]
p = SLVector(σ=10.0,ρ=26.0,β=8/3)
@test ff(a,p,0.0) isa SVector
@test typeof(ff(b,p,0.0)) <: SLArray
@test ff(c,p,0.0) isa Vector
@test ff(a,p,0.0) == ff(b,p,0.0)
@test ff(a,p,0.0) == ff(c,p,0.0)

@test ff.jac(a,p,0.0) isa SMatrix
@test typeof(ff.jac(b,p,0.0)) <: SMatrix
@test ff.jac(c,p,0.0) isa Matrix
@test ff.jac(a,p,0.0) == ff.jac(b,p,0.0)
@test ff.jac(a,p,0.0) == ff.jac(c,p,0.0)

# Test similar_type
@test ff(b,p,ForwardDiff.Dual(0.0,1.0)) isa SLArray
d = LVector(x=1.0,y=2.0,z=3.0)
@test ff(d,p,ForwardDiff.Dual(0.0,1.0)) isa LArray
@test ff.jac(b,p,ForwardDiff.Dual(0.0,1.0)) isa SArray
@test eltype(ff.jac(b,p,ForwardDiff.Dual(0.0,1.0))) <: ForwardDiff.Dual
@test ff.jac(d,p,ForwardDiff.Dual(0.0,1.0)) isa Array
@test eltype(ff.jac(d,p,ForwardDiff.Dual(0.0,1.0))) <: ForwardDiff.Dual


text:
@testset "interactions" begin
	fh = FloryInteraction()
	set_interaction!(fh, 1, 2, 0.5)
	set_interaction!(fh, 2, 1, 0.5)
	@test length(fh.chis) == 1
	@test_throws AssertionError set_interaction!(fh, 1, 1, 0.5)
end


text:
module projection

using FFTW, PaddedViews
using ..params, ..psf
using Base.Threads
using DSP
#import FFTW.cFFTWplan
export propagate, HMatrix

struct Delta
  x::Vector{Int64}
  y::Vector{Int64}
  function Delta(lf::LightFieldSimulation)
    obj = lf.obj
    sim = lf.par.sim
    offset = obj.center
    xidx = collect(repeat(1:obj.xlen, inner = obj.ylen))
    yidx = collect(repeat(obj.ylen:-1:1, outer = obj.xlen))
    x = (xidx .- offset) .* sim.osr
    y = (yidx .- offset) .* sim.osr
    new(x, y)
  end
end

struct HMatrix

    kernel::Array{ComplexF64,2}
    flayer::Array{ComplexF64,3}
    rlayer::Array{Float64,3}
    dsrlayer::Array{Float64,3}
    multiwdf::Array{Float64,6}
    perlayer::Int
    lenslets::Int
    N::Int
    plan #::FFTW.cFFTWplan{Complex{Float64},-1,true,3}
    function HMatrix(lf::LightFieldSimulation, originpsf::Array{ComplexF64,3}, samples::Int)
        p = lf.par
        obj = lf.obj
        N = p.sim.vpix
        lenslets = fld(samples,N)
        perlayer = N^2
        psfsize = size(originpsf,1)

        kernel = fresnelH(originpsf[:,:,1], p, p.opt.d)

        flayer = Array{ComplexF64,3}(undef,(psfsize, psfsize, perlayer))
        rlayer = Array{Float64,3}(undef,(psfsize, psfsize, perlayer))
        dsrlayer = Array{Float64,3}(undef,(samples, samples, perlayer))
        multiwdf = Array{Float64,6}(undef,(N, N, lenslets, lenslets, N, N))

		FFTW.set_num_threads(Threads.nthreads()>> 1)
        plan = plan_fft!(flayer,[1,2], flags=FFTW.MEASURE)

        new(kernel, flayer, rlayer, dsrlayer, multiwdf, perlayer, lenslets, N, plan)
    end
end

"Manually generated sampling--to ensure the center point remains fixed"
function samples(lf::LightFieldSimulation)
    img = lf.img
    sim = lf.par.sim

    half1 = img.center:-sim.osr:1
    half2 = (img.center + sim.osr):sim.osr:img.xlen
    s::Vector{Int} = vcat(half1, half2)
    sort!(s)
    # TODO: fix scope shit here
    while mod(length(s), sim.vpix) > 0
        s = s[2:end-1]
    end
    return s
end

"Makes 1D sinc kernel for separable LP img filter prior to 3x downsampling"
function sinc1d()
    # NOTE THIS ASSUMES OSR OF 3
    x = range(-2, stop=2, length=13)
    scale = sum(sinc.(x))
    return sinc.(x) ./ scale
end

"Image translation using PaddedViews"
function shift(img::Array{ComplexF64,2}, Δx::Int, Δy::Int)
    @views begin
        if  Δx > 0
            if Δy > 0
                return PaddedView(0, img[1+Δy:end, 1:end-Δx], size(img), (1, 1+Δx))
            elseif Δy < 0
                return PaddedView(0, img[1:end+Δy, 1:end-Δx], size(img), (1-Δy, 1+Δx))
            else
                return PaddedView(0, img[:, 1:end-Δx], size(img), (1, 1+Δx))
            end
        elseif Δx < 0
            if Δy > 0
                return PaddedView(0, img[1+Δy:end, 1-Δx:end], size(img))
            elseif Δy < 0
                return PaddedView(0, img[1:end+Δy, 1-Δx:end], size(img), (1-Δy, 1))
            else
                return PaddedView(0, img[:, 1-Δx:end], size(img))
            end
        else
            if Δy > 0
                return PaddedView(0, img[1+Δy:end, :], size(img))
            elseif Δy < 0
                return PaddedView(0, img[1:end+Δy, :], size(img), (1-Δy, 1))
            else
                return PaddedView(0, img[:,:], size(img))
            end
        end
    end
end

function shiftpsf!(H::HMatrix, img::Array{ComplexF64,2}, delta::Delta)
    Threads.@threads for i in 1:size(H.flayer,3)
        @inbounds H.flayer[:,:,i] .= shift(img, delta.x[i], delta.y[i])
    end
end

function stackmul!(H::HMatrix, tf::Array{ComplexF64,2})
    Threads.@threads for i in 1:size(H.flayer,3)
        @views @inbounds H.flayer[:,:,i] .= H.flayer[:,:,i] .* tf
    end
end

function psfmag!(H::HMatrix)
    Threads.@threads for i in 1:size(H.flayer,3)
        @views @inbounds H.rlayer[:,:,i] .= abs2.(H.flayer[:,:,i])
    end
end

function fresnelconv!(H::HMatrix)
    H.plan*H.flayer
    stackmul!(H, H.kernel)
    H.plan\H.flayer
    return
end

function downsample!(H::HMatrix, kernel1D::Vector{Float64}, samples::Vector{Int})
    Threads.@threads for i in 1:size(H.rlayer,3)
	    FFTW.set_num_threads(1)
		H.dsrlayer[:,:,i] .= view(@view(conv(kernel1D,kernel1D, H.rlayer[:,:,i])[6:end-6,6:end-6]),samples,samples)
    end
	FFTW.set_num_threads(8)
    return
end

function lfphase!(H::HMatrix, lf::LightFieldSimulation)
  N = lf.par.sim.vpix
  zlen = lf.obj.zlen
  samplelen = size(H.dsrlayer,1)
  lenslets = H.lenslets

  for i in 1:N, j in 1:N
    bylenslets = @view(H.dsrlayer[i:N:end, j:N:end, :])
    for a in 1:lenslets, b in 1:lenslets
      H.multiwdf[i, j, a, b, :, :] .= reshape(@view(bylenslets[a,b,:]), (N,N))
    end
  end
  layer4d = reshape(view(H.dsrlayer,:), samplelen, samplelen, N, N)
  for a in 1:lenslets, i in 1:N
    x = N * a + 1 - i;
    for b in 1:lenslets, j in 1:N
      y = N * b + 1 - j;
      @views layer4d[x, y, :, :] .= H.multiwdf[:, :, a, b, i, j];
    end
  end

  for i in 1:N, j in 1:N
    layer4d[:,:,i,j] .= rot180(layer4d[:,:,i,j])
  end

  return
end

function propagate(lf::LightFieldSimulation)
    # call to generate originpsf
    originpsf = originPSFproj(lf)
    obj = lf.obj
    sinckern = sinc1d()
    delta = Delta(lf)
    samplepts = samples(lf)
    Hlayer = HMatrix(lf, originpsf, length(samplepts))
    Himgs = Array{Float64,5}(undef,length(samplepts), length(samplepts), Hlayer.N, Hlayer.N, obj.zlen)
    # make H matrix
    for layer in 1:obj.zlen
        shiftpsf!(Hlayer, originpsf[:,:,layer], delta)
		stackmul!(Hlayer, lf.mla.array)
		fresnelconv!(Hlayer)
        psfmag!(Hlayer)
        downsample!(Hlayer, sinckern, samplepts)
		# I think missing reverse shift before phase space remapping, still wicked lol
        lfphase!(Hlayer, lf)
		Himgs[:,:,:,:,layer] .= reshape(view(Hlayer.dsrlayer,:),length(samplepts),length(samplepts), lf.par.sim.vpix, lf.par.sim.vpix)
	end
    return Himgs
end

end #module

#=

function postprocess(Himgs::Array{ComplexF64,3}, objspace::Space, Zidx::Array{Int,1})
    #TODO: Needs to be revised for new pipeline

    "Normalize the whole data set"
    Himgs = Himgs./maximum(Himgs)

    """Threshold all the images (grouped by Z-plane) by the tolerance
    value. Flooring all values below the cutoff to zero."""
    tol = 0.005

    for layer in 1:lf.obj.zlen
        H4Dslice = Himgs[:,:,Zidx .== layer]
        H4Dslice[H4Dslice .< (tol*maximum(H4Dslice))] .= 0
        Himgs[:,:,Zidx .== layer] = H4Dslice
    end
    return Himgs
end
=#


text:
using Oceananigans.AbstractOperations: UnaryOperation, Derivative, BinaryOperation, MultiaryOperation
using Oceananigans.AbstractOperations: KernelFunctionOperation
using Oceananigans.Operators: ℑxyᶜᶠᵃ, ℑxyᶠᶜᵃ
using Oceananigans.Fields: PressureField, compute_at!
using Oceananigans.BuoyancyModels: BuoyancyField

using Oceananigans.AbstractOperations: Δx, Δy, Δz, Ax, Ay, Az, volume

function simple_binary_operation(op, a, b, num1, num2)
    a_b = op(a, b)
    interior(a) .= num1
    interior(b) .= num2
    return a_b[2, 2, 2] == op(num1, num2)
end

function three_field_addition(a, b, c, num1, num2)
    a_b_c = a + b + c
    interior(a) .= num1
    interior(b) .= num2
    interior(c) .= num2
    return a_b_c[2, 2, 2] == num1 + num2 + num2
end

function x_derivative(a)
    dx_a = ∂x(a)

    arch = architecture(a)
    one_two_three = arch_array(arch, [1, 2, 3])

    for k in 1:3
        interior(a)[:, 1, k] .= one_two_three
        interior(a)[:, 2, k] .= one_two_three
        interior(a)[:, 3, k] .= one_two_three
    end

    return dx_a[2, 2, 2] == 1
end

function y_derivative(a)
    dy_a = ∂y(a)

    arch = architecture(a)
    one_three_five = arch_array(arch, [1, 3, 5])

    for k in 1:3
        interior(a)[1, :, k] .= one_three_five
        interior(a)[2, :, k] .= one_three_five
        interior(a)[3, :, k] .= one_three_five
    end

    return dy_a[2, 2, 2] == 2
end

function z_derivative(a)
    dz_a = ∂z(a)

    arch = architecture(a)
    one_four_seven = arch_array(arch, [1, 4, 7])

    for k in 1:3
        interior(a)[1, k, :] .= one_four_seven
        interior(a)[2, k, :] .= one_four_seven
        interior(a)[3, k, :] .= one_four_seven
    end

    return dz_a[2, 2, 2] == 3
end

function x_derivative_cell(arch)
    grid = RectilinearGrid(size=(3, 3, 3), extent=(3, 3, 3))
    a = Field(Center, Center, Center, arch, grid, nothing)
    dx_a = ∂x(a)

    one_four_four = arch_array(arch, [1, 4, 4])

    for k in 1:3
        interior(a)[:, 1, k] .= one_four_four
        interior(a)[:, 2, k] .= one_four_four
        interior(a)[:, 3, k] .= one_four_four
    end

    return dx_a[2, 2, 2] == 3
end

function times_x_derivative(a, b, location, i, j, k, answer)
    a∇b = @at location b * ∂x(a)
    return a∇b[i, j, k] == answer
end

for arch in archs
    @testset "Abstract operations [$(typeof(arch))]" begin
        @info "Testing abstract operations [$(typeof(arch))]..."

        grid = RectilinearGrid(size=(3, 3, 3), extent=(3, 3, 3))
        u, v, w = VelocityFields(arch, grid)
        c = Field(Center, Center, Center, arch, grid, nothing)

        @testset "Unary operations and derivatives [$(typeof(arch))]" begin
            for ψ in (u, v, w, c)
                for op in (sqrt, sin, cos, exp, tanh)
                    @test typeof(op(ψ)[2, 2, 2]) <: Number
                end

                for d_symbol in Oceananigans.AbstractOperations.derivative_operators
                    d = eval(d_symbol)
                    @test typeof(d(ψ)[2, 2, 2]) <: Number
                end
            end
        end

        @testset "Binary operations [$(typeof(arch))]" begin
            generic_function(x, y, z) = x + y + z
            for (ψ, ϕ) in ((u, v), (u, w), (v, w), (u, c), (generic_function, c), (u, generic_function))
                for op_symbol in Oceananigans.AbstractOperations.binary_operators
                    op = eval(op_symbol)
                    @test typeof(op(ψ, ϕ)[2, 2, 2]) <: Number
                end
            end
        end

        @testset "Multiary operations [$(typeof(arch))]" begin
            generic_function(x, y, z) = x + y + z
            for (ψ, ϕ, σ) in ((u, v, w), (u, v, c), (u, v, generic_function))
                for op_symbol in Oceananigans.AbstractOperations.multiary_operators
                    op = eval(op_symbol)
                    @test typeof(op((Center, Center, Center), ψ, ϕ, σ)[2, 2, 2]) <: Number
                end
            end
        end

        @testset "KernelFunctionOperations [$(typeof(arch))]" begin
            trivial_kernel_function(i, j, k, grid) = 1
            op = KernelFunctionOperation{Center, Center, Center}(trivial_kernel_function, grid)
            @test op isa KernelFunctionOperation

            less_trivial_kernel_function(i, j, k, grid, u, v) = @inbounds u[i, j, k] * ℑxyᶠᶜᵃ(i, j, k, grid, v)
            op = KernelFunctionOperation{Face, Center, Center}(less_trivial_kernel_function, grid, computed_dependencies=(u, v))
            @test op isa KernelFunctionOperation

            still_fairly_trivial_kernel_function(i, j, k, grid, u, v, μ) = @inbounds μ * ℑxyᶜᶠᵃ(i, j, k, grid, u) * v[i, j, k]
            op = KernelFunctionOperation{Center, Face, Center}(still_fairly_trivial_kernel_function, grid,
                                                               computed_dependencies=(u, v), parameters=0.1)
            @test op isa KernelFunctionOperation
        end

        @testset "Fidelity of simple binary operations" begin
            @info "  Testing simple binary operations..."
            num1 = Float64(π)
            num2 = Float64(42)
            grid = RectilinearGrid(size=(3, 3, 3), extent=(3, 3, 3))

            u, v, w = VelocityFields(arch, grid)
            T, S = TracerFields((:T, :S), arch, grid)

            for op in (+, *, -, /)
                @test simple_binary_operation(op, u, v, num1, num2)
                @test simple_binary_operation(op, u, w, num1, num2)
                @test simple_binary_operation(op, u, T, num1, num2)
                @test simple_binary_operation(op, T, S, num1, num2)
            end
            @test three_field_addition(u, v, w, num1, num2)
        end

        @testset "Derivatives" begin
            @info "  Testing derivatives..."
            grid = RectilinearGrid(size=(3, 3, 3), extent=(3, 3, 3),
                                          topology=(Periodic, Periodic, Periodic))

            u, v, w = VelocityFields(arch, grid)
            T, S = TracerFields((:T, :S), arch, grid)
            for a in (u, v, w, T)
                @test x_derivative(a)
                @test y_derivative(a)
                @test z_derivative(a)
            end

            @test x_derivative_cell(arch)
        end

        @testset "Combined binary operations and derivatives" begin
            @info "  Testing combined binary operations and derivatives..."
            arch = CPU()
            Nx = 3 # Δx=1, xC = 0.5, 1.5, 2.5
            grid = RectilinearGrid(size=(Nx, Nx, Nx), extent=(Nx, Nx, Nx))
            a, b = (Field(Center, Center, Center, arch, grid, nothing) for i in 1:2)

            set!(b, 2)
            set!(a, (x, y, z) -> x < 2 ? 3x : 6)

            #                            0   0.5   1   1.5   2   2.5   3
            # x -▶                  ∘ ~~~|--- * ---|--- * ---|--- * ---|~~~ ∘
            #        i Face:    0        1         2        3          4
            #        i Center:        0         1         2         3         4

            #              a = [    0,       1.5,      4.5,       6,        0    ]
            #              b = [    0,        2,        2,        2,        0    ]
            #          ∂x(a) = [        1.5,       3,       1.5,      -6         ]

            # x -▶                  ∘ ~~~|--- * ---|--- * ---|--- * ---|~~~ ∘
            #        i Face:    0        1         2         3         4
            #        i Center:        0         1         2         3         4

            # ccc: b * ∂x(a) = [             4.5,      4.5      -4.5,            ]
            # fcc: b * ∂x(a) = [        1.5,       6,        3,       -6         ]

            C = Center
            F = Face

            @test times_x_derivative(a, b, (C, C, C), 1, 2, 2, 4.5)
            @test times_x_derivative(a, b, (C, C, C), 2, 2, 2, 4.5)
            @test times_x_derivative(a, b, (C, C, C), 3, 2, 2, -4.5)

            @test times_x_derivative(a, b, (F, C, C), 1, 2, 2, 1.5)
            @test times_x_derivative(a, b, (F, C, C), 2, 2, 2, 6)
            @test times_x_derivative(a, b, (F, C, C), 3, 2, 2, 3)
            @test times_x_derivative(a, b, (F, C, C), 4, 2, 2, -6)
        end

        grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1),
                                      topology=(Periodic, Periodic, Bounded))

        buoyancy = SeawaterBuoyancy(gravitational_acceleration = 1,
                                             equation_of_state = LinearEquationOfState(α=1, β=1))

        model = NonhydrostaticModel(architecture = arch,
                                            grid = grid,
                                        buoyancy = buoyancy,
                                         tracers = (:T, :S))

        @testset "Construction of abstract operations [$(typeof(arch))]" begin
            @info "    Testing construction of abstract operations [$(typeof(arch))]..."

            u, v, w, T, S = fields(model)

            for ϕ in (u, v, w, T)
                for op in (sin, cos, sqrt, exp, tanh)
                    @test op(ϕ) isa UnaryOperation
                end

                for ∂ in (∂x, ∂y, ∂z)
                    @test ∂(ϕ) isa Derivative
                end

                for ψ in (u, v, w, T, S)
                    @test ψ ^ ϕ isa BinaryOperation
                    @test ψ * ϕ isa BinaryOperation
                    @test ψ + ϕ isa BinaryOperation
                    @test ψ - ϕ isa BinaryOperation
                    @test ψ / ϕ isa BinaryOperation

                    for χ in (u, v, w, T, S)
                        @test ψ * ϕ * χ isa MultiaryOperation
                        @test ψ + ϕ + χ isa MultiaryOperation
                    end
                end

                for metric in (Δx, Δy, Δz, Ax, Ay, Az, volume)
                    @test location(metric * ϕ) == location(ϕ)
                end
            end

            @test u ^ 2 isa BinaryOperation
            @test u * 2 isa BinaryOperation
            @test u + 2 isa BinaryOperation
            @test u - 2 isa BinaryOperation
            @test u / 2 isa BinaryOperation
        end

        @testset "BinaryOperations with GridMetricOperation [$(typeof(arch))]" begin
            grid = RectilinearGrid(size=(1, 1, 1), extent=(2, 3, 4))

            c = CenterField(arch, grid)
            c .= 1

            # Δx, Δy, Δz = 2, 3, 4
            # Ax, Ay, Az = 12, 8, 6
            # volume = 24
            op = c * Δx;     @test op[1, 1, 1] == 2
            op = c * Δy;     @test op[1, 1, 1] == 3
            op = c * Δz;     @test op[1, 1, 1] == 4
            op = c * Ax;     @test op[1, 1, 1] == 12
            op = c * Ay;     @test op[1, 1, 1] == 8
            op = c * Az;     @test op[1, 1, 1] == 6
            op = c * volume; @test op[1, 1, 1] == 24
        end
    end
end


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie
let
    x = 1:10
    lines(x, x.^2, color = :black, linewidth = 2, linestyle = ".-", label = "x²",
    figure = (resolution = (700,450), fontsize = 18, backgroundcolor = "#D0DFE6FF"),
    axis = (xlabel = "x", ylabel = "x²", backgroundcolor = :white))
    axislegend("legend", position = :lt)
    limits!(0,10,0,100)
    current_figure()
    #save("lines1.png"), current_figure(), px_per_unit = 2.0)
    save(joinpath(@__DIR__, "output", "lines1.png"), current_figure(), px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status("CairoMakie") # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
"""
    markov_knockoffs()

Generates knockoff of variables distributed as a discrete Markov Chain
with `K` states.

# Inputs
+ `Z`: Length `p` vector of `GenotypeState` where Z[i] = (ka, kb) is the 2
    haplotype motifs of state `i`
+ `Q`: `K × K × p - 1` array. `Q[:, :, j]` is a `K × K` matrix of transition
    probabilities for `j`th state, i.e. Q[l, k, j] = P(X_{j} = k | X_{j - 1} = l)
+ `q`: `K × 1` vector of initial probabilities

# Reference
Equations 4-5 of "Gene hunting with hidden Markov model knockoffs" by
Sesia, Sabatti, and Candes
"""
function markov_knockoffs(
    Z::Vector{Int},
    Q::Array{T, 3},
    q::Vector{T}
    ) where T <: AbstractFloat
    p = length(Z)
    statespace = size(Q, 1)
    # preallocated arrays
    Z̃ = zeros(Int, p)
    N = zeros(p, statespace)
    d = Categorical([1 / statespace for _ in 1:statespace])
    return markov_knockoffs!(Z̃, Z, N, d, Q, q) # algorithm 1 in Sesia et al
end

function markov_knockoffs!(
    Z̃::Vector{Int},
    Z::Vector{Int},
    N::AbstractMatrix,
    d::Categorical, # Categorical distribution from Distributions.jl
    Q::Array{T, 3},
    q::Vector{T}
    ) where T <: AbstractFloat
    fill!(N, 0)
    for j in 1:length(Z)
        update_normalizing_constants!(N, Z, Z̃, Q, q, j) # equation 5 in Sesia et al
        single_state_dmc_knockoff!(Z̃, Z, d, N, Q, q, j) # sample Z̃j
    end
    return Z̃
end

"""
    update_normalizing_constants!(Q::AbstractMatrix{T}, q::AbstractVector{T})

Computes normalizing constants recursively using equation (5).

# Inputs
+ `Q`: `K × K × p` array. `Q[:, :, j]` is a `K × K` matrix of transition
    probabilities for `j`th state, i.e. Q[l, k, j] = P(X_{j} = k | X_{j - 1} = l).
    The first transition matrix is not used.
+ `q`: `K × 1` vector of initial probabilities

# todo: efficiency
"""
function update_normalizing_constants!(
    N::AbstractMatrix{T},
    Z::AbstractVector{Int},
    Z̃::AbstractVector{Int},
    Q::Array{T, 3},
    q::AbstractVector{T},
    j::Int
    ) where T <: AbstractFloat
    statespace, p = size(Q, 1), size(Q, 3)
    if j == 1
        @inbounds mul!(@view(N[1, :]), Transpose(@view(Q[:, :, 2])), q)
    elseif j == p
        val = 0.0
        for l in 1:statespace
            @inbounds val += Q[Z[p-1], l, p] * Q[Z̃[p-1], l, p] / N[p-1, l]
        end
        N[j, :] .= val
    else
        for k in 1:statespace, l in 1:statespace
            @inbounds N[j, k] += Q[Z[j-1], l, j] * Q[Z̃[j-1], l, j] * Q[l, k, j + 1] / N[j - 1, l]
        end
    end
    return nothing
end

function single_state_dmc_knockoff!(
    Z̃::AbstractVector{Int},
    Z::AbstractVector{Int},
    d::Categorical, # Categorical distribution from Distributions.jl
    N::AbstractMatrix{T},
    Q::Array{T, 3},
    q::AbstractVector{T},
    j::Int
    ) where T
    statespace, p = size(Q, 1), size(Q, 3)
    if j == 1
        for z̃ in 1:statespace
            @inbounds d.p[z̃] = q[z̃] * Q[z̃, Z[2], 2] / N[1, Z[2]]
        end
    elseif j == p
        for z̃ in 1:statespace
            @inbounds d.p[z̃] = Q[Z[p-1], z̃, p] * Q[Z̃[p-1], z̃, p] / N[p-1, z̃] / N[p, 1]
        end
    else
        for z̃ in 1:statespace # todo: numerical error?
            @inbounds d.p[z̃] = Q[Z[j - 1], z̃, j] * Q[Z̃[j-1], z̃, j] * Q[z̃, Z[j+1], j+1] / N[j-1, z̃] / N[j, Z[j+1]]
        end
        # if j > 7158
        #     println(d.p)
        # end
    end
    @assert sum(d.p) ≈ 1 "single_state_dmc_knockoff!: probability should sum to 1 but was $(d.p) and j = $j"
    Z̃[j] = rand(d)
    return nothing
end


text:
function log_taylor(x::Number)
    if x <= 0
        return undef
    end

    e = float(Base.MathConstants.e)
    lb = 0.546
    ub = lb * e

    res = 0

    while x < lb
        x *= e
        res -= 1
    end

    while x > ub
        x /= e
        res += 1
    end

    s = 1

    for i in 1 : 40
        res += s * (x - 1)^i / i
        s *= -1
    end

    return res
end