text:
using Revise
module Capacity
using Revise

using Atom, Juno, UnPack
using Optim, Pkg, SparseArrays
using FixedEffectModels
using CSV, DataFrames, Serialization, Dates, Statistics, StatsBase, BenchmarkTools
using LinearAlgebra, Distributions, Distributed, VegaLite, Random
using CodecZlib, CodecZstd, GLM, XLSX, InteractiveUtils, Base.Iterators
using PGFPlotsX, Colors, LaTeXStrings
#, BlackBoxOptim
#using MCMCChains
using CUDA
using Transducers, BangBang, LineSearches, LoopVectorization, Zygote

using Finometrics

#https://github.com/clintonTE/Finometrics#master

#TODO: Find a way to test the Jacobian from Zygote
#maybe see if there is a way to do a hand calculation
using CUDA: CuVector, CuMatrix, CuArray
CUDA.allowscalar(false)

import Base: Symbol, show, length, iterate, broadcastable, *

abstract type AbstractCharacteristic end

const RANDOM_SEED_VALUE = 1111
Random.seed!(RANDOM_SEED_VALUE)

#########################CRITICAL MACROS####################
#NOTE: This macro allows for the parallelization routines to be shut down
#by flicking the above PARALLEL switch
#BROKEN!!!
#=macro mpar(expr)
  quote
    if PARALLEL[]
      :($(Threads.@threads($expr)))
    else
      :($($expr))
    end
  end
end

macro mspawn(expr)
  quote
    if PARALLEL[]
      :($(Threads.@spawn($expr)))
    else
      :($($expr))
    end
  end
end=#

###########################CONSTANTS##########################
const PARAMETER_FILE = pwd() * "\\Parameters.xlsx"
const PARALLEL = Ref{Bool}(true) #perhaps we can keep this the only global state var
const PARAM = Dict{Symbol,Any}()
const LARGE_VAL = prevfloat(Inf)^0.33
const LARGE_VAL32 = (prevfloat(Inf32)^0.5f0)/10f0
const SMALL_VAL = 1e-6
const SMALL_VAL32 = 1f-6

IN_SMALL_CSV_STREAM(p::String) = open(p)
IN_SMALL_CSV_STREAM(F::Function, p::String) = open(F, p)
const SMALL_CSV_EXTENSION = "csv"

IN_CSV_STREAM(p::String) = read(p) |> x->transcode(ZstdDecompressor,x)
IN_CSV_STREAM(F::Function, p::String) = open(F, ZstdDecompressorStream, p)
const CSV_EXTENSION = "csv.zstd"

function IN_BIN_STREAM(p::String)
  local obj
  open(ZstdDecompressorStream, p) do io
      obj = deserialize(io)
  end

  return obj
end

const BIN_EXTENSION = "bin.zstd"
function OUT_BIN_STREAM(p::String, obj, level=1)

  io = ZstdCompressorStream(open(p, "w"), level=level)
  try
    serialize(io, obj)
  finally
    close(io)
  end
end

#####abstract types
abstract type AbstractIterState{T<:Real} end
abstract type AbstractVolumeParts{TM<:AbstractMatrix, TV<:AbstractVector, T<:Real} end

#note- we ahve to be careful about using methods from AbstractVolumeParts since we ahve no V₀
abstract type AbstractVolumePartsIter{TM,TV,T} <: AbstractVolumeParts{TM,TV,T} end

abstract type AbstractAnnhilator{T<:Real} end

include("Utilities\\GeneralUtilities.jl")
#include("Measure\\DistributionUtilities.jl")  #disabled due to flux incompatibility
include("Utilities\\TestUtilities.jl")
include("Utilities\\Annihilator.jl")
include("Utilities\\HyperReg.jl")
include("Utilities\\AR1.jl")

include("Preliminary\\PrelimBeta.jl")
include("Preliminary\\PrelimFunds.jl")
#include("Preliminary\\PrelimPreprocessCRSP.jl")
include("Data\\TrailingCRSPMeasure.jl")
include("Data\\Inputs\\Parameters.jl")
include("Data\\Inputs\\IBES.jl")
include("Data\\Inputs\\CCM.jl")
include("Data\\Inputs\\CRSP.jl")
include("Data\\Inputs\\COMP.jl")
include("Data\\Inputs\\FF.jl")

include("Data\\MergeCRSPCOMP.jl")
include("Data\\CRSPweekly.jl")
include("Data\\CRSPmonthly.jl")
include("Measure\\Setup\\WeightSpec.jl")
include("Measure\\Setup\\Characteristics.jl")
include("Measure\\Setup\\MeasureSpec.jl")
include("Measure\\Setup\\Standardization.jl")
include("Measure\\Verification\\Returns.jl")
include("Measure\\Verification\\Funds.jl")
#include("VolumeParts\\VolumeParts.jl")

#include("VolumeParts\\VolumePartsLSq.jl")
include("Measure\\VolumeParts\\XYIter.jl")
include("Measure\\VolumeParts\\TestXv.jl")
include("Measure\\VolumeParts\\VolumePartsIter.jl")

#include("VolumeParts\\VolumePartsFlux.jl") #disabled due to flux incompatibility
#NOTE: below may be resurrected at some point
#include("Measure\\MCMC\\CapacityChains.jl")
#include("VolumeParts\\VolumePartsMCMC.jl")

include("Measure\\Iter\\IterModels.jl")
include("Measure\\VolumeParts\\VolumePartsIterLevelDemean.jl") #mostly for testing purposes
include("Measure\\VolumeParts\\VolumePartsIterTest.jl")
include("Measure\\Iter\\IterEstimate.jl")
include("Measure\\Iter\\IterPanelAnalysis.jl")
include("Measure\\Iter\\Iter.jl")

#include("Measure\\IterSolveFP.jl")
include("Measure\\Iter\\IterSolveOptim.jl")
include("Measure\\Iter\\IterBootstrap.jl")
#include("Measure\\Iter\\IterSolveBB.jl")

#include("Measure\\MCMC\\Posteriors.jl") #NOTE: these may be resurrected for use later
#include("Measure\\MCMC\\MCMCEstimate.jl")
#include("Measure\\MCMC\\MCMCModels.jl")
#include("Measure\\MCMC\\TestMCMC.jl")

include("Measure\\Verification\\Comomentum.jl")
#include("Measure\\Verification\\Comomentumbb.jl")

include("Measure\\Setup\\Measures.jl")
include("Data\\MomentumStrategies.jl")
include("Measure\\Setup\\Controls.jl")
include("Data\\FormPanel.jl")
include("Measure\\ComputeMeasure.jl")

include("Analysis\\Analyze.jl")
include("Analysis\\LeadLag.jl")
include("Analysis\\LPFilter.jl")
include("Analysis\\MovingBlockBootstrap.jl")
include("Analysis\\Verify.jl")
include("Analysis\\IO\\Tables.jl")
include("Analysis\\IO\\Graphs.jl")
include("Analysis\\IO\\LLTables.jl")
include("Analysis\\IO\\IO.jl")

REGRESSION_METHOD = Dict{Symbol, Any}(
#  :flux=>formfluxregression,
#  :opt=>formoptregression,
  :iter=>formiterregression,
  #:mcmc=>formmcmcregression
)

end # module

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using HOODESolver
using LinearAlgebra
using Plots
using Plots.PlotMeasures

function fctmain(n_tau, prec)
    setprecision(prec)
    u0 = BigFloat.([-34//100, 78//100, 67//100, -56//10])
    B = BigFloat.([12//100 -78//100 91//100 34//100
        -45//100 56//100 3//100 54//100
        -67//100 09//100 18//100 89//100
        -91//100 -56//100 11//100 -56//100])
    alpha =  BigFloat.([12//100, -98//100, 45//100, 26//100])
    beta =  BigFloat.([-4//100, 48//100, 23//100, -87//100])
    println("u0=$u0")
    println("B=$B")
    println("alpha=$alpha")
    println("beta=$beta")
    fct = (u,p,t) -> B*u + t*p[1] +p[2]
    t_max = big"1.0"
    epsilon=big"0.015"
    A = [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
    prob = HOODEProblem(fct,u0, (big"0.0",t_max), (alpha, beta), A, epsilon, B)
#    nbmaxtest=6
    nbmaxtest=12
    ordmax=17
    debord=3
    pasord=2
    y = ones(Float64, nbmaxtest, div(ordmax-debord,pasord)+1 )
    x=zeros(Float64,nbmaxtest)
    ind=1
    for order=debord:pasord:ordmax
        nb = 100
        indc = 1
        labels=Array{String,2}(undef, 1, div(order-debord,pasord)+1)
        ordprep = min(order+2,10)
        ordprep = order+2
        println("preparation ordre $ordprep")
        par_u0=missing
        while indc <= nbmaxtest
            println("u0=$u0")
            println("B=$B")
            println("alpha=$alpha")
            println("beta=$beta")
            @time res = solve(prob, nb_tau=n_tau, order=order, order_prep=ordprep, nb_t=nb,par_u0=par_u0, dense=false)
            par_u0=res.par_u0
            sol = res[end]
            solref=getexactsol(res.par_u0.parphi, u0, t_max)
            println("solref=$solref")
            println("nb=$nb sol=$sol")
            diff=solref-sol
            x[indc] = 1.0/nb
            println("nb=$nb dt=$(1.0/nb) normInf=$(norm(diff,Inf)) norm2=$(norm(diff))")
            resnorm = Float64(norm(diff,Inf))
            y[indc,ind] = (resnorm < 1) ? resnorm : NaN
            println("result=$y")
            println("log2(y)=$(log2.(y))")
            nb *= 2
            indc += 1
        end
        for i=debord:pasord:order
            labels[1,div((i-debord),pasord)+1] = " order=$i "
        end
        yv = y[:,1:ind]
        p=Plots.plot(
    x,
    yv,
    xlabel="Δt",
    xaxis=:log,
    ylabel="error",
    yaxis=:log,
    legend=:bottomright,
    label=labels,
    marker=2,
)
        pp=Plots.plot(
    x,
    yv,
    xlabel="Δt",
    xaxis=:log,
    ylabel="error",
    yaxis=:log,
    legend=:bottomright,
    label=labels,
    marker=2,
    bottom_margin=30px,
    left_margin=60px
)
        prec_v = precision(BigFloat)
        eps_v = convert(Float32,epsilon)
        Plots.savefig(p,"out/res9p_$(prec_v)_$(eps_v)_$(order)_$(ordprep)_$(n_tau)_exact.pdf")
        Plots.savefig(p,"out/res9p_$(prec_v)_$(eps_v)_$(order)_$(ordprep)_$(n_tau)_exact.png")
        Plots.savefig(pp,"out/res9pp_$(prec_v)_$(eps_v)_$(order)_$(ordprep)_$(n_tau)_exact.pdf")
        Plots.savefig(pp,"out/res9pp_$(prec_v)_$(eps_v)_$(order)_$(ordprep)_$(n_tau)_exact.png")
        ind+= 1
    end
end

# testODESolver()

# for i=3:9
#     fctMain(2^i)
# end
fctmain(32, 512)

ext: jl


text:
function create_displacement(mol, A::Int, i::Int, h)

    new_mol = deepcopy(mol)
    new_mol.atoms[A].xyz[i] += h

    return new_mol
end

function apply_gradient(mol, g, d = 0.001)
    Zvals = [A.Z for A = mol.atoms]
    Svals = [A.AtomicSymbol for A = mol.atoms]

    new_atoms = Fermi.Atom[]

    for i = eachindex(Zvals)
        x = mol.atoms[i].xyz[1] - d*g[i, 1]
        y = mol.atoms[i].xyz[2] - d*g[i, 2]
        z = mol.atoms[i].xyz[3] - d*g[i, 3]
        push!(new_atoms, Fermi.Atom(Svals[i], Zvals[i], (x,y,z)))
    end

    return Fermi.Molecule(new_atoms, mol.charge, mol.multiplicity)
end

function geom_rms(mol1, mol2)
    out = 0.0
    N = length(mol1.atoms)
    for a in 1:N
        out += sum((mol1.atoms[a].xyz .- mol2.atoms[a].xyz).^2)
    end
    return √(out / 3*N)
end

function findif_intgrad(X::String, mol, A, i, h=0.005)
    mol_disp = create_displacement(mol, A, i, h)
    I = Fermi.Integrals.IntegralHelper(molecule=mol_disp)
    Xplus = I[X]
    mol_disp = create_displacement(mol, A, i, -h)
    I = Fermi.Integrals.IntegralHelper(molecule=mol_disp)
    Xminus = I[X]
    g = (Xplus - Xminus) ./ (2*h)
    return g * PhysicalConstants.bohr_to_angstrom
end

function gradient_test(mol, energy_function, h=0.005)
    N = length(mol.atoms)
    Eplus = zeros(N,3)
    Eminus = zeros(N,3)

    # Plus
    for A = 1:N, i = 1:3
        mol_disp = create_displacement(mol, A, i, h)
        wfn = eval(Expr(:call, energy_function, mol_disp))
        Eplus[A,i] = wfn.energy
    end
    # Minus
    for A = 1:N, i = 1:3
        mol_disp = create_displacement(mol, A, i, -h)
        wfn = eval(Expr(:call, energy_function, mol_disp))
        Eminus[A,i] = wfn.energy
    end

    g = (Eplus - Eminus) ./ (2*h)

    return g
end

function opt_test(energy_function; h=0.005, d=0.01)
    @set printstyle none

    scf_Etol  = Options.get("scf_e_conv")
    scf_Dtol  = Options.get("scf_max_rms")
    old_mol = Fermi.Molecule()

    # Central
    wfn = eval(Expr(:call, energy_function, old_mol))
    oldE = wfn.energy
    println(format("Initial Energy: {:15.10f}", oldE))

    ite = 1
    dE = 1
    while abs(dE) > 1e-8
        g = gradient_test(old_mol, energy_function, h)
        new_mol = apply_gradient(old_mol, g, d)
        wfn = eval(Expr(:call, energy_function, new_mol))
        if abs(wfn.e_conv) > scf_Etol || abs(wfn.d_conv) > scf_Dtol
            display(new_mol)
            error("SCF not converged")
        end
        dE = oldE - wfn.energy
        RMS = √(sum(g.^2) / length(g))
        GRMS = geom_rms(old_mol, new_mol)
        println(format("Iter {:3}   Energy: {:15.10f}   ΔE: {:15.10f}   RMS: {:15.10f}   GRMS: {:15.10f}", ite, wfn.energy, dE, RMS, GRMS))
        oldE = wfn.energy
        old_mol = new_mol
        ite += 1
    end

    display(old_mol)
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
@recipe function f(e::Edge)
    linecolor --> :black
    label --> "Edge"
    xguide --> "x"
    yguide --> "y"
    zguide --> "z"
    unitformat --> :slash
    if occursin("GRBackend", string(typeof(plotattributes[:plot_object].backend)))
        aspect_ratio --> 1.0
    end
    @series begin
        seriestype --> :path3d
        [e.a[1], e.b[1]]*internal_length_unit, [e.a[2], e.b[2]]*internal_length_unit, [e.a[3], e.b[3]]*internal_length_unit
    end
end


text:
# # Custom structs

#=
In this example, we demonstrate implicit differentiation through functions that manipulate `NamedTuple`s, as a first step towards compatibility with general structs.
=#

using ComponentArrays
using ImplicitDifferentiation
using Zygote

using ChainRulesCore  #src
using ChainRulesTestUtils  #src
using Test  #src

# ## Implicit function wrapper

#=
We replicate a componentwise square function with `NamedTuple`s, taking `a=(x,y)` as input and returning `b=(u,v)`.
=#

forward(a::ComponentVector) = ComponentVector(u=a.x .^ 2, v=a.y .^ 2)

function conditions(a::ComponentVector, b::ComponentVector)
    return vcat(b.u .- a.x .^ 2, b.v .- a.y .^ 2)
end

implicit = ImplicitFunction(forward, conditions);

#=
In order to be able to call `Zygote.gradient`, we use `implicit` to define a convoluted version of the squared Euclidean norm, which takes a `ComponentVector` as input and returns a real number.
=#

function fullnorm(a::ComponentVector)
    b = implicit(a)
    return sum(b)
end

function halfnorm(a::ComponentVector)  #src
    b = implicit(a)  #src
    return sum(b.u)  #src
end;  #src

# ## Testing

a = ComponentVector(x=rand(3), y=rand(3))

# Let us first check that our weird squared norm returns the correct result.

fullnorm(a) ≈ sum(abs2, a)
halfnorm(a) ≈ sum(abs2, a.x)  #src

# Now we go one step further and compute its gradient, which involves the reverse rule for `implicit`.

Zygote.gradient(fullnorm, a)[1] ≈ 2a

#-

Zygote.gradient(halfnorm, a)[1] ≈ vcat(2a.x, zero(a.y))  #src

# The following tests are not included in the docs.  #src

@test fullnorm(a) ≈ sum(abs2, a)  #src
@test halfnorm(a) ≈ sum(abs2, a.x)  #src

@test Zygote.gradient(fullnorm, a)[1] ≈ 2a  #src
@test Zygote.gradient(halfnorm, a)[1] ≈ vcat(2a.x, zero(a.y))  #src


text:
using LinearAlgebra
using Test

using SetProg, SetProg.Sets
using Polyhedra
using MultivariatePolynomials
const MP = MultivariatePolynomials

import DynamicPolynomials
import MultivariateBases
const MB = MultivariateBases
const MonoBasis = MB.MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}

using JuMP
const MOIT = MOI.Test

function ci_square_test(optimizer, config::MOIT.Config,
                        inner::Bool, variable::SetProg.AbstractVariable,
                        metric::Function, objective_value, set_test, nvars=nothing)
    model = _model(optimizer)

    □ = polyhedron(HalfSpace([1, 0], 1.0) ∩ HalfSpace([-1, 0], 1) ∩ HalfSpace([0, 1], 1) ∩ HalfSpace([0, -1], 1))

    @variable(model, ◯, variable)
    if inner
        cref = @constraint(model, ◯ ⊆ □)
    else
        cref = @constraint(model, □ ⊆ ◯)
    end

    Δt = 0.5
    A = [1.0 Δt]
    E = [1.0 0.0]

    if variable.symmetric
        @constraint(model, A * ◯ ⊆ E * ◯)
    else
        @constraint(model, A * ◯ ⊆ E * ◯, S_procedure_scaling = 1.0)
    end

    @objective(model, inner ? MOI.MAX_SENSE : MOI.MIN_SENSE,
               metric(volume(◯)))

    SetProg.optimize!(model)
    if nvars !== nothing
        @test nvars == num_variables(model)
    end
    @test JuMP.termination_status(model) == MOI.OPTIMAL
    @test JuMP.objective_sense(model) == MOI.MAX_SENSE
    @test JuMP.objective_value(model) ≈ objective_value atol=config.atol rtol=config.rtol
    set_test(JuMP.value(◯))
end

function ci_ell_homogeneous_test(optimizer, config)
    ci_square_test(optimizer, config, true,
                   Ellipsoid(symmetric=true, dimension=2), nth_root, √15/4,
                   ◯ -> begin
                       @test ◯ isa Sets.Polar{Float64, Sets.Ellipsoid{Float64}}
                       @test Sets.polar(◯).Q ≈ Symmetric([1.0 -1/4; -1/4 1.0]) atol=config.atol rtol=config.rtol
                   end)
end

function ci_ell_nonhomogeneous_test(optimizer, config)
    ci_square_test(optimizer, config, true,
                   Ellipsoid(point=SetProg.InteriorPoint([0.0, 0.0])),
                   nth_root, √15/4,
                   ◯ -> begin
                       @test ◯ isa Sets.PerspectiveDual{Float64, Sets.Householder{Float64, Sets.ShiftedEllipsoid{Float64}, Float64}}
                       z = Sets.perspective_variable(◯)
                       x, y = Sets.space_variables(◯)
                       ◯_dual = Sets.perspective_dual(◯)
                       @test ◯_dual.p ≈ -z^2 + x^2 - x*y/2 + y^2 atol=config.atol rtol=config.rtol
                       @test ◯_dual.set.Q ≈ Symmetric([1.0 -1/4; -1/4 1.0]) atol=config.atol rtol=config.rtol
                       @test ◯_dual.set.b ≈ [0.0, 0.0] atol=config.atol rtol=config.rtol
                       @test ◯_dual.set.β ≈ -1.0 atol=config.atol rtol=config.rtol
                       @test Sets._householder(◯_dual.h) ≈ [-1.0 0.0 0.0
                                                             0.0 1.0 0.0
                                                             0.0 0.0 1.0] atol=config.atol rtol=config.rtol
                   end)
end

function ci_piecewise_semiell_homogeneous_test(optimizer, config)
    ci_square_test(
        optimizer, config, true,
        Ellipsoid(symmetric=true, piecewise=◇),
        set -> L1_heuristic(set), 19 / 6,
        p◯ -> begin
            @test p◯ isa Sets.Polar
            ◯ = p◯.set
            @test ◯ isa Sets.Piecewise{Float64, Sets.Ellipsoid{Float64}}
            @test length(◯.sets) == 4
            Q1 = Symmetric([ 1.0 -0.25
                            -0.25 1.0])
            Q2 = Symmetric([ 1.0 -1.0
                            -1.0  1.0])
            @test ◯.sets[1].Q ≈ Q1 atol=config.atol rtol=config.rtol
            @test ◯.sets[2].Q ≈ Q2 atol=config.atol rtol=config.rtol
            @test ◯.sets[3].Q ≈ Q2 atol=config.atol rtol=config.rtol
            @test ◯.sets[4].Q ≈ Q1 atol=config.atol rtol=config.rtol
        end,
        24,
    )
end

function ci_piecewise_semiell_mci_homogeneous_test(optimizer, config)
    polar_mci = polyhedron(convexhull([1.0, 0.0], [-1.0, 0.0], [0.0, 1.0], [0.0, -1.0], [1.0, 0.5], [-1.0, -0.5]), lib)
    ci_square_test(
        optimizer, config, true,
        Ellipsoid(symmetric=true, piecewise=polar_mci),
        set -> L1_heuristic(set), 0.9349838949990186,
        p◯ -> begin
            @test p◯ isa Sets.Polar
            ◯ = p◯.set
            @test ◯ isa Sets.Piecewise{Float64, Sets.Ellipsoid{Float64}}
            @test length(◯.sets) == 6
            Q1 = Symmetric([ 1.0 -1.0
                            -1.0  1.0])
            Q2 = Symmetric([ 1.0  0.0
                             0.0  0.0])
            Q3 = Symmetric([ 0.25 0.5
                             0.5  1.0])
            @test ◯.sets[1].Q ≈ Q1 atol=config.atol rtol=config.rtol
            @test ◯.sets[2].Q ≈ Q1 atol=config.atol rtol=config.rtol
            @test ◯.sets[3].Q ≈ Q2 atol=config.atol rtol=config.rtol
            @test ◯.sets[4].Q ≈ Q3 atol=config.atol rtol=config.rtol
            @test ◯.sets[5].Q ≈ Q2 atol=config.atol rtol=config.rtol
            @test ◯.sets[6].Q ≈ Q3 atol=config.atol rtol=config.rtol
        end,
        109,
    )
end

function ci_quad_nonhomogeneous_test(optimizer, config)
    ci_square_test(optimizer, config, true,
                   PolySet(degree=2, convex=true, point=SetProg.InteriorPoint([0.0, 0.0])),
                   set -> L1_heuristic(set, [1.0, 1.0]), 8/3,
                   ◯ -> begin
                       @test ◯ isa Sets.PerspectiveDual{Float64, Sets.Householder{Float64, Sets.ConvexPolynomialSet{Float64, MonoBasis, Float64}, Float64}}
                       z = Sets.perspective_variable(◯)
                       x, y = Sets.space_variables(◯)
                       ◯_dual = Sets.perspective_dual(◯)
                       # The coefficient of `x*y` does not influence the volume
                       # and with the values of the other parameters, it should
                       # simply be in the interval [-2, -0.5].
                       α = MP.coefficient(◯_dual.p, x*y)
                       @test α ≥ -2 - 2config.atol - config.rtol
                       @test α ≤ -0.5 + 0.5config.atol + config.rtol
                       @test ◯_dual.p ≈ -z^2 + x^2 + α*x*y + y^2 atol=config.atol rtol=config.rtol
                   end)
end

function ci_quartic_homogeneous_test(optimizer, config)
    ci_square_test(optimizer, config, true,
                   PolySet(symmetric=true, degree=4, convex=true),
                   set -> L1_heuristic(set, [1.0, 1.0]),
                   0.4,
                   ◯ -> begin
                       @test ◯ isa Sets.Polar{Float64, Sets.ConvexPolySet{Float64, MonoBasis, Float64}}
                       @test Sets.polar(◯).degree == 4
                       x, y = variables(Sets.polar(◯).p)
                       α = MP.coefficient(Sets.polar(◯).p, x^3*y) / 2
                       q = x^4 + 2α*x^3*y + 6x^2*y^2 + 2α*x*y^3 + y^4
                       @test all(eigvals(Matrix(Sets.polar(◯).p.Q)) .≥ -config.atol)
                       @test polynomial(Sets.polar(◯).p) ≈ q atol=config.atol rtol=config.rtol
                       convexity_proof = Sets.convexity_proof(◯)
                       @test convexity_proof.n == 4

                       hess = 6 * [2.0, α, 2.0, α, 2.0,
                                   2.0, 2.0, α, α, 2.0]
                       Hess = SetProg.SumOfSquares.MultivariateMoments.SymMatrix(hess, 4)
                       @test all(eigvals(Matrix(Hess)) .≥ -config.atol)
                       @test convexity_proof.Q ≈ hess atol=config.atol rtol=config.rtol
                   end)
end

const ci_tests = Dict(
    "ci_ell_homogeneous" =>
     ci_ell_homogeneous_test,
    "ci_ell_nonhomogeneous" =>
     ci_ell_nonhomogeneous_test,
    "ci_piecewise_semiell_homogeneous" =>
     ci_piecewise_semiell_homogeneous_test,
    "ci_piecewise_semiell_mci_homogeneous" =>
     ci_piecewise_semiell_mci_homogeneous_test,
    "ci_quad_nonhomogeneous" =>
     ci_quad_nonhomogeneous_test,
    "ci_quartic_homogeneous" =>
     ci_quartic_homogeneous_test
)

@test_suite ci


text:
https://en.wikipedia.org/wiki/J._K._Rowling
https://en.wikipedia.org/wiki/Magician_(fantasy)
https://en.wikipedia.org/wiki/Harry_Potter_(character)
https://en.wikipedia.org/wiki/Hermione_Granger
https://en.wikipedia.org/wiki/Ron_Weasley
https://en.wikipedia.org/wiki/Hogwarts_School_of_Witchcraft_and_Wizardry
https://en.wikipedia.org/wiki/Lord_Voldemort
https://en.wikipedia.org/wiki/Ministry_of_Magic
https://en.wikipedia.org/wiki/Muggle

max_stars_repo_name: pncnmnp/search-engine

max_issues_repo_name: pncnmnp/search-engine

max_forks_repo_name: pncnmnp/search-engine


text:
#####
##### `reshape`
#####

function rrule(::typeof(reshape), A::AbstractArray, dims::Tuple{Vararg{Union{Colon,Int}}})
    A_dims = size(A)
    function reshape_pullback(Ȳ)
        return (NoTangent(), reshape(Ȳ, A_dims), NoTangent())
    end
    return reshape(A, dims), reshape_pullback
end

function rrule(::typeof(reshape), A::AbstractArray, dims::Union{Colon,Int}...)
    A_dims = size(A)
    function reshape_pullback(Ȳ)
        ∂A = reshape(Ȳ, A_dims)
        ∂dims = broadcast(_ -> NoTangent(), dims)
        return (NoTangent(), ∂A, ∂dims...)
    end
    return reshape(A, dims...), reshape_pullback
end

#####
##### `repeat`
#####

function rrule(::typeof(repeat), xs::AbstractArray; inner=ntuple(_->1, ndims(xs)), outer=ntuple(_->1, ndims(xs)))

    function repeat_pullback(ȳ)
        dY = unthunk(ȳ)
        Δ′ = zero(xs)
        S = size(xs)

        # Loop through each element of Δ, calculate source dimensions, accumulate into Δ′
        for (dest_idx, val) in pairs(IndexCartesian(), dY)
            # First, round dest_idx[dim] to nearest gridpoint defined by inner[dim], then
            # wrap around based on original size S.
            src_idx = [mod1(div(dest_idx[dim] - 1, inner[dim]) + 1, S[dim]) for dim in 1:length(S)]
            Δ′[src_idx...] += val
        end
        return (NoTangent(), Δ′)
    end

    return repeat(xs; inner = inner, outer = outer), repeat_pullback
end

function rrule(::typeof(repeat), xs::AbstractVector, m::Integer)

    d1 = size(xs, 1)
    function repeat_pullback(ȳ)
        Δ′ = dropdims(sum(reshape(ȳ, d1, :); dims=2); dims=2)
        return (NoTangent(), Δ′, NoTangent())
    end

    return repeat(xs, m), repeat_pullback
end

function rrule(::typeof(repeat), xs::AbstractVecOrMat, m::Integer, n::Integer)
    d1, d2 = size(xs, 1), size(xs, 2)
    function repeat_pullback(ȳ)
        ȳ′ = reshape(ȳ, d1, m, d2, n)
        return NoTangent(), reshape(sum(ȳ′; dims=(2,4)), (d1, d2)), NoTangent(), NoTangent()
    end

    return repeat(xs, m, n), repeat_pullback
end

function rrule(T::typeof(repeat), xs::AbstractVecOrMat, m::Integer)

    # Workaround use of positional default (i.e. repeat(xs, m, n = 1)))
    y, full_pb = rrule(T, xs, m, 1)
    repeat_pullback(ȳ) = full_pb(ȳ)[1:3]
    return y, repeat_pullback
end

#####
##### `hcat`
#####

function rrule(::typeof(hcat), Xs::Union{AbstractArray, Number}...)
    Y = hcat(Xs...)  # note that Y always has 1-based indexing, even if X isa OffsetArray
    ndimsY = Val(ndims(Y))  # this avoids closing over Y, Val() is essential for type-stability
    sizes = map(size, Xs)   # this avoids closing over Xs
    function hcat_pullback(ȳ)
        dY = unthunk(ȳ)
        hi = Ref(0)  # Ref avoids hi::Core.Box
        dXs = map(sizes) do sizeX
            ndimsX = length(sizeX)
            lo = hi[] + 1
            hi[] += get(sizeX, 2, 1)
            ind = ntuple(ndimsY) do d
                if d==2
                    d > ndimsX ? lo : lo:hi[]
                else
                    d > ndimsX ? 1 : (:)
                end
            end
            if ndimsX > 0
                # Here InplaceableThunk breaks @inferred, removed for now
                # InplaceableThunk(@thunk(dY[ind...]), dX -> dX .+= view(dY, ind...))
                dY[ind...]
            else
                # This is a hack to perhaps avoid GPU scalar indexing
                sum(view(dY, ind...))
            end
        end
        return (NoTangent(), dXs...)
    end
    return Y, hcat_pullback
end

function rrule(::typeof(reduce), ::typeof(hcat), As::AbstractVector{<:AbstractVecOrMat})
    widths = map(A -> size(A,2), As)
    function reduce_hcat_pullback_2(dY)
        hi = Ref(0)
        dAs = map(widths) do w
            lo = hi[]+1
            hi[] += w
            dY[:, lo:hi[]]
        end
        return (NoTangent(), NoTangent(), dAs)
    end
    return reduce(hcat, As), reduce_hcat_pullback_2
end

function rrule(::typeof(reduce), ::typeof(hcat), As::AbstractVector{<:AbstractVector})
    axe = axes(As,1)
    function reduce_hcat_pullback_1(dY)
        hi = Ref(0)
        dAs = map(_ -> dY[:, hi[]+=1], axe)
        return (NoTangent(), NoTangent(), dAs)
    end
    return reduce(hcat, As), reduce_hcat_pullback_1
end

#####
##### `vcat`
#####

function rrule(::typeof(vcat), Xs::Union{AbstractArray, Number}...)
    Y = vcat(Xs...)
    ndimsY = Val(ndims(Y))
    sizes = map(size, Xs)
    function vcat_pullback(ȳ)
        dY = unthunk(ȳ)
        hi = Ref(0)
        dXs = map(sizes) do sizeX
            ndimsX = length(sizeX)
            lo = hi[] + 1
            hi[] += get(sizeX, 1, 1)
            ind = ntuple(ndimsY) do d
                if d==1
                    d > ndimsX ? lo : lo:hi[]
                else
                    d > ndimsX ? 1 : (:)
                end
            end
            if ndimsX > 0
                # InplaceableThunk(@thunk(dY[ind...]), dX -> dX .+= view(dY, ind...))
                dY[ind...]
            else
                sum(view(dY, ind...))
            end
        end
        return (NoTangent(), dXs...)
    end
    return Y, vcat_pullback
end

function rrule(::typeof(reduce), ::typeof(vcat), As::AbstractVector{<:AbstractVecOrMat})
    Y = reduce(vcat, As)
    ndimsY = Val(ndims(Y))
    heights = map(A -> size(A,1), As)
    function reduce_vcat_pullback(dY)
        hi = Ref(0)
        dAs = map(heights) do z
            lo = hi[]+1
            hi[] += z
            ind = ntuple(d -> d==1 ? (lo:hi[]) : (:), ndimsY)
            dY[ind...]
        end
        return (NoTangent(), NoTangent(), dAs)
    end
    return Y, reduce_vcat_pullback
end

#####
##### `cat`
#####

_val(::Val{x}) where {x} = x

function rrule(::typeof(cat), Xs::Union{AbstractArray, Number}...; dims)
    Y = cat(Xs...; dims=dims)
    cdims = dims isa Val ? Int(_val(dims)) : dims isa Integer ? Int(dims) : Tuple(dims)
    ndimsY = Val(ndims(Y))
    sizes = map(size, Xs)
    function cat_pullback(ȳ)
        dY = unthunk(ȳ)
        prev = fill(0, _val(ndimsY))  # note that Y always has 1-based indexing, even if X isa OffsetArray
        dXs = map(sizes) do sizeX
            ndimsX = length(sizeX)
            index = ntuple(ndimsY) do d
                if d in cdims
                    d > ndimsX ? (prev[d]+1) : (prev[d]+1:prev[d]+sizeX[d])
                else
                    d > ndimsX ? 1 : (:)
                end
            end
            for d in cdims
                prev[d] += get(sizeX, d, 1)
            end
            if ndimsX > 0
                # InplaceableThunk(@thunk(dY[index...]), dX -> dX .+= view(dY, index...))
                dY[index...]
            else
                sum(view(dY, index...))
            end
        end
        return (NoTangent(), dXs...)
    end
    return Y, cat_pullback
end

#####
##### `hvcat`
#####

function rrule(::typeof(hvcat), rows, values::Union{AbstractArray, Number}...)
    Y = hvcat(rows, values...)
    cols = size(Y,2)
    ndimsY = Val(ndims(Y))
    sizes = map(size, values)
    function hvcat_pullback(dY)
        prev = fill(0, 2)
        dXs = map(sizes) do sizeX
            ndimsX = length(sizeX)
            index = ntuple(ndimsY) do d
                if d in (1, 2)
                    d > ndimsX ? (prev[d]+1) : (prev[d]+1:prev[d]+sizeX[d])
                else
                    d > ndimsX ? 1 : (:)
                end
            end
            prev[2] += get(sizeX, 2, 1)
            if prev[2] == cols
                prev[2] = 0
                prev[1] += get(sizeX, 1, 1)
            end
            dY[index...]
        end
        return (NoTangent(), NoTangent(), dXs...)
    end
    return Y, hvcat_pullback
end

#####
##### `fill`
#####

function rrule(::typeof(fill), value::Any, dims::Tuple{Vararg{Int}})
    function fill_pullback(Ȳ)
        return (NoTangent(), sum(Ȳ), NoTangent())
    end
    return fill(value, dims), fill_pullback
end

function rrule(::typeof(fill), value::Any, dims::Int...)
    function fill_pullback(Ȳ)
        return (NoTangent(), sum(Ȳ), ntuple(_->NoTangent(), length(dims))...)
    end
    return fill(value, dims), fill_pullback
end


text:
include("abstraction.jl")
include("../utilities/tDATAsets.jl")
using DataFrames
using SparseArrays

function diam(vertices)
    submatrix = M_dist[vertices, vertices]
    if (submatrix .|> iszero |> sum) > length(vertices)
        return Inf
    else
        return maximum(M_dist[vertices, vertices])
    end
end

Random.seed!(0)
M2n = rand(Sphere(2), 10)
d, n = size(M2n)
ε = 2.0
# dims = n-1
dims = 8

# Initialization
kdtree = KDTree(M2n)
vrtx, dist = knn(kdtree, M2n, dims+1, true)
M_dist = spzeros(Float64, n,n)
for i in 1:n
  for (idxj, j) in enumerate(vrtx[i])
    M_dist[i,j] = dist[i][idxj]
  end
end

df = DataFrame(ID = UInt64[], Dimension = Int64[],
               Appearance = Float64[], Simplex = Set{Int64}[])
Δ = Dict()
for k ∈ 0:dims
    push!(Δ, k::Int64 => Set{Set{Int64}}())
    if k == 0
        for j ∈ 1:n
            splx = Set([j])
            push!(Δ[k], splx)
            push!(df, [hash(splx), 0, 0, splx])
        end
        continue
    end
    for face in Δ[k-1]
        for v ∈ 1:n
            if v ∉ face
                candy = [collect(face) ; v]
            else
                continue
            end
            δ = diam(candy)
            if 0 < δ ≤ ε
                splx = Set(candy)
                push!(Δ[k], splx)
                push!(df, [hash(splx), k, δ, splx])
            end
        end
    end
    # println(k)
    # # k = 6
    # δ = last.(dist)
    # id_VR = findall(δ .≤ ε)
    # if !(id_VR |> isempty)
    #     for id_splx ∈ id_VR
    #         splx = Set(vrtx[id_splx])
    #         push!(Δ[k], Set(splx))
    #     end
    # end
    # map(pop!, vrtx)
    # map(pop!, dist)
end
Δ = sort(Δ)
sort!(df, :Appearance, rev = true)
unique!(df, :ID)
show(stdout, "text/plain", df)
sort!(df, :Dimension, rev = true)

print(1)
# for k in dims:-1:0
# # k = 0
#     vert = view.(prevert, Ref((0:k+1).+1))
#     dist = getindex.(predist, k+1)
#     bit_VR = dist.≤2ε
#     preΔ[k+1] = vert[bit_VR] .|> sort |> unique
#     for β in preΔ[k]
#         # flag_Δ = false
#         # for α in preΔ[k+1]
#         #     # Skipping if β already in Δ
#         #     # If we have {[1], [2], [3], [1,2]},
#         #     # then we want {[3], [1,2]}
#         #     if β.V ⊆ α.V
#         #         flag_Δ = true
#         #         break
#         #     end
#         # end
#         # if flag_Δ continue end
#         push!(Δ[k], β)
#     end

#     # Setting next dimension or escape
#     if preΔ[k+1] |> isempty
#         break
#     else
#         Δ[k+1] = Int64[]
#     end
# end
#    Δ = sort(Δ)
# preΔ = sort(preΔ)

x, y = eachrow(M2n);
p1 = scatter(x, y, text = 1:n, size = (600, 600), color = :white)
for (idx, δ) ∈ enumerate(reduce(∪, dist) |> sort)
    p2 = plot(p1)
    for t = 1:10
        p2 = circle!(p2, (x[t], y[t]), δ/2, label = :none)
    end
    png(p2, "D:/TDAexmaple/$idx-δ=$δ.png")
end


text:
@testset "Test ranking" begin
    parameters1 = BAT.NamedTupleDist(
        p1 = -0.2..0.2,
        p2 = -0.1..0.1,
    )

    function testfunc1(params)
    c = [20.12, 5.56, 325.556]
    return c[1] * params.p1^2 + c[2] * params.p1 * params.p2 + c[3] * params.p2^2
    end

    measurements1 = (
    meas1 = Measurement(testfunc1, 0.2,
    uncertainties = (unc1=10.1, unc2=12.2, unc3=13.3), active=true),

    meas2 = Measurement(Observable(testfunc1, min=0, max=1000), 0.3,
    uncertainties = (unc1=20.1, unc2=20.2, unc3=23.3), active=false),

    meas3 = Measurement(Observable(testfunc1, min=0, max=1000), 0.3,
    uncertainties = (unc1=30.1, unc2=30.2, unc3=30.3), active=true),

    meas4 = MeasurementDistribution(Function[testfunc1, testfunc1, testfunc1],
    [0.28, 0.27, 0.29], uncertainties = (unc1=[0.11, 0.12, 0.13], unc2=[0.21, 0.22, 0.23], unc3=[0.31, 0.32, 0.33]),
    active = [true, false, true])
    )

    corr_matrix = to_correlation_matrix(measurements1,
    (:meas1, :meas2, 0.12), # will be ignored later in EFTfitterModel
    (:meas1, :meas3, 0.13),
    (:meas1, :meas4_bin1, 0.141),
    (:meas1, :meas4_bin2, 0.142), # will be ignored later in EFTfitterModel
    (:meas1, :meas4_bin3, 0.143), # will be ignored later in EFTfitterModel
    (:meas4_bin1, :meas4_bin2, 0.412), # will be ignored later in EFTfitterModel
    (:meas4_bin1, :meas4_bin2, 0.413),
    (:meas4_bin1, :meas4_bin1, 0.9) # will be ignored later in EFTfitterModel
    )

    @test corr_matrix ≈ [ 1.0    0.12  0.13  0.141  0.142  0.143;
                     0.12   1.0   0.0   0.0    0.0    0.0;
                     0.13   0.0   1.0   0.0    0.0    0.0;
                     0.141  0.0   0.0   1.0    0.413  0.0;
                     0.142  0.0   0.0   0.413  1.0    0.0;
                     0.143  0.0   0.0   0.0    0.0    1.0] rtol = 0.01

    correlations1 = (
    unc1 = NoCorrelation(active=true),

    # wrong matrix size for number of measurements, will be ignored if active=false:
    unc2 = Correlation([1.0 0.5 0.7;
                        0.5 1.0 0.6;
                        0.7 0.6 1.0], active=false),

    unc3 = Correlation(corr_matrix)
    )

    model1 = EFTfitterModel(parameters1, measurements1, correlations1)

    measurement_models, measurement_keys = EFTfitter.measurement_models(model1)
    unc_models, unc_keys = EFTfitter.uncertainty_models(model1)

    @testset "Test model creation" begin
        measurement_models, measurement_keys = EFTfitter.measurement_models(model1)
        @test length(measurement_models) == 4
        @test measurement_keys == [:meas1, :meas3, :meas4_bin1, :meas4_bin3]
        @test keys(measurement_models[1].measurements) == (:meas3, :meas4_bin1, :meas4_bin3)
        @test keys(measurement_models[2].measurements) == (:meas1, :meas4_bin1, :meas4_bin3)

        unc_models, unc_keys = EFTfitter.uncertainty_models(model1)
        @test length(unc_models) == 2
        @test unc_keys == [:unc1, :unc3]
        @test keys(unc_models[1].correlations) == (:unc3,)
        @test keys(unc_models[2].correlations) == (:unc1,)
    end

    @testset "Test criterion_value" begin
        @test EFTfitter.criterion_value(SmallestInterval(key=:p1, p=0.9), unc_models[1],) ≈ 0.264 rtol = 0.05
        @test EFTfitter.criterion_value(SmallestInterval(key=:p2, p=0.9), unc_models[1],) ≈ 0.066 rtol = 0.05

        @test EFTfitter.criterion_value(SumOfSmallestIntervals(p=0.9), unc_models[1],) ≈ [0.264, 0.066] rtol = 0.05
    end

    @testset "Test ranking function" begin
        meas_ranks = EFTfitter.rank_measurements(model1, sampling_algorithm=SobolSampler(nsamples = 10^5))
        @test meas_ranks.names == [:meas4_bin1, :meas4_bin3, :meas1, :meas3]
        @test meas_ranks.values ≈ [0.179104478, 0.11940298, 0.0, 0.0] rtol=0.1

        unc_ranks = EFTfitter.rank_uncertainties(model1, sampling_algorithm=SobolSampler(nsamples = 10^5), order=:names, rev=false)
        @test unc_ranks.names == [:unc1, :unc3]
        @test unc_ranks.values ≈ [0.02985, 0.24627] rtol=0.1
    end
end


text:
module Kerrs

end


text:
using COSMO
using MarketData
using PortfolioOpt
using PortfolioOpt.TestUtils
using Test

DEFAULT_SOLVER = optimizer_with_attributes(
    COSMO.Optimizer, "verbose" => false, "max_iter" => 900000
)

@testset "PortfolioOpt.jl" begin
    # Write your tests here.
end


text:
module QuDirac

    using Compat
    using Iterators.product

    if !(v"0.3-" <= VERSION < v"0.4-")
        warn("QuDirac v0.1 only officially supports the v0.3 release of Julia. Your version of Julia is $VERSION.")
    end

    ####################
    # String Constants #
    ####################
    const lang = "\u27E8"
    const rang = "\u27E9"
    const otimes = "\u2297"
    const vdots ="\u205E"

    ##################
    # Abstract Types #
    ##################
    abstract AbstractInner

    immutable UndefinedInner <: AbstractInner end
    immutable KroneckerDelta <: AbstractInner end

    abstract AbstractDirac{P<:AbstractInner,N}
    abstract DiracOp{P,N} <: AbstractDirac{P,N}
    abstract DiracState{P,N} <: AbstractDirac{P,N}

    abstract DiracScalar <: Number

    #############
    # Functions #
    #############
    # These functions will be in
    # QuBase when it releases
    tensor() = error("Cannot call tensor function without arguments")
    tensor(s...) = reduce(tensor, s)

    ######################
    # Include Statements #
    ######################
    include("labels.jl")
    include("innerexpr.jl")
    include("state.jl")
    include("opsum.jl")
    include("outerproduct.jl")
    include("printfuncs.jl")
    include("dictfuncs.jl")
    include("mapfuncs.jl")
    include("str_macros.jl")

    #################
    # default_inner #
    #################
    # Julia doesn't recompile functions within other
    # previously compiled functions, so we can't have a
    # get_default_inner() or something like that.
    #
    # Also, global optimization is poor, so we don't want
    # to use that either. Thus, we go for a function that
    # straight-up redefines the default constructors for
    # the relevant objects. This is hacky, but works for now,
    # seeing as how only a few functions actually "use"
    # the default ptype.
    function default_inner(ptype::AbstractInner)
        QuDirac.OpSum(dict::Dict) = OpSum(ptype, dict)
        QuDirac.Ket(dict::Dict) = Ket(ptype, dict)
        QuDirac.ket(label::StateLabel) = ket(ptype, label)
        QuDirac.ket(items...) = ket(ptype, StateLabel(items))
        info("QuDirac's default inner product type is currently $ptype")
    end

    default_inner(KroneckerDelta());

    export AbstractInner,
        UndefinedInner,
        KroneckerDelta,
        default_inner,
        AbstractDirac,
        DiracState,
        DiracOp,
        # All functions that conflict
        # with QuBase should be exported
        # below:
        tensor,
        commute,
        anticommute,
        normalize,
        normalize!

end # module QuDirac


text:
using Javis

function object(r, p = O, color = "black")
    sethue(color)
    circle(p, r, :fill)
    return p
end

function connector!(connection, p1, p2, color)
    sethue(color)
    push!(connection, [p1, p2])
    map(x -> line(x[1], x[2], :stroke), connection)
end

function ground(args...)
    background("black") # canvas background
    sethue("white") # pen color
end

function circ(p = O, color = "black", action = :fill, radius = 25, edge = "solid")
    sethue(color)
    setdash(edge)
    circle(p, radius, action)
end

function make_animation()

    # to store the connectors
    connection = []

    frames = 1000

    # setup the video
    myvideo = Video(500, 500)
    Background(1:frames, ground)

    # add the objects
    earth = Object(1:frames, (args...) -> object(5, O, "blue"), Point(200, 0))
    venus = Object(1:frames, (args...) -> object(4, O, "red"), Point(144, 0))

    # draw the orbits
    earth_orbit = Object((args...) -> circ(O, "white", :stroke, 200))
    venus_orbit = Object((args...) -> circ(O, "white", :stroke, 144))

    # move the planets
    # We need the planets to revolve according to their time periods.
    # Earth completes its one revolution in 365 days and Venus does that in 224.7 days.
    # Hence, we need to multiply (224.7/365) so that the time period matches properly i.e.,
    # when earth completes its full revolution, Venus has done (224.7/365) th of its revolution.
    act!(earth, Action(anim_rotate_around(12.5 * 2π, O)))
    act!(venus, Action(anim_rotate_around(12.5 * 2π * (224.7 / 365), O)))

    # draw the connectors
    Object(1:frames, (args...) -> connector!(connection, pos(earth), pos(venus), "#f05a4f"))

    # render
    render(myvideo, pathname = "cosmic_dance.gif")
end

make_animation()


text:
include("precompile_funcs.jl")

using FluxDash

using .PrecompileFuncs
PrecompileFuncs.precompile_nn_train()

app = FluxDash.make_app()

port = haskey(ENV, "PORT") ? parse(Int64, ENV["PORT"]) : 8050

FluxDash.Dash.run_server(app, "0.0.0.0", port)

max_stars_repo_name: DEICIDIA/fluxdash

max_issues_repo_name: DEICIDIA/fluxdash

max_forks_repo_name: DEICIDIA/fluxdash


text:
using Plots

x = -pi/2:0.1:pi/2

plot(x, sin.(x))
plot!(x, cos.(x))
savefig("trigonometric.png")


text:
module ConfigurationInteraction

using Fermi
using Fermi.Geometry: Molecule
using Fermi.Integrals: IntegralHelper
using Fermi.Output
using TensorOperations
using LinearAlgebra

function print_header()
    banner =
raw"""
================================================================================
     ____                __  _                            _    _
    / ___| ___   _ __   / _|(_)  __ _  _   _  _ __  __ _ | |_ (_)  ___   _ __
   | |    / _ \ | '_ \ | |_ | | / _` || | | || '__|/ _` || __|| | / _ \ | '_ \
   | |___| (_) || | | ||  _|| || (_| || |_| || |  | (_| || |_ | || (_) || | | |
    \____|\___/ |_| |_||_|  |_| \__, | \__,_||_|   \__,_| \__||_| \___/ |_| |_|
    ___         _               |___/       _    _
   |_ _| _ __  | |_  ___  _ __  __ _   ___ | |_ (_)  ___   _ __
    | | | '_ \ | __|/ _ \| '__|/ _` | / __|| __|| | / _ \ | '_ \
    | | | | | || |_|  __/| |  | (_| || (__ | |_ | || (_) || | | |
   |___||_| |_| \__|\___||_|   \__,_| \___| \__||_| \___/ |_| |_|

             Module by G.J.R. Aroeira and M.M. Davis
================================================================================
"""
    @output "\n{}\n" banner
end

abstract type AbstractCIWavefunction <: Fermi.AbstractCorrelatedWavefunction end

# Struct symbolizing the type of implementation for different CI methods
abstract type CIAlgorithm end
struct SparseHamiltonian <: CIAlgorithm end
struct ACI <: CIAlgorithm end

include("CASCI/CASCI.jl")

end #module CI


text:
using BuildingGeometry
using Test

@testset "BuildingGeometry.jl" begin
    # Write your tests here.
end


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

JLLWrappers.@generate_main_file_header("MariaDB_Connector_ODBC")
JLLWrappers.@generate_main_file("MariaDB_Connector_ODBC", UUID("1124d457-72a3-5205-9757-36a127c7d0f8"))
end  # module MariaDB_Connector_ODBC_jll


text:
using Test, Random, LinearAlgebra, YaoArrayRegister, YaoBlocks, YaoBase

A = rand(ComplexF64, 4, 4)
mg = matblock(A)
mg2 = copy(mg)
@test mg2 == mg
@test mat(ComplexF64, mg) ≈ A
@test_logs (:warn, "converting Complex{Float64} to eltype Complex{Float32}, consider create another matblock with eltype Complex{Float32}") mat(ComplexF32, mg)

mg2.mat[:, 2] .= 10
@test mg2 != mg
@test nqubits(mg) == 2
@test_throws DimensionMismatch matblock(randn(3,3))

reg = rand_state(2)
@test apply!(copy(reg), mg) |> statevec == mat(mg) * reg.state |> vec

@test matblock(X) == GeneralMatrixBlock(mat(X))

a = rand_unitary(2)
@test mat(matblock(a))' == mat(matblock(a)')


text:
@doc raw"""
    approx_besselroots(ν::Real, n::Integer) -> Vector{Float64}

Return the first ``n`` roots of [Bessel function](https://en.wikipedia.org/wiki/Bessel_function).
Note that this function is only 12-digits of precision.

```math
J_{\nu}(x) = \sum_{m=0}^{\infty}\frac{(-1)^j}{\Gamma(\nu+j+1)j!} \left(\frac{x}{2}\right)^{2j+\nu}
```

# Examples
```jldoctest
julia> ν = 0.3;

julia> roots = approx_besselroots(ν, 10);

julia> zeros = (x -> besselj(ν, x)).(roots);

julia> all(zeros .< 1e-12)
true
```
"""
function approx_besselroots(ν::Real, n::Integer)
    # FIXME (related issue #22 and #80)
    return approx_besselroots(Float64(ν), n)
end

function approx_besselroots(ν::Float64, n::Integer)
# DEVELOPERS NOTES:
#   ν = 0 --> Full Float64 precision for n ≤ 20 (Wolfram Alpha), and very
#     accurate approximations for n > 20 (McMahon's expansion)
#   -1 ≤ ν ≤ 5 : ν ~= 0 --> 12 decimal figures for the 6 first zeros
#     (Piessens's Chebyshev series approximations), and very accurate
#     approximations for the others (McMahon's expansion)
#   ν > 5 --> moderately accurate for the 6 first zeros and good
#     approximations for the others (McMahon's expansion)

# This code was originally written by L. L. Peixoto in MATLAB.
# Later modified by A. Townsend to work in Julia

    if n < 0
        throw(DomainError(n, "Input n must be a non-negative integer"))
    end

    x = zeros(n)
    if ν == 0
        for k in 1:min(n,20)
            x[k] = BESSELJ0_ROOTS[k]
        end
        for k in min(n,20)+1:n
            x[k] = McMahon(ν, k)
        end
    elseif -1 ≤ ν ≤ 5
        correctFirstFew = piessens(ν)
        for k in 1:min(n,6)
            x[k] = correctFirstFew[k]
        end
        for k in min(n,6)+1:n
            x[k] = McMahon(ν, k)
        end
    elseif 5 < ν
        for k in 1:n
            x[k] = McMahon(ν, k)
        end
    end
    return x
end

function McMahon(ν::Real, k::Integer)
    # FIXME (related issue #22 and #80)
    return McMahon(Float64(ν), k)
end

function McMahon(ν::Float64, k::Integer)
    # McMahon's expansion. This expansion gives very accurate approximation
    # for the sth zero (s ≥ 7) in the whole region ν ≥ -1, and moderate
    # approximation in other cases.
    μ = 4ν^2
    a1 = 1 / 8
    a3 = (7μ-31) / 384
    a5 = 4*(3779+μ*(-982+83μ)) / 61440 # Evaluate via Horner's method.
    a7 = 6*(-6277237+μ*(1585743+μ*(-153855+6949μ))) / 20643840
    a9 = 144*(2092163573+μ*(-512062548+μ*(48010494+μ*(-2479316+70197μ)))) / 11890851840
    a11 = 720 *(-8249725736393+μ*(1982611456181+μ*(-179289628602+μ*(8903961290 +
          μ*(-287149133 + 5592657μ))))) / 10463949619200
    a13 = 576 *(423748443625564327 + μ*(-100847472093088506 + μ*(8929489333108377 +
        μ*(-426353946885548+μ*(13172003634537+μ*(-291245357370 + 4148944183μ)))))) / 13059009124761600
    b = 0.25 * (2ν+4k-1)*π
    # Evaluate using Horner's scheme:
    x = b - (μ-1)*( ((((((a13/b^2 + a11)/b^2 + a9)/b^2 + a7)/b^2 + a5)/b^2 + a3)/b^2 + a1)/b)
    return x
end

function _piessens_chebyshev30(ν::Float64)
    # Piessens's Chebyshev series approximations (1984). Calculates the 6 first
    # zeros to at least 12 decimal figures in region -1 ≤ ν ≤ 5:
    pt = (ν-2)/3

    T1 = 1.0
    T2 = pt
    T3 = 2pt*T2 - T1
    T4 = 2pt*T3 - T2
    T5 = 2pt*T4 - T3
    T6 = 2pt*T5 - T4
    T7 = 2pt*T6 - T5
    T8 = 2pt*T7 - T6
    T9 = 2pt*T8 - T7
    T10 = 2pt*T9 - T8
    T11 = 2pt*T10 - T9
    T12 = 2pt*T11 - T10
    T13 = 2pt*T12 - T11
    T14 = 2pt*T13 - T12
    T15 = 2pt*T14 - T13
    T16 = 2pt*T15 - T14
    T17 = 2pt*T16 - T15
    T18 = 2pt*T17 - T16
    T19 = 2pt*T18 - T17
    T20 = 2pt*T19 - T18
    T21 = 2pt*T20 - T19
    T22 = 2pt*T21 - T20
    T23 = 2pt*T22 - T21
    T24 = 2pt*T23 - T22
    T25 = 2pt*T24 - T23
    T26 = 2pt*T25 - T24
    T27 = 2pt*T26 - T25
    T28 = 2pt*T27 - T26
    T29 = 2pt*T28 - T27
    T30 = 2pt*T29 - T28

    T = SVector(T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15,T16,T17,T18,T19,T20,T21,T22,T23,T24,T25,T26,T27,T28,T29,T30)
    return T
end

function piessens(ν::Float64)
    # Piessens's Chebyshev series approximations (1984). Calculates the 6 first
    # zeros to at least 12 decimal figures in region -1 ≤ ν ≤ 5:
    C = PIESSENS_C
    T = _piessens_chebyshev30(ν)
    y = C'*T
    _y = collect(y)
    _y[1] *= sqrt(ν+1)  # Scale the first root.
    return _y
end

function besselZeroRoots(m)
    # BESSEL0ROOTS ROOTS OF BESSELJ(0,x). USE ASYMPTOTICS.
    # Use McMahon's expansion for the remainder (NIST, 10.21.19):
    jk = Array{Float64}(undef, m)
    p = (1071187749376 / 315, 0.0, -401743168 / 105, 0.0, 120928 / 15,
         0.0, -124 / 3, 0.0, 1.0, 0.0)
    # First 20 are precomputed:
    @inbounds for jj = 1:min(m, 20)
        jk[jj] = BESSELJ0_ROOTS[jj]
    end
    @inbounds for jj = 21:min(m, 47)
        ak = π * (jj - .25)
        ak82 = (.125 / ak)^2
        jk[jj] = ak + .125 / ak * @evalpoly(ak82, 1.0, p[7], p[5], p[3])
    end
    @inbounds for jj = 48:min(m, 344)
        ak = π * (jj - .25)
        ak82 = (.125 / ak)^2
        jk[jj] = ak + .125 / ak * @evalpoly(ak82, 1.0, p[7], p[5])
    end
    @inbounds for jj = 345:min(m,13191)
        ak = π * (jj - .25)
        ak82 = (.125 / ak)^2
        jk[jj] = ak + .125 / ak * muladd(ak82, p[7], 1.0)
    end
    @inbounds for jj = 13192:m
        ak = π * (jj - .25)
        jk[jj] = ak + .125 / ak
    end
    return jk
end

function besselJ1(m)
    # BESSELJ1 EVALUATE BESSELJ(1,x)^2 AT ROOTS OF BESSELJ(0,x).
    # USE ASYMPTOTICS. Use Taylor series of (NIST, 10.17.3) and McMahon's
    # expansion (NIST, 10.21.19):
    Jk2 = Array{Float64}(undef, m)
    c = (-171497088497 / 15206400, 461797 / 1152, -172913 / 8064,
         151 / 80, -7 / 24, 0.0, 2.0)
    # First 10 are precomputed:
    @inbounds for jj = 1:min(m, 10)
        Jk2[jj] = BESSELJ1_ON_BESSELJ0_ROOTS[jj]
    end
    @inbounds for jj = 11:min(m, 15)
        ak = π * (jj - .25)
        ak2 = (1 / ak)^2
        Jk2[jj] = 1 / (π * ak) * muladd(@evalpoly(ak2, c[5], c[4], c[3],
                                                  c[2], c[1]), ak2^2, c[7])
    end
    @inbounds for jj = 16:min(m, 21)
        ak = π * (jj - .25)
        ak2 = (1 / ak)^2
        Jk2[jj] = 1 / (π * ak) * muladd(@evalpoly(ak2, c[5], c[4], c[3], c[2]),
                                        ak2^2, c[7])
    end
    @inbounds for jj = 22:min(m,55)
        ak = π * (jj - .25)
        ak2 = (1 / ak)^2
        Jk2[jj] = 1 / (π * ak) * muladd(@evalpoly(ak2, c[5], c[4], c[3]),
                                        ak2^2, c[7])
    end
    @inbounds for jj = 56:min(m,279)
        ak = π * (jj - .25)
        ak2 = (1 / ak)^2
        Jk2[jj] = 1 / (π * ak) * muladd(muladd(ak2, c[4], c[5]), ak2^2, c[7])
    end
    @inbounds for jj = 280:min(m,2279)
        ak = π * (jj - .25)
        ak2 = (1 / ak)^2
        Jk2[jj] = 1 / (π * ak) * muladd(ak2^2, c[5], c[7])
    end
    @inbounds for jj = 2280:m
        ak = π * (jj - .25)
        Jk2[jj] = 1 / (π * ak) * c[7]
    end
    return Jk2
end

function besselroots(ν::Real, n::Integer)
    @warn "`besselroots` has been renamed to `approx_besselroots`, and `besselroots` will be removed in the next major release."
    return approx_besselroots(ν, n)
end


text:
module TranscriptionOpt

using ..InfiniteOpt

# Load in the module files
include("model.jl")
include("measures.jl")
include("transcribe.jl")
include("optimize.jl")

# Export transcription datatypes
export TranscriptionData, TranscriptionModel

# Export transcription methods
export is_transcription_model, transcription_data, transcription_variable,
transcription_constraint, transcription_model, transcription_expression

end # end module


text:
## topological sort

"""
For each variable in graph, calculate all variables that depend on it.
This is essentially the opposite of `dependencies(nd::ExNode)`, but
operates on variable names rather than nodes.
"""
function dependents(g::AbstractExGraph)
    dpts = Dict{Symbol, Vector{Symbol}}()
    # prepare all dependent lists
    for nd in g.tape
        dpts[varname(nd)] = Symbol[]
    end
    for nd in g.tape
        vname = varname(nd)
        for dep in dependencies(nd)
            if haskey(dpts, dep)          # don't include global constants
                push!(dpts[dep], vname)
            end
        end
    end
    return dpts
end

function topsort_visit!(g::AbstractExGraph, dpts::Dict{Symbol, Vector{Symbol}},
                        temp_marked::Set{Symbol}, perm_marked::Set{Symbol},
                        sorted::Vector{Symbol}, vname::Symbol)
    if vname in temp_marked
        error("Expression graph isn't a DAG!")
    end
    if !in(vname, temp_marked) && !in(vname, perm_marked)
        push!(temp_marked, vname)
        for dpt in dpts[vname]
            topsort_visit!(g, dpts, temp_marked, perm_marked, sorted, dpt)
        end
        push!(perm_marked, vname)
        delete!(temp_marked, vname)
        push!(sorted, vname)
    end
end

"""Sort graph topologically"""
function topsort(g::AbstractExGraph)
    dpts = dependents(g)
    sorted = Symbol[]
    temp_marked = Set{Symbol}()
    perm_marked = Set{Symbol}()
    for vname in keys(dpts)
        topsort_visit!(g, dpts, temp_marked, perm_marked, sorted, vname)
    end
    sg  = Espresso.reset_tape(g)
    for vname in reverse(sorted)
        push!(sg, g[vname])
    end
    return sg
end

## expand const

"""Expand all constant vars in a given expression"""
function expand_const(g::AbstractExGraph, ex)
    st = Dict{Symbol, Any}()
    vnames = get_var_names(ex)
    for vname in vnames
        if haskey(g, vname) && isa(g[vname], ExNode{:constant})
            # st[vname] = g[vname].val
            st[vname] = getvalue(g[vname])
        end
    end
    return subs(ex, st)
end

reindex_from_beginning(g::ExGraph) = g

## inline / inline unknown

function subgraph_interm_subs_table(sub_g::ExGraph, dont_subs; prefix="")
    interm_vars = [varname(sub_nd) for sub_nd in sub_g.tape
                   if !in(varname(sub_nd), dont_subs)]
    new_names = [genname("$(prefix)_$(v)_") for v in interm_vars]
    return Dict(zip(interm_vars, new_names))
end

"""
Find definition of a called function and build its subgraph ready for inlining
"""
function make_subgraph(g::ExGraph, nd::ExNode{:call})
    mod = @get(g.ctx, :mod, Main)  # TODO: Main or current_graph()?
    fname = getexpr(nd).args[1]
    f = Core.eval(mod, fname)
    args = dependencies(nd)
    arg_types = ([typeof(getvalue(g[arg])) for arg in args]...,)
    params, sub_ex = funexpr(f, arg_types)
    sub_g = ExGraph(sub_ex; ctx=g.ctx)
    st = Dict(zip(params, args))           # rename internal params to actual arguments
    st[varname(sub_g[end])] = varname(nd)  # rename output var to this node's
    dont_subs = Set(keys(st))
    st = merge(st, subgraph_interm_subs_table(sub_g, dont_subs; prefix=fname))
    rename!(sub_g, st)
    return sub_g
end

function inline_subgraphs(g::ExGraph, inline_subs::Dict{Symbol, ExGraph})
    new_g = reset_tape(g)
    for nd in g
        vname = varname(nd)
        if haskey(inline_subs, vname)
            for sub_nd in inline_subs[vname]
                push!(new_g, sub_nd)
            end
        else
            push!(new_g, nd)
        end
    end
    return new_g
end

function inline_nodes(g::ExGraph, vnames::Set{Symbol})
    inline_subs = Dict(vname => make_subgraph(g, g[vname]) for vname in vnames)
    return inline_subgraphs(g, inline_subs)
end

iscall(x) = isa(x, Expr) && x.head == :call

function convert_call(g::AbstractExGraph, nd::Union{ExNode{:call}, ExNode{:bcast}})
    new_ex = expand_const(g, getexpr(nd)) |> simplify
    if isa(new_ex, Symbol) || (isa(new_ex, Expr) && new_ex.head == :ref)
        # convert to assignment
        return copy(nd; category=:(=), ex=new_ex)
    elseif isa(new_ex, Number) || isa(new_ex, AbstractArray)
        # convert to constant
        if haskey(g.ctx, :bitness)
            new_ex = force_bitness(new_ex, Val(g.ctx[:bitness]))
        end
        return copy(nd; category=:constant, ex=new_ex)
    elseif isa(new_ex, Tuple)
        return copy(nd; category=:tuple, ex=new_ex)
    else
        error("Call node $nd is simplified to an unknown non-call $new_ex")
    end
end

"""
Given a symbolic name, either adds `2` to the end
or increment existing number. Example:

    inc_var_name(:x)   # ==> x2
    inc_var_name(:x2)  # ==> x3

"""
function inc_var_name(var::Symbol)
    svar = string(var)
    r = match(r"(.*)(\d+)", svar)
    if r != nothing
        base = r.captures[1]
        num = parse(Int, r.captures[2]) + 1
        return Symbol("$base$num")
    else
        return Symbol("$(svar)2")
    end
end

function rename_from_beginning(g::ExGraph)
    new_g = deepcopy(g)
    vars = [getvar(nd) for nd in g]
    new_vars = [Symbol("var$i") for i=1:length(vars)]
    st = Dict(zip(vars, new_vars))
    rename!(new_g, st)
    return new_g
end

"""
A single number to represent a graph.
Insensitive to variable names.
"""
function graph_hash(g::ExGraph)
    return rename_from_beginning(g) |> to_expr |> hash
end


text:
stream = IOBuffer("1\n2\n3\n4\n\n6")

while !eof(stream)
    line = readline(stream)
    println(line)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# ----------------------------------------------------------
# ---------------  Build computational tree  ---------------
# ----------------------------------------------------------
# -- Inner nodes are FunctionOperatorComposite objects, ----
# --------- leaves are FunctionOperator objects. -----------
# ----------------------------------------------------------

# +, -, and * operators combines the mappings on the left and right side into a
# FunctionOperatorComposite object. Assertions on type and size made when applicable.

Base.:+(FO1::FunOp, FO2::FunOp) = begin
    assertType(FO1, FO2)
    assertAddDim(FO1, FO2)
    FunctionOperatorComposite(FO1, FO2,  :+)
end

Base.:+(FO::FunOp, S::LinearAlgebra.UniformScaling) = begin
    assertAddDimScaling(FO, S)
    FunctionOperatorComposite(FO, createScalingForAddSub(FO, S),  :+)
end

Base.:+(S::LinearAlgebra.UniformScaling, FO::FunOp) = begin
    assertAddDimScaling(FO, S)
    FunctionOperatorComposite(createScalingForAddSub(FO, S), FO,  :+)
end

Base.:-(FO1::FunOp, FO2::FunOp) = begin
    assertType(FO1, FO2)
    assertAddDim(FO1, FO2)
    FunctionOperatorComposite(FO1, FO2,  :-)
end

Base.:-(FO::FunOp, S::LinearAlgebra.UniformScaling) = begin
    assertAddDimScaling(FO, S)
    FunctionOperatorComposite(FO, createScalingForAddSub(FO, S),  :-)
end

Base.:-(S::LinearAlgebra.UniformScaling, FO::FunOp) =  begin
    assertAddDimScaling(FO, S)
    FunctionOperatorComposite(createScalingForAddSub(FO, S), FO,  :-)
end

Base.:*(FO1::FunOp, FO2::FunOp) = begin
    assertType(FO1,FO2)
    assertMultDim(FO1, FO2)
    FunctionOperatorComposite(FO1, FO2,  :*)
end

Base.:*(FO::FunctionOperator, S::LinearAlgebra.UniformScaling{Bool}) =
    FunctionOperator(FO, name = getName(FO) * " * I")

Base.:*(FO::FunctionOperatorComposite, S::LinearAlgebra.UniformScaling{Bool}) =
    FunctionOperatorComposite(FO, name = getName(FO) * " * I")

Base.:*(S::LinearAlgebra.UniformScaling{Bool}, FO::FunctionOperator) =
    FunctionOperator(FO, name = "I * " * getName(FO))

Base.:*(S::LinearAlgebra.UniformScaling{Bool}, FO::FunctionOperatorComposite) =
    FunctionOperatorComposite(FO, name = "I * " * getName(FO))

Base.:*(FO::FunOp, S::LinearAlgebra.UniformScaling) =
    FunctionOperatorComposite(FO, createScalingForMult(FO, S, FO.inDims),  :*)

Base.:*(S::LinearAlgebra.UniformScaling, FO::FunOp) =
    FunctionOperatorComposite(createScalingForMult(FO, S, FO.outDims), FO,  :*)

Base.:*(FO::FunOp, λ::Number) = FO * (λ*I)

Base.:*(λ::Number, FO::FunOp) = (λ*I) * FO

# Adjoint operator creates a new FunctionOperatorComposite object, toggles the adjoint field and
# switches the input and output dimension constraints (and also voids plan for FunctionOperatorComposite)

Base.:adjoint(FO::FunctionOperator) =
    FunctionOperator(FO, adjoint = !FO.adjoint, inDims = FO.outDims, outDims = FO.inDims)

Base.:adjoint(FO::FunctionOperatorComposite{T}) where {T} = begin
    FunctionOperatorComposite(FO, name = "("*getName(FO)*")'", adjoint = !FO.adjoint,
        inDims = FO.outDims, outDims = FO.inDims, plan_function = noplan)
end


text:
#!/usr/bin/env julia

using Base.Test
using TestSetExtensions
using PortAudio
using SampledSignals
using RingBuffers

# pull in some extra stuff we need to test the callback directly
using PortAudio: notifyhandle, notifycb_c, shim_processcb_c
using PortAudio: pa_shim_errmsg_t, pa_shim_info_t
using PortAudio: PA_SHIM_ERRMSG_ERR_OVERFLOW, PA_SHIM_ERRMSG_UNDERFLOW, PA_SHIM_ERRMSG_OVERFLOW

"Setup buffers to test callback behavior"
function setup_callback(inchans, outchans, nframes, synced)
    sourcebuf = RingBuffer{Float32}(inchans, nframes*2) # the microphone input should end up here
    sinkbuf = RingBuffer{Float32}(outchans, nframes*2) # the callback should copy this to cb_output
    errbuf = RingBuffer{pa_shim_errmsg_t}(1, 8)

    # pass NULL for i/o we're not using
    info = pa_shim_info_t(
        inchans > 0 ? pointer(sourcebuf) : C_NULL,
        outchans > 0 ? pointer(sinkbuf) : C_NULL,
        pointer(errbuf),
        synced, notifycb_c,
        inchans > 0 ? notifyhandle(sourcebuf) : C_NULL,
        outchans > 0 ? notifyhandle(sinkbuf) : C_NULL,
        notifyhandle(errbuf)
    )
    flags = Culong(0)

    cb_input = rand(Float32, inchans, nframes) # simulate microphone input
    cb_output = rand(Float32, outchans, nframes) # this is where the output should go

    function processfunc()
        ccall(shim_processcb_c, Cint,
            (Ptr{Float32}, Ptr{Float32}, Culong, Ptr{Void}, Culong, Ptr{Void}),
            cb_input, cb_output, nframes, C_NULL, flags, pointer_from_objref(info))
    end

    (sourcebuf, sinkbuf, errbuf, cb_input, cb_output, processfunc)
end

function test_callback(inchans, outchans, synced)
    nframes = 8
    (sourcebuf, sinkbuf, errbuf,
     cb_input, cb_output, process) = setup_callback(inchans, outchans,
                                                    nframes, synced)
    if outchans > 0
        testout = rand(Float32, outchans, nframes) # generate some test data to play
        write(sinkbuf, testout) # fill the output ringbuffer
    end
    @test process() == PortAudio.paContinue
    if outchans > 0
        # testout -> sinkbuf -> cb_output
        @test cb_output == testout
    end
    if inchans > 0
        # cb_input -> sourcebuf
        @test read(sourcebuf, nframes) == cb_input
    end
    @test framesreadable(errbuf) == 0
end

"""
    test_callback_underflow(inchans, outchans; nframes=8, underfill=3, synced=false)

Test that the callback works on underflow conditions. underfill is the numer of
frames we feed in, which should be less than nframes.
"""
function test_callback_underflow(inchans, outchans, synced)
    nframes = 8
    underfill = 3 # must be less than nframes
    (sourcebuf, sinkbuf, errbuf,
     cb_input, cb_output, process) = setup_callback(inchans, outchans,
                                                    nframes, synced)
    outchans > 0 || error("Can't test underflow with no output")
    testout = rand(Float32, outchans, underfill)
    write(sinkbuf, testout) # underfill the output ringbuffer
    # call callback (partial underflow)
    @test process() == PortAudio.paContinue
    @test cb_output[:, 1:underfill] == testout
    @test cb_output[:, (underfill+1):nframes] == zeros(Float32, outchans, (nframes-underfill))
    errs = readavailable(errbuf)
    if inchans > 0
        received = readavailable(sourcebuf)
        if synced
            @test size(received, 2) == underfill
            @test received == cb_input[:, 1:underfill]
            @test length(errs) == 2
            @test Set(errs) == Set([PA_SHIM_ERRMSG_UNDERFLOW, PA_SHIM_ERRMSG_OVERFLOW])
        else
            @test size(received, 2) == nframes
            @test received == cb_input
            @test length(errs) == 1
            @test errs[1] == PA_SHIM_ERRMSG_UNDERFLOW
        end
    else
        @test length(errs) == 1
        @test errs[1] == PA_SHIM_ERRMSG_UNDERFLOW
    end

    # call again (total underflow)
    @test process() == PortAudio.paContinue
    @test cb_output == zeros(Float32, outchans, nframes)
    errs = readavailable(errbuf)
    if inchans > 0
        received = readavailable(sourcebuf)
        if synced
            @test size(received, 2) == 0
            @test length(errs) == 2
            @test Set(errs) == Set([PA_SHIM_ERRMSG_UNDERFLOW, PA_SHIM_ERRMSG_OVERFLOW])
        else
            @test size(received, 2) == nframes
            @test received == cb_input
            @test length(errs) == 1
            @test errs[1] == PA_SHIM_ERRMSG_UNDERFLOW
        end
    else
        @test length(errs) == 1
        @test errs[1] == PA_SHIM_ERRMSG_UNDERFLOW
    end
end

function test_callback_overflow(inchans, outchans, synced)
    nframes = 8
    (sourcebuf, sinkbuf, errbuf,
     cb_input, cb_output, process) = setup_callback(inchans, outchans,
                                                    nframes, synced)
    inchans > 0 || error("Can't test overflow with no input")
    @test frameswritable(sinkbuf) == nframes*2

    # the first time it should half-fill the input ring buffer
    if outchans > 0
        testout = rand(Float32, outchans, nframes)
        write(sinkbuf, testout)
    end
    @test framesreadable(sourcebuf) == 0
    outchans > 0 && @test frameswritable(sinkbuf) == nframes
    @test process() == PortAudio.paContinue
    @test framesreadable(errbuf) == 0
    @test framesreadable(sourcebuf) == nframes
    outchans > 0 && @test frameswritable(sinkbuf) == nframes*2

    # now run the process func again to completely fill the input ring buffer
    outchans > 0 && write(sinkbuf, testout)
    @test framesreadable(sourcebuf) == nframes
    outchans > 0 && @test frameswritable(sinkbuf) == nframes
    @test process() == PortAudio.paContinue
    @test framesreadable(errbuf) == 0
    @test framesreadable(sourcebuf) == nframes*2
    outchans > 0 && @test frameswritable(sinkbuf) == nframes*2

    # now this time the process func should overflow the input buffer
    outchans > 0 && write(sinkbuf, testout)
    @test framesreadable(sourcebuf) == nframes*2
    outchans > 0 && @test frameswritable(sinkbuf) == nframes
    @test process() == PortAudio.paContinue
    @test framesreadable(sourcebuf) == nframes*2
    errs = readavailable(errbuf)
    if outchans > 0
        if synced
            # if input and output are synced, thec callback didn't pull from
            # the output ringbuf
            @test frameswritable(sinkbuf) == nframes
            @test cb_output == zeros(Float32, outchans, nframes)
            @test length(errs) == 2
            @test Set(errs) == Set([PA_SHIM_ERRMSG_UNDERFLOW, PA_SHIM_ERRMSG_OVERFLOW])
        else
            @test frameswritable(sinkbuf) == nframes*2
            @test length(errs) == 1
            @test errs[1] == PA_SHIM_ERRMSG_OVERFLOW
        end
    else
        @test length(errs) == 1
        @test errs[1] == PA_SHIM_ERRMSG_OVERFLOW
    end
end

@testset ExtendedTestSet "PortAudio Tests" begin
    @testset "Reports version" begin
        io = IOBuffer()
        PortAudio.versioninfo(io)
        result = split(String(take!((io))), "\n")
        # make sure this is the same version I tested with
        @test startswith(result[1], "PortAudio V19")
    end

    @testset "using correct shim version" begin
        @test PortAudio.shimhash() == "87021557a9f999545828eb11e4ebad2cd278b734dd91a8bd3faf05c89912cf80"
    end

    @testset "Basic callback functionality" begin
        @testset "basic duplex (no sync)" begin
            test_callback(2, 3, false)
        end
        @testset "basic input-only (no sync)" begin
            test_callback(2, 0, false)
        end
        @testset "basic output-only (no sync)" begin
            test_callback(0, 2, false)
        end
        @testset "basic no input or output (no sync)" begin
            test_callback(0, 0, false)
        end
        @testset "basic duplex (sync)" begin
            test_callback(2, 3, true)
        end
        @testset "basic input-only (sync)" begin
            test_callback(2, 0, true)
        end
        @testset "basic output-only (sync)" begin
            test_callback(0, 2, true)
        end
        @testset "basic no input or output (sync)" begin
            test_callback(0, 0, true)
        end
    end

    @testset "Ouput underflow" begin
        @testset "underflow duplex (nosync)" begin
            test_callback_underflow(2, 3, false)
        end
        @testset "underflow output-only (nosync)" begin
            test_callback_underflow(0, 3, false)
        end
        @testset "underflow duplex (sync)" begin
            test_callback_underflow(2, 3, true)
        end
        @testset "underflow output-only (sync)" begin
            test_callback_underflow(0, 3, true)
        end
    end

    @testset "Input overflow" begin
        @testset "overflow duplex (nosync)" begin
            test_callback_overflow(2, 3, false)
        end
        @testset "overflow input-only (nosync)" begin
            test_callback_overflow(2, 0, false)
        end
        @testset "overflow duplex (sync)" begin
            test_callback_overflow(2, 3, true)
        end
        @testset "overflow input-only (sync)" begin
            test_callback_overflow(2, 0, true)
        end
    end
end


text:
# Dunavant, D. A. (1985), High degree efficient symmetrical Gaussian quadrature
# rules for the triangle. Int. J. Numer. Meth. Engng., 21: 1129–1148. doi:
# 10.1002/nme.1620210612
function _get_gauss_tridata(n::Int)
    if (n == 1)
    xw=[0.33333333333333 0.33333333333333 1.00000000000000 / 2.0];
    elseif (n == 2)
    xw=[0.16666666666667 0.16666666666667 0.33333333333333 / 2.0
    0.16666666666667 0.66666666666667 0.33333333333333 / 2.0
    0.66666666666667 0.16666666666667 0.33333333333333 / 2.0];
    elseif (n == 3)
    xw=[0.33333333333333 0.33333333333333 -0.56250000000000 / 2.0
    0.20000000000000 0.20000000000000 0.52083333333333 / 2.0
    0.20000000000000 0.60000000000000 0.52083333333333 / 2.0
    0.60000000000000 0.20000000000000 0.52083333333333 / 2.0];
    elseif (n == 4)
    xw=[0.44594849091597 0.44594849091597 0.22338158967801 / 2.0
    0.44594849091597 0.10810301816807 0.22338158967801 / 2.0
    0.10810301816807 0.44594849091597 0.22338158967801 / 2.0
    0.09157621350977 0.09157621350977 0.10995174365532 / 2.0
    0.09157621350977 0.81684757298046 0.10995174365532 / 2.0
    0.81684757298046 0.09157621350977 0.10995174365532 / 2.0];
    elseif (n == 5)
    xw=[0.33333333333333 0.33333333333333 0.22500000000000 / 2.0
    0.47014206410511 0.47014206410511 0.13239415278851 / 2.0
    0.47014206410511 0.05971587178977 0.13239415278851 / 2.0
    0.05971587178977 0.47014206410511 0.13239415278851 / 2.0
    0.10128650732346 0.10128650732346 0.12593918054483 / 2.0
    0.10128650732346 0.79742698535309 0.12593918054483 / 2.0
    0.79742698535309 0.10128650732346 0.12593918054483 / 2.0];
    elseif (n == 6)
    xw=[0.24928674517091 0.24928674517091 0.11678627572638 / 2.0
    0.24928674517091 0.50142650965818 0.11678627572638 / 2.0
    0.50142650965818 0.24928674517091 0.11678627572638 / 2.0
    0.06308901449150 0.06308901449150 0.05084490637021 / 2.0
    0.06308901449150 0.87382197101700 0.05084490637021 / 2.0
    0.87382197101700 0.06308901449150 0.05084490637021 / 2.0
    0.31035245103378 0.63650249912140 0.08285107561837 / 2.0
    0.63650249912140 0.05314504984482 0.08285107561837 / 2.0
    0.05314504984482 0.31035245103378 0.08285107561837 / 2.0
    0.63650249912140 0.31035245103378 0.08285107561837 / 2.0
    0.31035245103378 0.05314504984482 0.08285107561837 / 2.0
    0.05314504984482 0.63650249912140 0.08285107561837 / 2.0];
    elseif (n == 7)
    xw=[0.33333333333333 0.33333333333333 -0.14957004446768 / 2.0
    0.26034596607904 0.26034596607904 0.17561525743321 / 2.0
    0.26034596607904 0.47930806784192 0.17561525743321 / 2.0
    0.47930806784192 0.26034596607904 0.17561525743321 / 2.0
    0.06513010290222 0.06513010290222 0.05334723560884 / 2.0
    0.06513010290222 0.86973979419557 0.05334723560884 / 2.0
    0.86973979419557 0.06513010290222 0.05334723560884 / 2.0
    0.31286549600487 0.63844418856981 0.07711376089026 / 2.0
    0.63844418856981 0.04869031542532 0.07711376089026 / 2.0
    0.04869031542532 0.31286549600487 0.07711376089026 / 2.0
    0.63844418856981 0.31286549600487 0.07711376089026 / 2.0
    0.31286549600487 0.04869031542532 0.07711376089026 / 2.0
    0.04869031542532 0.63844418856981 0.07711376089026 / 2.0];
    elseif (n == 8)
    xw=[0.33333333333333 0.33333333333333 0.14431560767779 / 2.0
    0.45929258829272 0.45929258829272 0.09509163426728 / 2.0
    0.45929258829272 0.08141482341455 0.09509163426728 / 2.0
    0.08141482341455 0.45929258829272 0.09509163426728 / 2.0
    0.17056930775176 0.17056930775176 0.10321737053472 / 2.0
    0.17056930775176 0.65886138449648 0.10321737053472 / 2.0
    0.65886138449648 0.17056930775176 0.10321737053472 / 2.0
    0.05054722831703 0.05054722831703 0.03245849762320 / 2.0
    0.05054722831703 0.89890554336594 0.03245849762320 / 2.0
    0.89890554336594 0.05054722831703 0.03245849762320 / 2.0
    0.26311282963464 0.72849239295540 0.02723031417443 / 2.0
    0.72849239295540 0.00839477740996 0.02723031417443 / 2.0
    0.00839477740996 0.26311282963464 0.02723031417443 / 2.0
    0.72849239295540 0.26311282963464 0.02723031417443 / 2.0
    0.26311282963464 0.00839477740996 0.02723031417443 / 2.0
    0.00839477740996 0.72849239295540 0.02723031417443 / 2.0];
    else
        throw(ArgumentError("unsupported order for triangle gauss-legendre integration"))
    end
    return xw
end


text:
# # Group Isomorphism

# ## Set up Basic Functions

using TightBindingLattice
using Formatting
function display_matrix(io::IO, matrix::AbstractMatrix; prefix::AbstractString="")
    width = ceil(Int, maximum(length("$item") for item in matrix)/4)*4
    for row in eachrow(matrix)
        for (icol, col) in enumerate(row)
            print(io, icol == 1 ? prefix : " ")
            printfmt(io, "{:>$(width)s}", "$col")
        end
        println(io)
    end
end

# ## Generate Two Isomorphic Groups

group1 = FiniteGroup([
    1 2 3 4;
    2 1 4 3;
    3 4 2 1;
    4 3 1 2;
]);
group2 = let mtab1 = group1.multiplication_table,
             mtab2 = zeros(Int, (4, 4)),
             mapping = [1, 3, 2, 4]
             for i in 1:4, j in 1:4
                 mtab2[mapping[i], mapping[j]] = mapping[mtab1[i,j]]
             end
             FiniteGroup(mtab2)
         end;

# ## Group Multiplication Tables

println("Multiplication Table of G₁")
println("--------------------------")
display_matrix(stdout, group_multiplication_table(group1))
println()

println("Multiplication Table of G₂")
println("--------------------------")
display_matrix(stdout, group_multiplication_table(group2))
println()

# ## Group Isomorphism

println("Group Isomorphism  ϕ: G₁ → G₂")
println("-----------------------------")
ϕ = group_isomorphism(group1, group2)
for g in 1:group_order(group1)
    println("  ϕ($g) = $(ϕ[g])")
end
println()

mtab2 = zeros(Int, (group_order(group1), group_order(group1)))
for g in 1:group_order(group1), h in 1:group_order(group1)
    mtab2[ϕ[g], ϕ[h]] = ϕ[ group_product(group1, g, h) ]
end

println("Multiplication Table of ϕ(G₁)")
println("-----------------------------")
display_matrix(stdout, mtab2)
println()
println("  ϕ(g)⋅ϕ(h) = ϕ(g⋅h)")
println("  ϕ(G₁) ≡ G₂ should hold.")


text:
module FefferyAntdCharts
using Dash

const resources_path = realpath(joinpath( @__DIR__, "..", "deps"))
const version = "0.0.1-rc4"

include("jl/''_antdarea.jl")
include("jl/''_antdbar.jl")
include("jl/''_antdchord.jl")
include("jl/''_antdcolumn.jl")
include("jl/''_antdline.jl")
include("jl/''_antdpie.jl")
include("jl/''_antdradar.jl")
include("jl/''_antdscatter.jl")
include("jl/''_antdstock.jl")
include("jl/''_antdsunburst.jl")
include("jl/''_antdwordcloud.jl")

function __init__()
    DashBase.register_package(
        DashBase.ResourcePkg(
            "feffery_antd_charts",
            resources_path,
            version = version,
            [
                DashBase.Resource(
    relative_package_path = "feffery_antd_charts.min.js",
    external_url = "https://unpkg.com/feffery_antd_charts@0.0.1-rc4/feffery_antd_charts/feffery_antd_charts.min.js",
    dynamic = nothing,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "feffery_antd_charts.min.js.map",
    external_url = "https://unpkg.com/feffery_antd_charts@0.0.1-rc4/feffery_antd_charts/feffery_antd_charts.min.js.map",
    dynamic = true,
    async = nothing,
    type = :js
)
            ]
        )

    )
end
end

max_stars_repo_name: CNFeffery/feffery-antd-charts

max_issues_repo_name: CNFeffery/feffery-antd-charts

max_forks_repo_name: CNFeffery/feffery-antd-charts


text:
# Autogenerated wrapper script for Clipper_jll for x86_64-apple-darwin
export libcclipper

JLLWrappers.@generate_wrapper_header("Clipper")
JLLWrappers.@declare_library_product(libcclipper, "@rpath/libcclipper.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcclipper,
        "lib/libcclipper.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using StatsPlots

function plotdistributions(names, xs; kwargs...)
    names_ = [string(i, ": ", name) for (i, name) in enumerate(names)]
    labels = map((names, g) -> fill(names, size(g)), names, xs);
    violin(labels, xs; lw=0, xrotation=20, alpha=.8, legend=false, kwargs...)
    boxplot!(labels, xs, alpha=.25)
end

function plotlayers(model, lossfn, batch; filter_fn = l -> true, sz = (900, 600), kwargs...)
    names, acts, grads = layerstats(model, lossfn, batch, filter_fn = filter_fn)
    gradsplot = plotdistributions(names, grads, title = "gradients")
    actsplot = plotdistributions(names, acts, title = "activations")
    l = @layout [a; b]
    plot(actsplot, gradsplot, layout = l, size = sz, kwargs...)
end


text:
using Oceananigans.Fields: ZeroField

#####
##### Momentum advection operators
#####

# Alternate names for advective fluxes
@inline _advective_momentum_flux_Uu(args...) = advective_momentum_flux_Uu(args...)
@inline _advective_momentum_flux_Vu(args...) = advective_momentum_flux_Vu(args...)
@inline _advective_momentum_flux_Wu(args...) = advective_momentum_flux_Wu(args...)

@inline _advective_momentum_flux_Uv(args...) = advective_momentum_flux_Uv(args...)
@inline _advective_momentum_flux_Vv(args...) = advective_momentum_flux_Vv(args...)
@inline _advective_momentum_flux_Wv(args...) = advective_momentum_flux_Wv(args...)

@inline _advective_momentum_flux_Uw(args...) = advective_momentum_flux_Uw(args...)
@inline _advective_momentum_flux_Vw(args...) = advective_momentum_flux_Vw(args...)
@inline _advective_momentum_flux_Ww(args...) = advective_momentum_flux_Ww(args...)

@inline _advective_tracer_flux_x(args...) = _advective_tracer_flux_x(args...)
@inline _advective_tracer_flux_y(args...) = _advective_tracer_flux_y(args...)
@inline _advective_tracer_flux_z(args...) = _advective_tracer_flux_z(args...)

const ZeroU = NamedTuple{(:u, :v, :w), Tuple{ZeroField, ZeroField, ZeroField}}

# Compiler hints
@inline div_𝐯u(i, j, k, grid, advection, ::ZeroU, u) = zero(eltype(grid))
@inline div_𝐯v(i, j, k, grid, advection, ::ZeroU, v) = zero(eltype(grid))
@inline div_𝐯w(i, j, k, grid, advection, ::ZeroU, w) = zero(eltype(grid))

@inline div_𝐯u(i, j, k, grid, advection, U, ::ZeroField) = zero(eltype(grid))
@inline div_𝐯v(i, j, k, grid, advection, U, ::ZeroField) = zero(eltype(grid))
@inline div_𝐯w(i, j, k, grid, advection, U, ::ZeroField) = zero(eltype(grid))

@inline div_𝐯u(i, j, k, grid, ::Nothing, U, u) = zero(eltype(grid))
@inline div_𝐯v(i, j, k, grid, ::Nothing, U, v) = zero(eltype(grid))
@inline div_𝐯w(i, j, k, grid, ::Nothing, U, w) = zero(eltype(grid))

@inline div_𝐯u(i, j, k, grid, ::Nothing, ::ZeroU, u) = zero(eltype(grid))
@inline div_𝐯v(i, j, k, grid, ::Nothing, ::ZeroU, v) = zero(eltype(grid))
@inline div_𝐯w(i, j, k, grid, ::Nothing, ::ZeroU, w) = zero(eltype(grid))

@inline div_𝐯u(i, j, k, grid, ::Nothing, U, ::ZeroField) = zero(eltype(grid))
@inline div_𝐯v(i, j, k, grid, ::Nothing, U, ::ZeroField) = zero(eltype(grid))
@inline div_𝐯w(i, j, k, grid, ::Nothing, U, ::ZeroField) = zero(eltype(grid))

"""
    div_𝐯u(i, j, k, grid, advection, U, u)

Calculate the advection of momentum in the ``x``-direction using the conservative form, ``𝛁⋅(𝐯 u)``,

    1/Vᵘ * [δxᶠᵃᵃ(ℑxᶜᵃᵃ(Ax * u) * ℑxᶜᵃᵃ(u)) + δy_fca(ℑxᶠᵃᵃ(Ay * v) * ℑyᵃᶠᵃ(u)) + δz_fac(ℑxᶠᵃᵃ(Az * w) * ℑzᵃᵃᶠ(u))]

which ends up at the location `fcc`.
"""
@inline function div_𝐯u(i, j, k, grid, advection, U, u)
    return 1/Vᶠᶜᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, _advective_momentum_flux_Uu, advection, U[1], u) +
                                    δyᵃᶜᵃ(i, j, k, grid, _advective_momentum_flux_Vu, advection, U[2], u) +
                                    δzᵃᵃᶜ(i, j, k, grid, _advective_momentum_flux_Wu, advection, U[3], u))
end

"""
    div_𝐯v(i, j, k, grid, advection, U, v)

Calculate the advection of momentum in the ``y``-direction using the conservative form, ``𝛁⋅(𝐯 v)``,

    1/Vʸ * [δx_cfa(ℑyᵃᶠᵃ(Ax * u) * ℑxᶠᵃᵃ(v)) + δyᵃᶠᵃ(ℑyᵃᶜᵃ(Ay * v) * ℑyᵃᶜᵃ(v)) + δz_afc(ℑxᶠᵃᵃ(Az * w) * ℑzᵃᵃᶠ(w))]

which ends up at the location `cfc`.
"""
@inline function div_𝐯v(i, j, k, grid, advection, U, v)
    return 1/Vᶜᶠᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, _advective_momentum_flux_Uv, advection, U[1], v) +
                                    δyᵃᶠᵃ(i, j, k, grid, _advective_momentum_flux_Vv, advection, U[2], v)    +
                                    δzᵃᵃᶜ(i, j, k, grid, _advective_momentum_flux_Wv, advection, U[3], v))
end

"""
    div_𝐯w(i, j, k, grid, advection, U, w)

Calculate the advection of momentum in the ``z``-direction using the conservative form, ``𝛁⋅(𝐯 w)``,

    1/Vʷ * [δx_caf(ℑzᵃᵃᶠ(Ax * u) * ℑxᶠᵃᵃ(w)) + δy_acf(ℑzᵃᵃᶠ(Ay * v) * ℑyᵃᶠᵃ(w)) + δzᵃᵃᶠ(ℑzᵃᵃᶜ(Az * w) * ℑzᵃᵃᶜ(w))]

which ends up at the location `ccf`.
"""
@inline function div_𝐯w(i, j, k, grid, advection, U, w)
    return 1/Vᶜᶜᶠ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, _advective_momentum_flux_Uw, advection, U[1], w) +
                                    δyᵃᶜᵃ(i, j, k, grid, _advective_momentum_flux_Vw, advection, U[2], w) +
                                    δzᵃᵃᶠ(i, j, k, grid, _advective_momentum_flux_Ww, advection, U[3], w))
end


text:
#!/usr/bin/env julia

ProjDir = dirname(@__FILE__)

# Create rectilinear grid VTK file.

using WriteVTK
import Compat.UTF8String

FloatType = Float32
vtk_filename_noext = "rl"
outfiles = UTF8String[]

dim = 2
x = range(0.0, stop=1.0, length=30)
y = range(0.0, stop=1.0, length=30)

Ni, Nj, Nk = length(x), length(y), 1
z = ones(FloatType, Nk)

# Create some scalar and vectorial data.
p = zeros(FloatType, Ni, Nj, Nk)
q = zeros(FloatType, Ni, Nj, Nk)

for i = 1:Ni, j=1:Nj
  p[i, j, 1] = 10*i
  q[i, j, 1] = 10*i*i
end

cdata = zeros(FloatType, Ni-1, Nj-1)
for j = 1:Nj-1, i = 1:Ni-1
  cdata[i, j] = 2i + 20 * sin(3*pi * (j-1) / (Nj-2))
end

# Test extents (this is optional!!)
ext = [0, Ni-1, 0, Nj-1, 0, Nk-1] + 42

# Initialise new vtr file (rectilinear grid).
vtk = vtk_grid(ProjDir*"/"*vtk_filename_noext*"_$(dim)D", x, y; extent=ext)

# Add data.
vtk_point_data(vtk, p, "p_values")
vtk_point_data(vtk, q, "q_values")
vtk_cell_data(vtk, cdata, "myCellData")

# Save and close vtk file.
append!(outfiles, vtk_save(vtk))


text:
struct HeightDiameter <: Function
formula::Function
b
end

#individual parameters
function calculate_height(HD::HeightDiameter,dbh,spp)
    return(HD.formula(dbh,get(HD.b,spp,0)))
end

#separate arrays for dbh and spp
function calculate_height(HD::HeightDiameter,dbh::AbstractArray,spp::AbstractArray)
    par=map(x->get(HD.b,x,0),spp) #lookup species coeffs from dict
    map(HD.formula,dbh,par) #calculate
end

##Named Equations
#2 parameter equations, mainly from the LMFOR R package
#sorted alphabetically
Curtis=(x,b)->4.5+b[1]*((x/1+x)^b[2])
Michailoff=(x,b)->4.5+b[1]*exp(-b[2]*x^(-1))
Meyer=(x,b)->4.5+b[1]*(1-exp(-b[2]*x))
Micment=(x,b)->4.5+b[1]*x/(b[2]+x)
Micment2=(x,b)-> 4.5+x/(b[1]+b[2]*x)
Naslund=(x,b)->4.5+(x^2/(b[1]+b[2]*x)^2)
Naslund2=(x,b)-> 4.5+x^2/(b[1]+exp(b[2])*x)^2
Naslund3=(x,b)-> 4.5+x^2/(exp(b[1]) + b[2]*x)^2
Naslund4=(x,b)-> 4.5+x^2/(exp(b[1]) +exp(b[2])*x)^2
Power=(x,b)-> 4.5+b[1]*x^b[2]
Wyckoff=(x,b)->4.5+exp(b[1]+(b[2]/(x+1)))

#3 parameter equations, mainly from LMFOR R package
#sorted alphabetically
Chapman=(x,b)->4.5+b[1]*(1-exp(-b[2]*x))^b[3]
Gompertz=(x,b)->4.5+b[1]*exp(-b[2]*exp(-b[3]*x))
HossfeldIV=(x,b)->4.5+b[1]/(1+1/(b[2]*x^b[3]))
Korf=(x,b)->4.5+b[1]*exp(-b[2]*x^-b[3])
Logistic=(x,b)->4.5+b[1]/(1+b[2]*exp(-b[3]*x))
Monserud=(x,b)->4.5+exp(b[1]+(b[2]*x*b[3]))
Prodan=(x,b)->4.5+x^2/(b[1]+b[2]*x+b[3]*x^2)
Ratkowsky=(x,b)->4.5+b[1]*exp(-b[2]/(x+b[3]))
Sibbesen=(x,b)->4.5+b[1]*x^(b[2]*x^(-b[3]))
Weibull=(x,b)->4.5+b[1]*(1-exp(-b[2]*x^b[3]))
###


text:
using Test
using NQCDynamics
using Unitful
using UnitfulAtomic
using RecursiveArrayTools
using LinearAlgebra: diag
using StatsBase
using StochasticDiffEq
using ComponentArrays
using NQCDynamics: DynamicsMethods, DynamicsUtils
using NQCDynamics.DynamicsMethods: ClassicalMethods, IntegrationAlgorithms

atoms = Atoms([:H, :C])
model = CompositeFrictionModel(Free(3), ConstantFriction(3, 1))
sim = RingPolymerSimulation{MDEF}(atoms, model, 3; temperature=10u"K")

v = RingPolymerArray(randn(size(sim)))
r = RingPolymerArray(randn(size(sim)))

@testset "friction!" begin
    gtmp = zeros(length(r), length(r))
    gtmp = zeros(ndofs(sim)*natoms(sim),ndofs(sim)*natoms(sim),nbeads(sim))
    F = ClassicalMethods.friction!(gtmp, r, sim, 0.0)
    hmass = sim.calculator.model.friction_model.γ/atoms.masses[1]
    cmass = sim.calculator.model.friction_model.γ/atoms.masses[2]
    for i=1:length(sim.beads)
        @test diag(F[:,:,i]) ≈ [hmass, hmass, hmass, cmass, cmass, cmass]
    end
end

prob = DynamicsMethods.create_problem(ComponentVector(v=v,r=r), (0.0, 0.5), sim)

@testset "step_C!" begin
    dt = 0.5
    c = RingPolymers.cayley_propagator(sim.beads, dt; half=true)
    ω_k = RingPolymers.get_matsubara_frequencies(length(sim.beads), sim.beads.ω_n)
    vbefore = copy(v)
    rbefore = copy(r)
    IntegrationAlgorithms.step_C!(v,r,c)
    for I in CartesianIndices(r)
        A = [0 1; -ω_k[I[3]] 0]
        a = exp(A*dt/2) * [rbefore[I], vbefore[I]]
        @test a ≈ [r[I], v[I]] rtol=1e-4
    end
end

@testset "step_O!" begin
    cache = StochasticDiffEq.alg_cache(IntegrationAlgorithms.BCOCB(), prob, ArrayPartition(v,r),0,0,sim,0,0,0,Float64,0,0,0,0,0,0,Val{true})
    integrator = init(prob,IntegrationAlgorithms.BCOCB();dt=0.5)
    IntegrationAlgorithms.step_O!(cache.friction, integrator, v, r, 0.0)
end

@testset "ThermalLangevin" begin
    atoms = Atoms([:H, :C])
    sim = RingPolymerSimulation{ThermalLangevin}(atoms, NQCModels.Free(), 3; temperature=100u"K", γ=0.01)

    v = RingPolymerArray(zeros(size(sim)))
    r = RingPolymerArray(zeros(size(sim)))

    sol = run_trajectory(ArrayPartition(v,r), (0.0, 1e4), sim; dt=1, output=(:kinetic))
    # @test mean(sol.kinetic) ≈ austrip(100u"K") * length(sim.beads) rtol=5e-1
end

@testset "MDEF" begin
    atoms = Atoms([:H, :C])
    model = CompositeFrictionModel(Harmonic(), RandomFriction(1))
    sim = RingPolymerSimulation{MDEF}(atoms, model, 3; temperature=100u"K")

    v = RingPolymerArray(zeros(size(sim)))
    r = RingPolymerArray(zeros(size(sim)))

    sol = run_trajectory(ArrayPartition(v,r), (0.0, 1e2), sim; dt=1, output=(:kinetic))
end


text: const ASSET_FINGERPRINT = "d2be8332b8999e6db5fd"

max_stars_repo_name: ky-s/MyGenieApp

max_issues_repo_name: ky-s/MyGenieApp

max_forks_repo_name: ky-s/MyGenieApp


text:
__precompile__()

module Dexter
    import MathProgBase

    include("solver.jl")

    include("dd/ordering.jl")
    include("dd/core.jl")
    include("dd/construction.jl")
    include("dd/computation.jl")

    include("dd/bb/bb_tree.jl")
    include("dd/mathprog/decomposition.jl")
    include("dd/mathprog/function_analysis.jl")
    include("ip/oa/ecp.jl")
end

max_stars_repo_name: ctjandra/ddopt-ecp

max_issues_repo_name: ctjandra/ddopt-ecp

max_forks_repo_name: ctjandra/ddopt-ecp


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ d16d2b70-fe21-11ea-2b51-bf319474c8a9
using Plots

# ╔═╡ 256a2d40-fe22-11ea-1a5d-bb46637d959b
md"""
Solution to the Circle Area Derivation Projects:

1. To be revealed.
"""

# ╔═╡ 6b6bdc10-fe24-11ea-2da7-b708c0948f86
@bind n html"<input type=range>"

# ╔═╡ d1bae400-fe26-11ea-079c-c3eca8a09c6d
"n=$(n-1)"

# ╔═╡ 1a147440-fe23-11ea-0c95-9bb75706304c
begin
	tmin = 0
	tmax = 4π
	tvec = range(tmin, tmax, length = 100)

	plot(sin.(tvec), cos.(tvec),legend=false)
	#n=10
	svec = range(tmin,tmax/2,length=n)

	plot!(sin.(svec),cos.(svec))
	trig = [(0,1),(0,0),(sin(svec[2]),cos(svec[2]))]

	plot!(trig)
end

# ╔═╡ 2d17afb0-fe2a-11ea-10bf-b94e000267e8

# ╔═╡ be9ee910-fe22-11ea-1448-0d5ec4b08821
#md"""

#$$h=r\sin(\frac{\theta}{2})$$
#$$b=r\cos(\frac{\theta}{2})$$
#"""

# ╔═╡ Cell order:
# ╟─d16d2b70-fe21-11ea-2b51-bf319474c8a9
# ╟─256a2d40-fe22-11ea-1a5d-bb46637d959b
# ╟─d1bae400-fe26-11ea-079c-c3eca8a09c6d
# ╟─6b6bdc10-fe24-11ea-2da7-b708c0948f86
# ╠═1a147440-fe23-11ea-0c95-9bb75706304c
# ╠═2d17afb0-fe2a-11ea-10bf-b94e000267e8
# ╟─be9ee910-fe22-11ea-1448-0d5ec4b08821


text:
#############
# VtreeNode
#############

"Root of the vtree node hiearchy"
abstract type VtreeNode <: TreeNode end

const Vtree{VN<:VtreeNode} = AbstractVector{<:VN}

#############
# Constructors
#############

"Construct a balanced vtree"
function balanced_vtree(::Type{VN}, num_vars::Int)::Vtree{VN} where {VN <: VtreeNode}
    return node2dag(balanced_vtree_root(VN, 1, num_vars))
end

"Construct a balanced vtree root node"
function balanced_vtree_root(::Type{VN}, first::Int, last::Int)::VN where VN
    @assert last >= first "Must have $last >= $first"
    if last == first
        return VN(Var(first))
    else
        return VN(balanced_vtree_root(VN, first, first+(last-first+1)÷2-1),
                  balanced_vtree_root(VN, first+(last-first+1)÷2, last))
    end
end

function random_vtree(::Type{VN}, num_variables::Integer; vtree_mode::String="balanced")::Vtree{VN} where {VN <: VtreeNode}
    @assert vtree_mode in ["linear", "balanced", "rand"]
    vars = Var.(Random.randperm(num_variables))

    leaves = map(vars) do v
        VN(v)
    end

    vtree = [Vector{VN}(); leaves]
    right = popfirst!(vtree)
    while !isempty(vtree)
        left = popfirst!(vtree)
        v = VN(left, right)
        if vtree_mode == "linear"
            pushfirst!(vtree, v)
        elseif vtree_mode == "balanced"
            push!(vtree, v)
        elseif vtree_mode == "rand"
            pushrand!(vtree, v)
        end
        right = popfirst!(vtree)
    end
    node2dag(right)
end

"""
Construct PlainVtree top town, using method specified by split_method.
"""
function top_down_vtree(::Type{VN}, vars::Vector{Var}, split_method::Function)::Vtree{VN} where {VN <: VtreeNode}
    node2dag(top_down_root(VN, vars, split_method))
end

function top_down_root(::Type{VN}, vars::Vector{Var}, split_method::Function)::VN  where {VN <: VtreeNode}
    @assert !isempty(vars) "Cannot construct a vtree with zero variables"
    if length(vars) == 1
        VN(vars[1])
    else
        (X, Y) = split_method(vars)
        prime = top_down_root(VN, X, split_method)
        sub = top_down_root(VN, Y, split_method)
        VN(prime, sub)
    end
end

"""
Construct PlainVtree bottom up, using method specified by combine_method!.
"""
function bottom_up_vtree(::Type{VN}, vars::Vector{Var}, combine_method!::Function)::Vtree{VN} where {VN <: VtreeNode}
    vars = copy(vars)
    ln = Vector{VN}()
    node_cache = Dict{Var, VN}() # map from variable to *highest* level node

    "1. construct leaf node"
    for var in vars
        n = VN(var)
        node_cache[var] = n
        push!(ln, n)
    end

    "2. construct inner node"
    while length(vars) > 1
        matches = combine_method!(vars) # vars are mutable
        for (left, right) in matches
            n = VN(node_cache[left], node_cache[right])
            node_cache[left] = node_cache[right] = n
            push!(ln, n)
        end
    end

    "3. clean up"
    node2dag(ln[end])
end

#############
# Methods
#############

variables(v::Vtree) = variables(v[end])
num_variables(v::Vtree) = num_variables(v[end])

num_variables(n::VtreeNode) = num_variables(NodeType(n), n::VtreeNode)
num_variables(::Leaf, ::VtreeNode) = 1
num_variables(::Inner, n::VtreeNode) = length(variables(n))

"""
Check vtree validation, variables(parent) = variables(left) + variables(right)
"""
isvalid(n::VtreeNode)::Bool = isvalid(NodeType(n), n)

function isvalid(::Inner, n::VtreeNode)::Bool
    return num_variables(n) == num_variables(n.left) + num_variables(n.right) &&
        isequal(variables(n), union(variables(n.left), variables(n.right))) &&
        isvalid(n.right) &&
        isvalid(n.left)
end

function isvalid(::Leaf, n::VtreeNode)::Bool true; end

function isvalid(vtree::Vtree)::Bool
    return isvalid(vtree[end])
end

"""
Return the length from vtree inner node `n` to leaf node which contains `var`
"""
path_length(n::VtreeNode, var::Var)::Int = path_length(NodeType(n), n, var)

function path_length(::Inner, n::VtreeNode, var::Var)::Int
    @assert var in variables(n)
    if var in variables(n.left)
        return 1 + path_length(n.left, var)
    else
        return 1 + path_length(n.right, var)
    end
end

function path_length(::Leaf, n::VtreeNode, var::Var)::Int
    @assert variables(n) == [var]
    return 0
end


text:
module TestLinkTree

using Phylo
using DataFrames
using JuliaDB

using Test
using IterableTables: getiterator

species = ["Dog", "Cat", "Human"]
ntips = 10
df = DataFrame(species = species, count = [10, 20, 3])
observations = ["Dog", "Cat", "Dog", "Dog"]
jdb = table((species = observations, count = 1:4))

@testset "LinkTree()" begin
    lts = RootedTree(species)
    @test nodedatatype(typeof(lts)) ≡ Dict{String, Any}
    @test branchdatatype(typeof(lts)) ≡ Dict{String, Any}
    @test leafinfotype(typeof(lts)) ≡ Dict{String, Any}

    ltdf = Phylo.LT{OneRoot, DataFrame, Float64}(df)
    @test nodedatatype(typeof(ltdf)) ≡ Dict{String, Any}
    @test branchdatatype(typeof(ltdf)) ≡ Dict{String, Any}
    @test leafinfotype(typeof(ltdf)) ≡ DataFrame

    rt = Unrooted
    lb = LinkBranch{rt, String, Nothing}
    ln = LinkNode{rt, String, Vector{Int}, lb}
    ltjdb = LinkTree{rt, String, ln, lb, typeof(jdb)}(jdb)
    @test nodedatatype(typeof(ltjdb)) ≡ Vector{Int}
    @test branchdatatype(typeof(ltjdb)) ≡ Nothing
    @test leafinfotype(typeof(ltjdb)) ≡ typeof(jdb)
end

end