text:
include("Lanczos_Ritz_Utilities.jl")
include("../../src/iterative_lanczos.jl")

# ------------------------------------------------------------------------------
# Demonstrating the loss of orthogonality of the lanczos vectors.
# And the linear dependence of eigenvalues of the exterior of the spectrum.
# ------------------------------------------------------------------------------
n = 64
A = Diagonal(LinRange(-1, 1, n).^3)
il = IterativeLanczos(A, rand(n))
for _ in 1: n - 1
    il()
end
Q = GetQMatrix(il)
T = GetTMatrix(il)

fig = heatmap(Q'*Q, size=(722, 512))
fig2 = heatmap(Q'*A*Q, size=(722, 512))
fig3 = heatmap(
    (A*Q[:, 1:end - 1] - Q*T[:, 1:end - 1]).|> abs,
    size=(1200, 768), dpi=250
)
savefig(fig, "$(@__DIR__)/plots/fig3.png")
savefig(fig2, "$(@__DIR__)/plots/fig4.png")

# ------------------------------------------------------------------------------
# All the ritz values during the computations process. and plotting it.
# ------------------------------------------------------------------------------

function RiztTrajectoryPlot(filename, n=64; itr_offset=20, full_ortho=false, title="")
    A = Diagonal(LinRange(-1, 1, n).^3)
    il = IterativeLanczos(A, ones(n))
    il.reorthogonalize = full_ortho
    FoundRitzValues = Vector{Vector{Float64}}()
    # push!(FoundRitzValues, [GetTMatrix(il)])
    TrueEigenValues = diag(A)

    for II in 1: n - 1
        il()
        T = GetTMatrix(il)
        λs, _ = eigen(T)
        if II >= itr_offset
            push!(FoundRitzValues, λs)
        end
    end

    for RitzValue in FoundRitzValues
        sort!(RitzValue, rev=true)
    end

    fig4 = scatter(title=title, legend=false)
    for Idx in 1: div(itr_offset,2)
        RitzTrajectory = Vector{Float64}()
        for RitzValues in FoundRitzValues
            if Idx <= length(RitzValues)
                push!(RitzTrajectory, RitzValues[Idx])
            end
        end
        plot!(
            fig4,
            # Idx: length(RitzTrajectory) + Idx - 1,
            itr_offset:(length(RitzTrajectory) + itr_offset - 1),
            RitzTrajectory, size=(750,750),
            dpi=250,
            markershape=:cross,
            linestyle=:solid
        )
    end
    for Idx in 1:div(itr_offset, 2)
        RitzTrajectory = Vector{Float64}()
        for RitzValues in FoundRitzValues
            if Idx <= length(RitzValues)
                push!(RitzTrajectory, RitzValues[end - Idx + 1])
            end
        end
        plot!(
            fig4,
            # Idx: length(RitzTrajectory) + Idx - 1,
            itr_offset:(length(RitzTrajectory) + itr_offset - 1),
            RitzTrajectory,
            dpi=250,
            markershape=:xcross,
            linestyle=:solid
        )
    end
    xlabel!(fig4, "iterations")
    ylabel!(fig4, "ritz values")
    display(fig4)
    savefig(fig4, "$(@__DIR__)/plots/$(filename).png")

return end

RiztTrajectoryPlot("ritz_trajectory_plot_floats", title="ritz trajectory no-reorthgonalization")
RiztTrajectoryPlot("ritz_trajectory_plot_exact", full_ortho=true, title="ritz trajectory with full-reorthogonalizatoin")


text:
# Buffer
# ======

struct _Buffer
    content::Cstring
    use::Cuint
    size::Cuint
    alloc::Ptr{Cvoid}
    contentIO::Cstring
end

mutable struct Buffer
    ptr::Ptr{_Buffer}

    function Buffer()
        buf_ptr = ccall(
            (:xmlBufferCreate, libxml2),
            Ptr{_Buffer},
            ())
        buf = new(buf_ptr)
        finalizer(finalize_buffer, buf)
        return buf
    end
end

function finalize_buffer(buf)
    ccall(
        (:xmlBufferFree, libxml2),
        Cvoid,
        (Ptr{Cvoid},),
        buf.ptr)
end

ext: jl


text:
# Setup
testpath = pwd()
scratch = tempdir()
cd(scratch)

## --- Times table, file IO, mallocarray

# Attempt to compile
# We have to start a new Julia process to get around the fact that Pkg.test
# disables `@inbounds`, but ironically we can use `--compile=min` to make that
# faster.
status = run(`julia --compile=min $testpath/scripts/times_table.jl`)
@test isa(status, Base.Process)
@test status.exitcode == 0

# Attempt to run
println("5x5 times table:")
status = run(`./times_table 5 5`)
@test isa(status, Base.Process)
@test status.exitcode == 0
@test parsedlm(Int64, c"table.tsv", '\t') == (1:5)*(1:5)'

## --- Random number generation

# Compile...
# We have to start a new Julia process to get around the fact that Pkg.test
# disables `@inbounds`, but ironically we can use `--compile=min` to make that
# faster.
status = run(`julia --compile=min $testpath/scripts/rand_matrix.jl`)
@test isa(status, Base.Process)
@test status.exitcode == 0

# Run...
println("5x5 random matrix:")
status = run(`./rand_matrix 5 5`)
@test isa(status, Base.Process)
@test status.exitcode == 0

## --- Test LoopVectorization integration

if Bool(LoopVectorization.VectorizationBase.has_feature(Val{:x86_64_avx2}))
    # Compile...
    status = run(`julia --compile=min $testpath/scripts/loopvec_product.jl`)
    @test isa(status, Base.Process)
    @test status.exitcode == 0

    # Run...
    println("10x10 table sum:")
    status = run(`./loopvec_product 10 10`)
    @test isa(status, Base.Process)
    @test status.exitcode == 0
    @test parsedlm(c"product.tsv",'\t')[] == 3025
end

# Compile...
status = run(`julia --compile=min $testpath/scripts/loopvec_matrix.jl`)
@test isa(status, Base.Process)
@test status.exitcode == 0

# Run...
println("10x5 matrix product:")
status = run(`./loopvec_matrix 10 5`)
@test isa(status, Base.Process)
@test status.exitcode == 0
A = (1:10) * (1:5)'
@test parsedlm(c"table.tsv",'\t') == A' * A

## --- Test string handling

# Compile...
status = run(`julia --compile=min $testpath/scripts/print_args.jl`)
@test isa(status, Base.Process)
@test status.exitcode == 0

# Run...
println("String indexing and handling:")
status = run(`./print_args foo bar`)
@test isa(status, Base.Process)
@test status.exitcode == 0

## --- Clean up

cd(testpath)

## ---


text:
using Plots
using IsingModels

function mag_callback!(model,m)
  m[] += calc_avg_magnetization(model)
end

h = 0.0
J = 1.0
Tc = 2.269185*J
Tog = 0.1*Tc
n = 100

model = IsingModel2D(h,J,Tog,n,seed=convert(UInt64,42))
randomdata!(model)
anim = @animate for i in 1:100
  mcsweep(model)
  heatmap!(model.data, c = :cool)
end
gif(anim, "low_T_ising.gif")

allup!(model)
mcsweeps(1000,model)

Tfacs = range(0.1,1.;step=0.1)
ms = zeros(length(Tfacs))
for i in 1:length(Tfacs)
  T = Tfacs[i]*Tc
  model.T = T
  mcsweeps(100,model)
  m = Ref{Float64}(ms[i])
  mcsweeps(1000,model;cb=mag_callback!,cbargs=m)
  ms[i] = abs(m[] / 1000.0)
  println(ms[i])
end

plt = plot(Tfacs,ms)
savefig("avg_mag.png")

anim = @animate for i in 1:100
  mcsweep(model)
  heatmap!(model.data, c = :cool)
end
gif(anim, "critical_T_ising.gif")


text:
module FailureRepresentation

include("SystemUnderTest.jl")
include("Adversary.jl")
include("Autoencoder.jl")
include("utils.jl")

using .Adversary
using .Autoencoder
using .SystemUnderTest
using BSON: @save
using CUDA
using Flux
using Flux: onecold
using Parameters
using Random
using Statistics

export
    SystemUnderTest,
    load,
    predict,
    getdata,
    Adversary,
    Autoencoder,
    load,
    autoencoder,
    sampled_validation_iteration

@with_kw mutable struct Metrics
    precision = [] # of examples recognized as failures, what percent ARE failures?
    recall = [] # of all examples that are failures, what percent are correctly classified?
end

"""
Load models for black-box system under test `sut` and dataset autoencoder `autoencoder`.
"""
function loadmodels()
    Core.eval(Main, :(import NNlib)) # required for .load
    Core.eval(Main, :(import Flux)) # required for .load
    sut = SystemUnderTest.load("models/sut.bson")
    autoencoder = Autoencoder.load("models/autoencoder.bson")
    encoder = Autoencoder.load("models/encoder.bson", :encoder)
    return sut, autoencoder, encoder
end

"""
Calculate the failure rate of the samples given the `model`, the inputs `X`, and targets `Y`.
"""
failure_rate(model, X, Y) = round(1 - mean(onecold(model(X)) .== onecold(Y)), digits=3)

"""
Main sampled validation iteration loop. See diagram for details.
"""
function sampled_validation_iteration(; T=1, iters=1, m=500, seedoffset=1)
    sut, autoencoder, encoder = loadmodels()
    _, testdata = SystemUnderTest.getdata()

    metrics = Metrics()
    metrics_rand = Metrics()
    local adversary = nothing # initialize to "no model"
    local mapping
    local X
    local Y

    # sampled validation iteration: loop t
    for t in 1:T
        @info "Model design interation, loop $t"
        Random.seed!(t*seedoffset)
        X, Y = rand(testdata, m) # `m` samples of testdata
        @info string("Sampled failure rate: ", failure_rate(sut, X, Y))

        # model evaluation iteration: loop i (UNUSED in non-continual version)
        for i in 1:iters
            ## Adversarial step:
            # train adversary using samples of failures from low dimensional representation (linked to true samples)
            # adversary is a classifier that determines if input is a failure or not
            adversary, encoded_testdata, true_testdata = Adversary.train(adversary, X, Y, sut, encoder)
            mapping = getmapping(true_testdata) # maps i -> true (x,y)

            # Save out adversary for analysis
            @save "models/adversary_$t.bson" adversary

            ## SUT evaluation
            # pass samples into SUT to evaluate if they're true failures
            last_num_cands = 0
            num_actual_failures = 0
            for selection_mode in [:debug, :adversary, :random]
                if selection_mode == :debug
                    X_candidates, Y_candidates = Adversary.select_candidates(adversary, encoded_testdata, mapping; debug=true)
                    num_actual_failures = length(X_candidates)
                elseif selection_mode == :adversary
                    X_candidates, Y_candidates = Adversary.select_candidates(adversary, encoded_testdata, mapping)
                elseif selection_mode == :random
                    X_candidates, Y_candidates = Adversary.rand_select_candidates(encoded_testdata, mapping, last_num_cands)
                else
                    error("No selection mode: $selection_mode")
                end

                @info "Number of candidate failures: $(length(X_candidates))"

                ## Run candidate failures on SUT
                failures = 0
                for (x,y) in zip(X_candidates, Y_candidates)
                    sut_prediction::Int = SystemUnderTest.predict(sut, x)[1]
                    true_target::Int = onecold(y)[1]
                    failure::Bool = sut_prediction != true_target
                    @info sut_prediction, true_target, failure, adversary(encoder(x))[1]
                    failures += failure
                end
                precision = failures/length(X_candidates)
                precision = isnan(precision) ? 0 : precision
                recall = failures/num_actual_failures
                recall = isnan(recall) ? 0 : recall
                if selection_mode == :adversary
                    push!(metrics.precision, precision)
                    push!(metrics.recall, recall)
                    last_num_cands = length(X_candidates)
                elseif selection_mode == :random
                    # Random sampling of encoded_testdata, using the same number of selected candidates as the adversary
                    push!(metrics_rand.precision, precision)
                    push!(metrics_rand.recall, recall)
                end
                @info "Precision: $precision"
                @info "Recall: $recall"
                println("—"^40)
            end
        end
    end

    println()
    println()
    @info "Random selection precision: $(round(mean(metrics_rand.precision), digits=4)) ±$(round(std(metrics_rand.precision), digits=2))"
    @info "Adversary selection precision: $(round(mean(metrics.precision), digits=4)) ±$(round(std(metrics.precision), digits=2))"
    println("—"^40)
    @info "Random selection recall: $(round(mean(metrics_rand.recall), digits=4)) ±$(round(std(metrics_rand.recall), digits=2))"
    @info "Adversary selection recall: $(round(mean(metrics.recall), digits=4)) ±$(round(std(metrics.recall), digits=2))"

    # use adversary to generate failures.
    # TEST: sample from adversary to get failures
    return X, Y, mapping, adversary, metrics, metrics_rand
end

end # module


text:
## Exercise 6-8
## The greatest common divisor (GCD) of a and b is the largest number that divides both of them with no remainder.

## One way to find the GCD of two numbers is based on the observation that if r is the remainder when a is divided by b, then gcd(a, b) = gcd(b, r). As a base case, we can use gcd(a, 0) = a.

## Write a function called gcd that takes parameters a and b and returns their greatest common divisor.
println("Ans: ")

function gcd(a, b)
    remainder = a % b
    if remainder == 0
        return b
    end
    return gcd(b, remainder)
end

println(gcd(210, 55))
println(gcd(12576, 4052))
println(gcd(4, 3))

println("End.")


text:
using NBodyIPs
using JuLIP, Test, StaticArrays, ForwardDiff, Combinatorics
using BenchmarkTools

using JuLIP.Potentials: evaluate, evaluate_d

using LinearAlgebra: norm
using Printf

const blinvariants = NBodyIPs.BLInvariants.invariants
const blinvariants_d = NBodyIPs.BLInvariants.invariants_d
const blinvariants_ed = NBodyIPs.BLInvariants.invariants_ed

include("aux_testing.jl")

all_blinvariants(r) = vcat(blinvariants(r)...)  # [I1; I2]
ad_blinvariants(r) = ForwardDiff.jacobian(all_blinvariants, r)

println("-------------------------------------------")
println("   Testing implementation of `blinvariants`")
println("-------------------------------------------")

# TODO: test correctness of the blinvariants implementation
#       against the MAGMA output

println("[1] Correctness of gradients")
for dim in [3, 6 ] #, 10]
   r = 1.0 + SVector(rand(dim)...)
   println("---------------")
   println("dim = $dim")
   println("---------------")

   I = all_blinvariants(r)
   dI1, dI2 = blinvariants_d(r)
   dI = [hcat(dI1...)'; hcat(dI2...)']
   dIh = zeros(size(dI))
   r0 = Vector(r)
   errs = []
   for p = 2:9
      h = .1^p
      dIh = zeros(size(dI))
      for j = 1:length(r)
         r0[j] += h
         Ih = all_blinvariants(SVector(r0...))
         dIh[:, j] = (Ih - I) / h
         r0[j] -= h
      end
      push!(errs, norm(dIh - dI, Inf))
      @printf(" %d | %.2e \n", p, errs[end])
   end
   println("---------------")
   @test minimum(errs) <= 1e-3 * maximum(errs)
end
println()

println("[2] Symmetry")
for dim in [3, 6 ]# , 10]
   for n = 1:3
      r = 1.0 + SVector(rand(dim)...)
      I = all_blinvariants(r)
      for rπ in simplex_permutations(r)
         @test I ≈ all_blinvariants(SVector(rπ...))
      end
      print(".")
   end
end
println()

println("[3] blinvariants_ed")
for dim in [3, 6 ]# , 10]
   for n = 1:3
      r = 1.0 + SVector(rand(dim)...)
      I1, I2 = blinvariants(r)
      dI1, dI2 = blinvariants_d(r)
      J1, J2, dJ1, dJ2 = blinvariants_ed(r)
      @test all(i ≈ j for (i,j) in zip(I1, J1))
      @test all(i ≈ j for (i,j) in zip(I2, J2))
      @test all(di ≈ dj for (di,dj) in zip(dI1, dJ1))
      @test all(di ≈ dj for (di,dj) in zip(dI2, dJ2))
      print(".")
   end
end
println()

println("-------------------------------------------")
println("   Testing implementation of `bainvariants`")
println("-------------------------------------------")

const BAI = NBodyIPs.BAInvariants
all_bainvariants(x) = vcat(BAI.invariants(x)...)
ad_bainvariants(x) = ForwardDiff.jacobian(all_bainvariants, x)
const bainvariants = BAI.invariants
const bainvariants_d = BAI.invariants_d
const bainvariants_ed = BAI.invariants_ed

println("[1] Correctness of gradients")
for dim in [3, 6 ] #, 10]
   # r ≡ x
   r = 1.0 + SVector(rand(dim)...)
   println("---------------")
   println("dim = $dim")
   println("---------------")

   I = all_bainvariants(r)
   dI1, dI2 = bainvariants_d(r)
   dI = [hcat(dI1...)'; hcat(dI2...)']
   dIh = zeros(size(dI))
   r0 = Vector(r)
   errs = []
   for p = 2:11
      h = .1^p
      dIh = zeros(size(dI))
      for j = 1:length(r)
         r0[j] += h
         Ih = all_bainvariants(SVector(r0...))
         dIh[:, j] = (Ih - I) / h
         r0[j] -= h
      end
      push!(errs, norm(dIh - dI, Inf))
      @printf(" %d | %.2e \n", p, errs[end])
   end
   println("---------------")
   @test minimum(errs) <= 1e-3 * maximum(errs)
end
println()

println("[3] all_bainvariants_ed")
for dim in [3, 6 ] # , 10]
   for n = 1:3
      r = 1.0 + SVector(rand(dim)...)
      I1, I2 = bainvariants(r)
      dI1, dI2 = bainvariants_d(r)
      J1, J2, dJ1, dJ2 = bainvariants_ed(r)
      @test all(i ≈ j for (i,j) in zip(I1, J1))
      @test all(i ≈ j for (i,j) in zip(I2, J2))
      @test all(di ≈ dj for (di,dj) in zip(dI1, dJ1))
      @test all(di ≈ dj for (di,dj) in zip(dI2, dJ2))
      print(".")
   end
end
println()

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
# Autogenerated wrapper script for Catalyst_jll for x86_64-apple-darwin
export catalyst_replay_exe, libcatalyst

JLLWrappers.@generate_wrapper_header("Catalyst")
JLLWrappers.@declare_library_product(libcatalyst, "@rpath/libcatalyst.3.dylib")
JLLWrappers.@declare_executable_product(catalyst_replay_exe)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcatalyst,
        "lib/libcatalyst.3.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        catalyst_replay_exe,
        "bin/catalyst_replay",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# additional tests for config
fs2()

foofig(s) = (write(joinpath(td, "config.md"), s); F.process_config())

@testset "config" begin
    # ================================
    # asssignments go to GLOBAL
    foofig(raw"""
        @def var = 5
        """)
    @test haskey(F.GLOBAL_VARS, "var")
    @test F.GLOBAL_VARS["var"][1] == 5

    # lxdefs go to GLOBAL
    foofig(raw"""
        \newcommand{\hello}{goodbye}
        """)
    @test haskey(F.GLOBAL_LXDEFS, "\\hello")
    @test F.GLOBAL_LXDEFS["\\hello"].def == "goodbye"

    # combination of lxdefs
    foofig(raw"""
        \newcommand{\hello}{goodbye}
        \newcommand{\hellob}{\hello}
        \newcommand{\helloc}{\hellob}
        """)

    @test F.GLOBAL_LXDEFS["\\hello"].from < F.GLOBAL_LXDEFS["\\hello"].to <
          F.GLOBAL_LXDEFS["\\hellob"].from < F.GLOBAL_LXDEFS["\\hellob"].to <
          F.GLOBAL_LXDEFS["\\helloc"].from < F.GLOBAL_LXDEFS["\\helloc"].to

    @test fd2html(raw"""\helloc"""; dir=td, internal=true) == "goodbye"
end

@testset "i381" begin
    foofig(raw"""
        @def tags = []
        """)

    s = """
        @def tags = ["tag1", "tag2"]

        ~~~
        {{for tag in tags}}
        {{fill tag}}
        {{end}}
        ~~~
        """ |> fd2html_td
    @test isapproxstr(s, "tag1 tag2")
end


text:
module ImarisWriter

using Libdl
#using Lz4_jll, Zlib_jll, HDF5_jll

# Need the following DLLs
# 1. bpImarisWriter96.dll
# 2. lz4.dll
# 3. zlib.dll
# 4. hdf5.dll

# Write your package code here.

include("imariswriter_common.jl")
include("imariswriter_api.jl")

const bpImarisWriter96_path = joinpath( @__DIR__(), "..", "bpImarisWriter96.dll")
const bpImarisWriter96 = :bpImarisWriter96

function __init__()
    global bpImarisWriter96_h
    bpImarisWriter96_h =  dlopen( bpImarisWriter96_path )
end

end


text:
using BasicFWI
using jInv.Mesh
using jInv.Utils
using jInv.ForwardShare

if nworkers()==1
	addprocs(2)
end

@everywhere begin
	using BasicFWI
	using jInv.Mesh
	using jInv.Utils
	using jInv.ForwardShare
end

# setup model and attenuation function
nx = 96; nz = 48
domain = [-1.,1,0,1]
Mr = getRegularMesh(domain,[nx,nz])

# the slowness model
m = ones(nx,nz);
m[20:70,10:32] .= 3;

# attenuation for BC
padx = 6; padz = 6
a    = 20;
xc = getCellCenteredGrid(Mr)
gammaxL = a*(xc[:,1] .- xc[padx,1]).^2;
gammaxL[padx+1:end,:] .= 0
gammaxR = a*(xc[:,1] .- xc[end-padx+1,1]).^2
gammaxR[1:end-padx,:] .= 0

gammax = gammaxL + gammaxR

gammaz = a*(xc[:,2] .- xc[end-padz+1,2]).^2
gammaz[:,1:end-padz] .= 0

gamma = gammax + gammaz

# parameters for the Helmholtz (units in km)
h = [50.0; 50.0;]/1000
n = [nx; nz;]
omega = 2*pi*[0.2;0.5;1;1.2;1.4;1.5;2.0;3.0]
nfreq = length(omega)

# generate sources
q = zeros(tuple(n.+1...)); q[padx+1:2:end-padx-1,1] .= 1e4
Q = sdiag(vec(q))
Q = Q[:,(LinearIndices(sum(Q,dims=2) .!= 0))[findall(sum(Q,dims=2) .!= 0)]]
nsrc = size(Q,2)

# receivers
p = zeros(tuple(n.+1...)); p[padx+1:end-padx-1,1] .= 1
P = sdiag(vec(p))
P = P[:,(LinearIndices(sum(P,dims=2) .!= 0))[findall(sum(P,dims=2) .!= 0)]]
nrec = size(P,2)

pFor  = getFWIparam(omega,gamma,Q,P,Mr)
pForp = getFWIparam(omega,gamma,Q,P,Mr,true)

# inversion mesh and forward mesh are the same here
M2Mp = ones(length(pForp))


text:
"""
Generates a non-trivial random MOI linear program by adding variables
and constraints to MOI-compatible Optimizer `optimizer`

minimize `c' * x`
subject to `Ax <= b, x >= 0`
where `x in R^{n}, A in R^{m*n}, b in R^{m}, c in R^{n}`

Note: Mutates the `optimizer` object
"""
function generate_lp(optimizer, n, m)
    s = rand(m)
    s = 2*s .- 1
    λ = max.(-s, 0)
    s = max.(s, 0)
    x̂ = rand(n)
    A = rand(m, n)
    b = A * x̂ .+ s
    c = -A' * λ

    x = MOI.add_variables(optimizer, n)

    # define objective
    objective_function = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(c, x), 0.0)
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objective_function)
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    # set constraints
    for i in 1:m
        MOI.add_constraint(
            optimizer,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(A[i,:], x), 0.), MOI.LessThan(b[i]),
        )
    end

    for i in 1:n
        MOI.add_constraint(optimizer, MOI.SingleVariable(x[i]), MOI.GreaterThan(0.0))
    end
end

"""
Generates a non-trivial random MOI convex quadratic program
by adding variables and constraints to MOI compatible Optimizer `optimizer`

minimize `0.5 * x' * Q * x  + q' * x`
subject to `Gx <= h, Ax == b`
where `x in R^{n}, Q in R^{n*n}, q in R^{n}, G in R^{m,n}, h in R^{m}, A in R^{p*n}, b in R^{p}`

Note: (1) Mutates the `optimizer` object
      (2) Matrix `Q` is Positive Semidefinite
"""
function generate_qp(optimizer,n,m,p)
    x̂ = rand(n)
    Q = rand(n, n)
    Q = Q' * Q # ensure PSD
    q = rand(n)
    G = rand(m, n)
    h = G * x̂ + rand(m)
    A = rand(p, n)
    b = A * x̂

    x = MOI.add_variables(optimizer, n)

    # define objective
    quadratic_terms = MOI.ScalarQuadraticTerm{Float64}[]
    for i in 1:n
        for j in i:n # indexes (i,j), (j,i) will be mirrored. specify only one kind
            push!(quadratic_terms, MOI.ScalarQuadraticTerm(Q[i,j],x[i],x[j]))
        end
    end

    objective_function = MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm.(q, x), quadratic_terms, 0.0)
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), objective_function)
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    # set constraints
    for i in 1:m
        MOI.add_constraint(
            optimizer,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(G[i,:], x), 0.),MOI.LessThan(h[i])
        )
    end

    for i in 1:p
        MOI.add_constraint(
            optimizer,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(A[i,:], x), 0.),MOI.EqualTo(b[i])
        )
    end
end


text:
# TODO
# - add an interface via points & distance metrics
# - get rid of `nobs` in `cutree()`

module Fastcluster

    # # Load in `deps.jl`, complaining if it does not exist
    # const depsjl_path = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
    # if !isfile(depsjl_path)
    #     error("Fastcluster.jl not installed properly, run Pkg.build(\"Fastcluster\"), restart Julia and try again")
    # end
    # include(depsjl_path)

    ##############################################################################
    ##
    ## Dependencies
    ##
    ##############################################################################

    # no dependencies

    ##############################################################################
    ##
    ## Exported methods and types
    ##
    ##############################################################################

    export linkage, cutree

    # global const __pkgsrcdir = @__DIR__
    # global const __libfastcluster = joinpath(__pkgsrcdir, "../lib/libfastcluster.dylib")
    global const __pkgsrcdir = @__DIR__
    global const libfastcluster = joinpath(__pkgsrcdir, "../deps/usr/libfastcluster.so")

    ##############################################################################
    ##
    ## Load files
    ##
    ##############################################################################

    include("linkage.jl")

    ##############################################################################
    ##
    ## Initialize libfastcluster
    ##
    ##############################################################################

    function __init__()

        # Always check your dependencies from `deps.jl`
        #check_deps()

        # make sure that libfastcluster.so exists and can be called
        try
            d = @__DIR__
            t = ccall((:validate, libfastcluster),
                Int32,
                ()
                )
            if t != 42
                error("Failed to load libfastcluster.")
            end
        catch
            error("Failed to load libfastcluster.")
        end

    end

end


text:
struct D2Q13{
    T <: Real,
    Abscissaes <: AbstractMatrix{<:Integer},
    Weights <: AbstractVector{<:Real},
} <: Quadrature
    abscissae::Abscissaes
    weights::Weights
    speed_of_sound_squared::T
end
function D2Q13()
    # High Knudsen Number Thermal Flows with the D2Q13 Lattice Boltzmann Model
    # P. Lopez & Y. Bayazitoglu
    w_0 = 3 / 8
    w_1 = 1 / 12
    w_2 = 1 / 16
    w_3 = 1 / 96
    cs = 2.0

    return D2Q13(
        [   #A #B          #C          #D
            0 1 -1 0 0 1 -1 1 -1 2 -2 0 0
            0 0 0 1 -1 1 -1 -1 1 0 0 -2 2
        ],
        [w_0, w_1, w_1, w_1, w_1, w_2, w_2, w_2, w_2, w_3, w_3, w_3, w_3],
        cs,
    )
end

order(q::D2Q13) = 5
Base.show(io::IO, q::D2Q13) = show(io, "D2Q13")
Base.string(q::D2Q13) = "D2Q13"


text:
"""
    ArcSwift()

Arc-Swift transition system for dependency parsing.

Described in [Qi & Manning 2017](https://nlp.stanford.edu/pubs/qi2017arcswift.pdf).
"""
struct ArcSwift <: AbstractTransitionSystem end

initconfig(::ArcSwift, graph::DependencyTree) = ArcSwiftConfig(graph)
initconfig(::ArcSwift, words) = ArcSwiftConfig(words)

transition_space(::ArcSwift, labels=[]; max_k=5) =
    isempty(labels) ? [LeftArc.(1:max_k)..., RightArc.(1:max_k)..., Shift()] :
    [[LeftArc(k, l) for k in 1:max_k for l in labels]...,
     [RightArc(k, l) for k in 1:max_k for l in labels]...,
     Shift()]

projective_only(::ArcSwift) = true

struct ArcSwiftConfig <: AbstractParserConfiguration
    σ::Vector{Int}
    β::Vector{Int}
    A::Vector{Token}
end

function ArcSwiftConfig(words)
    σ = [0]
    β = collect(1:length(words))
    A = Token.(words)
    ArcSwiftConfig(σ, β, A)
end

function ArcSwiftConfig(gold::DependencyTree)
    σ = [0]
    β = collect(1:length(gold))
    A = [Token(t, head=-1) for t in gold]
    ArcSwiftConfig(σ, β, A)
end

stack(cfg::ArcSwiftConfig)  = cfg.σ
buffer(cfg::ArcSwiftConfig) = cfg.β

tokens(cfg::ArcSwiftConfig) = cfg.A
tokens(cfg::ArcSwiftConfig, is) = [token(cfg, i) for i in is if 0 <= i <= length(cfg.A)]

function leftarc(cfg::ArcSwiftConfig, k::Int, args...; kwargs...)
    # Assert a head-dependent relation between the word at the front
    # of the input buffer and the word at the top of the stack; pop
    # the stack.
    i = length(cfg.σ) - k + 1
    s, b = cfg.σ[i], cfg.β[1]
    A = copy(cfg.A)
    if s > 0
        A[s] = Token(A[s]; head=b, kwargs...)
    end
    ArcSwiftConfig(cfg.σ[1:i-1], cfg.β, A)
end

function rightarc(cfg::ArcSwiftConfig, k::Int, args...; kwargs...)
    # Assert a head-dependent relation between the word on the top of
    # the σ and the word at front of the input buffer; shift the
    # word at the front of the input buffer to the stack.
    i = length(cfg.σ) - k + 1
    s, b = cfg.σ[i], cfg.β[1]
    A = copy(cfg.A)
    A[b] = Token(A[b]; head=s, kwargs...)
    ArcSwiftConfig([cfg.σ[1:i] ; b], cfg.β[2:end], A)
end

function shift(cfg::ArcSwiftConfig)
    # Remove the word from the front of the input buffer and push it
    # onto the stack.
    ArcSwiftConfig([cfg.σ ; cfg.β[1]], cfg.β[2:end], cfg.A)
end

isfinal(cfg::ArcSwiftConfig) =
    all(t -> all(h >= 0 for h in t.head), cfg.A) #&& length(cfg.σ) > 0 && length(cfg.β) > 0

"""
    static_oracle(cfg::ArcSwiftConfig, tree, arc)

Oracle function for arc-swift dependency parsing.

Described in [Qi & Manning 2017](https://nlp.stanford.edu/pubs/qi2017arcswift.pdf).
"""
function static_oracle(cfg::ArcSwiftConfig, gold, arc=untyped)
    S = length(cfg.σ)
    if S >= 1 && length(cfg.β) >= 1
        b = cfg.β[1]
        for k in 1:S
            i = S - k + 1
            s = cfg.σ[i]
            has_arc(gold, b, s) && return LeftArc(k, arc(gold[s])...)
            has_arc(gold, s, b) && return RightArc(k, arc(gold[b])...)
        end
    end
    return Shift()
end

==(cfg1::ArcSwiftConfig, cfg2::ArcSwiftConfig) =
    cfg1.σ == cfg2.σ && cfg1.β == cfg2.β && cfg1.A == cfg2.A

# TODO
function possible_transitions(cfg::ArcSwiftConfig, graph::DependencyTree, arc=untyped)
    TransitionOperator[static_oracle(cfg, graph, arc)]
end

function Base.show(io::IO, c::ArcSwiftConfig)
    A = join(["$i $(t.form) $(t.head)" for (i,t) in enumerate(tokens(c))], ", ")
    print(io, "ArcSwiftConfig($(c.σ),$(c.β),$A)")
end


text:
let
    mdp = GridWorld()

    pomdp = FullyObservablePOMDP(mdp)

    @test observations(pomdp) == states(pomdp)
    @test n_observations(pomdp) == n_states(pomdp)
    @test statetype(pomdp) == obstype(pomdp)

    @test observations(pomdp) == states(pomdp)
    @test n_observations(pomdp) == n_states(pomdp)
    @test statetype(pomdp) == obstype(pomdp)

    s_po = initialstate(pomdp, MersenneTwister(1))
    s_mdp = initialstate(mdp, MersenneTwister(1))
    @test s_po == s_mdp

    solver = ValueIterationSolver(max_iterations = 100)
    mdp_policy = solve(solver, mdp)
    pomdp_policy = solve(solver, pomdp)
    @test mdp_policy.util == pomdp_policy.util
end


text:
using Compat
using BenchmarkTools
using LinearAlgebra
using SparseArrays
using LuxurySparse

n=16
Id = IMatrix(1<<n)
V = randn(ComplexF64, 1<<n)
Dv = Diagonal(V)

Pm = pmrand(ComplexF64, 1<<n)
Sp = SparseMatrixCSC(Pm)

Ds = randn(ComplexF64, 1<<10,1<<10)

mats = [Id, Sp, Pm, Dv]
for A in mats
    bg = BenchmarkGroup()
    for B in mats
        if !(A===B===Ds)
            println("============ $(typeof(A))-$(typeof(B)) =============")
            display(@benchmark $A * $B)
            println()
        end
    end
    display(run(bg, verbose=true))
end


text:
# AutoDiff Factory
abstract type AbstractAutoDiffFactory end

struct AutoDiffFactory <: AbstractAutoDiffFactory
    Jgₓ::AutoDiff.StateJacobian
    Jgᵤ::AutoDiff.ControlJacobian
    ∇f::AdjointStackObjective
end

# Counters
abstract type AbstractCounter end

mutable struct NLPCounter <: AbstractCounter
    objective::Int
    gradient::Int
    hessian::Int
    jacobian::Int
    jtprod::Int
    hprod::Int
end
NLPCounter() = NLPCounter(0, 0, 0, 0, 0, 0)

function Base.empty!(c::NLPCounter)
    for attr in fieldnames(NLPCounter)
        setfield!(c, attr, 0)
    end
end

# Active set utils
function _check(val, val_min, val_max)
    violated_inf = findall(val .< val_min)
    violated_sup = findall(val .> val_max)
    n_inf = length(violated_inf)
    n_sup = length(violated_sup)
    err_inf = norm(val_min[violated_inf] .- val[violated_inf], Inf)
    err_sup = norm(val[violated_sup] .- val_max[violated_sup] , Inf)
    return (n_inf, err_inf, n_sup, err_sup)
end

function _inf_pr(nlp::AbstractNLPEvaluator, cons)
    (n_inf, err_inf, n_sup, err_sup) = _check(cons, nlp.g_min, nlp.g_max)
    return max(err_inf, err_sup)
end

function active_set(c, c♭, c♯; tol=1e-8)
    @assert length(c) == length(c♭) == length(c♯)
    active_lb = findall(c .< c♭ .+ tol)
    active_ub = findall(c .> c♯ .- tol)
    return active_lb, active_ub
end

# Scaler utils
abstract type AbstractScaler end

scale_factor(h, tol, η) = max(tol, η / max(1.0, h))

struct MaxScaler{T} <: AbstractScaler
    scale_obj::T
    scale_cons::AbstractVector{T}
    g_min::AbstractVector{T}
    g_max::AbstractVector{T}
end
function MaxScaler(g_min, g_max)
    @assert length(g_min) == length(g_max)
    m = length(g_min)
    sc = similar(g_min) ; fill!(sc, 1.0)
    return MaxScaler(1.0, sc, g_min, g_max)
end
function MaxScaler(nlp::AbstractNLPEvaluator, u0::AbstractVector;
                   η=100.0, tol=1e-8)
    n = n_variables(nlp)
    m = n_constraints(nlp)
    update!(nlp, u0)
    ∇g = similar(u0) ; fill!(∇g, 0)
    gradient!(nlp, ∇g, u0)
    s_obj = scale_factor(norm(∇g, Inf), tol, η)

    # TODO: avoid forming whole Jacobian
    jac = similar(u0, (m, n))
    s_cons = similar(u0, m)
    jacobian!(nlp, jac, u0)
    for i in eachindex(s_cons)
        ∇c = @view jac[i, :]
        s_cons[i] = scale_factor(norm(∇c, Inf), tol, η)
    end

    return MaxScaler(s_obj, s_cons, s_cons .* nlp.g_min, s_cons .* nlp.g_max)
end


text:
@testset "emulate" begin
    # try emulating one cycle
    c8 = Chip()
    CHIP8.loadApplication(c8, "pong2.c8")
    @test_nowarn CHIP8.emulateCycle(c8)

    # now run entire emulation
    @test_nowarn emulate("pong2.c8")
end


text:
# Autogenerated wrapper script for PPL_jll for x86_64-w64-mingw32-cxx03
export libppl, libppl_c, ppl_config, ppl_lcdd, ppl_pips

using GMP_jll
JLLWrappers.@generate_wrapper_header("PPL")
JLLWrappers.@declare_library_product(libppl, "libppl-14.dll")
JLLWrappers.@declare_library_product(libppl_c, "libppl_c-4.dll")
JLLWrappers.@declare_executable_product(ppl_config)
JLLWrappers.@declare_executable_product(ppl_lcdd)
JLLWrappers.@declare_executable_product(ppl_pips)
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll)
    JLLWrappers.@init_library_product(
        libppl,
        "bin\\libppl-14.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libppl_c,
        "bin\\libppl_c-4.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        ppl_config,
        "bin\\ppl-config.exe",
    )

    JLLWrappers.@init_executable_product(
        ppl_lcdd,
        "bin\\ppl_lcdd.exe",
    )

    JLLWrappers.@init_executable_product(
        ppl_pips,
        "bin\\ppl_pips.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using ISA #References: Sandoval, Steven, and Phillip L. De Leon. "The Instantaneous Spectrum: A General Framework for Time-Frequency Analysis." IEEE Transactions on Signal Processing 66.21 (2018): 5679-5693.

t = Array(0.0:0.005:2.0)

using Interact
@manipulate for a= 0:0.05:1, ω = -5:0.1:20, φ = -pi:pi/50:pi
    a₀(t) = a
    ω₀(t) = ω
    φ₀ = φ
    𝐶₀ = Tuple([a₀,ω₀,φ₀])

    isaPlot3d(𝐶₀, t)
end

max_stars_repo_name: kricheso/ISA

max_issues_repo_name: kricheso/ISA

max_forks_repo_name: kricheso/ISA


text:
module AdditionTest

using Test
using ForneyLab
using ForneyLab: outboundType, isApplicable
using ForneyLab: SPAdditionOutNGG, SPAdditionOutNGP, SPAdditionOutNPG, SPAdditionIn1GNG, SPAdditionIn1PNG, SPAdditionIn2GGN, SPAdditionIn2PGN, SPAdditionIn1GNP, SPAdditionIn2GPN, SPAdditionOutNPP, SPAdditionIn1PNP, SPAdditionIn2PPN, MAdditionNGG

@testset "Addition node construction through + syntax" begin
    g = FactorGraph()
    @RV x ~ GaussianMeanVariance(constant(0.0), constant(1.0))
    @RV y ~ GaussianMeanVariance(constant(0.0), constant(1.0))
    @RV z = x + y
    @test isa(z, Variable)
    @test isa(g.nodes[:addition_1], Addition)

    g = FactorGraph()
    @RV x ~ GaussianMeanVariance(constant(0.0), constant(1.0))
    @RV z = x + 1.0
    @test isa(z, Variable)
    @test isa(g.nodes[:addition_1], Addition)
end

@testset "Addition node construction through - syntax" begin
    g = FactorGraph()
    @RV x ~ GaussianMeanVariance(constant(0.0), constant(1.0))
    @RV y ~ GaussianMeanVariance(constant(0.0), constant(1.0))
    @RV z = x - y
    @test isa(z, Variable)
    @test isa(g.nodes[:addition_1], Addition)
    # This syntax works by changing the order interfaces are attached in. Below checks whether the pairings are succesful by
    # matching out and in2 interfaces of the corresponding nodes
    @test g.nodes[:addition_1].i[:out] == g.nodes[:gaussianmeanvariance_1].i[:out].partner
    @test g.nodes[:addition_1].i[:in2] == g.nodes[:gaussianmeanvariance_2].i[:out].partner

    g = FactorGraph()
    @RV x ~ GaussianMeanVariance(constant(0.0), constant(1.0))
    @RV z = x - 1.0
    @test isa(z, Variable)
    @test isa(g.nodes[:addition_1], Addition)
end

#-------------
# Update rules
#-------------

@testset "SPAdditionOutNGG" begin
    @test SPAdditionOutNGG <: SumProductRule{Addition}
    @test outboundType(SPAdditionOutNGG) == Message{GaussianMeanVariance}
    @test isApplicable(SPAdditionOutNGG, [Nothing, Message{Gaussian}, Message{Gaussian}])
    @test !isApplicable(SPAdditionOutNGG, [Nothing, Message{PointMass}, Message{PointMass}])
    @test !isApplicable(SPAdditionOutNGG, [Message{Gaussian}, Message{Gaussian}, Nothing])

    @test ruleSPAdditionOutNGG(nothing, Message(Univariate, GaussianMeanVariance, m=1.0, v=2.0), Message(Univariate, GaussianMeanVariance, m=3.0, v=4.0)) == Message(Univariate, GaussianMeanVariance, m=4.0, v=6.0)
    @test ruleSPAdditionOutNGG(nothing, Message(Multivariate, GaussianMeanVariance, m=[1.0], v=mat(2.0)), Message(Multivariate, GaussianMeanVariance, m=[3.0], v=mat(4.0))) == Message(Multivariate, GaussianMeanVariance, m=[4.0], v=mat(6.0))
end

@testset "SPAdditionIn2GGN" begin
    @test SPAdditionIn2GGN <: SumProductRule{Addition}
    @test outboundType(SPAdditionIn2GGN) == Message{GaussianMeanVariance}
    @test isApplicable(SPAdditionIn2GGN, [Message{Gaussian}, Message{Gaussian}, Nothing])

    @test ruleSPAdditionIn2GGN(Message(Univariate, GaussianMeanVariance, m=3.0, v=4.0), Message(Univariate, GaussianMeanVariance, m=1.0, v=2.0), nothing) == Message(Univariate, GaussianMeanVariance, m=2.0, v=6.0)
    @test ruleSPAdditionIn2GGN(Message(Multivariate, GaussianMeanVariance, m=[3.0], v=mat(4.0)), Message(Multivariate, GaussianMeanVariance, m=[1.0], v=mat(2.0)), nothing) == Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(6.0))
end

@testset "SPAdditionIn1GNG" begin
    @test SPAdditionIn1GNG <: SumProductRule{Addition}
    @test outboundType(SPAdditionIn1GNG) == Message{GaussianMeanVariance}
    @test isApplicable(SPAdditionIn1GNG, [Message{Gaussian}, Nothing, Message{Gaussian}])

    @test ruleSPAdditionIn1GNG(Message(Univariate, GaussianMeanVariance, m=3.0, v=4.0), nothing, Message(Univariate, GaussianMeanVariance, m=1.0, v=2.0)) == Message(Univariate, GaussianMeanVariance, m=2.0, v=6.0)
    @test ruleSPAdditionIn1GNG(Message(Multivariate, GaussianMeanVariance, m=[3.0], v=mat(4.0)), nothing, Message(Multivariate, GaussianMeanVariance, m=[1.0], v=mat(2.0))) == Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(6.0))
end

@testset "SPAdditionOutNGP" begin
    @test SPAdditionOutNGP <: SumProductRule{Addition}
    @test outboundType(SPAdditionOutNGP) == Message{GaussianMeanVariance}
    @test isApplicable(SPAdditionOutNGP, [Nothing, Message{Gaussian}, Message{PointMass}])

    @test ruleSPAdditionOutNGP(nothing, Message(Univariate, GaussianMeanVariance, m=1.0, v=2.0), Message(Univariate, PointMass, m=3.0)) == Message(Univariate, GaussianMeanVariance, m=4.0, v=2.0)
    @test ruleSPAdditionOutNGP(nothing, Message(Multivariate, GaussianMeanVariance, m=[1.0], v=mat(2.0)), Message(Multivariate, PointMass, m=[3.0])) == Message(Multivariate, GaussianMeanVariance, m=[4.0], v=mat(2.0))
end

@testset "SPAdditionOutNPG" begin
    @test SPAdditionOutNPG <: SumProductRule{Addition}
    @test outboundType(SPAdditionOutNPG) == Message{GaussianMeanVariance}
    @test isApplicable(SPAdditionOutNPG, [Nothing, Message{PointMass}, Message{Gaussian}])

    @test ruleSPAdditionOutNPG(nothing, Message(Univariate, PointMass, m=3.0), Message(Univariate, GaussianMeanVariance, m=1.0, v=2.0)) == Message(Univariate, GaussianMeanVariance, m=4.0, v=2.0)
    @test ruleSPAdditionOutNPG(nothing, Message(Multivariate, PointMass, m=[3.0]), Message(Multivariate, GaussianMeanVariance, m=[1.0], v=mat(2.0))) == Message(Multivariate, GaussianMeanVariance, m=[4.0], v=mat(2.0))
end

@testset "SPAdditionIn1PNG" begin
    @test SPAdditionIn1PNG <: SumProductRule{Addition}
    @test outboundType(SPAdditionIn1PNG) == Message{GaussianMeanVariance}
    @test isApplicable(SPAdditionIn1PNG, [Message{PointMass}, Nothing, Message{Gaussian}])

    @test ruleSPAdditionIn1PNG(Message(Univariate, PointMass, m=3.0), nothing, Message(Univariate, GaussianMeanVariance, m=1.0, v=2.0)) == Message(Univariate, GaussianMeanVariance, m=2.0, v=2.0)
    @test ruleSPAdditionIn1PNG(Message(Multivariate, PointMass, m=[3.0]), nothing, Message(Multivariate, GaussianMeanVariance, m=[1.0], v=mat(2.0))) == Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(2.0))
end

@testset "SPAdditionIn2PGN" begin
    @test SPAdditionIn2PGN <: SumProductRule{Addition}
    @test outboundType(SPAdditionIn2PGN) == Message{GaussianMeanVariance}
    @test isApplicable(SPAdditionIn2PGN, [Message{PointMass}, Message{Gaussian}, Nothing])

    @test ruleSPAdditionIn2PGN(Message(Univariate, PointMass, m=3.0), Message(Univariate, GaussianMeanVariance, m=1.0, v=2.0), nothing) == Message(Univariate, GaussianMeanVariance, m=2.0, v=2.0)
    @test ruleSPAdditionIn2PGN(Message(Multivariate, PointMass, m=[3.0]), Message(Multivariate, GaussianMeanVariance, m=[1.0], v=mat(2.0)), nothing) == Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(2.0))
end

@testset "SPAdditionIn1GNP" begin
    @test SPAdditionIn1GNP <: SumProductRule{Addition}
    @test outboundType(SPAdditionIn1GNP) == Message{GaussianMeanVariance}
    @test isApplicable(SPAdditionIn1GNP, [Message{Gaussian}, Nothing, Message{PointMass}])

    @test ruleSPAdditionIn1GNP(Message(Univariate, GaussianMeanVariance, m=3.0, v=4.0), nothing, Message(Univariate, PointMass, m=1.0)) == Message(Univariate, GaussianMeanVariance, m=2.0, v=4.0)
    @test ruleSPAdditionIn1GNP(Message(Multivariate, GaussianMeanVariance, m=[3.0], v=mat(4.0)), nothing, Message(Multivariate, PointMass, m=[1.0])) == Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(4.0))
end

@testset "SPAdditionIn2GPN" begin
    @test SPAdditionIn2GPN <: SumProductRule{Addition}
    @test outboundType(SPAdditionIn2GPN) == Message{GaussianMeanVariance}
    @test isApplicable(SPAdditionIn2GPN, [Message{Gaussian}, Message{PointMass}, Nothing])

    @test ruleSPAdditionIn2GPN(Message(Univariate, GaussianMeanVariance, m=3.0, v=4.0), Message(Univariate, PointMass, m=1.0), nothing) == Message(Univariate, GaussianMeanVariance, m=2.0, v=4.0)
    @test ruleSPAdditionIn2GPN(Message(Multivariate, GaussianMeanVariance, m=[3.0], v=mat(4.0)), Message(Multivariate, PointMass, m=[1.0]), nothing) == Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(4.0))
end

@testset "SPAdditionOutNPP" begin
    @test SPAdditionOutNPP <: SumProductRule{Addition}
    @test outboundType(SPAdditionOutNPP) == Message{PointMass}
    @test isApplicable(SPAdditionOutNPP, [Nothing, Message{PointMass}, Message{PointMass}])

    @test ruleSPAdditionOutNPP(nothing, Message(Univariate, PointMass, m=1.0), Message(Univariate, PointMass, m=3.0)) == Message(Univariate, PointMass, m=4.0)
    @test ruleSPAdditionOutNPP(nothing, Message(Multivariate, PointMass, m=[1.0]), Message(Multivariate, PointMass, m=[3.0])) == Message(Multivariate, PointMass, m=[4.0])
end

@testset "SPAdditionIn2PPN" begin
    @test SPAdditionIn2PPN <: SumProductRule{Addition}
    @test outboundType(SPAdditionIn2PPN) == Message{PointMass}
    @test isApplicable(SPAdditionIn2PPN, [Message{PointMass}, Message{PointMass}, Nothing])

    @test ruleSPAdditionIn2PPN(Message(Univariate, PointMass, m=3.0), Message(Univariate, PointMass, m=1.0), nothing) == Message(Univariate, PointMass, m=2.0)
    @test ruleSPAdditionIn2PPN(Message(Multivariate, PointMass, m=[3.0]), Message(Multivariate, PointMass, m=[1.0]), nothing) == Message(Multivariate, PointMass, m=[2.0])
end

@testset "SPAdditionIn1PNP" begin
    @test SPAdditionIn1PNP <: SumProductRule{Addition}
    @test outboundType(SPAdditionIn1PNP) == Message{PointMass}
    @test isApplicable(SPAdditionIn1PNP, [Message{PointMass}, Nothing, Message{PointMass}])

    @test ruleSPAdditionIn1PNP(Message(Univariate, PointMass, m=3.0), nothing, Message(Univariate, PointMass, m=1.0)) == Message(Univariate, PointMass, m=2.0)
    @test ruleSPAdditionIn1PNP(Message(Multivariate, PointMass, m=[3.0]), nothing, Message(Multivariate, PointMass, m=[1.0])) == Message(Multivariate, PointMass, m=[2.0])
end

@testset "MAdditionNGG" begin
    @test MAdditionNGG <: MarginalRule{Addition}
    @test isApplicable(MAdditionNGG, [Nothing, Message{Gaussian}, Message{Gaussian}])

    @test ruleMAdditionNGG(Message(Univariate, GaussianWeightedMeanPrecision, xi=1.0, w=2.0), Message(Univariate, GaussianWeightedMeanPrecision, xi=3.0, w=4.0), Message(Univariate, GaussianWeightedMeanPrecision, xi=5.0, w=6.0)) == ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[4.0, 6.0], w=[6.0 2.0; 2.0 8.0])
    @test ruleMAdditionNGG(Message(Multivariate, GaussianWeightedMeanPrecision, xi=[1.0], w=mat(2.0)), Message(Multivariate, GaussianWeightedMeanPrecision, xi=[3.0], w=mat(4.0)), Message(Multivariate, GaussianWeightedMeanPrecision, xi=[5.0], w=mat(6.0))) == ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[4.0, 6.0], w=[6.0 2.0; 2.0 8.0])
end

end # module


text:
module CommunalHelperDreamFloatySpaceBlock

using ..Ahorn, Maple
using Ahorn.CommunalHelper

@mapdef Entity "CommunalHelper/DreamFloatySpaceBlock" DreamFloatySpaceBlock(
    x::Integer,
    y::Integer,
    width::Integer=Maple.defaultBlockWidth,
    height::Integer=Maple.defaultBlockHeight,
    featherMode::Bool=false,
    oneUse::Bool=false,
    refillCount::Integer=-1,
    below::Bool=false,
)

const placements = Ahorn.PlacementDict(
    "Dream Floaty Space Block (Communal Helper)" => Ahorn.EntityPlacement(
        DreamFloatySpaceBlock,
        "rectangle"
    ),
)

Ahorn.minimumSize(entity::DreamFloatySpaceBlock) = 8, 8
Ahorn.resizable(entity::DreamFloatySpaceBlock) = true, true

Ahorn.selection(entity::DreamFloatySpaceBlock) = Ahorn.getEntityRectangle(entity)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::DreamFloatySpaceBlock, room::Maple.Room)
    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    renderDreamBlock(ctx, 0, 0, width, height, entity.data)
end

end

max_stars_repo_name: vitellaryjr/CommunalHelper

max_issues_repo_name: vitellaryjr/CommunalHelper

max_forks_repo_name: vitellaryjr/CommunalHelper


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    LWR(var₁=>param₁, var₂=>param₂, ...)

Locally weighted regression estimation solver.

## Parameters

* `neighbors` - Number of neighbors (default to all the data)
* `distance`  - A distance from Distances.jl (default to `Euclidean()`)
* `weightfun` - Weighting function (default to `exp(-3*h^2/2)`)

### References

* Stone 1977. *Consistent non-parametric regression.*
* Cleveland 1979. *Robust locally weighted regression and smoothing scatterplots.*
* Cleveland & Grosse 1991. *Computational methods for local regression.*
"""
@estimsolver LWR begin
  @param neighbors = nothing
  @param distance = Euclidean()
  @param weightfun = h -> exp(-3*h^2)
end

function solve(problem::EstimationProblem, solver::LWR)
  # retrieve problem info
  pdata = data(problem)
  pdomain = domain(problem)

  mactypeof = Dict(name(v) => mactype(v) for v in variables(problem))

  # result for each variable
  μs = []; σs = []

  for covars in covariables(problem, solver)
    for var in covars.names
      # get user parameters
      varparams = covars.params[(var,)]

      # determine value type
      V = mactypeof[var]

      # retrieve non-missing data
      locs = findall(!ismissing, pdata[var])
      𝒟 = view(pdata, locs)
      n = nelements(𝒟)

      # determine number of nearest neighbors to use
      k = isnothing(varparams.neighbors) ? n : varparams.neighbors

      # determine distance type
      D = varparams.distance

      # weight function
      w = varparams.weightfun

      @assert n > 0 "estimation requires data"

      @assert k ≤ n "invalid number of neighbors"

      # fit search tree
      X = [coordinates(centroid(𝒟, i)) for i in 1:nelements(𝒟)]
      if D isa NearestNeighbors.MinkowskiMetric
        tree = KDTree(X, D)
      else
        tree = BallTree(X, D)
      end

      # lookup non-missing values
      z = 𝒟[var]

      # pre-allocate memory for results
      varμ = Vector{V}(undef, nelements(pdomain))
      varσ = Vector{V}(undef, nelements(pdomain))

      # estimation loop
      for loc in traverse(pdomain, LinearPath())
        x = coordinates(centroid(pdomain, loc))

        # find neighbors
        is, ds = knn(tree, x, k)
        δs = ds ./ maximum(ds)

        # weighted least-squares
        Wₗ = Diagonal(w.(δs))
        Xₗ = [ones(eltype(x), k) reduce(hcat, X[is])']
        zₗ = view(z, is)
        θₗ = Xₗ'*Wₗ*Xₗ \ Xₗ'*Wₗ*zₗ

        # linear combination of response values
        xₒ = [one(eltype(x)); x]
        ẑₒ = θₗ ⋅ xₒ
        rₗ = Wₗ*Xₗ*(Xₗ'*Wₗ*Xₗ\xₒ)
        r̂ₒ = norm(rₗ)

        varμ[loc] = ẑₒ
        varσ[loc] = r̂ₒ
      end

      push!(μs, var => varμ)
      push!(σs, Symbol(var,"_variance") => varσ)
    end
  end

  georef((; μs..., σs...), pdomain)
end


text:
## CHOLMOD
const CHOLMOD_TRUE  = int32(1)
const CHOLMOD_FALSE = int32(0)

## itype defines the types of integer used:
const CHOLMOD_INT  = int32(0)  # all integer arrays are int
const CHOLMOD_LONG = int32(2)  # all integer arrays are UF_long
ityp(::Type{Int32}) = CHOLMOD_INT
ityp(::Type{Int64}) = CHOLMOD_LONG

## dtype defines what the numerical type is (double or float):
const CHOLMOD_DOUBLE = int32(0)        # all numerical values are double
const CHOLMOD_SINGLE = int32(1)        # all numerical values are float
dtyp(::Type{Float32}) = CHOLMOD_SINGLE
dtyp(::Type{Float64}) = CHOLMOD_DOUBLE
dtyp(::Type{Complex64}) = CHOLMOD_SINGLE
dtyp(::Type{Complex128}) = CHOLMOD_DOUBLE

## xtype defines the kind of numerical values used:
const CHOLMOD_PATTERN = int32(0)       # pattern only, no numerical values
const CHOLMOD_REAL    = int32(1)       # a real matrix
const CHOLMOD_COMPLEX = int32(2)       # a complex matrix (ANSI C99 compatible)
const CHOLMOD_ZOMPLEX = int32(3)       # a complex matrix (MATLAB compatible)
xtyp(::Type{Float32}) = CHOLMOD_REAL
xtyp(::Type{Float64}) = CHOLMOD_REAL
xtyp(::Type{Complex64}) = CHOLMOD_COMPLEX
xtyp(::Type{Complex128}) = CHOLMOD_COMPLEX

## Scaling modes, selected by the scale input parameter:
const CHOLMOD_SCALAR = int32(0)        # A = s*A
const CHOLMOD_ROW    = int32(1)        # A = diag(s)*A
const CHOLMOD_COL    = int32(2)        # A = A*diag(s)
const CHOLMOD_SYM    = int32(3)        # A = diag(s)*A*diag(s)

## Types of systems to solve
const CHOLMOD_A    = int32(0)          # solve Ax=b
const CHOLMOD_LDLt = int32(1)          # solve LDL'x=b
const CHOLMOD_LD   = int32(2)          # solve LDx=b
const CHOLMOD_DLt  = int32(3)          # solve DL'x=b
const CHOLMOD_L    = int32(4)          # solve Lx=b
const CHOLMOD_Lt   = int32(5)          # solve L'x=b
const CHOLMOD_D    = int32(6)          # solve Dx=b
const CHOLMOD_P    = int32(7)          # permute x=Px
const CHOLMOD_Pt   = int32(8)          # permute x=P'x

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
module LatticeVis

using Quaternions, GeometryTypes, Reactive, Colors
using GLFW, GLWindow, GLAbstraction, GLVisualize

export RGBA, Signal
export Point3f0

import Base: merge, expand, push!, ==
import GLAbstraction: center!

include("Lattice.jl")
export sc, bcc, fcc, diamond
export merge, expand, center!

include("Window.jl")
export init_window, camera

include("KeyBinds.jl")
export Key, HotKeys, push!
export init_screenshots!, init_QWEASD_controls!, init_angles

include("Visualize.jl")
export glsphere, draw

end # module


text:
function t2i(x::Int, sample_rate::Real) # Time to index
    ceil(Int64, x * sample_rate + 1)
end

function rs_bnd_both(
    x::Neuro1DRealSignal,
    side::Symbol;
    position::Tuple{Int, Int})

    p = t2i.(position, x.sample_rate)

    Neuro1DRealSignal(
        signal = x.signal[minimum(p):maximum(p), :],
        sample_rate = x.sample_rate
    )
end

function rs_bnd_single(
    x::Neuro1DRealSignal;
    position::Int)

    p = t2i(position, x.sample_rate)

    Neuro1DRealSignal(
        signal = side == :right ? x.signal[1:p, :] : x.signal[p:end, :],
        sample_rate = x.sample_rate
    )
end

function rs_bnd(
    x::Neuro1DRealSignal,
    side = :both;
    kwargs...)::Neuro1DRealSignal

    if side == :both
        rs_bnd_both(x, side; kwargs...)
    elseif method in (:left, :right)
        rs_bnd_single(x; kwargs...)
    else
        throw(ArgumentError("invalid side"))
    end
end


text:
using FastTransforms, Compat
using Compat.Test

import FastTransforms: normalizecolumns!, maxcolnorm

@testset "Spherical harmonic API" begin
    let n = 512
        A = sphrandn(Float64, n, n);
        normalizecolumns!(A);

        B = sph2fourier(A)
        C = fourier2sph(B)
        println("The backward difference between slow plan and original: ", maxcolnorm(A-C))

        P = plan_sph2fourier(A)
        B = P*A
        C = P\B

        println("The backward difference between slow plan and original: ", maxcolnorm(A-C))

        n = 1024
        A = sphrandn(Float64, n, n);
        normalizecolumns!(A);

        B = sph2fourier(A; sketch = :none)
        C = fourier2sph(B; sketch = :none)
        println("The backward difference between thin plan and original: ", maxcolnorm(A-C))

        P = plan_sph2fourier(A; sketch = :none)
        B = P*A
        C = P\B

        println("The backward difference between thin plan and original: ", maxcolnorm(A-C))
    end
end


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

mutable struct PatchSettings <: SwaggerModel
    patchMode::Any # spec type: Union{ Nothing, String } # spec name: patchMode

    function PatchSettings(;patchMode=nothing)
        o = new()
        validate_property(PatchSettings, Symbol("patchMode"), patchMode)
        setfield!(o, Symbol("patchMode"), patchMode)
        o
    end
end # type PatchSettings

const _property_map_PatchSettings = Dict{Symbol,Symbol}(Symbol("patchMode")=>Symbol("patchMode"))
const _property_types_PatchSettings = Dict{Symbol,String}(Symbol("patchMode")=>"String")
Base.propertynames(::Type{ PatchSettings }) = collect(keys(_property_map_PatchSettings))
Swagger.property_type(::Type{ PatchSettings }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_PatchSettings[name]))}
Swagger.field_name(::Type{ PatchSettings }, property_name::Symbol) =  _property_map_PatchSettings[property_name]

const _allowed_PatchSettings_patchMode = ["Manual", "AutomaticByOS", "AutomaticByPlatform"]

function check_required(o::PatchSettings)
    true
end

function validate_property(::Type{ PatchSettings }, name::Symbol, val)
    if name === Symbol("patchMode")
        Swagger.validate_param(name, "PatchSettings", :enum, val, _allowed_PatchSettings_patchMode)
    end
end


text:
using Graphics: @mustimplement

export Motor, StepperMotor, MotorDirection

abstract type Motor <: Device end
abstract type StepperMotor <: Motor end

@enum MotorDirection begin
  MOTOR_FORWARD
  MOTOR_BACKWARD
  MOTOR_STILL
end

Base.close(motor::Motor) = nothing

@mustimplement direction(motor::Motor)
@mustimplement direction(motor::Motor, dir::MotorDirection)

getMotors(scanner::MPIScanner) = getDevices(scanner, Motor)
function getMotor(scanner::MPIScanner)
  motors = getMotors(scanner)
  if length(motors) > 1
    error("The scanner has more than one motor device. Therefore, a single motor cannot be retrieved unambiguously.")
  else
    return motors[1]
  end
end

"""
Returns x,y, and z values and applies a coordinate transformation
"""
function getXYZValues(gauss::GaussMeter)
  values = [getXValue(gauss), getYValue(gauss), getZValue(gauss)]
  return gauss.params.coordinateTransformation*values
end

include("SimulatedMotor.jl")
#include("TinkerforgeSilentStepper.jl")


text:
GTR(α::F, β::F, γ::F, δ::F, ϵ::F, η::F,
    πA::F, πC::F, πG::F, πT::F;
    safe::Bool=true) where F <: Float64 =
  GTRabs(α, β, γ, δ, ϵ, η, πA, πC, πG, πT, safe=safe)

GTR(α::F, β::F, γ::F, δ::F, ϵ::F,
    πA::F, πC::F, πG::F, πT::F;
    safe::Bool=true) where F <: Float64 =
  GTRrel(α, β, γ, δ, ϵ, πA, πC, πG, πT, safe=safe)

function GTR(θ_vec::A,
               π_vec::A;
               safe::Bool=true) where A <: AbstractArray
  if safe && length(π_vec) != 4
    error("Incorrect base frequency vector length")
  end
  if length(θ_vec) == 5
    return GTRrel(θ_vec[1], θ_vec[2], θ_vec[3], θ_vec[4], θ_vec[5],
                  π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T],
                  safe=safe)
  elseif length(θ_vec) == 6
    return GTRabs(θ_vec[1], θ_vec[2], θ_vec[3], θ_vec[4], θ_vec[5], θ_vec[6],
                  π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T],
                  safe=safe)
  else
    error("Parameter vector length incompatiable with absolute or relative rate form of substitution model")
  end
end

function GTRrel(θ_vec::A,
                π_vec::A;
                safe::Bool=true) where A <: AbstractArray
  if safe
    if length(θ_vec) != 5
      error("Incorrect parameter vector length")
    elseif length(π_vec) != 4
      error("Incorrect base frequency vector length")
    end
  end
  return GTRrel(θ_vec[1], θ_vec[2], θ_vec[3], θ_vec[4], θ_vec[5],
                π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T],
                safe=safe)
end

function GTRabs(θ_vec::A,
                π_vec::A;
                safe::Bool=true) where A <: AbstractArray
  if safe
    if length(θ_vec) != 6
      error("Incorrect parameter vector length")
    elseif length(π_vec) != 4
      error("Incorrect base frequency vector length")
    end
  end
  return GTRabs(θ_vec[1], θ_vec[2], θ_vec[3], θ_vec[4], θ_vec[5], θ_vec[6],
                π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T],
                safe=safe)
end


text:
###############################################################################
#
#   Absolute Power Series Ring
#
###############################################################################

doc"""
    *(a::fmpz, b::AbsSeriesElem)
> Return $a\times b$.
"""
function *(a::fmpz, b::AbsSeriesElem)
   len = length(b)
   z = parent(b)()
   fit!(z, len)
   set_prec!(z, precision(b))
   for i = 1:len
      z = setcoeff!(z, i - 1, a*coeff(b, i - 1))
   end
   set_length!(z, normalise(z, len))
   return z
end

doc"""
    *(a::AbsSeriesElem, b::fmpz)
> Return $a\times b$.
"""
*(a::AbsSeriesElem, b::fmpz) = b*a

doc"""
    ==(x::AbsSeriesElem, y::fmpz)
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
==(x::AbsSeriesElem, y::fmpz) = precision(x) == 0 || ((length(x) == 0 && iszero(y))
                                       || (length(x) == 1 && coeff(x, 0) == y))

doc"""
    ==(x::fmpz, y::AbsSeriesElem)
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
==(x::fmpz, y::AbsSeriesElem) = y == x

doc"""
    divexact(a::AbsSeriesElem, b::fmpz)
> Return $a/b$ where the quotient is expected to be exact.
"""
function divexact(x::AbsSeriesElem, y::fmpz)
   iszero(y) && throw(DivideError())
   lenx = length(x)
   z = parent(x)()
   fit!(z, lenx)
   set_prec!(z, precision(x))
   for i = 1:lenx
      z = setcoeff!(z, i - 1, divexact(coeff(x, i - 1), y))
   end
   return z
end

function (a::Generic.AbsSeriesRing{T})(b::fmpz) where {T <: RingElement}
   if b == 0
      z = Generic.AbsSeries{T}(Array{T}(0), 0, a.prec_max)
   else
      z = Generic.AbsSeries{T}([base_ring(a)(b)], 1, a.prec_max)
   end
   z.parent = a
   return z
end

###############################################################################
#
#   Relative Power Series Ring
#
###############################################################################

doc"""
    *(a::fmpz, b::RelSeriesElem)
> Return $a\times b$.
"""
function *(a::fmpz, b::RelSeriesElem)
   len = pol_length(b)
   z = parent(b)()
   fit!(z, len)
   set_prec!(z, precision(b))
   set_val!(z, valuation(b))
   for i = 1:len
      z = setcoeff!(z, i - 1, a*polcoeff(b, i - 1))
   end
   set_length!(z, normalise(z, len))
   renormalize!(z)
   return z
end

doc"""
    *(a::RelSeriesElem, b::fmpz)
> Return $a\times b$.
"""
*(a::RelSeriesElem, b::fmpz) = b*a

doc"""
    ==(x::RelSeriesElem, y::fmpz)
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
==(x::RelSeriesElem, y::fmpz) = precision(x) == 0 ||
                  ((pol_length(x) == 0 && iszero(y)) || (pol_length(x) == 1 &&
                    valuation(x) == 0 && polcoeff(x, 0) == y))

doc"""
    ==(x::fmpz, y::RelSeriesElem)
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
==(x::fmpz, y::RelSeriesElem) = y == x

doc"""
    divexact(a::RelSeriesElem, b::fmpz)
> Return $a/b$ where the quotient is expected to be exact.
"""
function divexact(x::RelSeriesElem, y::fmpz)
   iszero(y) && throw(DivideError())
   lenx = pol_length(x)
   z = parent(x)()
   fit!(z, lenx)
   set_prec!(z, precision(x))
   set_val!(z, valuation(x))
   for i = 1:lenx
      z = setcoeff!(z, i - 1, divexact(polcoeff(x, i - 1), y))
   end
   return z
end

function (a::Generic.RelSeriesRing{T})(b::fmpz) where {T <: RingElement}
   if iszero(b)
      z = Generic.RelSeries{T}(Array{T}(0), 0, a.prec_max, a.prec_max)
   else
      z = Generic.RelSeries{T}([base_ring(a)(b)], 1, a.prec_max, 0)
   end
   z.parent = a
   return z
end

###############################################################################
#
#   Polynomial Ring
#
###############################################################################

doc"""
    *(a::fmpz, b::PolyElem)
> Return $a\times b$.
"""
function *(a::fmpz, b::PolyElem)
   len = length(b)
   z = parent(b)()
   fit!(z, len)
   for i = 1:len
      z = setcoeff!(z, i - 1, a*coeff(b, i - 1))
   end
   set_length!(z, normalise(z, len))
   return z
end

doc"""
    *(a::PolyElem, b::fmpz)
> Return $a\times b$.
"""
*(a::PolyElem, b::fmpz) = b*a

doc"""
    ==(x::PolyElem, y::fmpz)
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
==(x::PolyElem, y::fmpz) = ((length(x) == 0 && iszero(y))
                        || (length(x) == 1 && coeff(x, 0) == y))

doc"""
    ==(x::fmpz, y::PolyElem)
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
==(x::fmpz, y::PolyElem) = y == x

doc"""
    divexact(a::PolyElem, b::fmpz)
> Return $a/b$ where the quotient is expected to be exact.
"""
function divexact(a::PolyElem, b::fmpz)
   iszero(b) && throw(DivideError())
   z = parent(a)()
   fit!(z, length(a))
   for i = 1:length(a)
      z = setcoeff!(z, i - 1, divexact(coeff(a, i - 1), b))
   end
   set_length!(z, length(a))
   return z
end

# ambiguities

function *(a::fmpz, b::PolyElem{fmpz})
   len = length(b)
   z = parent(b)()
   fit!(z, len)
   for i = 1:len
      z = setcoeff!(z, i - 1, a*coeff(b, i - 1))
   end
   set_length!(z, normalise(z, len))
   return z
end

*(a::PolyElem{fmpz}, b::fmpz) = b*a

==(x::PolyElem{fmpz}, y::fmpz) = ((length(x) == 0 && iszero(y))
                        || (length(x) == 1 && coeff(x, 0) == y))

==(x::fmpz, y::PolyElem{fmpz}) = y == x

function divexact(a::PolyElem{fmpz}, b::fmpz)
   iszero(b) && throw(DivideError())
   z = parent(a)()
   fit!(z, length(a))
   for i = 1:length(a)
      z = setcoeff!(z, i - 1, divexact(coeff(a, i - 1), b))
   end
   set_length!(z, length(a))
   return z
end

###############################################################################
#
#   Residue Ring
#
###############################################################################

doc"""
    *(a::ResElem, b::fmpz)
> Return $a\times b$.
"""
*(a::ResElem, b::fmpz) = parent(a)(data(a) * b)

doc"""
    *(a::fmpz, b::ResElem)
> Return $a\times b$.
"""
*(a::fmpz, b::ResElem) = parent(b)(a * data(b))

doc"""
    +(a::ResElem, b::fmpz)
> Return $a + b$.
"""
+(a::ResElem, b::fmpz) = parent(a)(data(a) + b)

doc"""
    +(a::fmpz, b::ResElem)
> Return $a + b$.
"""
+(a::fmpz, b::ResElem) = parent(b)(a + data(b))

doc"""
    -(a::ResElem, b::fmpz)
> Return $a - b$.
"""
-(a::ResElem, b::fmpz) = parent(a)(data(a) - b)

doc"""
    -(a::fmpz, b::ResElem)
> Return $a - b$.
"""
-(a::fmpz, b::ResElem) = parent(b)(a - data(b))

doc"""
    ==(x::ResElem, y::fmpz)
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
function ==(a::ResElem, b::fmpz)
   z = base_ring(a)(b)
   return data(a) == mod(z, modulus(a))
end

doc"""
    ==(x::fmpz, y::ResElem)
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
function ==(a::fmpz, b::ResElem)
   z = base_ring(b)(a)
   return data(b) == mod(z, modulus(b))
end

function (a::Generic.ResRing{T})(b::fmpz) where {T <: RingElement}
   z = Generic.Res{T}(mod(base_ring(a)(b), modulus(a)))
   z.parent = a
   return z
end

function (a::Generic.ResRing{fmpz})(b::fmpz)
   z = Generic.Res{fmpz}(mod(base_ring(a)(b), modulus(a)))
   z.parent = a
   return z
end

# ambiguities

*(a::ResElem{fmpz}, b::fmpz) = parent(a)(data(a) * b)

*(a::fmpz, b::ResElem{fmpz}) = b*a

+(a::ResElem{fmpz}, b::fmpz) = parent(a)(data(a) + b)

+(a::fmpz, b::ResElem{fmpz}) = b + a

-(a::ResElem{fmpz}, b::fmpz) = parent(a)(data(a) - b)

-(a::fmpz, b::ResElem{fmpz}) = parent(b)(a - data(b))

function ==(a::ResElem{fmpz}, b::fmpz)
   z = base_ring(a)(b)
   return data(a) == mod(z, modulus(a))
end

==(a::fmpz, b::ResElem{fmpz}) = b == a

###############################################################################
#
#   Multivariate Polynomial Ring
#
###############################################################################

function *(a::Generic.MPoly{T}, n::fmpz) where T <: RingElem
   N = size(a.exps, 1)
   r = parent(a)()
   fit!(r, length(a))
   j = 1
   for i = 1:length(a)
      c = a.coeffs[i]*n
      if c != 0
         r.coeffs[j] = c
         monomial_set!(r.exps, j, a.exps, i, N)
         j += 1
      end
   end
   r.length = j - 1
   resize!(r.coeffs, r.length)
   return r
end

*(n::fmpz, a::Generic.MPoly{T}) where T <: RingElem = a*n

function ==(a::Generic.MPoly{T}, n::fmpz) where T <: RingElem
   N = size(a.exps, 1)
   if n == 0
      return a.length == 0
   elseif a.length == 1
      return a.coeffs[1] == n && monomial_iszero(a.exps, 1, N)
   end
   return false
end

function evaluate(a::Generic.MPoly{T}, A::Array{fmpz, 1}) where {T <: RingElement}
   if iszero(a)
      return base_ring(a)()
   end
   N = size(a.exps, 1)
   ord = parent(a).ord
   if ord == :lex
      start_var = 1
   else
      start_var = 2
   end
   while a.length > 1 || (a.length == 1 && !monomial_iszero(a.exps, a.length, N))
      k = main_variable(a, start_var)
      p = main_variable_extract(a, k)
      a = evaluate(p, A[k - start_var + 1])
   end
   if a.length == 0
      return base_ring(a)()
   else
      return a.coeffs[1]
   end
end

###############################################################################
#
#   Sparse Polynomial Ring
#
###############################################################################

function *(a::Generic.SparsePoly{T}, n::fmpz) where T <: RingElem
   r = parent(a)()
   fit!(r, length(a))
   j = 1
   for i = 1:length(a)
      c = a.coeffs[i]*n
      if c != 0
         r.coeffs[j] = c
         r.exps[j] = a.exps[i]
         j += 1
      end
   end
   r.length = j - 1
   return r
end

*(n::fmpz, a::Generic.SparsePoly{T}) where T <: RingElem = a*n

function ==(a::Generic.SparsePoly{T}, b::fmpz) where T <: RingElem
   return length(a) == 0 ? b == 0 : a.length == 1 &
          a.exps[1] == 0 && a.coeffs[1] == b
end

==(a::fmpz, b::Generic.SparsePoly{T}) where T <: RingElem = b == a

function divexact(a::Generic.SparsePoly{T}, b::fmpz) where T <: RingElem
   len = length(a)
   exps = deepcopy(a.exps)
   coeffs = [divexact(a.coeffs[i], b) for i in 1:len]
   return parent(a)(coeffs, exps)
end

###############################################################################
#
#   Matrix Space
#
###############################################################################

doc"""
    *(x::fmpz, y::MatElem)
> Return $x\times y$.
"""
function *(x::fmpz, y::MatElem)
   z = similar(y)
   for i = 1:rows(y)
      for j = 1:cols(y)
         z[i, j] = x*y[i, j]
      end
   end
   return z
end

doc"""
    *(x::MatElem, y::fmpz)
> Return $x\times y$.
"""
*(x::MatElem, y::fmpz) = y*x

doc"""
    +(x::fmpz, y::MatElem)
> Return $S(x) + y$ where $S$ is the parent of $y$.
"""
function +(x::fmpz, y::MatElem)
   z = similar(y)
   R = base_ring(y)
   for i = 1:rows(y)
      for j = 1:cols(y)
         if i != j
            z[i, j] = deepcopy(y[i, j])
         else
            z[i, j] = y[i, j] + R(x)
         end
      end
   end
   return z
end

doc"""
    +(x::MatElem, y::fmpz)
> Return $x + S(y)$ where $S$ is the parent of $x$.
"""
+(x::MatElem, y::fmpz) = y + x

doc"""
    -(x::fmpz, y::MatElem)
> Return $S(x) - y$ where $S$ is the parent of $y$.
"""
function -(x::fmpz, y::MatElem)
   z = similar(y)
   R = base_ring(y)
   for i = 1:rows(y)
      for j = 1:cols(y)
         if i != j
            z[i, j] = -y[i, j]
         else
            z[i, j] = x - y[i, j]
         end
      end
   end
   return z
end

doc"""
    -(x::MatElem, y::fmpz)
> Return $x - S(y)$, where $S$ is the parent of $x$
"""
function -(x::MatElem, y::fmpz)
   z = similar(x)
   R = base_ring(x)
   for i = 1:rows(x)
      for j = 1:cols(x)
         if i != j
            z[i, j] = deepcopy(x[i, j])
         else
            z[i, j] = x[i, j] - y
         end
      end
   end
   return z
end

doc"""
    ==(x::MatElem, y::fmpz)
> Return `true` if $x == S(y)$ arithmetically, where $S$ is the parent of $x$,
> otherwise return `false`.
"""
function ==(x::MatElem, y::fmpz)
   for i = 1:min(rows(x), cols(x))
      if x[i, i] != y
         return false
      end
   end
   for i = 1:rows(x)
      for j = 1:cols(x)
         if i != j && !iszero(x[i, j])
            return false
         end
      end
   end
   return true
end

doc"""
    ==(x::fmpz, y::MatElem)
> Return `true` if $S(x) == y$ arithmetically, where $S$ is the parent of $y$,
> otherwise return `false`.
"""
==(x::fmpz, y::MatElem) = y == x

doc"""
    divexact(x::MatElem, y::fmpz)
> Return $x/y$, i.e. the matrix where each of the entries has been divided by
> $y$. Each division is expected to be exact.
"""
function divexact(x::MatElem, y::fmpz)
   z = similar(x)
   for i = 1:rows(x)
      for j = 1:cols(x)
         z[i, j] = divexact(x[i, j], y)
      end
   end
   return z
end

function (a::Generic.MatSpace{T})(b::fmpz_mat) where {T <: RingElement}
  if a.rows != rows(b) || a.cols != cols(b)
    error("incompatible matrix dimensions")
  end
  A = a()
  R = base_ring(a)
  for i=1:a.rows
    for j=1:a.cols
      A[i,j] = R(b[i,j])
    end
  end
  return A
end

###############################################################################
#
#   Residue Ring
#
###############################################################################

###############################################################################
#
#   Fraction Field
#
###############################################################################

//(x::T, y::fmpz) where {T <: RingElem} = x//parent(x)(y)

//(x::fmpz, y::T) where {T <: RingElem} = parent(y)(x)//y

doc"""
    *(a::FracElem, b::fmpz)
> Return $a\times b$.
"""
function *(a::FracElem, b::fmpz)
   c = base_ring(a)(b)
   g = gcd(denominator(a), c)
   n = numerator(a)*divexact(c, g)
   d = divexact(denominator(a), g)
   return parent(a)(n, d)
end

doc"""
    *(a::fmpz, b::FracElem)
> Return $a\times b$.
"""
function *(a::fmpz, b::FracElem)
   c = base_ring(b)(a)
   g = gcd(denominator(b), c)
   n = numerator(b)*divexact(c, g)
   d = divexact(denominator(b), g)
   return parent(b)(n, d)
end

doc"""
    +(a::FracElem, b::fmpz)
> Return $a + b$.
"""
function +(a::FracElem, b::fmpz)
   n = numerator(a) + denominator(a)*b
   d = denominator(a)
   g = gcd(n, d)
   return parent(a)(divexact(n, g), divexact(d, g))
end

doc"""
    -(a::FracElem, b::fmpz)
> Return $a - b$.
"""
function -(a::FracElem, b::fmpz)
   n = numerator(a) - denominator(a)*b
   d = denominator(a)
   g = gcd(n, d)
   return parent(a)(divexact(n, g), divexact(d, g))
end

doc"""
    +(a::fmpz, b::FracElem)
> Return $a + b$.
"""
+(a::fmpz, b::FracElem) = b + a

doc"""
    -(a::fmpz, b::FracElem)
> Return $a - b$.
"""
function -(a::fmpz, b::FracElem)
   n = a*denominator(b) - numerator(b)
   d = denominator(b)
   g = gcd(n, d)
   return parent(b)(divexact(n, g), divexact(d, g))
end

doc"""
    ==(x::FracElem, y::fmpz)
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
function ==(x::FracElem, y::fmpz)
   return (isone(denominator(x)) && numerator(x) == y) || (numerator(x) == denominator(x)*y)
end

doc"""
    ==(x::fmpz, y::FracElem)
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
==(x::fmpz, y::FracElem) = y == x

doc"""
    divexact(a::FracElem, b::fmpz)
> Return $a/b$.
"""
function divexact(a::FracElem, b::fmpz)
   iszero(b) && throw(DivideError())
   c = base_ring(a)(b)
   g = gcd(numerator(a), c)
   n = divexact(numerator(a), g)
   d = denominator(a)*divexact(c, g)
   return parent(a)(n, d)
end

doc"""
    divexact(a::fmpz, b::FracElem)
> Return $a/b$.
"""
function divexact(a::fmpz, b::FracElem)
   iszero(b) && throw(DivideError())
   c = base_ring(b)(a)
   g = gcd(numerator(b), c)
   n = denominator(b)*divexact(c, g)
   d = divexact(numerator(b), g)
   return parent(b)(n, d)
end

function (a::Generic.FracField{T})(b::fmpz) where {T <: RingElement}
   z = Generic.Frac{T}(base_ring(a)(b), one(base_ring(a)))
   z.parent = a
   return z
end


text:
function fill_halos!(f, domain, f0)
    nx, ny, nz = size(domain)
    hx, hy, hz = halo(domain)
    for k = 1-hx:nz+hx
        for j = 1-hy:ny+hy
            for i = 1-hx:0
                x, y, z = domain[i,j,k]
                f[i,j,k] = f0(z)
            end
            for i = nx+1:nx+hx
                x, y, z = domain[i,j,k]
                f[i,j,k] = f0(z)
            end
        end
        for j = 1-hy:0
            for i = 1:nx
                x, y, z = domain[i,j,k]
                f[i,j,k] = f0(z)
            end
        end
        for j = ny+1:ny+hy
            for i = 1:nx
                x, y, z = domain[i,j,k]
                f[i,j,k] = f0(z)
            end
        end
    end
    for k = 1-hz:0
        for j = 1:ny
            for i = 1:nx
                x′, y′, z′ = domain[i,j,1]
                x, y, z = domain[i,j,k]
                f[i,j,k] = f[i,j,1] + f0(z) - f0(z′)
            end
        end
    end
    for k = nz+1:nz+hz
        for j = 1:ny
            for i = 1:nx
                x′, y′, z′ = domain[i,j,nz]
                x, y, z = domain[i,j,k]
                f[i,j,k] = f[i,j,nz] + f0(z) - f0(z′)
            end
        end
    end
    return f
end

function fill_ψ_halos!(ψ, d::Domain, p::Params)
    return fill_halos!(ψ, d, (z) -> background_ψ(z, p))
end

function fill_ϕ_halos!(ϕ, d::Domain, p::Params)
    return fill_halos!(ϕ, d, (z) -> background_ϕ(z, p))
end

function fill_ψ′_halos!(ψ′, d::Domain)
    return fill_halos!(ψ′, d, (z) -> 0.0)
end

function fill_ϕ′_halos!(ϕ′, d::Domain)
    return fill_halos!(ϕ′, d, (z) -> 0.0)
end


text:
################################################################################
# Schemas for terms

# step 1: extract all Term symbols
# step 2: create empty Schema (Dict)
# step 3: for each term, create schema entry based on column from data store

# TODO: handle streaming (Data.RowTable) by iterating over rows and updating
# schemas in place

terms(t::FormulaTerm) = union(terms(t.lhs), terms(t.rhs))
terms(t::InteractionTerm) = terms(t.terms)
terms(t::FunctionTerm{Fo,Fa,names}) where {Fo,Fa,names} = Term.(names)
terms(t::AbstractTerm) = [t]
terms(t::MatrixTerm) = terms(t.terms)
terms(t::TupleTerm) = mapreduce(terms, union, t)

needs_schema(::AbstractTerm) = true
needs_schema(::ConstantTerm) = false
needs_schema(t) = false
# first possible fix for #97
needs_schema(::Union{CategoricalTerm, ContinuousTerm, InterceptTerm}) = false

"""
    StatsModels.Schema

Struct that wraps a `Dict` mapping `Term`s to their concrete forms.  This exists
mainly for dispatch purposes and to support possibly more sophisticated behavior
in the future.

A `Schema` behaves for all intents and purposes like an immutable `Dict`, and
delegates the constructor as well as `getindex`, `get`, `merge!`, `merge`,
`keys`, and `haskey` to the wrapped `Dict`.
"""
struct Schema
    schema::Dict{Term,AbstractTerm}
    Schema(x...) = new(Dict{Term,AbstractTerm}(x...))
end

function Base.show(io::IO, schema::Schema)
    n = length(schema.schema)
    println(io, "StatsModels.Schema with $n ", n==1 ? "entry:" : "entries:")
    for (k,v) in schema.schema
        println(io, "  ", k, " => ", v)
    end
end

Base.getindex(schema::Schema, key) = getindex(schema.schema, key)
Base.get(schema::Schema, key, default) = get(schema.schema, key, default)
Base.merge(a::Schema, b::Schema) = Schema(merge(a.schema, b.schema))
Base.merge!(a::Schema, b::Schema) = (merge!(a.schema, b.schema); a)

Base.keys(schema::Schema) = keys(schema.schema)
Base.haskey(schema::Schema, key) = haskey(schema.schema, key)

"""
    schema([terms::AbstractVector{<:AbstractTerm}, ]data, hints::Dict{Symbol})
    schema(term::AbstractTerm, data, hints::Dict{Symbol})

Compute all the invariants necessary to fit a model with `terms`.  A schema is a dict that
maps `Term`s to their concrete instantiations (either `CategoricalTerm`s or
`ContinuousTerm`s.  "Hints" may optionally be supplied in the form of a `Dict` mapping term
names (as `Symbol`s) to term or contrast types.  If a hint is not provided for a variable,
the appropriate term type will be guessed based on the data type from the data column: any
numeric data is assumed to be continuous, and any non-numeric data is assumed to be
categorical.

Returns a [`StatsModels.Schema`](@ref), which is a wrapper around a `Dict`
mapping `Term`s to their concrete instantiations (`ContinuousTerm` or
`CategoricalTerm`).

# Example

```jldoctest 1
julia> d = (x=sample([:a, :b, :c], 10), y=rand(10));

julia> ts = [Term(:x), Term(:y)];

julia> schema(ts, d)
StatsModels.Schema with 2 entries:
  y => y
  x => x

julia> schema(ts, d, Dict(:x => HelmertCoding()))
StatsModels.Schema with 2 entries:
  y => y
  x => x

julia> schema(term(:y), d, Dict(:y => CategoricalTerm))
StatsModels.Schema with 1 entry:
  y => y
```

Note that concrete `ContinuousTerm` and `CategoricalTerm` and un-typed `Term`s print the
same in a container, but when printed alone are different:

```jldoctest 1
julia> sch = schema(ts, d)
StatsModels.Schema with 2 entries:
  y => y
  x => x

julia> term(:x)
x(unknown)

julia> sch[term(:x)]
x(DummyCoding:3→2)

julia> sch[term(:y)]
y(continuous)
```
"""
schema(data, hints=Dict{Symbol,Any}()) = schema(columntable(data), hints)
schema(dt::D, hints=Dict{Symbol,Any}()) where {D<:ColumnTable} =
    schema(Term.(collect(fieldnames(D))), dt, hints)
schema(ts::AbstractVector{<:AbstractTerm}, data, hints::Dict{Symbol}) =
    schema(ts, columntable(data), hints)

# handle hints:
schema(ts::AbstractVector{<:AbstractTerm}, dt::ColumnTable,
      hints::Dict{Symbol}=Dict{Symbol,Any}()) =
    sch = Schema(t=>concrete_term(t, dt, hints) for t in ts)

schema(f::TermOrTerms, data, hints::Dict{Symbol}) =
    schema(filter(needs_schema, terms(f)), data, hints)

schema(f::TermOrTerms, data) = schema(f, data, Dict{Symbol,Any}())

"""
    concrete_term(t::Term, data[, hint])

Create concrete term from the placeholder `t` based on a data source and
optional hint.  If `data` is a table, the `getproperty` is used to extract the
appropriate column.

The `hint` can be a `Dict{Symbol}` of hints, or a specific hint, a concrete term
type (`ContinuousTerm` or `CategoricalTerm`), or an instance of some
`<:AbstractContrasts`, in which case a `CategoricalTerm` will be created using
those contrasts.

If no hint is provided (or `hint==nothing`), the `eltype` of the data is used:
`Number`s are assumed to be continuous, and all others are assumed to be
categorical.

# Example

```jldoctest
julia> concrete_term(term(:a), [1, 2, 3])
a(continuous)

julia> concrete_term(term(:a), [1, 2, 3], nothing)
a(continuous)

julia> concrete_term(term(:a), [1, 2, 3], CategoricalTerm)
a(DummyCoding:3→2)

julia> concrete_term(term(:a), [1, 2, 3], EffectsCoding())
a(EffectsCoding:3→2)

julia> concrete_term(term(:a), [1, 2, 3], Dict(:a=>EffectsCoding()))
a(EffectsCoding:3→2)

julia> concrete_term(term(:a), (a = [1, 2, 3], b = rand(3)))
a(continuous)
```
"""
concrete_term(t::Term, d, hints::Dict{Symbol}) =
    concrete_term(t, d, get(hints, t.sym, nothing))
concrete_term(t::Term, dt::ColumnTable, hint) =
    concrete_term(t, getproperty(dt, t.sym), hint)
concrete_term(t::Term, dt::ColumnTable, hints::Dict{Symbol}) =
    concrete_term(t, getproperty(dt, t.sym), get(hints, t.sym, nothing))
concrete_term(t::Term, d) = concrete_term(t, d, nothing)

# second possible fix for #97
concrete_term(t, d, hint) = t

concrete_term(t::Term, xs::AbstractVector{<:Number}, ::Nothing) = concrete_term(t, xs, ContinuousTerm)
function concrete_term(t::Term, xs::AbstractVector, ::Type{ContinuousTerm})
    μ, σ2 = StatsBase.mean_and_var(xs)
    min, max = extrema(xs)
    ContinuousTerm(t.sym, promote(μ, σ2, min, max)...)
end
# default contrasts: dummy coding
concrete_term(t::Term, xs::AbstractVector, ::Nothing) = concrete_term(t, xs, CategoricalTerm)
concrete_term(t::Term, xs::AbstractArray, ::Type{CategoricalTerm}) = concrete_term(t, xs, DummyCoding())

function concrete_term(t::Term, xs::AbstractArray, contrasts::AbstractContrasts)
    contrmat = ContrastsMatrix(contrasts, intersect(levels(xs), unique(xs)))
    CategoricalTerm(t.sym, contrmat)
end

"""
    apply_schema(t, schema::StatsModels.Schema[, Mod::Type = Nothing])

Return a new term that is the result of applying `schema` to term `t` with
destination model (type) `Mod`.  If `Mod` is omitted, `Nothing` will be used.

When `t` is a `ContinuousTerm` or `CategoricalTerm` already, the term will be returned
unchanged _unless_ a matching term is found in the schema.  This allows
selective re-setting of a schema to change the contrast coding or levels of a
categorical term, or to change a continuous term to categorical or vice versa.

When defining behavior for custom term types, it's best to dispatch on
[`StatsModels.Schema`](@ref) for the second argument.  Leaving it as `::Any` will work
in _most_ cases, but cause method ambiguity in some.
"""
apply_schema(t, schema) = apply_schema(t, schema, Nothing)
apply_schema(t, schema, Mod::Type) = t
apply_schema(terms::TupleTerm, schema, Mod::Type) = sum(apply_schema.(terms, Ref(schema), Mod))

apply_schema(t::Term, schema::Schema, Mod::Type) = schema[t]
apply_schema(ft::FormulaTerm, schema::Schema, Mod::Type) =
    FormulaTerm(apply_schema(ft.lhs, schema, Mod),
                collect_matrix_terms(apply_schema(ft.rhs, schema, Mod)))
apply_schema(it::InteractionTerm, schema::Schema, Mod::Type) =
    InteractionTerm(apply_schema(it.terms, schema, Mod))

# for re-setting schema (in setcontrasts!)
apply_schema(t::Union{ContinuousTerm, CategoricalTerm}, schema::Schema, Mod::Type) =
    get(schema, term(t.sym), t)
apply_schema(t::MatrixTerm, sch::Schema, Mod::Type) =
    MatrixTerm(apply_schema.(t.terms, Ref(sch), Mod))

# TODO: special case this for <:RegressionModel ?
function apply_schema(t::ConstantTerm, schema::Schema, Mod::Type)
    t.n ∈ (-1, 0, 1) ||
        throw(ArgumentError("can't create InterceptTerm from $(t.n) " *
                            "(only -1, 0, and 1 allowed)"))
    InterceptTerm{t.n==1}()
end

"""
    has_schema(t::T) where {T<:AbstractTerm}

Return `true` if `t` has a schema, meaning that `apply_schema` would be a no-op.
"""
has_schema(t::AbstractTerm) = true
has_schema(t::ConstantTerm) = false
has_schema(t::Term) = false
has_schema(t::Union{ContinuousTerm,CategoricalTerm}) = true
has_schema(t::InteractionTerm) = all(has_schema(tt) for tt in t.terms)
has_schema(t::TupleTerm) = all(has_schema(tt) for tt in t)
has_schema(t::FormulaTerm) = has_schema(t.lhs) && has_schema(t.rhs)

struct FullRank
    schema::Schema
    already::Set{AbstractTerm}
end

FullRank(schema) = FullRank(schema, Set{AbstractTerm}())

Base.get(schema::FullRank, key, default) = get(schema.schema, key, default)
Base.merge(a::FullRank, b::FullRank) = FullRank(merge(a.schema, b.schema),
                                                union(a.already, b.already))

function apply_schema(t::FormulaTerm, schema::Schema, Mod::Type{<:StatisticalModel})
    schema = FullRank(schema)

    # Models with the drop_intercept trait do not support intercept terms,
    # usually because they include one implicitly.
    if drop_intercept(Mod)
        if hasintercept(t)
            throw(ArgumentError("Model type $Mod doesn't support intercept " *
                                "specified in formula $t"))
        end
        # start parsing as if we already had the intercept
        push!(schema.already, InterceptTerm{true}())
    elseif implicit_intercept(Mod) && !hasintercept(t) && !omitsintercept(t)
        t = FormulaTerm(t.lhs, InterceptTerm{true}() + t.rhs)
    end

    # only apply rank-promoting logic to RIGHT hand side
    FormulaTerm(apply_schema(t.lhs, schema.schema, Mod),
                collect_matrix_terms(apply_schema(t.rhs, schema, Mod)))
end

# strategy is: apply schema, then "repair" if necessary (promote to full rank
# contrasts).
#
# to know whether to repair, need to know context a term appears in.  main
# effects occur in "own" context.

"""
    apply_schema(t::AbstractTerm, schema::StatsModels.FullRank, Mod::Type)

Apply a schema, under the assumption that when a less-than-full rank model
matrix would be produced, categorical terms should be "promoted" to full rank
(where a categorical variable with ``k`` levels would produce ``k`` columns,
instead of ``k-1`` in the standard contrast coding schemes).  This step is
applied automatically when `Mod <: StatisticalModel`, but other types of models
can opt-in by adding a method like

```
StatsModels.apply_schema(t::FormulaTerm, schema::StatsModels.Schema, Mod::Type{<:MyModelType}) =
    apply_schema(t, StatsModels.FullRank(schema), mod)
```

See the section on [Modeling categorical data](@ref) in the docs for more
information on how promotion of categorical variables works.
"""
function apply_schema(t::ConstantTerm, schema::FullRank, Mod::Type)
    push!(schema.already, t)
    apply_schema(t, schema.schema, Mod)
end

apply_schema(t::InterceptTerm, schema::FullRank, Mod::Type) = (push!(schema.already, t); t)

# TODO: maybe change this to t::Any of ::AbstractTerm to catch
# categorical/continuous terms?
function apply_schema(t::AbstractTerm, schema::FullRank, Mod::Type)
    push!(schema.already, t)
    t = apply_schema(t, schema.schema, Mod) # continuous or categorical now
    apply_schema(t, schema, Mod, t) # repair if necessary
end

function apply_schema(t::InteractionTerm, schema::FullRank, Mod::Type)
    push!(schema.already, t)
    terms = apply_schema.(t.terms, Ref(schema.schema), Mod)
    terms = apply_schema.(terms, Ref(schema), Mod, Ref(t))
    InteractionTerm(terms)
end

# context doesn't matter for non-categorical terms
apply_schema(t, schema::FullRank, Mod::Type, context::AbstractTerm) = t
# when there's a context, check to see if any of the terms already seen would be
# aliased by this term _if_ it were full rank.
function apply_schema(t::CategoricalTerm, schema::FullRank, Mod::Type, context::AbstractTerm)
    aliased = drop_term(context, t)
    @debug "$t in context of $context: aliases $aliased\n  seen already: $(schema.already)"
    for seen in schema.already
        if symequal(aliased, seen)
            @debug "  aliased term already present: $seen"
            return t
        end
    end
    # aliased term not seen already:
    # add aliased term to already seen:
    push!(schema.already, aliased)
    # repair:
    new_contrasts = ContrastsMatrix(FullDummyCoding(), t.contrasts.levels)
    t = CategoricalTerm(t.sym, new_contrasts)
    @debug "  aliased term absent, repairing: $t"
    t
end

drop_term(from, to) = symequal(from, to) ? ConstantTerm(1) : from
drop_term(from::FormulaTerm, to) = FormulaTerm(from.lhs, drop_term(from.rhs, to))
drop_term(from::MatrixTerm, to) = MatrixTerm(drop_term(from.terms, to))
drop_term(from::TupleTerm, to) =
    tuple((t for t = from if !symequal(t, to))...)
function drop_term(from::InteractionTerm, t)
    terms = drop_term(from.terms, t)
    length(terms) > 1 ? InteractionTerm(terms) : terms[1]
end

"""
    termsyms(t::Terms.Term)

Extract the set of symbols referenced in this term.

This is needed in order to determine when a categorical term should have
standard (reduced rank) or full rank contrasts, based on the context it occurs
in and the other terms that have already been encountered.
"""
termsyms(t::AbstractTerm) = Set()
termsyms(t::InterceptTerm{true}) = Set(1)
termsyms(t::ConstantTerm) = Set((t.n,))
termsyms(t::Union{Term, CategoricalTerm, ContinuousTerm}) = Set([t.sym])
termsyms(t::InteractionTerm) = mapreduce(termsyms, union, t.terms)
termsyms(t::FunctionTerm) = Set([t.exorig])

symequal(t1::AbstractTerm, t2::AbstractTerm) = issetequal(termsyms(t1), termsyms(t2))

"""
    termvars(t::AbstractTerm)

The data variables that this term refers to.
"""
termvars(::AbstractTerm) = Symbol[]
termvars(t::Union{Term, CategoricalTerm, ContinuousTerm}) = [t.sym]
termvars(t::InteractionTerm) = mapreduce(termvars, union, t.terms)
termvars(t::TupleTerm) = mapreduce(termvars, union, t, init=Symbol[])
termvars(t::MatrixTerm) = termvars(t.terms)
termvars(t::FormulaTerm) = union(termvars(t.lhs), termvars(t.rhs))
termvars(t::FunctionTerm{Fo,Fa,names}) where {Fo,Fa,names} = collect(names)


text:
#--------------------------------------------------------------------
# DNSS.jl
# Soham 03-2022
# Simulate the collapse of a self-gravitating scalar field
# in Vaidya spacetime
#--------------------------------------------------------------------

using ForwardDiff, QuadGK

function psi(x::T)::T where {T<:Real}
    (a, b) = (-0.5, 0.5)
    # (a, b) = (-1, 1)
    if (a <= x <= b)
        return x * (1-(x/a))^8  * (1 - (x/b))^8
    else
        return 0.0
    end
end

function psi_spherical(u::T, v::T)::T where {T<:Number}
    t = v + u
    r = v - u
    α = π/2 - k*t
    if  r == 0.0
        return cos(k*t)
    else
        return cos(k*t) * (sin(k*r) / (k*r))
    end
end

function checkinitialdata(a::Field{S}) where {S}
    PS = derivative(a.space)
    dψ = Field(PS, dpsi)
    r  = Field(PS, grr)
    D  = derivative(PS)
    @show norm(-(2 / a) * (D * a) + 4 * π * r * (dψ)^2)
end

function dpsi(x::T)::T where {T<:Real}
    return ForwardDiff.derivative(psi, x)
end

function guv(v::T)::T where {T<:Real}
    integral, error = quadgk(x->(x * dpsi(x)^2), -1, v, rtol=1e-12)
    return exp(2π * integral)
end

function resolve(U::NTuple{3, Field{S}}, var::Symbol)::NTuple{3, Field{S}} where {S}
    PS = first(U).space
    @show PS
    D = derivative(PS)
    I = identity(PS)
    B = incomingboundary(PS)

    if var == :a
        (a, η, ψ) = U
        L = (D*η)*D + (4π * η * (D*ψ)^2) * I
        lna = linsolve(L + B, B*log(a))
        a   = exp(lna)
        return (a, η, ψ)
    end
end

# TODO: Test the simplest initial data for spherically symmetric wave
A  = ChebyshevGL{U, 20, Float64}(0.0, 1.0)
B  = ChebyshevGL{V, 20, Float64}(0.0, 1.0)
PS = ProductSpace(A, B)

k  = 0.01
ψ  = extractUboundary(Field(PS, psi_spherical), :incoming)
r  = extractUboundary(Field(PS, (u,v)->v-u), :incoming)
a  = extractUboundary(Field(PS, (u,v)->2), :incoming)
@show norm(lineconstraints((a, r, ψ)))


text:
using TargetedLearning.TMLEs
import TargetedLearning.TMLEs.TMLE

facts("Test TMLEs") do
    n,p = 100, 10
    w = [ones(n) rand(n, p)]
    a = round(rand(n))
    y = round(rand(n))

    qfit = lreg([w a], y)
    logitQnA1 = linpred(qfit, [w ones(n)])
    logitQnA0 = linpred(qfit, [w zeros(n)])
    gn1 = predict(lreg(w,a),w)

    est=tmle(logitQnA1, logitQnA0, gn1, a, y, param=ATE(1,0))

    @fact est --> is_a(TMLE)

    hA1 = 1./gn1
    hA0 = -1./(1 - gn1)
    hAA = ifelse(a.==1, hA1, hA0)
    fluc = lreg(hAA, y, offset = ifelse(a.==1, logitQnA1, logitQnA0))
    QstarA1 = predict(fluc, hA1, offset=logitQnA1)
    QstarA0 = predict(fluc, hA0, offset=logitQnA0)
    psi = mean(QstarA1 .- QstarA0)
    ic = hAA .* (y .- ifelse(a.==1, QstarA1, QstarA0)) .+ QstarA1 .- QstarA0 .- psi

    @fact est.psi --> roughly(psi)
    @fact est.ic --> roughly(ic)

    context("With weighted fluctuation") do
        @fact tmle(logitQnA1, logitQnA0, gn1, a, y, weightedfluc=true) --> is_a(TMLE)
        @pending "more weighted fluc tests" --> nothing
    end
end


text:
struct PerfectElectricConductor <: AbstractBoundaryCondition end
const PEC = PerfectElectricConductor()


text:
using Millboard
using Base.Test

include(joinpath(dirname(@__FILE__),"numbers.jl"))
include(joinpath(dirname(@__FILE__),"strings.jl"))
include(joinpath(dirname(@__FILE__),"dicts.jl"))
include(joinpath(dirname(@__FILE__),"enums.jl"))
include(joinpath(dirname(@__FILE__),"functions.jl"))
include(joinpath(dirname(@__FILE__),"tuples.jl"))
include(joinpath(dirname(@__FILE__),"types.jl"))
include(joinpath(dirname(@__FILE__),"unicode.jl"))
include(joinpath(dirname(@__FILE__),"regex.jl"))

if Base.have_color
  include(joinpath(dirname(@__FILE__),"colored.jl"))
end


text:
using Printf
using Oceananigans.Grids: size_summary, scalar_summary

location_str(::Type{Face})    = "Face"
location_str(::Type{Center})  = "Center"
location_str(::Type{Nothing}) = "⋅"
show_location(LX, LY, LZ) = "($(location_str(LX)), $(location_str(LY)), $(location_str(LZ)))"
show_location(field::AbstractField) = show_location(location(field)...)

function Base.summary(field::Field)
    LX, LY, LZ = location(field)
    prefix = string(size_summary(size(field)), " Field{$LX, $LY, $LZ}")

    grid_name = typeof(field.grid).name.wrapper
    reduced_dims = reduced_dimensions(field)

    suffix = reduced_dims === () ?
        string(" on ", grid_name, " on ", summary(architecture(field))) :
        string(" reduced over dims = ", reduced_dims,
               " on ", grid_name, " on ", summary(architecture(field)))

    return string(prefix, suffix)
end

data_summary(field) = string("max=", scalar_summary(maximum(field)), ", ",
                             "min=", scalar_summary(minimum(field)), ", ",
                             "mean=", scalar_summary(mean(field)))

function Base.show(io::IO, field::Field)

    prefix =
        string("$(summary(field))\n",
               "├── grid: ", summary(field.grid), '\n',
               "├── boundary conditions: ", summary(field.boundary_conditions), '\n')

    middle = isnothing(field.operand) ? "" :
        string("├── operand: ", summary(field.operand), '\n',
               "├── status: ", summary(field.status), '\n')

    suffix = string("└── data: ", summary(field.data), '\n',
                    "    └── ", data_summary(field))

    print(io, prefix, middle, suffix)
end

Base.summary(status::FieldStatus) = "time=$(status.time)"

Base.summary(::ZeroField{N}) where N = "ZeroField{$N}"
Base.show(io::IO, z::ZeroField) = print(io, summary(z))

Base.show(io::IO, ::MIME"text/plain", f::AbstractField) = show(io, f)