text:
Random.seed!(1351)
n, p = 100, 5
X  = randn(n, p)
X1 = hcat(X, ones(100))
θ  = randn(p+1)
y  = X1 * θ .+ 0.1 .* randn(n)

θ_ls = X1 \ y

@testset "AdaBoostReg" begin
    m = AdaBoostRegressor(random_state=0, n_estimators=100)
    f, _, _ = fit(m, 1, X, y)
    @test 0.9 ≤ f[1].score(X, y) ≤ 0.999
    @test 0.2 ≤ norm(predict(m, f, X) .- y)/norm(y) ≤ 0.25
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:estimators, :estimator_weights, :estimator_errors, :feature_importances_)
end

@testset "BaggingReg" begin
    m = BaggingRegressor(random_state=0)
    f, _, _ = fit(m, 1, X, y)
    @test 0.9 ≤ f[1].score(X, y) ≤ 0.999
    @test 0.15 ≤ norm(predict(m, f, X) .- y)/norm(y) ≤ 0.25
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:estimators, :estimators_samples, :estimators_features, :oob_score, :oob_prediction)
end

@testset "XTreeReg" begin
    m = ExtraTreesRegressor(random_state=0)
    f, _, _ = fit(m, 1, X, y)
    @test 0.9 ≤ f[1].score(X, y) ≤ 0.999
    @test 0.15 ≤ norm(predict(m, f, X) .- y)/norm(y) ≤ 0.25
    fp = fitted_params(m, f)
    @test keys(fp) == (:estimators, :feature_importances, :n_features, :n_outputs, :oob_score, :oob_prediction)
end

@testset "GBReg" begin
    m = GradientBoostingRegressor(random_state=0)
    f, _, _ = fit(m, 1, X, y)
    @test 0.9 ≤ f[1].score(X, y) ≤ 0.999
    @test 0.03 ≤ norm(predict(m, f, X) .- y)/norm(y) ≤ 0.05
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:feature_importances, :train_score, :loss, :init, :estimators, :oob_improvement)
end

@testset "RFReg" begin
    m = RandomForestRegressor(random_state=0, oob_score=true)
    f, _, _ = fit(m, 1, X, y)
    @test 0.9 ≤ f[1].score(X, y) ≤ 0.999
    @test 0.15 ≤ norm(predict(m, f, X) .- y)/norm(y) ≤ 0.25
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:estimators, :feature_importances, :n_features, :n_outputs, :oob_score, :oob_prediction)
end

# =================
# classifiers
# =================

Xc2, yc2 = gen_classif(classes=["A", "B", "C"])

@testset "AdaboostClf" begin
    m, f = simple_test_classif_prob(AdaBoostClassifier(), Xc2, yc2)
    fp = fitted_params(m, f)
    @test keys(fp) == (:estimators, :estimator_weights, :estimator_errors, :classes, :n_classes)
    infos = info_dict(m)
    @test infos[:input_scitype] == MLJBase.Table(MLJBase.Continuous)
    @test infos[:target_scitype] == AbstractVector{<:MLJBase.Finite}
    @test !isempty(infos[:docstring])
end

@testset "BaggingClf" begin
    m, f = simple_test_classif_prob(BaggingClassifier(), Xc2, yc2)
    fp = fitted_params(m, f)
    @test keys(fp) == (:base_estimator, :estimators, :estimators_samples, :estimators_features, :classes, :n_classes, :oob_score, :oob_decision_function)
    infos = info_dict(m)
    @test infos[:input_scitype] == MLJBase.Table(MLJBase.Continuous)
    @test infos[:target_scitype] == AbstractVector{<:MLJBase.Finite}
    @test !isempty(infos[:docstring])
end

@testset "GradBoostClf" begin
    m, f = simple_test_classif_prob(GradientBoostingClassifier(), Xc2, yc2)
    fp = fitted_params(m, f)
    @test keys(fp) == (:n_estimators, :feature_importances, :train_score, :loss, :init, :estimators, :oob_improvement)
    infos = info_dict(m)
    @test infos[:input_scitype] == MLJBase.Table(MLJBase.Continuous)
    @test infos[:target_scitype] == AbstractVector{<:MLJBase.Finite}
    @test !isempty(infos[:docstring])
end

@testset "RForestClf" begin
    m, f = simple_test_classif_prob(RandomForestClassifier(), Xc2, yc2)
    fp = fitted_params(m, f)
    @test keys(fp) == (:estimators, :classes, :n_classes, :n_features, :n_outputs, :feature_importances, :oob_score, :oob_decision_function)
    infos = info_dict(m)
    @test infos[:input_scitype] == MLJBase.Table(MLJBase.Continuous)
    @test infos[:target_scitype] == AbstractVector{<:MLJBase.Finite}
    @test !isempty(infos[:docstring])
end

@testset "XTreeClf" begin
    m, f = simple_test_classif_prob(ExtraTreesClassifier(), Xc2, yc2)
    fp = fitted_params(m, f)
    @test keys(fp) == (:estimators, :classes, :n_classes, :feature_importances, :n_features, :n_outputs, :oob_score, :oob_decision_function)
    infos = info_dict(m)
    @test infos[:input_scitype] == MLJBase.Table(MLJBase.Continuous)
    @test infos[:target_scitype] == AbstractVector{<:MLJBase.Finite}
    @test !isempty(infos[:docstring])
end


text:
"""
    test_bqg_rossbywave(; kwargs...)

Evolves a Rossby wave and compares with the analytic solution.
"""
function test_bqg_rossbywave(stepper, dt, nsteps)
    nx = 64
  beta = 2.0
    Lx = 2π
    mu = 0.0
    nu = 0.0

  gr  = TwoDGrid(nx, Lx)
  x, y = gridpoints(gr)

  # the following if statement is called so that all the cases of
  # Problem() fuction are tested
  if stepper=="ForwardEuler"
    eta = zeros(nx, nx)
  else
    eta(x, y) = 0*x
  end

  prob = BarotropicQG.InitialValueProblem(nx=nx, Lx=Lx, eta=eta, beta=beta, mu=mu, nu=nu, stepper=stepper, dt=dt)
  sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid

  x, y = gridpoints(g)

  # the Rossby wave initial condition
   ampl = 1e-2
  kwave = 3.0*2π/g.Lx
  lwave = 2.0*2π/g.Ly
      ω = -p.beta*kwave/(kwave^2.0 + lwave^2.0)
     ζ0 = @. ampl*cos(kwave*x)*cos(lwave*y)
    ζ0h = rfft(ζ0)

  BarotropicQG.set_zeta!(prob, ζ0)

  stepforward!(prob, nsteps)
  dealias!(sol, g)
  BarotropicQG.updatevars!(prob)

  ζ_theory = @. ampl*cos(kwave*(x - ω/kwave*cl.t)) * cos(lwave*y)

  isapprox(ζ_theory, v.zeta, rtol=g.nx*g.ny*nsteps*1e-12)
end

"""
    test_stochasticforcingbudgets(; kwargs...)

Tests if the energy budgets are closed for BarotropicQG with stochastic forcing.
"""
function test_bqg_stochasticforcingbudgets(; n=256, dt=0.01, L=2π, nu=1e-7, nnu=2, mu=1e-1, message=false)
  n, L  = 256, 2π
  nu, nnu = 1e-7, 2
  mu = 1e-1
  dt, tf = 0.005, 0.1/mu
  nt = round(Int, tf/dt)
  ns = 1

  # Forcing
  kf, dkf = 12.0, 2.0
  ε = 0.1

  gr  = TwoDGrid(n, L)
  x, y = gridpoints(gr)

  Kr = [ gr.kr[i] for i=1:gr.nkr, j=1:gr.nl]

  force2k = zero(gr.Krsq)
  @. force2k = exp.(-(sqrt(gr.Krsq)-kf)^2/(2*dkf^2))
  @. force2k[gr.Krsq .< 2.0^2 ] = 0
  @. force2k[gr.Krsq .> 20.0^2 ] = 0
  @. force2k[Kr .< 2π/L] = 0
  ε0 = parsevalsum(force2k.*gr.invKrsq/2.0, gr)/(gr.Lx*gr.Ly)
  force2k .= ε/ε0 * force2k

  Random.seed!(1234)

  function calcFq!(Fqh, sol, t, cl, v, p, g)
    eta = exp.(2π*im*rand(Float64, size(sol)))/sqrt(cl.dt)
    eta[1, 1] = 0
    @. Fqh = eta * sqrt(force2k)
    nothing
  end

  prob = BarotropicQG.ForcedProblem(nx=n, Lx=L, nu=nu, nnu=nnu, mu=mu, dt=dt,
   stepper="RK4", calcFq=calcFq!, stochastic=true)

  sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid

  BarotropicQG.set_zeta!(prob, 0*x)
  E = Diagnostic(BarotropicQG.energy,      prob, nsteps=nt)
  D = Diagnostic(BarotropicQG.dissipation, prob, nsteps=nt)
  R = Diagnostic(BarotropicQG.drag,        prob, nsteps=nt)
  W = Diagnostic(BarotropicQG.work,        prob, nsteps=nt)
  diags = [E, D, W, R]

  # Step forward

  stepforward!(prob, diags, round(Int, nt))

  BarotropicQG.updatevars!(prob)

  cfl = cl.dt*maximum([maximum(v.v)/g.dx, maximum(v.u)/g.dy])

  E, D, W, R = diags

  t = round(mu*cl.t, digits=2)

  i₀ = 1
  dEdt = (E[(i₀+1):E.i] - E[i₀:E.i-1])/cl.dt
  ii = (i₀):E.i-1
  ii2 = (i₀+1):E.i

  # dEdt = W - D - R?
  # If the Ito interpretation was used for the work
  # then we need to add the drift term
  # total = W[ii2]+ε - D[ii] - R[ii]      # Ito
  total = W[ii2] - D[ii] - R[ii]        # Stratonovich

  residual = dEdt - total

  if message
    println("step: %04d, t: %.1f, cfl: %.3f, time: %.2f s\n", cl.step, cl.t, cfl, t)
  end
  # println(mean(abs.(residual)))
  isapprox(mean(abs.(residual)), 0, atol=1e-4)
end

"""
    test_stochasticforcingbudgets(; kwargs...)

Tests if the energy budgets are closed for BarotropicQG with stochastic forcing.
"""
function test_bqg_deterministicforcingbudgets(; n=256, dt=0.01, L=2π, nu=1e-7, nnu=2, mu=1e-1, message=false)
  n, L  = 256, 2π
  nu, nnu = 1e-7, 2
  mu = 1e-1
  dt, tf = 0.005, 0.1/mu
  nt = round(Int, tf/dt)
  ns = 1

  # Forcing = 0.01cos(4x)cos(5y)cos(2t)
  gr  = TwoDGrid(n, L)
  x, y = gridpoints(gr)
  k0, l0 = gr.kr[2], gr.l[2]

  f = @. 0.01*cos(4*k0*x)*cos(5*l0*y)
  fh = rfft(f)
  function calcFq!(Fqh, sol, t, cl, v, p, g)
    @. Fqh = fh*cos(2*t)
    nothing
  end

  prob = BarotropicQG.ForcedProblem(nx=n, Lx=L, nu=nu, nnu=nnu, mu=mu, dt=dt,
   stepper="RK4", calcFq=calcFq!, stochastic=false)

  sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid

  BarotropicQG.set_zeta!(prob, 0*x)
  E = Diagnostic(BarotropicQG.energy,      prob, nsteps=nt)
  D = Diagnostic(BarotropicQG.dissipation, prob, nsteps=nt)
  R = Diagnostic(BarotropicQG.drag,        prob, nsteps=nt)
  W = Diagnostic(BarotropicQG.work,        prob, nsteps=nt)
  diags = [E, D, W, R]

  # Step forward

  stepforward!(prob, diags, round(Int, nt))

  BarotropicQG.updatevars!(prob)

  cfl = cl.dt*maximum([maximum(v.v)/g.dx, maximum(v.u)/g.dy])

  E, D, W, R = diags

  t = round(mu*cl.t, digits=2)

  i₀ = 1
  dEdt = (E[(i₀+1):E.i] - E[i₀:E.i-1])/cl.dt
  ii = (i₀):E.i-1
  ii2 = (i₀+1):E.i

  # dEdt = W - D - R?
  total = W[ii2] - D[ii] - R[ii]

  residual = dEdt - total

  if message
    println("step: %04d, t: %.1f, cfl: %.3f, time: %.2f s\n", prob.step, cl.t, cfl, tc)
  end
  # println(mean(abs.(residual)))
  isapprox(mean(abs.(residual)), 0, atol=1e-8)
end

"""
    test_bqg_nonlinearadvection(dt, stepper; kwargs...)

Tests the advection term in the twodturb module by timestepping a
test problem with timestep dt and timestepper identified by the string stepper.
The test problem is derived by picking a solution ζf (with associated
streamfunction ψf) for which the advection term J(ψf, ζf) is non-zero. Next, a
forcing Ff is derived according to Ff = ∂ζf/∂t + J(ψf, ζf) - nuΔζf. One solution
to the vorticity equation forced by this Ff is then ζf. (This solution may not
be realized, at least at long times, if it is unstable.)
"""
function test_bqg_advection(dt, stepper; n=128, L=2π, nu=1e-2, nnu=1, mu=0.0, message=false)
  n, L  = 128, 2π
  nu, nnu = 1e-2, 1
  mu = 0.0
  tf = 1.0
  nt = round(Int, tf/dt)

  gr  = TwoDGrid(n, L)
  x, y = gridpoints(gr)

  psif = @. sin(2x)*cos(2y) + 2sin(x)*cos(3y)
  qf = @. -8sin(2x)*cos(2y) - 20sin(x)*cos(3y)

  Ff = @. -(
    nu*( 64sin(2x)*cos(2y) + 200sin(x)*cos(3y) )
    + 8*( cos(x)*cos(3y)*sin(2x)*sin(2y) - 3cos(2x)*cos(2y)*sin(x)*sin(3y) )
  )

  Ffh = rfft(Ff)

  # Forcing
  function calcFq!(Fqh, sol, t, cl, v, p, g)
    Fqh .= Ffh
    nothing
  end

  prob = BarotropicQG.ForcedProblem(nx=n, Lx=L, nu=nu, nnu=nnu, mu=mu, dt=dt, stepper=stepper, calcFq=calcFq!)
  sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
  BarotropicQG.set_zeta!(prob, qf)

  # Step forward
  stepforward!(prob, round(Int, nt))
  BarotropicQG.updatevars!(prob)
  isapprox(v.q, qf, rtol=1e-13)
end

"""
    test_bqg_formstress(dt, stepper; kwargs...)

Tests the form stress term that forces the domain-averaged zonal flow U(t).
"""
function test_bqg_formstress(dt, stepper; n=128, L=2π, nu=0.0, nnu=1, mu=0.0, message=false)
  n, L  = 128, 2π
  nu, nnu = 1e-2, 1
  mu = 0.0
  tf = 1
  nt = 1

  gr  = TwoDGrid(n, L)
  x, y = gridpoints(gr)

  zetai = @. -20*sin(10*x)*cos(10*y)
  topoPV(x, y) = @. cos(10x)*cos(10y)
  F(t) = 0 #no forcing

  answer = 0.25 # this is what <v*eta> should be

  prob = BarotropicQG.ForcedProblem(nx=n, Lx=L, nu=nu, nnu=nnu, mu=mu, dt=dt, stepper=stepper, eta=topoPV, calcFU = F)
  BarotropicQG.set_zeta!(prob, zetai)
  BarotropicQG.updatevars!(prob)

  # Step forward
  stepforward!(prob, nt)
  isapprox(prob.timestepper.N[1, 1], answer, rtol=1e-13)
end

function test_bqg_energyenstrophy()
  nx, Lx  = 64, 2π
  ny, Ly  = 64, 3π
  g  = TwoDGrid(nx, Lx, ny, Ly)
  k0, l0 = g.k[2], g.l[2] # fundamental wavenumbers
  x, y = gridpoints(g)

  energy_calc = 29/9
  enstrophy_calc = 2701/162

    eta = @. cos(10k0*x)*cos(10l0*y)
   psi0 = @. sin(2k0*x)*cos(2l0*y) + 2sin(k0*x)*cos(3l0*y)
  zeta0 = @. -((2k0)^2+(2l0)^2)*sin(2k0*x)*cos(2l0*y) - (k0^2+(3l0)^2)*2sin(k0*x)*cos(3l0*y)

  prob = BarotropicQG.InitialValueProblem(nx=nx, Lx=Lx, ny=ny, Ly=Ly, eta = eta, stepper="ForwardEuler")
  sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
  BarotropicQG.set_zeta!(prob, zeta0)
  BarotropicQG.updatevars!(prob)

  energyzeta0 = BarotropicQG.energy(prob)
  enstrophyzeta0 = BarotropicQG.enstrophy(prob)

  isapprox(energyzeta0, energy_calc, rtol=1e-13) && isapprox(enstrophyzeta0, enstrophy_calc, rtol=1e-13) &&
  BarotropicQG.addforcing!(prob.timestepper.N, sol, cl.t, cl, v, p, g)==nothing
end

function test_bqg_meanenergyenstrophy()
  nx, Lx  = 64, 2π
  ny, Ly  = 96, 3π
  g  = TwoDGrid(nx, Lx, ny, Ly)
  k0, l0 = g.k[2], g.l[2] # fundamental wavenumbers
  x, y = gridpoints(g)

  calcFU(t) = 0.0
  eta(x, y) = @. cos(10x)*cos(10y)
  psi0 = @. sin(2k0*x)*cos(2l0*y) + 2sin(k0*x)*cos(3l0*y)
 zeta0 = @. -((2k0)^2+(2l0)^2)*sin(2k0*x)*cos(2l0*y) - (k0^2+(3l0)^2)*2sin(k0*x)*cos(3l0*y)
  beta = 10.0
  U = 1.2

  energy_calc = 29/9
  enstrophy_calc = 2701/162

  prob = BarotropicQG.ForcedProblem(nx=nx, Lx=Lx, ny=ny, Ly=Ly, beta=beta, eta=eta, calcFU = calcFU,
                                    stepper="ForwardEuler")

  BarotropicQG.set_zeta!(prob, zeta0)
  BarotropicQG.set_U!(prob, U)
  BarotropicQG.updatevars!(prob)

  energyU = BarotropicQG.meanenergy(prob)
  enstrophyU = BarotropicQG.meanenstrophy(prob)

  energyzeta0 = BarotropicQG.energy(prob)
  enstrophyzeta0 = BarotropicQG.enstrophy(prob)

  (isapprox(energyU, 0.5*U^2, rtol=1e-13) &&
    isapprox(enstrophyU, beta*U, rtol=1e-13) &&
    isapprox(energyzeta0, energy_calc, rtol=1e-13) &&
    isapprox(enstrophyzeta0, enstrophy_calc, rtol=1e-13)
  )
end

function test_bqg_problemtype(T=Float32)
  prob = BarotropicQG.Problem(T=T)

  (typeof(prob.sol)==Array{Complex{T},2} && typeof(prob.grid.Lx)==T && typeof(prob.grid.x)==Array{T,2} && typeof(prob.vars.u)==Array{T,2})
end

ext: jl


text:
using BinaryBuilder

include("../common.jl")

# Collection of sources required to build OpenBLAS
name = "OpenBLAS"
version = v"0.3.10"

sources = openblas_sources(version)
script = openblas_script()
platforms = openblas_platforms(;experimental=true)
products = openblas_products()
dependencies = openblas_dependencies()

# Build the tarballs
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"6", lock_microarchitecture=false)


text:
## --- Weighted mean of an array

    """
    ```julia
    (wx, wσ, mswd) = awmean(x, σ)
    ```
    Weighted mean, absent the geochonologist's MSWD correction to uncertainty.
    """
    function awmean(x, σ)
        n = length(x)

        sum_of_values = sum_of_weights = χ2 = 0.0
        @inbounds @simd for i=1:n
            sum_of_values += x[i] / (σ[i]*σ[i])
            sum_of_weights += 1 / (σ[i]*σ[i])
        end
        wx = sum_of_values / sum_of_weights

        @inbounds @simd for i=1:n
            χ2 += (x[i] - wx) * (x[i] - wx) / (σ[i] * σ[i])
        end
        mswd = χ2 / (n-1)
        wσ = sqrt(1.0 / sum_of_weights)
        return wx, wσ, mswd
    end
    # function awmean(x::Array{<:Number}, σ::Array{<:Number})
    #     n = length(x)
    #
    #     sum_of_values = sum_of_weights = χ2 = 0.0
    #     @turbo for i=1:n
    #         sum_of_values += x[i] / (σ[i]*σ[i])
    #         sum_of_weights += 1 / (σ[i]*σ[i])
    #     end
    #     wx = sum_of_values / sum_of_weights
    #
    #     @turbo for i=1:n
    #         χ2 += (x[i] - wx) * (x[i] - wx) / (σ[i] * σ[i])
    #     end
    #     mswd = χ2 / (n-1)
    #     wσ = sqrt(1.0 / sum_of_weights)
    #     return wx, wσ, mswd
    # end
    export awmean

    """
    ```julia
    (wx, wσ, mswd) = gwmean(x, σ)
    ```
    Geochronologist's weighted mean, with "MSWD correction" to uncertainty,
    i.e., wσ is increased by a factor of sqrt(mswd)
    """
    function gwmean(x, σ)
        n = length(x)

        sum_of_values = sum_of_weights = χ2 = 0.0
        @inbounds @simd for i=1:n
            sum_of_values += x[i] / (σ[i]*σ[i])
            sum_of_weights += 1 / (σ[i]*σ[i])
        end
        wx = sum_of_values / sum_of_weights

        @inbounds @simd for i=1:n
            χ2 += (x[i] - wx) * (x[i] - wx) / (σ[i] * σ[i])
        end
        mswd = χ2 / (n-1)
        wσ = sqrt(mswd / sum_of_weights)
        return wx, wσ, mswd
    end
    # function gwmean(x::Array{<:Number}, σ::Array{<:Number})
    #     n = length(x)
    #     sum_of_values = sum_of_weights = χ2 = 0.0
    #     @turbo for i=1:n
    #         sum_of_values += x[i] / (σ[i]*σ[i])
    #         sum_of_weights += 1 / (σ[i]*σ[i])
    #     end
    #     wx = sum_of_values / sum_of_weights
    #
    #     @turbo for i=1:n
    #         χ2 += (x[i] - wx) * (x[i] - wx) / (σ[i] * σ[i])
    #     end
    #     mswd = χ2 / (n-1)
    #     wσ = sqrt(mswd / sum_of_weights)
    #     return wx, wσ, mswd
    # end
    export gwmean

    """
    ```julia
    MSWD(x, σ)
    ```
    Return the Mean Square of Weighted Deviates (AKA the reduced chi-squared
    statistic) of a dataset with values `x` and one-sigma uncertainties `σ`
    """
    function MSWD(x, σ)
        sum_of_values = sum_of_weights = χ2 = 0.0
        n = length(x)

        @inbounds @simd for i=1:n
            w = 1 / (σ[i]*σ[i])
            sum_of_values += w * x[i]
            sum_of_weights += w
        end
        wx = sum_of_values / sum_of_weights

        @inbounds @simd for i=1:n
            χ2 += (x[i] - wx) * (x[i] - wx) / (σ[i] * σ[i])
        end

        return χ2 / (n-1)
    end
    # function MSWD(x::Array{<:Number}, σ::Array{<:Number})
    #     sum_of_values = sum_of_weights = χ2 = 0.0
    #     n = length(x)
    #
    #     @turbo for i=1:n
    #         w = 1 / (σ[i]*σ[i])
    #         sum_of_values += w * x[i]
    #         sum_of_weights += w
    #     end
    #     wx = sum_of_values / sum_of_weights
    #
    #     @turbo for i=1:n
    #         χ2 += (x[i] - wx) * (x[i] - wx) / (σ[i] * σ[i])
    #     end
    #
    #     return χ2 / (n-1)
    # end
    export MSWD

## -- The York (1968) two-dimensional linear regression with x and y uncertainties
    # as commonly used in isochrons

    # Custom type to hold York fit resutls
    struct YorkFit{T}
        intercept::T
        intercept_sigma::T
        slope::T
        slope_sigma::T
        mswd::T
    end

    """
    ```julia
    yorkfit(x, σx, y, σy)
    ```
    Uses the York (1968) two-dimensional least-squares fit to calculate `a`, `b`,
    and uncertanties `σa`, `σb` for the equation `y = a + bx`, given `x`, `y` and
    uncertaintes `σx`, ``σy`.
    """
    function yorkfit(x, σx, y, σy; niterations=10)

        ## 1. Ordinary linear regression (to get a first estimate of slope and intercept)

        # Check for missing data
        t = (x.==x) .& (y.==y) .& (σx.==σx) .& (σy.==σy)
        x = x[t]
        y = y[t]
        σx = σx[t]
        σy = σy[t]

        # Calculate the ordinary least-squares fit
        # For the equation y=a+bx, m(1)=a, m(2)=b
        g = [ones(length(x)) x]
        m = (g'*g)\g'*y
        b = m[2]
        a = m[1]

        ## 2. Now, let's define parameters needed by the York fit

        # Weighting factors
        ωx = 1.0 ./ σx.^2
        ωy = 1.0 ./ σy.^2

        # terms that don't depend on a or b
        α = sqrt.(ωx .* ωy)

        x̄ = _vmean(x, :)
        ȳ = _vmean(y, :)
        r = sum((x .- x̄).*(y .- ȳ)) ./ (sqrt(sum((x .- x̄).^2)) * sqrt(sum((y .- ȳ).^2)))

        ## 3. Perform the York fit (must iterate)
        W = ωx.*ωy ./ (b^2*ωy + ωx - 2*b*r.*α)

        X̄ = sum(W.*x) / sum(W)
        Ȳ = sum(W.*y) / sum(W)

        U = x .- X̄
        V = y .- Ȳ

        sV = W.^2 .* V .* (U./ωy + b.*V./ωx - r.*V./α)
        sU = W.^2 .* U .* (U./ωy + b.*V./ωx - b.*r.*U./α)
        b = sum(sV) ./ sum(sU)

        a = Ȳ - b .* X̄
        for i = 2:niterations
            W .= ωx.*ωy ./ (b^2*ωy + ωx - 2*b*r.*α)

            X̄ = sum(W.*x) / sum(W)
            Ȳ = sum(W.*y) / sum(W)

            U .= x .- X̄
            V .= y .- Ȳ

            sV .= W.^2 .* V .* (U./ωy + b.*V./ωx - r.*V./α)
            sU .= W.^2 .* U .* (U./ωy + b.*V./ωx - b.*r.*U./α)
            b = sum(sV) ./ sum(sU)

            a = Ȳ - b .* X̄
        end

        ## 4. Calculate uncertainties and MSWD
        β = W .* (U./ωy + b.*V./ωx - (b.*U+V).*r./α)

        u = X̄ .+ β
        v = Ȳ .+ b.*β

        xm = sum(W.*u)./sum(W)
        ym = sum(W.*v)./sum(W)

        σb = sqrt(1.0 ./ sum(W .* (u .- xm).^2))
        σa = sqrt(1.0 ./ sum(W) + xm.^2 .* σb.^2)

        # MSWD (reduced chi-squared) of the fit
        mswd = 1.0 ./ length(x) .* sum( (y .- a.-b.* x).^2 ./ (σy.^2 + b.^2 .* σx.^2) )

        ## Results
        return YorkFit(a, σa, b, σb, mswd)
    end
    export yorkfit

max_stars_repo_name: brenhinkeller/StatGeochemBase

max_issues_repo_name: brenhinkeller/StatGeochemBase

max_forks_repo_name: brenhinkeller/StatGeochemBase


text:
function run_powermodels_pf(json_path)
    pm = load_pm_from_json(json_path)
    active_powermodels_silence!(pm)
    pm = check_powermodels_data!(pm)

    # calculate branch power flows
    if pm["pm_model"] == "ACNative"
        result = _PM.compute_ac_pf(pm)
    elseif pm["pm_model"] == "DCNative"
        result = _PM.compute_dc_pf(pm)
    else
        model = get_model(pm["pm_model"])
        solver = get_solver(pm)
        result = _PM.run_pf(
            pm,
            model,
            solver,
            setting = Dict("output" => Dict("branch_flows" => true)),
        )
    end

    # add result to net data
    _PM.update_data!(pm, result["solution"])
    # calculate branch power flows
    if pm["ac"]
        flows = _PM.calc_branch_flow_ac(pm)
    else
        flows = _PM.calc_branch_flow_dc(pm)
    end
    # add flow to net and result
    _PM.update_data!(result["solution"], flows)
    # _PM.update_data!(pm, result["solution"])
    # _PM.update_data!(pm, flows)
    return result
end

function run_powermodels_opf(json_path)
    pm = load_pm_from_json(json_path)
    active_powermodels_silence!(pm)

    model = get_model(pm["pm_model"])
    solver = get_solver(pm)

    cl = check_current_limit!(pm)

    if cl == 0
        pm = check_powermodels_data!(pm)
        result = _PM.run_opf(
            pm,
            model,
            solver,
            setting = Dict("output" => Dict("branch_flows" => true)),
        )
    else

        for (key, value) in pm["gen"]
           value["pmin"] /= pm["baseMVA"]
           value["pmax"] /= pm["baseMVA"]
           value["qmax"] /= pm["baseMVA"]
           value["qmin"] /= pm["baseMVA"]
           value["pg"] /= pm["baseMVA"]
           value["qg"] /= pm["baseMVA"]
           value["cost"] *= pm["baseMVA"]
        end

        for (key, value) in pm["branch"]
           value["c_rating_a"] /= pm["baseMVA"]
        end

        for (key, value) in pm["load"]
           value["pd"] /= pm["baseMVA"]
           value["qd"] /= pm["baseMVA"]
        end

        result = _PM._run_opf_cl(
            pm,
            model,
            solver,
            setting = Dict("output" => Dict("branch_flows" => true)),
        )
    end

    return result
end

# function run_powermodels_custom(json_path)
#
#     pm = load_pm_from_json(json_path)
#
#     active_powermodels_silence!(pm)
#     pm = check_powermodels_data!(pm)
#
#     model = get_model(pm["pm_model"])
#     solver = get_solver(pm)
#
#     result = _PM.run_pf(pm, model, solver)
#     # add branch flows
#     _PM.update_data!(pm, result["solution"])
#     flows = _PM.calc_branch_flow_ac(pm)
#     _PM.update_data!(result["solution"], flows)
#     return result
# end

function run_powermodels_tnep(json_path)
    pm = load_pm_from_json(json_path)
    active_powermodels_silence!(pm)
    pm = check_powermodels_data!(pm)
    model = get_model(pm["pm_model"])
    solver = get_solver(pm)

    result = _PM.run_tnep(
        pm,
        model,
        solver,
        setting = Dict("output" => Dict("branch_flows" => true)),
    )
    return result
end

function run_powermodels_ots(json_path)
    pm = load_pm_from_json(json_path)
    active_powermodels_silence!(pm)
    pm = check_powermodels_data!(pm)
    model = get_model(pm["pm_model"])
    solver = get_solver(pm)

    result = _PM.run_ots(
        pm,
        model,
        solver,
        setting = Dict("output" => Dict("branch_flows" => true)),
    )
    return result
end


text:
abstract type AbstractJoint{T,N} end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, constraint::AbstractJoint{T,N}) where {T,N}
    summary(io, constraint)
end

### General functions
getT(joint::AbstractJoint{T}) where T = T
Base.length(joint::AbstractJoint{T,N}) where {T,N} = N

### Constraints and derivatives
## Position level constraint wrappers
@inline g(joint::AbstractJoint, body1::Body, body2::Body, Δt) = constraintmat(joint) * g(joint, body1.state, body2.state, Δt)
@inline g(joint::AbstractJoint, body1::Origin, body2::Body, Δt) = constraintmat(joint) * g(joint, body2.state, Δt)
@inline g(joint::AbstractJoint, body1::Body, body2::Body) = constraintmat(joint) * g(joint, body1.state, body2.state)
@inline g(joint::AbstractJoint, body1::Origin, body2::Body) = constraintmat(joint) * g(joint, body2.state)

## Discrete-time position derivatives (for dynamics)
# Wrappers 1
@inline function ∂g∂ʳposa(joint::AbstractJoint, body1::Body, body2::Body)
    if body2.id == joint.childid
        return constraintmat(joint) * ∂g∂ʳposa(joint, body1.state, body2.state)
    else
        return ∂g∂ʳposa(joint)
    end
end
@inline function ∂g∂ʳposb(joint::AbstractJoint, body1::Body, body2::Body)
    if body2.id == joint.childid
        return constraintmat(joint) * ∂g∂ʳposb(joint, body1.state, body2.state)
    else
        return ∂g∂ʳposb(joint)
    end
end
@inline function ∂g∂ʳposb(joint::AbstractJoint, body1::Origin, body2::Body)
    if body2.id == joint.childid
        return constraintmat(joint) * ∂g∂ʳposb(joint, body2.state)
    else
        return ∂g∂ʳposb(joint)
    end
end

# Wrappers 2
@inline ∂g∂ʳposa(joint::AbstractJoint{T,N}) where {T,N} = szeros(T, N, 6)
@inline ∂g∂ʳposb(joint::AbstractJoint{T,N}) where {T,N} = szeros(T, N, 6)
∂g∂ʳposa(joint::AbstractJoint, statea::State, stateb::State) = ∂g∂ʳposa(joint, posargsk(statea)..., posargsk(stateb)...)
∂g∂ʳposb(joint::AbstractJoint, statea::State, stateb::State) = ∂g∂ʳposb(joint, posargsk(statea)..., posargsk(stateb)...)
∂g∂ʳposb(joint::AbstractJoint, stateb::State) = ∂g∂ʳposb(joint, posargsk(stateb)...)

# Derivatives accounting for quaternion specialness
@inline function ∂g∂ʳposa(joint::AbstractJoint, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion)
    X, Q = ∂g∂posa(joint, xa, qa, xb, qb)
    Q = Q * LVᵀmat(qa)

    return [X Q]
end
@inline function ∂g∂ʳposb(joint::AbstractJoint, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion)
    X, Q = ∂g∂posb(joint, xa, qa, xb, qb)
    Q = Q * LVᵀmat(qb)

    return [X Q]
end
@inline function ∂g∂ʳposb(joint::AbstractJoint, xb::AbstractVector, qb::UnitQuaternion)
    X, Q = ∂g∂posb(joint, xb, qb)
    Q = Q * LVᵀmat(qb)

    return [X Q]
end

## Discrete-time velocity derivatives (for dynamics)
# Wrappers 1
@inline function ∂g∂ʳvela(joint::AbstractJoint, body1::Body, body2::Body, Δt)
    if body2.id == joint.childid
        return constraintmat(joint) * ∂g∂ʳvela(joint, body1.state, body2.state, Δt)
    else
        return ∂g∂ʳvela(joint)
    end
end
@inline function ∂g∂ʳvelb(joint::AbstractJoint, body1::Body, body2::Body, Δt)
    if body2.id == joint.childid
        return constraintmat(joint) * ∂g∂ʳvelb(joint, body1.state, body2.state, Δt)
    else
        return ∂g∂ʳvelb(joint)
    end
end
@inline function ∂g∂ʳvelb(joint::AbstractJoint, body1::Origin, body2::Body, Δt)
    if body2.id == joint.childid
        return constraintmat(joint) * ∂g∂ʳvelb(joint, body2.state, Δt)
    else
        return ∂g∂ʳvelb(joint)
    end
end

# Wrappers 2
@inline ∂g∂ʳvela(joint::AbstractJoint{T,N}) where {T,N} = szeros(T, N, 6)
@inline ∂g∂ʳvelb(joint::AbstractJoint{T,N}) where {T,N} = szeros(T, N, 6)
∂g∂ʳvela(joint::AbstractJoint, statea::State, stateb::State, Δt) = ∂g∂ʳvela(joint, posargsnext(statea, Δt)..., posargsnext(stateb, Δt)..., fullargssol(statea)..., Δt)
∂g∂ʳvelb(joint::AbstractJoint, statea::State, stateb::State, Δt) = ∂g∂ʳvelb(joint, posargsnext(statea, Δt)..., posargsnext(stateb, Δt)..., fullargssol(stateb)..., Δt)
∂g∂ʳvelb(joint::AbstractJoint, stateb::State, Δt) = ∂g∂ʳvelb(joint, posargsnext(stateb, Δt)..., fullargssol(stateb)..., Δt)

# Derivatives accounting for quaternion specialness
@inline function ∂g∂ʳvela(joint::AbstractJoint, x2a::AbstractVector, q2a::UnitQuaternion, x2b::AbstractVector, q2b::UnitQuaternion,
        x1a::AbstractVector, v1a::AbstractVector, q1a::UnitQuaternion, ω1a::AbstractVector, Δt
    )

    X, Q = ∂g∂posa(joint, x2a, q2a, x2b, q2b)
    V = X * Δt
    Ω = Q * Lmat(q1a) * derivωbar(ω1a, Δt)

    return [V Ω]
end
@inline function ∂g∂ʳvelb(joint::AbstractJoint, x2a::AbstractVector, q2a::UnitQuaternion, x2b::AbstractVector, q2b::UnitQuaternion,
        x1b::AbstractVector, v1b::AbstractVector, q1b::UnitQuaternion, ω1b::AbstractVector, Δt
    )

    X, Q = ∂g∂posb(joint, x2a, q2a, x2b, q2b)
    V = X * Δt
    Ω = Q * Lmat(q1b) * derivωbar(ω1b, Δt)

    return [V Ω]
end
@inline function ∂g∂ʳvelb(joint::AbstractJoint, x2b::AbstractVector, q2b::UnitQuaternion,
        x1b::AbstractVector, v1b::AbstractVector, q1b::UnitQuaternion, ω1b::AbstractVector, Δt
    )

    X, Q = ∂g∂posb(joint, x2b, q2b)
    V = X * Δt
    Ω = Q * Lmat(q1b) * derivωbar(ω1b, Δt)

    return [V Ω]
end

## Continuous-time position derivatives (for constraint solver)
# Wrappers 1
@inline function ∂g∂posac(joint::AbstractJoint, body1::Body, body2::Body)
    if body2.id == joint.childid
        return constraintmat(joint) * ∂g∂posac(joint, body1.state, body2.state)
    else
        return ∂g∂posac(joint)
    end
end
@inline function ∂g∂posbc(joint::AbstractJoint, body1::Body, body2::Body)
    if body2.id == joint.childid
        return constraintmat(joint) * ∂g∂posbc(joint, body1.state, body2.state)
    else
        return ∂g∂posbc(joint)
    end
end
@inline function ∂g∂posbc(joint::AbstractJoint, body1::Origin, body2::Body)
    if body2.id == joint.childid
        return constraintmat(joint) * ∂g∂posbc(joint, body2.state)
    else
        return ∂g∂posbc(joint)
    end
end

# Wrappers 2
@inline ∂g∂posac(joint::AbstractJoint{T,N}) where {T,N} = szeros(T, N, 7)
@inline ∂g∂posbc(joint::AbstractJoint{T,N}) where {T,N} = szeros(T, N, 7)
∂g∂posac(joint::AbstractJoint, statea::State, stateb::State) = hcat(∂g∂posa(joint, posargsc(statea)..., posargsc(stateb)...)...)
∂g∂posbc(joint::AbstractJoint, statea::State, stateb::State) = hcat(∂g∂posb(joint, posargsc(statea)..., posargsc(stateb)...)...)
∂g∂posbc(joint::AbstractJoint, stateb::State) = hcat(∂g∂posb(joint, posargsc(stateb)...)...)

### Forcing (for dynamics)
## Wrappers
@inline function applyFτ!(joint::AbstractJoint, body1::Body, body2::Body, clear::Bool)
    applyFτ!(joint, body1.state, body2.state, clear)
    return
end
@inline function applyFτ!(joint::AbstractJoint, ::Origin, body2::Body, clear::Bool)
    applyFτ!(joint, body2.state, clear)
    return
end


text:
import Solve
import JSON

function main()
  println("Listening at port "*ARGS[1])
  server = listen(ARGS[1])

  begin
    run = true
    while run
      sock = accept(server)
      println("Accepted connection.")
      while isopen(sock) && run
        line = readline(sock)
        if !isempty(line)
          command = JSON.parse(line)
          if command["type"] == "close"
            run = false
          else
            println("Running command...")
            result = run_command(command)
            println("Sending result...")
            write(sock,"[")
            JSON.print(sock,result)
            write(sock,"]\n")
          end
        end
      end
    end
  end
end

function run_command(command)
  try
    if command["type"] == "solve"
      Solve.solve_problem(command["problem"],command["prefs"])
    elseif command["type"] == "solvefile"
      prob = command["prob_file"]
      pref = command["pref_file"]
      newcommand = Dict("type" => "solve",
                     "problem" => JSON.parsefile(prob),
                     "prefs" => JSON.parsefile(pref))
      open(command["solution_file"],"w") do sol
        JSON.print(sol,run_command(newcommand))
      end
      Dict("result" => "saved")
    else
      Dict("result" => "error",
           "message" => "Unkonwn command type: "*command["type"])
    end
  catch e
    buffer = IOBuffer()
    showerror(buffer,e,catch_backtrace())
	estr = takebuf_string(buffer)
	println(estr)
    Dict("result" => "error",
         "message" => "Exception thrown: \n"*estr)
  end
end

prob = "/Users/davidlittle/Desktop/TA2016_17.json"
pref = "data/prefs.json"

function run_solver(prob,pref,solutionfile)
  open(solutionfile,"w") do sol
    command = Dict("type" => "solve",
                   "problem" => JSON.parsefile(prob),
                   "prefs" => JSON.parsefile(pref))
    JSON.print(sol,run_command(command))
  end
end

max_stars_repo_name: haberdashPI/csd-ta-assign

max_issues_repo_name: haberdashPI/csd-ta-assign

max_forks_repo_name: haberdashPI/csd-ta-assign


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

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

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

const _property_map_ContainerServiceOrchestratorProfile = Dict{Symbol,Symbol}(Symbol("orchestratorType")=>Symbol("orchestratorType"))
const _property_types_ContainerServiceOrchestratorProfile = Dict{Symbol,String}(Symbol("orchestratorType")=>"String")
Base.propertynames(::Type{ ContainerServiceOrchestratorProfile }) = collect(keys(_property_map_ContainerServiceOrchestratorProfile))
Swagger.property_type(::Type{ ContainerServiceOrchestratorProfile }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ContainerServiceOrchestratorProfile[name]))}
Swagger.field_name(::Type{ ContainerServiceOrchestratorProfile }, property_name::Symbol) =  _property_map_ContainerServiceOrchestratorProfile[property_name]

const _allowed_ContainerServiceOrchestratorProfile_orchestratorType = ["Swarm", "DCOS", "Custom", "Kubernetes"]

function check_required(o::ContainerServiceOrchestratorProfile)
    true
end

function validate_property(::Type{ ContainerServiceOrchestratorProfile }, name::Symbol, val)
    if name === Symbol("orchestratorType")
        Swagger.validate_param(name, "ContainerServiceOrchestratorProfile", :enum, val, _allowed_ContainerServiceOrchestratorProfile_orchestratorType)
    end
end


text:
type Mem
    lastq::Array{Int64,1}
    seed_save::Int64
    v::Matrix{Int64}
    recipd::Float64
    poly::Vector{Int64}
    Mem(dim::Int64) = new(Array(Int64,dim), 0, zeros(Int64,40,30), 0.0, Vector{Int64}(40))
end

"""
# % Usage:
# %   X = sobol(dim, N)
# %
# % Inputs:
# %    dim                number of variables, the MAX number of variables is 40
# %    N                  number of samples
# %
# % Output:
# %     X                matrix [N x dim] with the quasti-random samples
# %
# % ------------------------------------------------------------------------
"""
function sobol(dim::Int64, N::Int64)
    X = Array(Float64,(N,dim))
    nextseed, MeM = sobol!(X, dim, N)

    return X, nextseed, MeM
end

#Create new sobol generator and store values in X
function sobol!(X, dim::Int64, N::Int64)
    nextseed::Int64 = 2^floor(log2(N)+1)

    assert(size(X) == (N,dim))

    MeM = InitSobol(dim)

    nextseed  = getNewSobolArray!(dim, nextseed, MeM, X, N)

    return nextseed, MeM
end

#Get next set of numbers and store in X
function sobol!(X, nextseed, MeM::Mem)
    (N,dim) = size(X)
    nextseed = getNewSobolArray!(dim, nextseed, MeM, X, N)

    return nextseed
end

function getNewSobolArray!(dim, seed, MeM, X, N)
    for j = 1:N
        seed::Int64 = max(floor(seed),0)
        if  seed == 0
            l = 1
            MeM.lastq = zeros(1,dim)
        elseif  seed == MeM.seed_save + 1
            l = smallest0bit(seed)
        elseif  seed <= MeM.seed_save
            MeM.seed_save = 0
            l = 1
            MeM.lastq[1:dim] = 0

            for seed_temp = MeM.seed_save:seed-1
                l = smallest0bit(seed_temp)
                for i = 1:dim
                    MeM.lastq[i] = MeM.lastq[i] $ MeM.v[i,l] # bitxor
               	end
            end

            l = smallest0bit(seed)

        elseif (MeM.seed_save + 1) < seed

            for seed_temp = (MeM.seed_save + 1):(seed - 1)
               	l = smallest0bit(seed_temp)
               	for i = 1:dim
                    MeM.lastq[i] = MeM.lastq[i] $ MeM.v[i,l] # bitxor
               	end
            end

            l = smallest0bit(seed)

        end

        #qrv = Array(Float64,dim)

        for i = 1 : dim
            X[j,i] = MeM.lastq[i] * MeM.recipd; #qrv[i] = MeM.lastq[i] * MeM.recipd;
            MeM.lastq[i] =  MeM.lastq[i] $ MeM.v[i,l]  # bitxor
        end

        MeM.seed_save = seed;
        seed = seed + 1;
    end
    return seed
end

function smallest0bit(b)
    i = 0
    b = floor(b)
    while true
        i += 1
        halfb = floor(b/2)
        if b == 2halfb
            break
       	end
        b = halfb
    end
    return i
end

function InitSobol(dim::Int64)
    MeM = Mem(dim)
    MeM.seed_save = -1
    #     MeM.v = zeros(40,30) # done in constructor instead

    MeM.v[1:40,1] = 1
    MeM.v[3:40,2] = [1, 3, 1, 3, 1, 3, 3, 1,
                     3, 1, 3, 1, 3, 1, 1, 3, 1, 3,
                     1, 3, 1, 3, 3, 1, 3, 1, 3, 1,
                     3, 1, 1, 3, 1, 3, 1, 3, 1, 3 ]'

    MeM.v[4:40,3] = [7, 5, 1, 3, 3, 7, 5,
                     5, 7, 7, 1, 3, 3, 7, 5, 1, 1,
                     5, 3, 3, 1, 7, 5, 1, 3, 3, 7,
                     5, 1, 1, 5, 7, 7, 5, 1, 3, 3 ]'

    MeM.v[6:40,4] = [1, 7, 9,13,11,
                     1, 3, 7, 9, 5,13,13,11, 3,15,
                     5, 3,15, 7, 9,13, 9, 1,11, 7,
                     5,15, 1,15,11, 5, 3, 1, 7, 9 ]';

    MeM.v[8:40,5] = [9, 3,27,
                     15,29,21,23,19,11,25, 7,13,17,
                     1,25,29, 3,31,11, 5,23,27,19,
                     21, 5, 1,17,13, 7,15, 9,31, 9 ]'

    MeM.v[14:40,6] = [37,33, 7, 5,11,39,63,
                      27,17,15,23,29, 3,21,13,31,25,
                      9,49,33,19,29,11,19,27,15,25 ]'

    MeM.v[20:40,7] = [13,
                      33,115, 41, 79, 17, 29,119, 75, 73,105,
                      7, 59, 65, 21,  3,113, 61, 89, 45,107 ]';

    MeM.v[38:40,8] = [7, 23, 39 ]'
    MeM.poly[1:40]= [1,   3,   7,  11,  13,  19,  25,  37,  59,  47,
                     61,  55,  41,  67,  97,  91, 109, 103, 115, 131,
                     193, 137, 145, 143, 241, 157, 185, 167, 229, 171,
                     213, 191, 253, 203, 211, 239, 247, 285, 369, 299 ]
   MeM.v[1,:] = 1
   for i = 2 : dim
       j = MeM.poly[i]
       m = 0

      	while true
           j = floor( j / 2 )
           if  j <= 0
               break
           end
           m += 1
       end

       j = MeM.poly[i]
       includ = Array(Bool,m)
       for k = m : -1 : 1
           j2 = floor( j / 2 )
           includ[k] =  j != 2 * j2
           j = j2
       end

      	for j = m + 1 : 30
           newv = MeM.v[i,j-m]
           l = 1
           for k = 1 : m
               l = 2 * l
               if  includ[k] != 0
                   newv =  newv $ l * MeM.v[i,j-k] # bitxor
               end
           end
           MeM.v[i,j] = newv
       end
     end

     l = 1
     for j = 29 : -1 : 1
         l = 2 * l
         MeM.v[:,j] = MeM.v[:,j] * l
     end

     MeM.recipd = 1.0 / ( 2 * l )
     MeM.lastq[1:dim] = 0
     return MeM
 end

 """
 `using Winston`
 Plots the first 4 outputs of sobel(2,512)
 """
 function test_sobol()
     X, nextseed, MeM = sobol(2,512)
     figure(); pp = plot(X[:,1],X[:,2],".b"); hold(true)

     X, nextseed, MeM = sobol(X, nextseed, MeM )
     plot(X[:,1],X[:,2],".g")

     X, nextseed, MeM = sobol(X, nextseed, MeM )
     plot(X[:,1],X[:,2],".r")

     X, nextseed, MeM = sobol(X, nextseed, MeM )
     plot(X[:,1],X[:,2],".m")

     display(pp)
 end


text:
lattice_primitives(constants::NTuple{1}) = (Point(1.0),), SMatrix{1,1}(1.0)

Base.getindex(lat::Lattice{1},i::CartesianIndex) = getindex(lat,i[1])
Base.getindex(lat::Lattice{1},i) = map(i->lat[i],i)
function Base.getindex(lat::Lattice{1},i::Real)
    if lat.type==:Cartesian
        return lat.origin + i*lat.dx*lat.e1
    else
        r = lat.r0 + i*lat.dr
        return lat.origin + r
    end
end

# latticeindex(lat::Lattice{1},p::Point{1}) = (p-lat.origin).x/lat.dx
latticeindex(lat::Lattice{1},p::Point{1}) = (p[1]-lat.origin[1])/lat.dx

function Base.show(io::IO, lat::Lattice{1})
    if lat.type==:Cartesian
        if !get(io, :sub, false)
            print(io, "Cartesian Lattice: \n",
                "\tdx = ", fmt("3.2f",lat.dx), "\n",
                "\torigin: ", fmt("3.2f",lat.x0))
        elseif !get(io, :sub1, false)
            print(io,
                "\n\tdx = ", fmt("3.2f",lat.dx), "\n",
                "\torigin: ", fmt("3.2f",lat.x0))
        else
            print(io,
                "\n\t\tdx = ", fmt("3.2f",lat.dx),"\n",
                "\t\torigin: ", fmt("3.2f",lat.x0))
        end
    elseif lat.type==:Polar
        println(io,"not yet implemented")
    else
        println(io, "unknown lattice type")
    end
end


text:
## https://github.com/macd/BaryRational.jl/blob/main/src/aaa.jl
## https://arxiv.org/pdf/1612.00337.pdf for paper

import BaryRational: aaa, prz
export aaa, prz

"""
    BaryRationalFunction(xs, ys)
    fit(BaryRationalFunction, xs, ys)

Fit a rational function to the points `(xs,ys)` using the `aaa` method from [BaryRational](https://github.com/macd/BaryRational.jl/).

References: Nakatsukasa, Sete, and Trefethen; THE AAA ALGORITHM FOR RATIONAL APPROXIMATION [	arXiv:1612.00337](https://arxiv.org/pdf/1612.00337.pdf)

> Even on a disk or interval the algorithm may outperform existing methods, and on more complicated domains it is especially competitive. The core ideas are (1) representation of the rational approximant in barycentric form with interpolation at certain support points and (2) greedy selection of the support points to avoid exponential instabilities. The name AAA stands for "adaptive Antoulas--Anderson" in honor of the authors who introduced a scheme based on [Approximation of Large-Scale Dynamical Systems](https://doi.org/10.1137/1.9780898718713).

"""
struct BaryRationalFunction{T,X,P} <: AbstractRationalFunction{T,X,P}
    aaa::BaryRational.AAAapprox{Array{T,1}}
    function BaryRationalFunction(xs, ys; var=:x, kwargs...)
        o = BaryRational.aaa(collect(xs), collect(ys), kwargs...)
        T = eltype(o.x)
        X = var
        P = Polynomial{T, X}
        new{T,X,P}(o)
    end
end

## need a convert(BaryRationalFunction, p::RationalFunction) method...
BaryRationalFunction(p::AbstractPolynomial, q::AbstractPolynomial) =
    throw(ArgumentError("Method not defined for BaryRationalFunction"))

Polynomials.fit(::Type{<:BaryRationalFunction}, xs, ys; var=:x, kwargs...) =
    BaryRationalFunction(xs, ys; var=var, kwargs...)

Base.show(io::IO, F::BaryRationalFunction) = print(io, "BaryRational approximation")

Base.numerator(pq::BaryRationalFunction) = pqs(pq)[1]
Base.denominator(pq::BaryRationalFunction) = pqs(pq)[2]

# evaluation
(F::BaryRationalFunction)(x) = F.aaa(x)

function Base.convert(PQ::Type{RationalFunction}, F::BaryRationalFunction)
    P = eltype(PQ)
    z = variable(P) // one(P)
    F.aaa(z)
end

function pqs(F::BaryRationalFunction{T,X,P}) where {T,X,P}
    pq = convert(RationalFunction, F)
    pqs(pq)
end


text:
struct Stats{A,B,C,D,I}
    number_stats::OrderedDict{Symbol, A}
    string_stats::OrderedDict{Symbol, B}
    bool_stats::OrderedDict{Symbol, C}
    null_stats::OrderedDict{Symbol, D}
    interval::I

    function Stats(interval::I = Second(10)) where {I}
        A = typeof(init_stat(Number))
        B = typeof(init_stat(String))
        C = typeof(init_stat(Bool))
        D = typeof(init_stat(Nothing))
        new{A,B,C,D,I}(
            OrderedDict{Symbol, A}(),
            OrderedDict{Symbol, B}(),
            OrderedDict{Symbol, C}(),
            OrderedDict{Symbol, D}(),
            interval
        )
    end
end

init_stat(::Type{<:Number}) = GroupBy(DateTime, Series(Mean(), Extrema(), ExpandingHist(100)))
init_stat(::Type{String}) = GroupBy(DateTime, CountMap(String))
init_stat(::Type{Bool}) = GroupBy(DateTime, CountMap(Bool))
init_stat(::Type{Nothing}) = GroupBy(DateTime, Counter(Nothing))

get_stats(o::Stats, ::Type{<:Number}) = o.number_stats
get_stats(o::Stats, ::Type{String}) = o.string_stats
get_stats(o::Stats, ::Type{Bool}) = o.bool_stats
get_stats(o::Stats, ::Type{Nothing}) = o.null_stats

#-----------------------------------------------------------------------------# fit!
function StatsBase.fit!(o::Stats, x::Pair{Symbol, T}) where {T}
    stat = get(get_stats(o, T), x[1], init_stat(T))
    get_stats(o, T)[x[1]] = fit!(stat, round(now(), o.interval) => x[2])
    o
end

#-----------------------------------------------------------------------------# printing
Base.show(io::IO, o::Stats) = print_tree(io, o)
AbstractTrees.printnode(io::IO, o::Stats) = print(io, "Stats | Every: $(o.interval)")
function AbstractTrees.children(o::Stats)
    d = reduce((a,b) -> merge(vcat, a, b), [o.number_stats, o.string_stats, o.bool_stats, o.null_stats])
    collect(d)
end

#-----------------------------------------------------------------------------# Plot
function _xlim(o::Stats)
    a, b = typemax(DateTime), typemin(DateTime)
    for dict in [o.number_stats, o.string_stats, o.bool_stats, o.null_stats]
        for (id, gb) in dict
            a = min(a, minimum(keys(gb.value)))
            b = max(b, maximum(keys(gb.value)))
        end
    end
    Dates.value(a), Dates.value(b)
end

@recipe function f(o::Stats, type=3)
    i = 0
    layout --> sum(length, (o.number_stats, o.string_stats, o.bool_stats, o.null_stats))
    xlims --> _xlim(o)
    xrotation --> 15
    for (id, gb) in o.number_stats
        i += 1
        @series begin
            subplot --> i
            title --> string(id)
            [(t, t + o.interval) => Series(series.stats[type]) for (t, series) in value(gb)]
        end
    end
    for (id, gb) in Iterators.flatten((o.string_stats , o.bool_stats, o.null_stats))
        i += 1
        @series begin
            subplot --> i
            title --> string(id)
            [(t, t + o.interval) => stat for (t, stat) in value(gb)]
        end
    end
end

# @recipe function f(o::Stats)
#     dicts = [o.number_stats, o.string_stats, o.bool_stats, o.null_stats]
#     nplots = sum(length(keys(d)) for d in dicts)
#     a, b = typemax(DateTime), typemin(DateTime)
#     layout --> nplots
#     i = 0
#     for dict in dicts
#         i += 1
#         for (k, v) in dict
#             @series begin
#                 label --> string(k)
#                 subplot --> i
#                 a = min(a, minimum(keys(v.value)))
#                 b = max(b, maximum(keys(v.value)))
#                 [(t, t + o.interval) => stat for (t,stat) in v.value]
#             end
#         end
#     end
#     @info a, b
#     xlim --> (a, b)
#     nothing
# end


text:
export scan

import Base: show

"""
    scan(::Type{R}, scanFn::F, seed::R) where { R, F <: Function }
    scan(scanFn::F) where { F <: Function }

Creates a scan operator, which applies a given accumulator `scanFn` function to each value emmited by the source
Observable, and returns each intermediate result with an optional seed value. If a seed value is specified,
then that value will be used as the initial value for the accumulator.
If no seed value is specified, the first item of the source is used as the seed.

# Arguments
- `::Type{R}`: the type of data of transformed value, may be or may not be the same as type of input source
- `scanFn::Function`: accumulator function with `(data::T, current::R) -> R` signature
- `seed::R`: optional initial value for accumulator function

# Producing

Stream of type `<: Subscribable{R}`

# Examples
```jldoctest
using Rocket

source = from([ 1, 2, 3 ])
subscribe!(source |> scan(+), logger())
;

# output

[LogActor] Data: 1
[LogActor] Data: 3
[LogActor] Data: 6
[LogActor] Completed

```

```jldoctest
using Rocket

source = from([ 1, 2, 3 ])
subscribe!(source |> scan(Vector{Int}, (d, c) -> [ c..., d ], Int[]), logger())
;

# output

[LogActor] Data: [1]
[LogActor] Data: [1, 2]
[LogActor] Data: [1, 2, 3]
[LogActor] Completed

```

See also: [`AbstractOperator`](@ref), [`RightTypedOperator`](@ref), [`ProxyObservable`](@ref), [`reduce`](@ref), [`logger`](@ref)
"""
scan(::Type{R}, scanFn::F, seed::R) where { R, F <: Function } = ScanOperator{R, F}(scanFn, seed)

# ------------------------------------------------------------------------------------------------ #
# Seed version of scan operator (typed with R also)
# ------------------------------------------------------------------------------------------------ #

struct ScanOperator{R, F} <: RightTypedOperator{R}
    scanFn  :: F
    seed    :: R
end

function on_call!(::Type{L}, ::Type{R}, operator::ScanOperator{R, F}, source) where { L, R, F }
    return proxy(R, source, ScanProxy{L, R, F}(operator.scanFn, operator.seed))
end

struct ScanProxy{L, R, F} <: ActorProxy
    scanFn  :: F
    seed    :: R
end

actor_proxy!(::Type{R}, proxy::ScanProxy{L, R, F}, actor::A) where { L, R, A, F } = ScanActor{L, R, A, F}(proxy.scanFn, actor, proxy.seed)

mutable struct ScanActor{L, R, A, F} <: Actor{L}
    scanFn  :: F
    actor   :: A
    current :: R
end

getcurrent(actor::ScanActor)         = actor.current
setcurrent!(actor::ScanActor, value) = actor.current = value

function on_next!(actor::ScanActor{L, R}, data::L) where { L, R }
    update = actor.scanFn(data, getcurrent(actor))
    setcurrent!(actor, update)
    next!(actor.actor, update)
end

on_error!(actor::ScanActor, err) = error!(actor.actor, err)
on_complete!(actor::ScanActor)   = complete!(actor.actor)

Base.show(io::IO, ::ScanOperator{R}) where R = print(io, "ScanOperator( -> $R)")
Base.show(io::IO, ::ScanProxy{L})    where L = print(io, "ScanProxy($L)")
Base.show(io::IO, ::ScanActor{L})    where L = print(io, "ScanActor($L)")

# ------------------------------------------------------------------------------------------------ #
# No seed version of scan operator (output data stream type is inferred from input)
# ------------------------------------------------------------------------------------------------ #

scan(scanFn::F) where { F <: Function } = ScanNoSeedOperator{F}(scanFn)

struct ScanNoSeedOperator{F} <: InferableOperator
    scanFn  :: F
end

operator_right(operator::ScanNoSeedOperator, ::Type{L}) where L = L

function on_call!(::Type{L}, ::Type{L}, operator::ScanNoSeedOperator{F}, source) where { L, F }
    return proxy(L, source, ScanNoSeedProxy{F}(operator.scanFn))
end

struct ScanNoSeedProxy{F} <: ActorProxy
    scanFn :: F
end

actor_proxy!(::Type{L}, proxy::ScanNoSeedProxy{F}, actor::A) where { L, A, F } = ScanNoSeedActor{L, A, F}(proxy.scanFn, actor, nothing)

mutable struct ScanNoSeedActor{L, A, F} <: Actor{L}
    scanFn  :: F
    actor   :: A
    current :: Union{L, Nothing}
end

getcurrent(actor::ScanNoSeedActor)         = actor.current
setcurrent!(actor::ScanNoSeedActor, value) = actor.current = value

function on_next!(actor::ScanNoSeedActor{L}, data::L) where L
    current = getcurrent(actor)
    if current === nothing
        setcurrent!(actor, data)
    else
        setcurrent!(actor, actor.scanFn(data, current))
    end
    next!(actor.actor, getcurrent(actor))
end

on_error!(actor::ScanNoSeedActor, err) = error!(actor.actor, err)
on_complete!(actor::ScanNoSeedActor)   = complete!(actor.actor)

Base.show(io::IO, ::ScanNoSeedOperator)         = print(io, "ScanNoSeedOperator(L -> L)")
Base.show(io::IO, ::ScanNoSeedProxy)            = print(io, "ScanNoSeedProxy()")
Base.show(io::IO, ::ScanNoSeedActor{L}) where L = print(io, "ScanNoSeedActor($L)")


text:
# Tests for the functions library in deps/examples

using CxxWrap
using Base.Test
using Compat
cxx_available = true
try
  using Cxx
catch
  cxx_available = false
end

const functions_lib_path = CxxWrap._l_functions

# Wrap the functions defined in C++
wrap_modules(functions_lib_path)

# Test functions from the CppHalfFunctions module
@test CppHalfFunctions.half_d(3) == 1.5
@show methods(CppHalfFunctions.half_d)
@test CppHalfFunctions.half_i(-2) == -1
@test CppHalfFunctions.half_u(3) == 1
@test CppHalfFunctions.half_lambda(2.) == 1.

# Test functions from the CppTestFunctions module
@test CppTestFunctions.concatenate_numbers(4, 2.) == "42"
@test length(methods(CppTestFunctions.concatenate_numbers)) == (Sys.WORD_SIZE == 64 ? 4 : 2) # due to overloads
@test CppTestFunctions.concatenate_strings(2, "ho", "la") == "holahola"
@test CppTestFunctions.test_int32_array(Int32[1,2])
@test CppTestFunctions.test_int64_array(Int64[1,2])
@test CppTestFunctions.test_float_array(Float32[1.,2.])
@test CppTestFunctions.test_double_array([1.,2.])
if !(is_windows() && Sys.WORD_SIZE == 32)
  @test_throws ErrorException CppTestFunctions.test_exception()
end
ta = [1.,2.]
@test CppTestFunctions.test_array_len(ta) == 2
@test CppTestFunctions.test_array_get(ta, Int64(0)) == 1.
@test CppTestFunctions.test_array_get(ta, Int64(1)) == 2.
CppTestFunctions.test_array_set(ta, Int64(0), 3.)
CppTestFunctions.test_array_set(ta, Int64(1), 4.)
@test ta[1] == 3.
@test ta[2] == 4.
@test CppTestFunctions.test_type_name("IO") == "IO"

@test CppTestFunctions.test_long_long() == 42
@test CppTestFunctions.test_short() == 43

# Test GC protection array
a = "str1"
b = "str2"
c = "str3"
protect_arr = CxxWrap._gc_protected
start_len = length(protect_arr)
CppTestFunctions.test_protect_from_gc(a)
CppTestFunctions.test_protect_from_gc(b)
@test length(protect_arr) == start_len + 2
@test protect_arr[end-1] == a
@test protect_arr[end] == b
CppTestFunctions.test_unprotect_from_gc(a)
@test length(protect_arr) == start_len + 2
@test protect_arr[end-1] == nothing
@test protect_arr[end] == b
CppTestFunctions.test_protect_from_gc(c)
@test length(protect_arr) == start_len + 2
@test protect_arr[end-1] == c
@test protect_arr[end] == b
@test CppTestFunctions.test_julia_call(1,2) == 2

# Performance tests
const test_size = 50000000
const numbers = rand(test_size)
output = zeros(test_size)

# Build a function to loop over the test array
function make_loop_function(name)
    fname = Symbol(:half_loop_,name,:!)
    inner_name = Symbol(:half_,name)
    @eval begin
        function $(fname)(n::Array{Float64,1}, out_arr::Array{Float64,1})
            test_length = length(n)
          for i in 1:test_length
                out_arr[i] = $(inner_name)(n[i])
          end
        end
    end
end

# Julia version
half_julia(d::Float64) = d*0.5

# C version
half_c(d::Float64) = ccall((:half_c, functions_lib_path), Cdouble, (Cdouble,), d)

# Bring C++ versions into scope
using CppHalfFunctions.half_d, CppHalfFunctions.half_lambda, CppHalfFunctions.half_loop_cpp!, CppHalfFunctions.half_loop_jlcall!

@static if cxx_available
  # Cxx.jl version
  cxx"""
  double half_cxx(const double d)
  {
    return 0.5*d;
  }
  """
  half_cxxjl(d::Float64) = @cxx half_cxx(d)
end

# Make the looping functions
make_loop_function(:julia)
make_loop_function(:c)
make_loop_function(:d) # C++ with regular C++ function pointer
make_loop_function(:lambda) # C++ lambda, so using std::function
if cxx_available
  make_loop_function(:cxxjl) # Cxx.jl version
end

# test that a "half" function does what it should
function test_half_function(f)
  input = [2.]
  output = [0.]
  f(input, output)
  @test output[1] == 1.
end
test_half_function(half_loop_julia!)
test_half_function(half_loop_c!)
test_half_function(half_loop_d!)
test_half_function(half_loop_lambda!)
test_half_function(half_loop_cpp!)
if cxx_available
  test_half_function(half_loop_cxxjl!)
end

# Run timing tests
println("---- Half test timings ----")
println("Julia test:")
@time half_loop_julia!(numbers, output)
@time half_loop_julia!(numbers, output)
@time half_loop_julia!(numbers, output)

println("C test:")
@time half_loop_c!(numbers, output)
@time half_loop_c!(numbers, output)
@time half_loop_c!(numbers, output)

println("C++ test:")
@time half_loop_d!(numbers, output)
@time half_loop_d!(numbers, output)
@time half_loop_d!(numbers, output)

if cxx_available
  println("Cxx.jl test:")
  @time half_loop_cxxjl!(numbers, output)
  @time half_loop_cxxjl!(numbers, output)
  @time half_loop_cxxjl!(numbers, output)
end

println("C++ lambda test:")
@time half_loop_lambda!(numbers, output)
@time half_loop_lambda!(numbers, output)
@time half_loop_lambda!(numbers, output)

println("C++ test, loop in the C++ code:")
@time half_loop_cpp!(numbers, output)
@time half_loop_cpp!(numbers, output)
@time half_loop_cpp!(numbers, output)

const small_in = rand(test_size÷100)
small_out = zeros(test_size÷100)

println("jl_call inside C++ loop (array is 100 times smaller than other tests):")
@time half_loop_jlcall!(small_in, small_out)
@time half_loop_jlcall!(small_in, small_out)
@time half_loop_jlcall!(small_in, small_out)


text:
module BlockTrees

import ClusterTrees

struct BlockTree{T}
    test_cluster::T
    trial_cluster::T
end

testcluster(blocktree::BlockTree) = blocktree.test_cluster
trialcluster(blocktree::BlockTree) = blocktree.trial_cluster

ClusterTrees.root(t::BlockTree) = (ClusterTrees.root(testcluster(t)), ClusterTrees.root(trialcluster(t)))
function ClusterTrees.data(t::BlockTree, n)
    (
        ClusterTrees.data(testcluster(t), n[1]),
        ClusterTrees.data(trialcluster(t), n[2]))
end

function ClusterTrees.children(b::BlockTree, node)
    test_chds = ClusterTrees.children(testcluster(b), node[1])
    trial_chds = ClusterTrees.children(trialcluster(b), node[2])
    ((ch[1],ch[2]) for ch in Iterators.product(test_chds, trial_chds))
end

function ClusterTrees.haschildren(b::BlockTree, node)
    !ClusterTrees.haschildren(testcluster(b), node[1]) && return false
    !ClusterTrees.haschildren(trialcluster(b), node[2]) && return false
    return true
end

function ClusterTrees.LevelledTrees.numlevels(bt)
    num_test_levels = ClusterTrees.LevelledTrees.numlevels(testcluster(bt))
    num_trial_levels = ClusterTrees.LevelledTrees.numlevels(trialcluster(bt))
    @assert num_test_levels == num_trial_levels
    return num_test_levels
end

end # module BlockTrees


text:
module Schedulers

include("../../docstrings.jl")

import ..Rule
import ..EGraph
using ..Util
using ..Rules

include("./abstractscheduler.jl")
include("./backoffscheduler.jl")
include("./simplescheduler.jl")
include("./scoredscheduler.jl")

export AbstractScheduler
export SimpleScheduler
export BackoffScheduler
export ScoredScheduler

export cansaturate
export cansearch
export inform!
export setiter!

end


text:
"""
predict(state::EKFSlamState, vehicle::Vehicle, Q::AbstractMatrix, dt::AbstractFloat)

Predict SLAM pose and covariance matrix according to equations of motion for vehicle.
Q is the covariance matrix for speed and steering angle gamma.
dt is the timestep
"""
function predict(state::EKFSlamState, vehicle::Vehicle, Q::AbstractMatrix, dt::AbstractFloat)

    x = state.x
    P = state.cov
    phi = x[3]

    g = vehicle.measured_gamma
    v = vehicle.measured_speed
    w = vehicle.wheelbase

    s = sin(g + phi)
    c = cos(g + phi)
    vts = v*dt*s
    vtc = v*dt*c

    # Jacobians
    Gv = [1 0 -vts;
          0 1 vtc;
          0 0 1]
    Gu = [dt*c -vts;
          dt*s  vtc;
          dt*sin(g)/w v*dt*cos(g)/w]

    # Predict covariance
    P[1:3,1:3] = Gv*P[1:3,1:3]*Gv' + Gu*Q*Gu'
    if size(P, 1) > 3
        P[1:3,4:end] = Gv*P[1:3,4:end]
        P[4:end,1:3] = P[1:3,4:end]'
    end

    # Predict pose
    x[1:3] = [x[1] + vtc;
              x[2] + vts;
              mpi_to_pi(phi + v*dt*sin(g)/w)]
    x, P
end

function update(state::EKFSlamState, z, R, idf)
    x = state.x
    P = state.cov

    lenz = size(z,2)
    lenx = length(x)
    H = zeros(2*lenz, lenx)
    v = zeros(2*lenz)
    RR = zeros(2*lenz, 2*lenz)
    for i = 1:lenz
        ii = 2*i + (-1:0)
        zp, H[ii,:] = predict_observation(x, idf[i])

        v[ii] = [z[1,i] - zp[1]; mpi_to_pi(z[2,i] - zp[2])]
        RR[ii,ii] = R
    end

    # Calculate the KF (or EKF) update given the prior state [x,P]
    # the innovation [v,RR] and the (linearised) observation model H.
    # The result is calculated using Cholesky factorisation, which
    # is more numerically stable than a naive implementation.
    PHt = P*H'
    S = H*PHt + RR
    S = (S + S')*0.5 # make symmetric
    C = inv(chol(S)) # triangular matrix
    W1 = PHt*C
    W = W1*C'

    x += W*v
    P -= W1*W1'
    x, P
end

"""
Add new feature/landmark measurements z and return augmented state
vector x and covariance matrix P
"""
function add_features(state::EKFSlamState, z, R)

    x = state.x
    P = state.cov
    phi = x[3]

    for i = 1:size(z, 2)
        len = length(x)

        # Range and bearing measurement i
        r, b = z[1,i], z[2,i]

        s, c = sin(phi + b), cos(phi + b)

        # augment x
        x = [x; x[1] + r*c; x[2] + r*s]

        # jacobians
        Gv = [1 0 -r*s; 0 1 r*c]
        Gz = [c -r*s; s r*c]

        # Augment P
        # TODO: Implement resizing more efficiently(?) For ideas, see
        # https://groups.google.com/forum/#!topic/julia-users/B4OUYPFM5L8
        P = hcat(P, zeros(size(P, 1), 2))
        P = vcat(P, zeros(2, size(P, 2)))

        rng = len+1:len+2
        P[rng,rng] = Gv*P[1:3,1:3]*Gv' + Gz*R*Gz' # feature cov
        P[rng,1:3] = Gv*P[1:3,1:3] # vehicle to feature xcorr
        P[1:3,rng] = P[rng,1:3]'
        if len > 3
            rnm = 4:len
            P[rng,rnm] = Gv*P[1:3,rnm] # map to feature xcorr
            P[rnm,rng] = P[rng,rnm]'
        end
    end
    x, P
end


text:
abstract type FieldProperty end
struct ChargeDensity <: FieldProperty end
struct CurrentDensity{D} <: FieldProperty end

mutable struct Particles{N}
    positions::Array{Float,2}
    velocities::Array{Float,2}
    accelerations::Array{Float,2}
    mass::Float
    charge::Float
    dim::Dim{N}

    function Particles{N}(positions::Array{<:Any,2}, mass, charge) where {N}
        dims, count = size(positions)
        @checkeq(dims, N)
        velocities = zeros(N, count)
        accelerations = zeros(N, count)
        new(positions, velocities, accelerations, mass, charge, Dim{N}())
    end
end

count(particles::Particles) = size(particles.positions)[2]
mass(particles::Particles) = particles.mass
charge(particles::Particles) = particles.charge

function advance_positions!(particles::Particles, Δt)
    particles.positions .+= particles.velocities .* Δt
end

function advance_velocities!(particles::Particles, Δt)
    particles.velocities .+= particles.accelerations .* Δt
end

function update_accelerations!(
    particles::Particles{N},
    fields::Fields{N},
) where {N}
    E_interp = VectorFieldInterpolator(fields, ElectricField)
    B_interp = VectorFieldInterpolator(fields, MagneticField)

    q_over_m = charge(particles) / mass(particles)
    for p = 1:count(particles)
        r = particles.positions[:, p]
        v = particles.velocities[:, p]
        E = interpolate(E_interp, r)
        B = interpolate(B_interp, r)
        particles.accelerations[:, p] .=
            lorentz_force(q_over_m, v, E, B, Dim{N}())
    end
end

function lorentz_force(q_over_m, v, E, B, dim::Dim{N}) where {N}
    q_over_m .* (0 .* E[1:N] .+ v_cross_B(v, B, dim) .+ [1.0, 0])
end

v_cross_B(v, B, ::Dim{2}) = [v[2] * B[3], -v[1] * B[3]]

v_cross_B(v, B, ::Dim{3}) = cross(v, B)

function update_j!(fields::Fields{N}, particles::Particles{N}) where {N}
    for n = 1:N
        sample_current_density!(fields.j[n], fields.grid, particles, Dim{n}())
    end
end

function sample_charge_density!(
    field::Array{Float,N},
    grid::Grid{N},
    particles::Particles{N},
) where {N}
    sample_property!(field, grid, particles, fill(LowerEdges, N), ChargeDensity)
end

function sample_current_density!(
    field::Array{Float,N},
    grid::Grid{N},
    particles::Particles{N},
    ::Dim{D},
) where {N,D}
    coordinate_references = fill(LowerEdges, N)
    coordinate_references[D] = Centers
    sample_property!(
        field,
        grid,
        particles,
        coordinate_references,
        CurrentDensity{D},
    )
end

function sample_property!(
    field::Array{Float,N},
    grid::Grid{N},
    particles::Particles{N},
    coordinate_references,
    property_type::Type{<:FieldProperty},
) where {N}
    @dbgasserteq(size(field), size(grid))

    field .= 0

    shape = size(field)
    weights = zeros(2, N)

    for p = 1:count(particles)
        position = particles.positions[:, p]
        indices = find_indices_at(grid, position, coordinate_references)
        if any((indices .< 1) .| (indices .>= shape))
            continue
        end

        weights .= 0
        for n = 1:N
            coords = coordinates(grid, coordinate_references[n])
            weights[2, n] =
                (position[n] - coords[n][indices[n]]) * cell_extents⁻¹(grid)[n]
            weights[1, n] = 1 - weights[2, n]
        end

        distribute_property!(
            field,
            compute_property(particles, grid, property_type, p),
            indices,
            weights,
        )
    end
end

function compute_property(
    particles::Particles{N},
    grid::Grid{N},
    ::Type{ChargeDensity},
    ::Integer,
) where {N}
    charge(particles) * cell_volume⁻¹(grid)
end

function compute_property(
    particles::Particles{N},
    grid::Grid{N},
    ::Type{CurrentDensity{D}},
    index::Integer,
) where {N,D}
    charge(particles) * particles.velocities[D, index] * cell_volume⁻¹(grid)
end

function distribute_property!(field::Array{Float,2}, property, indices, weights)
    for i = 1:2, j = 1:2
        field[indices[1]+i-1, indices[2]+j-1] +=
            property * weights[i, 1] * weights[j, 2]
    end
end

function distribute_property!(field::Array{Float,3}, property, indices, weights)
    for i = 1:2, j = 1:2, k = 1:2
        field[indices[1]+i-1, indices[2]+j-1, indices[3]+k-1] +=
            property * weights[i, 1] * weights[j, 2] * weights[k, 3]
    end
end

max_stars_repo_name: lars-frogner/PIC

max_issues_repo_name: lars-frogner/PIC

max_forks_repo_name: lars-frogner/PIC


text:
using EllipsisNotation
using Base.Test

include("basic.jl")
include("more_generic.jl")


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

@recipe function f(v::EmpiricalVarioplane; showrange=true, model=GaussianVariogram)
  # retrieve info
  θs = v.θs
  γs = v.γs

  # abscissa is the same for all variograms
  rs, _, __ = values(γs[1])

  # variogram values for all variograms
  zs = map(γs) do γ
    _, zs, __ = values(γ)

    # handle NaN values (i.e. empty bins)
    isnan(zs[1]) && (zs[1] = 0)
    for i in 2:length(zs)
      isnan(zs[i]) && (zs[i] = zs[i-1])
    end

    zs
  end
  Z = reduce(hcat, zs)

  # exploit symmetry
  θs = range(0, 2π, length=2*length(θs))
  Z  = [Z Z]

  # plot settings
  projection --> :polar
  ylims --> (0, Inf)

  # plot varioplane
  @series begin
    seriestype --> :heatmap
    θs, rs, Z
  end

  # plot ranges
  if showrange
    ls = [range(fit(model, γ)) for γ in γs]
    @series begin
      seriestype --> :path
      linecolor --> :black
      label --> "range"
      θs, [ls; ls]
    end
  end
end


text:
# include the package
using LatticePhysics

# obtain the unitcell of the square octagon lattice (kitaev bonds)
unitcell = getUnitcellSquareOctagon(5)

# build a lattice for visual confirmation that this is the correct unitcell
lattice = getLatticeInBox(unitcell, [12.0, 6.0])
# plot the lattice for visual confirmation (including bond coloring)
plotLattice2D(
    lattice,
    openfile=true,
    colorcode_bonds=Dict("tx"=>[255,0,0], "ty"=>[175, 215, 175], "tz"=>[175, 185, 215]),
    site_radius=-0.25,
    bond_thickness=-0.25
)


text:
include("../model/SensorModel.jl")
using Base.Test

let cartesian_point = [4 2]

    homogeneous_point = cart2hom(cartesian_point)

    @test homogeneous_point[1] ≈ 4
    @test homogeneous_point[2] ≈ 2
    @test homogeneous_point[3] ≈ 1
end

let homogeneous_point = [4 2 1]

    cartesian_point = hom2cart(homogeneous_point)

    @test cartesian_point[1] ≈ 4
    @test cartesian_point[2] ≈ 2
end

let homogeneous_point = [4 2 0.5]

    cartesian_point = hom2cart(homogeneous_point)

    @test cartesian_point[1] ≈ 8
    @test cartesian_point[2] ≈ 4
end

let homogeneous_point = [4; 2; 0.5]

    cartesian_point = hom2cart(homogeneous_point)

    @test cartesian_point[1] ≈ 8
    @test cartesian_point[2] ≈ 4
end

let homogeneous_point = [4 2 0]

    cartesian_point = hom2cart(homogeneous_point)

    @test cartesian_point[1] ≈ Inf
    @test cartesian_point[2] ≈ Inf
end

let cartesian_point = [1.5 0.5]
    rotation = π/2

    homogeneous_point = cart2hom(cartesian_point)
    rotated_point = rotate(rotation) * transpose(homogeneous_point)
    rotated_cartesian_point = hom2cart(rotated_point)

    @test rotated_cartesian_point[1] ≈ -0.5
    @test rotated_cartesian_point[2] ≈ 1.5
end

let cartesian_point = [1.5 0.5]
    translation_x = 3
    translation_y = -4

    homogeneous_point = cart2hom(cartesian_point)
    translated_point =  translate(translation_x, translation_y) * transpose(homogeneous_point)
    translated_cartesian_point = hom2cart(translated_point)

    @test translated_cartesian_point[1] ≈ 4.5
    @test translated_cartesian_point[2] ≈ -3.5
end

let cartesian_point = [1.5 0.5]
    rotation = π/2
    translation_x = 3

    homogeneous_point = cart2hom(cartesian_point)
    first_rotated_then_translated_point = translate(translation_x, 0) * rotate(rotation) * transpose(homogeneous_point)
    transformed_point = hom2cart(first_rotated_then_translated_point)

    @test transformed_point[1] ≈ 2.5
    @test transformed_point[2] ≈ 1.5
end

max_stars_repo_name: grgomrton/robot-model

max_issues_repo_name: grgomrton/robot-model

max_forks_repo_name: grgomrton/robot-model


text:
rules = Dict(r[1]=>strip(r[2], [' ', '"']) for r in split.(Iterators.takewhile(!isempty, eachline("a19.txt")), ':'))
messages = [m for m in Iterators.dropwhile(!isempty, eachline("a19.txt")) if !isempty(m)]

function rulesasregexes(rules::Dict)
  regexes = Dict()
  function resolve(n)
    if (haskey(regexes, n))
      return regexes[n]
    end
    s = rules[n];
    if !isnothing(match(r"\d+", s))
      for m in eachmatch(r"\d+", s)
        s = replace(s, m.match=>resolve(m.match); count=1)
      end
      s = "(" * s * ")"
    end
    regexes[n] = replace(s, " "=>"")
  end
  for (n,s) in rules
    resolve(n)
  end
  regexes
end

regexes = rulesasregexes(rules)

function matchescompletely(r, s)
  findfirst(r, s) == 1:length(s)
end

#part 1
println("$(count(m -> matchescompletely(Regex(regexes["0"]), m), messages))")

#part2
#assume rule 31 does not depend on 8

function backwardsregex(r)
  function parenswitch(c)
    if c == ')'
      '('
    elseif c == '('
      ')'
    else
      c
    end
  end
  join(parenswitch.(c for c in reverse(r)))
end

thirtyone = backwardsregex(regexes["31"])
fourtytwo = backwardsregex(regexes["42"])

function matchesnewrulezero(m)
  count = 0
  m = reverse(m)
  r = findfirst(Regex(thirtyone), m)
  while !isnothing(r) && first(r) == 1
    m = SubString(m, last(r)+1)
    count += 1
    r = findfirst(Regex(thirtyone), m)
  end
  count != 0 && matchescompletely(Regex(fourtytwo * "{$count}" * fourtytwo * "+"), m)
end

println("$(count(matchesnewrulezero, messages))")


text:
@testset "GolubYe - geneigsolve full" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        @testset for orth in (cgs2, mgs2, cgsr, mgsr)
            A = rand(T,(n,n)) .- one(T)/2
            A = (A+A')/2
            B = rand(T,(n,n)) .- one(T)/2
            B = sqrt(B*B')
            v = rand(T,(n,))
            alg = GolubYe(orth = orth, krylovdim = n, maxiter = 1, tol = precision(T), verbosity = 1)
            n1 = div(n,2)
            D1, V1, info = @constinferred geneigsolve((wrapop(A), wrapop(B)), wrapvec(v), n1, :SR; orth = orth, krylovdim = n, maxiter = 1, tol = precision(T), ishermitian = true, isposdef = true, verbosity = 2)
            @test KrylovKit.geneigselector((A, B), eltype(v); orth = orth, krylovdim = n, maxiter = 1, tol = precision(T), ishermitian = true, isposdef = true) isa GolubYe
            n2 = n-n1
            D2, V2, info = @constinferred geneigsolve((wrapop(A), wrapop(B)), wrapvec(v), n2, :LR, alg)
            @test vcat(D1[1:n1],reverse(D2[1:n2])) ≈ eigvals(A, B)

            U1 = hcat(unwrapvec.(V1)...)
            U2 = hcat(unwrapvec.(V2)...)
            @test U1'*B*U1 ≈ I
            @test U2'*B*U2 ≈ I

            @test A*U1 ≈ B*U1*Diagonal(D1)
            @test A*U2 ≈ B*U2*Diagonal(D2)
        end
    end
end

@testset "GolubYe - geneigsolve iteratively" begin
    @testset for T in (Float64, ComplexF64)
        @testset for orth in (cgs2, mgs2, cgsr, mgsr)
            A = rand(T,(N,N)) .- one(T)/2
            A = (A+A')/2
            B = rand(T,(N,N)) .- one(T)/2
            B = sqrt(B*B')
            v = rand(T,(N,))
            alg = GolubYe(orth = orth, krylovdim = 3*n, maxiter = 100, tol = cond(B)*precision(T))
            D1, V1, info1 = @constinferred geneigsolve((wrapop(A), wrapop(B)), wrapvec(v), n, :SR, alg)
            D2, V2, info2 = geneigsolve((wrapop(A), wrapop(B)), wrapvec(v), n, :LR, alg)

            l1 = info1.converged
            l2 = info2.converged
            @test l1 > 0
            @test l2 > 0
            @test D1[1:l1] ≊ eigvals(A, B)[1:l1]
            @test D2[1:l2] ≊ eigvals(A, B)[N:-1:N-l2+1]

            U1 = hcat(unwrapvec.(V1)...)
            U2 = hcat(unwrapvec.(V2)...)
            @test U1'*B*U1 ≈ I
            @test U2'*B*U2 ≈ I

            R1 = hcat(unwrapvec.(info1.residual)...)
            R2 = hcat(unwrapvec.(info2.residual)...)
            @test A*U1 ≈ B*U1*Diagonal(D1) + R1
            @test A*U2 ≈ B*U2*Diagonal(D2) + R2
        end
    end
end


text:
using DelimitedFiles

function read_arrays_jl_c(filename)
  fmt = "%ld%c"
  c = Ref{Cchar}(0)
  x = Ref{Int}(0)
  n = 0
  set = UInt8[',', '#']
  f = ccall(:fopen, Ptr{Cvoid}, (Cstring, Cstring), filename, "r")
  while ccall(:fscanf, Cint, (Ptr{Cvoid}, Cstring, Ref{Int}, Ref{Cchar}), f, fmt, x, c) == 2
    if c[] in set
      n += 1
    end
  end
  ccall(:fclose, Ptr{Cvoid}, (Ptr{Cvoid},), f)

  keys = zeros(Int, n)
  indexes = zeros(Int, n)
  elements = zeros(Int, n)

  count = 1
  k = -1
  j = 0
  f = ccall(:fopen, Ptr{Cvoid}, (Cstring, Cstring), filename, "r")
  while ccall(:fscanf, Cint, (Ptr{Cvoid}, Cstring, Ref{Int}, Ref{Cchar}), f, fmt, x, c) == 2
    if c[] == Int8('#')
      k = x[]
      j = 0
    else
      keys[count] = k
      indexes[count] = j
      j += 1
      elements[count] = x[]
      count += 1
    end
  end
  ccall(:fclose, Ptr{Cvoid}, (Ptr{Cvoid},), f)

  return keys, indexes, elements
end


text:
abstract type Branch <: Device end

abstract type ACBranch <: Branch end

abstract type DCBranch <: Branch end

get_from_bus(b::T) where {T <: Branch} = b.arc.from
get_to_bus(b::T) where {T <: Branch} = b.arc.to


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

JLLWrappers.@generate_main_file_header("CFITSIO")
JLLWrappers.@generate_main_file("CFITSIO", UUID("b3e40c51-02ae-5482-8a39-3ace5868dcf4"))
end  # module CFITSIO_jll


text:
"""
Created in February, 2022 by
[chifi - an open source software dynasty.](https://github.com/orgs/ChifiSource)
by team
[toolips](https://github.com/orgs/ChifiSource/teams/toolips)
This software is MIT-licensed.
### Toolips
**Toolips.jl** is a **fast**, **asynchronous**, **low-memory**, **full-stack**,
and **reactive** web-development framework **always** written in **pure** Julia.
##### Module Composition
- [**Toolips**](https://github.com/ChifiSource/Toolips.jl)
"""
module Toolips
using Crayons
using Sockets, HTTP, Pkg, ParseNotEval, Dates
import Base: getindex, setindex!, push!, get, string
#==
SuperTypes
==#
"""
### abstract type Servable
Servables can be written to a Connection via thier f() function and the
interface. They can also be indexed with strings or symbols to change properties
##### Consistencies
- f::Function - Function whose output to be written to http().
"""
abstract type Servable <: Any end

"""
### abstract type StyleComponent <: Servable
No different from a normal Servable, simply an abstract type step for the
interface to separate working with Animations and Styles.
### Servable Consistencies
```
Servables can be written to a Connection via thier f() function and the
interface. They can also be indexed with strings or symbols to change properties
##### Consistencies
- f::Function - Function whose output to be written to http().
```
"""
abstract type StyleComponent <: Servable end

"""
### abstract type ToolipsServer
ToolipsServers are returned whenever the ServerTemplate.start() field is
called. If you are running your server as a module, it should be noted that
commonly a global start() method is used and returns this server, and dev is
where this module is loaded, served, and revised.
##### Consistencies
- routes::Dict - The server's route => function dictionary.
- extensions::Dict - The server's currently loaded extensions.
- server::Any - The server, whatever type it may be...
"""
abstract type ToolipsServer end

"""
### abstract type ServerExtension
Server extensions are loaded into the server on startup, and
can have a few different abilities according to their type
field's value. This value can be either a Symbol or a Vector of Symbols.
##### Consistencies
- type::T where T == Vector{Symbol}  || T == Symbol
"""
abstract type ServerExtension end

"""
### Connection
- routes::Dict
- http::HTTP.Stream
- extensions::Dict
The connection type is passed into route functions and pages as an argument.
This is both for functions, as well as Servable.f() methods. This constructor
    should not be called directly. Instead, it is called by the server and
    passed through the function pipeline. Indexing a Connection will return
        the extension named with that symbol.
##### example
```
                  #  v The Connection
home = route("/") do c::Connection
    c[Logger].log(1, "We can index extensions by type or symbol")
    c[:logger].log(1, "see?")
    c.routes["/"] = c::Connection -> write!(c, "rerouting!")
    httpstream = c.http
    write!(c, "Hello world!")
    myheading::Component = h("myheading", 1, text = "Whoa!")
    write!(c, myheading)
end
```
------------------
##### field info
- **routes::Dict** - A dictionary of routes where the keys
are the routed URL and the values are the functions to
those keys.
- **http::HTTP.Stream** - The stream for this current peer's connection.
- **extensions::Dict** - A dictionary of extensions to load with the
name to reference as keys and the extension as the pair.
------------------
##### constructors
- Connection(routes::Dict, http::HTTP.Stream, extensions::Dict)
"""
mutable struct Connection
    routes::Dict
    http::HTTP.Stream
    extensions::Dict
    function Connection(routes::Dict, http::HTTP.Stream, extensions::Dict)
        new(routes, http, extensions)::Connection
    end
end

#==
Includes/Exports
==#
include("interface/Servables.jl")
include("server/Core.jl")
include("interface/Interface.jl")

# Core Server
export ServerTemplate, Route, Connection, WebServer
# Server Extensions
export Logger, Files, Document
# Servables
export File, Component
export Animation, Style

export img, link, meta, input, a, p, h, button, ul, li, divider, form, br, i
export title, span, iframe, svg, element, label, script, nav, button, form
export element, label, script, nav, button, form
# High-level api
export push!, getindex, setindex!, properties!, components
export animate!, style!, keyframe!, delete_keyframe!, @keyframe!
export route, routes, route!, write!, stop!, unroute!, navigate!, stop!
export has_extension
export getargs, getarg, postargs, postarg, get, post, getip, getpost

#==
Project API
==#
"""
### create_serverdeps(name::String, inc::String) -> _
------------------
Creates the essential portions of the webapp file structure, where name is the
project's name and inc is any extensions or strings to incorporate at the top
of the file.
#### example
```
create_serverdeps("ToolipsApp")
```
"""
function create_serverdeps(name::String, inc::String = "")
    Pkg.generate(name)
    Pkg.activate(name)
    Pkg.add(url = "https://github.com/ChifiSource/Toolips.jl.git")
    Pkg.add("Revise")
    dir = pwd() * "/"
    src = dir * name * "/src"
    logs = dir * name * "/logs"
    mkdir(logs)
    touch(name * "/dev.jl")
    touch(name * "/prod.jl")
    touch(logs * "/log.txt")
    rm(src * "/$name.jl")
    touch(src * "/$name.jl")
    open(src * "/$name.jl", "w") do io
        write(io, """
module $name
using Toolips
$inc

function home(c::Connection)
    write!(c, p("helloworld", text = "hello world!"))
end

fourofour = route("404") do c
    write!(c, p("404message", text = "404, not found!"))
end

\"\"\"
start()
\"\"\"
function start(IP::String = "127.0.0.1", PORT::Integer = 8000,
    extensions::Dict = Dict(:logger => Logger()))
    rs = routes(route("/", home), fourofour)
    server = ServerTemplate(IP, PORT, rs, extensions = extensions)
    server.start()
end

end # - module
        """)
    end

end

"""
### new_app(::String) -> _
------------------
Creates a minimalistic app, usually used for creating endpoints -- but can
be used for anything. For an app with a real front-end, it might make sense to
add some extensions.
#### example
```
using Toolips
Toolips.new_app("ToolipsApp")
```
"""
function new_app(name::String = "ToolipsApp")
    create_serverdeps(name)
    servername = name * "Server"
    open(name * "/dev.jl", "w") do io
        write(io, """
        #==
        dev.jl is an environment file. This file loads and starts servers, and
        defines environmental variables, setting the scope a lexical step higher
        with modularity.
        ==#
        using Pkg; Pkg.activate(".")
        using Toolips
        using Revise

        IP = "127.0.0.1"
        PORT = 8000
        extensions = Dict(:logger => Logger())
        using $name
        $servername = $name.start(IP, PORT, extensions)
        """)
    end
    open(name * "/prod.jl", "w") do io
        write(io, """
        using Pkg; Pkg.activate(".")
        using Toolips

        IP = "127.0.0.1"
        PORT = 8000
        extensions = Dict(:logger => Logger())
        using $name
        $servername = $name.start(IP, PORT, extensions)
        """)
    end
end

"""
### new_webapp(::String) -> _
------------------
Creates a fully-featured web-app. Adds ToolipsModifier, ideal for full-stack
web-sites.
#### example
```
using Toolips
Toolips.new_webapp("ToolipsApp")
```
"""
function new_webapp(name::String = "ToolipsApp")
    servername = name * "Server"
    create_serverdeps(name, "using ToolipsModifier")
    open(name * "/dev.jl", "w") do io
        write(io, """
        #==
        dev.jl is an environment file. This file loads and starts servers, and
        defines environmental variables, setting the scope a lexical step higher
        with modularity.
        ==#
        using Pkg; Pkg.activate(".")
        using Toolips
        using Revise
        using ToolipsModifier

        IP = "127.0.0.1"
        PORT = 8000
        #==
        Extension description
        :logger -> Logs messages into both a file folder and the terminal.
        :public -> Routes the files from the public directory.
        :mod -> ToolipsModifier; allows us to make Servables reactive. See ?(on)
        ==#
        extensions = Dict(:logger => Logger(), :public => Files("public"),
        :mod => Modifier())
        using $name
        $servername = $name.start(IP, PORT, extensions)
        """)
    end
    open(name * "/prod.jl", "w") do io
        write(io, """
        using Pkg; Pkg.activate(".")
        using Toolips
        using ToolipsModifier

        IP = "127.0.0.1"
        PORT = 8000
        extensions = Dict(:logger => Logger(), :public => Files("public"))
        using $name
        $servername = $name.start(IP, PORT, extensions)
        """)
    end
    public = pwd() * "/$name/public"
    mkdir(public)
end
# --

end


text:
using Optim

# rosenbrock(x)=(1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
# o = optimize(rosenbrock, [0.0,0.0])
# println(o)

ndims = 10
A = randn(ndims,ndims)
H = A'*A
b = randn(ndims)
g(x)=sum(0.5*x'*H*x+x'*b)
println(optimize(g, zeros(ndims)))
# gradient: Hx+b
# hessian: H
# minimum: -H^-1 b
-H\b


text:
function advanced_input_ui()

    # Title of section
    title = node(:div, tachyons_css, "Advanced Mode Options") |> class"f4 lh-title"

    # Sources
    o1 = Observable("")
    sources = vbox(node(:div, "Current source file: ",
                      attributes = Dict(:style => "margin-top: 12px")) |> class"b",
                   filepicker(value = o1))

    # Grounds
    o2 = Observable("")
    grounds = vbox(node(:div, "Ground point file: ",
                      attributes = Dict(:style => "margin-top: 12px")) |> class"b",
                   filepicker(value = o2))

    adv = vbox(title,
             sources,
             grounds)

    # Define scope
    s = Scope()
    s.dom = adv
    @private s["source"] = o1
    @private s["ground"] = o2

    s
end


text:
module GNNGraphs

using SparseArrays
using Functors: @functor
using CUDA
import Graphs
using Graphs: AbstractGraph, outneighbors, inneighbors, adjacency_matrix, degree, has_self_loops, is_directed
import Flux
using Flux: batch
import NearestNeighbors
import NNlib
import LearnBase
import StatsBase
import KrylovKit
using ChainRulesCore
using LinearAlgebra, Random

include("gnngraph.jl")
export GNNGraph,
       node_features,
       edge_features,
       graph_features

include("query.jl")
export adjacency_list,
       edge_index,
       get_edge_weight,
       graph_indicator,
       has_multi_edges,
       is_directed,
       is_bidirected,
       normalized_laplacian,
       scaled_laplacian,
       # from Graphs
       adjacency_matrix,
       degree,
       has_self_loops,
       inneighbors,
       outneighbors

include("transform.jl")
export add_nodes,
       add_edges,
       add_self_loops,
       getgraph,
       negative_sample,
       rand_edge_split,
       remove_self_loops,
       remove_multi_edges,
       # from Flux
       batch,
       unbatch,
       # from SparseArrays
       blockdiag

include("generate.jl")
export rand_graph,
       knn_graph

include("sampling.jl")
export sample_neighbors

include("operators.jl")
# Base.intersect

include("convert.jl")
include("utils.jl")

end #module


text:
#=
    arguments to renormalize must be in nonincreasing magnitude order
        and adjacent argument magnitudes (abs(value)) must be separated by
        ?4 significant bits?
=#

renormalize{T<:StdFloat}(a0::T,a1::T) = eftSumGTE(a0,a1)
renormalizeAs1{T<:StdFloat}(a0::T,a1::T) = a0+a1

function renormalize{T<:StdFloat}(a0::T,a1::T,a2::T)
   s,t1 = eftSumGTE(a1,a3)
   b0,t0 = eftSumGTE(a0,s)

   b1,t0 = eftSumGTE(t0,s)
   b2 = t0+t1

   b0,b1,b2
end

function renormalizeAs2{T<:StdFloat}(a0::T,a1::T,a2::T)
   s,t1 = eftSumGTE(a1,a3)
   b0,t0 = eftSumGTE(a0,s)

   b1 = t0+s

   b0,b1
end

function renormalizeAs1{T<:StdFloat}(a0::T,a1::T,a2::T)
   b0  = a1+a3
   b0 += a0

   b0
end

#=
   follows Algorithm A.9
   in Mixed Precision Iterative Methods using High Precision Arithmetic
   by Hidehiko Hasegawa, 2013
=#

function renormalize{T<:StdFloat}(a0::T,a1::T,a2::T,a3::T)
   s,t2 = eftSumGTE(a2,a3)
   s,t1 = eftSumGTE(a1,s)
   b0,t0 = eftSumGTE(a0,s)

   s,t1 = eftSumGTE(t1,t2)
   b1,t0 = eftSumGTE(t0,s)
   b2,b3 = eftSumGTE(t0,t1)

   b0,b1,b2,b3
end

function renormalizeAs3{T<:StdFloat}(a0::T,a1::T,a2::T,a3::T)
   s,t2 = eftSumGTE(a2,a3)
   s,t1 = eftSumGTE(a1,s)
   b0,t0 = eftSumGTE(a0,s)

   s,t1 = eftSumGTE(t1,t2)
   b1,t0 = eftSumGTE(t0,s)
   b2 = t0+t1

   b0,b1,b2
end

function renormalizeAs2{T<:StdFloat}(a0::T,a1::T,a2::T,a3::T)
   s,t2 = eftSumGTE(a2,a3)
   s,t1 = eftSumGTE(a1,s)
   b0,t0 = eftSumGTE(a0,s)

   s = t1+t2
   b1 = t0+s

   b0,b1
end

function renormalizeAs1{T<:StdFloat}(a0::T,a1::T,a2::T,a3::T)
   b0 = a2+a3
   b0 += a1
   b0 += a0

   b0
end

#=
   follows Algorithm A.8
   in Mixed Precision Iterative Methods using High Precision Arithmetic
   by Hidehiko Hasegawa, 2013
=#

function renormalizeAs4{T<:StdFloat}(a0::T,a1::T,a2::T,a3::T,a4::T)
   s,t3 = eftSumGTE(a3,a4)
   s,t2 = eftSumGTE(a2,s)
   s,t1 = eftSumGTE(a1,s)
   b0,t0 = eftSumGTE(a0,s)

   s,t2 = eftSumGTE(t2,t3)
   s,t1 = eftSumGTE(t1,s)
   b1,t0 = eftSumGTE(t0,s)
   s,t1 = eftSumGTE(t1,t2)
   b2,t0 = eftSumGTE(t0,s)
   b3 = t0+t1

   b0,b1,b2,b3
end

function renormalizeAs3{T<:StdFloat}(a0::T,a1::T,a2::T,a3::T,a4::T)
   s,t3 = eftSumGTE(a3,a4)
   s,t2 = eftSumGTE(a2,s)
   s,t1 = eftSumGTE(a1,s)
   b0,t0 = eftSumGTE(a0,s)

   s,t2 = eftSumGTE(t2,t3)
   s,t1 = eftSumGTE(t1,s)
   b1,t0 = eftSumGTE(t0,s)
   s,t1 = eftSumGTE(t1,t2)
   b2 = t0+s

   b0,b1,b2
end

function renormalizeAs2{T<:StdFloat}(a0::T,a1::T,a2::T,a3::T,a4::T)
   s,t3 = eftSumGTE(a3,a4)
   s,t2 = eftSumGTE(a2,s)
   s,t1 = eftSumGTE(a1,s)
   b0,t0 = eftSumGTE(a0,s)

   s = t2+t3
   s = t1+s
   b1 = t0+s

   b0,b1
end


text:
module LazyArrays

# Use README as the docstring of the module:
@doc read(joinpath(dirname(@__DIR__), "README.md"), String) LazyArrays

using Base, Base.Broadcast, LinearAlgebra, FillArrays, StaticArrays, ArrayLayouts, MatrixFactorizations, SparseArrays
import LinearAlgebra.BLAS

import Base: AbstractArray, AbstractMatrix, AbstractVector,
        ReinterpretArray, ReshapedArray, AbstractCartesianIndex, Slice,
             RangeIndex, BroadcastStyle, copyto!, length, broadcastable, axes,
             getindex, eltype, tail, IndexStyle, IndexLinear, getproperty,
             *, +, -, /, \, ==, isinf, isfinite, sign, angle, show, isless,
         fld, cld, div, min, max, minimum, maximum, mod,
         <, ≤, >, ≥, promote_rule, convert, copy,
         size, step, isempty, length, first, last, ndims,
         getindex, setindex!, setindex, intersect, @_inline_meta, inv,
         sort, sort!, issorted, sortperm, diff, accumulate, cumsum, sum, in, broadcast,
         eltype, parent, real, imag,
         conj, transpose, adjoint, permutedims, vec,
         exp, log, sqrt, cos, sin, tan, csc, sec, cot,
                   cosh, sinh, tanh, csch, sech, coth,
                   acos, asin, atan, acsc, asec, acot,
                   acosh, asinh, atanh, acsch, asech, acoth, (:),
         AbstractMatrix, AbstractArray, checkindex, unsafe_length, OneTo, one, zero,
        to_shape, _sub2ind, print_matrix, print_matrix_row, print_matrix_vdots,
      checkindex, Slice, @propagate_inbounds, @_propagate_inbounds_meta,
      _in_range, _range, _rangestyle, Ordered,
      ArithmeticWraps, floatrange, reverse, unitrange_last,
      AbstractArray, AbstractVector, axes, (:), _sub2ind_recurse, broadcast, promote_eltypeof,
      similar, @_gc_preserve_end, @_gc_preserve_begin,
      @nexprs, @ncall, @ntuple, tuple_type_tail,
      all, any, isbitsunion, issubset, replace_with_centered_mark, replace_in_print_matrix,
      unsafe_convert, strides, union, map, searchsortedfirst, searchsortedlast, searchsorted

import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, Broadcasted, broadcasted,
                        combine_eltypes, DefaultArrayStyle, instantiate, materialize,
                        materialize!, eltypes

import LinearAlgebra: AbstractTriangular, AbstractQ, checksquare, pinv, fill!, tilebufsize, Abuf, Bbuf, Cbuf, dot, factorize, qr, lu, cholesky,
                        norm2, norm1, normInf, normp, normMinusInf, diag, det, logabsdet, tr, AdjOrTrans, triu, tril,
                        lmul!, rmul!

import LinearAlgebra.BLAS: BlasFloat, BlasReal, BlasComplex

import FillArrays: AbstractFill, getindex_value

import StaticArrays: StaticArrayStyle

import ArrayLayouts: MatMulVecAdd, MatMulMatAdd, MulAdd, Lmul, Rmul, Ldiv, Dot, Mul, _inv,
                        transposelayout, conjlayout, sublayout, triangularlayout, triangulardata,
                        reshapedlayout, diagonallayout, tridiagonallayout, symtridiagonallayout, bidiagonallayout, symmetriclayout, hermitianlayout,
                        adjointlayout, sub_materialize, mulreduce,
                        check_mul_axes, _mul_eltype, check_ldiv_axes, ldivaxes, colsupport, rowsupport,
                        _fill_lmul!, scalarone, scalarzero, fillzeros, zero!, layout_getindex, _copyto!,
                        AbstractQLayout, StridedLayout, layout_replace_in_print_matrix

import Base: require_one_based_indexing, oneto

export Mul, Applied, MulArray, MulVector, MulMatrix, InvMatrix, PInvMatrix,
        Hcat, Vcat, Kron, BroadcastArray, BroadcastMatrix, BroadcastVector, cache, Ldiv, Inv, PInv, Diff, Cumsum, Accumulate,
        applied, materialize, materialize!, ApplyArray, ApplyMatrix, ApplyVector, apply, @~, LazyArray,
        PaddedArray, PaddedVector, PaddedMatrix

include("lazyapplying.jl")
include("lazybroadcasting.jl")
include("linalg/linalg.jl")
include("cache.jl")
include("lazyconcat.jl")
include("lazysetoperations.jl")
include("lazyoperations.jl")
include("lazymacro.jl")

# support x^2
Base.to_power_type(x::LazyArray) = x

# Special broadcasting for BlockArrays.jl
map(::typeof(length), A::BroadcastArray{OneTo{Int},1,Type{OneTo}}) = A.args[1]
map(::typeof(length), A::BroadcastArray{<:Fill,1,Type{Fill},<:NTuple{2,Any}}) = A.args[2]
map(::typeof(length), A::BroadcastArray{<:Zeros,1,Type{Zeros}}) = A.args[1]
map(::typeof(length), A::BroadcastArray{<:Vcat,1,Type{Vcat}}) = broadcast(+,map.(length,A.args)...)
broadcasted(::LazyArrayStyle{1}, ::typeof(length), A::BroadcastArray{OneTo{Int},1,Type{OneTo}}) = A.args[1]
broadcasted(::LazyArrayStyle{1}, ::typeof(length), A::BroadcastArray{<:Fill,1,Type{Fill},<:NTuple{2,Any}}) = A.args[2]

end # module


text:
## Platforms for certain bases and frames

export FourierPlatform
"""
    struct FourierPlatform{T} <: BasisPlatform

A platform of Fourier series on [0,1]
"""
struct FourierPlatform{T} <: BasisPlatform
end

correctparamformat(::FourierPlatform, ::Int) = true
correctparamformat(::FourierPlatform, _) = false

FourierPlatform() = FourierPlatform{Float64}()

unsafe_dictionary(p::FourierPlatform{T}, n) where {T} = Fourier{T}(n)

SolverStyle(p::FourierPlatform, ::OversamplingStyle) = DualStyle()

measure(p::FourierPlatform{T}) where {T} = FourierWeight{T}()

export ChebyshevPlatform
"""
    struct ChebyshevPlatform{T} <: BasisPlatform

A platform of ChebyshevT dictionaries on [-1,1]
"""
struct ChebyshevPlatform{T} <: BasisPlatform
end

correctparamformat(::ChebyshevPlatform, ::Int) = true
correctparamformat(::ChebyshevPlatform, _) = false

ChebyshevPlatform() = ChebyshevPlatform{Float64}()

unsafe_dictionary(p::ChebyshevPlatform{T}, n) where {T} = ChebyshevT{T}(n)

measure(platform::ChebyshevPlatform{T}) where {T} = ChebyshevWeight{T}()

export FourierExtensionPlatform
"""
    struct FourierExtensionPlatform <: FramePlatform

A platform of fourier extensions with [0,1] as bouding box.
"""
struct FourierExtensionPlatform <: FramePlatform
    basisplatform   ::  FourierPlatform
    domain          ::  Domain

    function FourierExtensionPlatform(basisplatform, domain::Domain{T}) where {T}
        # Only supported for intervals
        # @assert issubset(domain, UnitInterval{T}())
        new(basisplatform, domain)
    end
end

correctparamformat(platform::FourierExtensionPlatform, param) =
    correctparamformat(platform.basisplatform, param)

function FourierExtensionPlatform(domain::Domain{T}) where {T}
    basisplatform = FourierPlatform{T}()
    FourierExtensionPlatform(basisplatform, domain)
end

unsafe_dictionary(p::FourierExtensionPlatform, n) =
    extensionframe(p.domain, unsafe_dictionary(p.basisplatform, n))

dualdictionary(platform::FourierExtensionPlatform, param, measure::Measure; options...) =
   extensionframe(dualdictionary(platform.basisplatform, param, supermeasure(measure); options...), platform.domain)

measure(platform::FourierExtensionPlatform) = measure(dictionary(platform, 1))

# discrete_gram_measure(p::FourierExtensionPlatform, n, L;
#             dict = dictionary(p, n), normalizedsampling = false, options...) =
#     normalizedsampling ?
#     restrict(BasisFunctions.NormalizedDiracComb(oversampling_grid(dictionary(p.basisplatform, n), L)), p.domain) :
#     restrict(discretemeasure(oversampling_grid(dictionary(p.basisplatform, n), L)), p.domain)

export OPSExtensionFramePlatform
"""
    struct OPSExtensionFramePlatform <: FramePlatform

A platform of classical orthogonal polynomial extensions. The bouding box is [-1,1]
"""
struct OPSExtensionFramePlatform <: FramePlatform
    basisplatform   ::Platform
    domain          ::Domain
    function OPSExtensionFramePlatform(dict::BasisFunctions.OPS, domain::Domain1d)
        new(platform(dict),domain)
    end
end

correctparamformat(p::OPSExtensionFramePlatform, param) =
    correctparamformat(p.basisplatform, param)

unsafe_dictionary(p::OPSExtensionFramePlatform, n) =
    ExtensionFrame(p.domain, unsafe_dictionary(p.basisplatform, n))

measure(platform::OPSExtensionFramePlatform) = restrict(measure(platform.basisplatform), platform.domain)

dualdictionary(platform::OPSExtensionFramePlatform, param, measure::Measure; options...) =
   extensionframe(dualdictionary(platform.basisplatform, param, supermeasure(measure); options...), platform.domain)

discretemeasure(::DiscreteStyle, p::OPSExtensionFramePlatform, n, L;
            dict = dictionary(p, n), options...) =
    _opsnodesmeasure(resize(dict,L))
_opsnodesmeasure(frame::ExtensionFrame) = restrict(gauss_rule(basis(frame)), support(frame))
_opsnodesmeasure(ops::BasisFunctions.OrthogonalPolynomials) = gauss_rule(ops)