text:
using Printf
using LinearAlgebra
using Random

include("../common/PWGrid_v01.jl")
include("../common/ortho_gram_schmidt.jl")
include("../common/wrappers_fft.jl")
include("../common/gen_lattice_pwscf.jl")

include("op_K.jl")
include("op_Vpot.jl")
include("op_H.jl")
include("calc_rho.jl")
include("calc_Etot.jl")
include("Kprec.jl")
include("diag_lobpcg.jl")
include("calc_grad.jl")
include("Sch_solve_Emin_sd.jl")
include("Sch_solve_Emin_cg.jl")
include("structure_factor.jl")

include("gen_rho.jl")
include("gen_dr.jl")
include("calc_ewald.jl")

function test_main( ns1::Int64,ns2::Int64,ns3::Int64 )

    Ns = [ns1,ns2,ns3]
    LatVecs = gen_lattice_sc(16.0)

    pw = PWGrid( Ns, LatVecs )

    Npoints = pw.Npoints
    Ω = pw.Ω
    r = pw.r
    G = pw.G
    G2 = pw.G2

    Xpos = reshape( [0.0, 0.0, 0.0], (3,1) )

    Sf = structure_factor( Xpos, G )

    E_nn = calc_ewald( pw, Xpos, Sf )

    Vg = zeros(ComplexF64,Npoints)
    prefactor = -4*pi/Ω
    for ig=2:Npoints
        Vg[ig] = prefactor/G2[ig]
    end
    Vpot = real( G_to_R(Ns, Vg .* Sf) ) * Npoints

    println("sum(Sf) = ", sum(Sf))
    println("sum(Vg) = ", sum(Vg))
    println("sum(Vpot) = ", sum(Vpot))
    for ip = 1:10
        @printf("%8d %18.10f\n", ip, Vpot[ip])
    end
    @printf("Ω = %f\n", Ω)
    @printf("maximum(Vpot) = %18.10f\n", maximum(Vpot))
    @printf("minimum(Vpot) = %18.10f\n", minimum(Vpot))

    #
    Nstates = 1
    Random.seed!(2222)
    psi = randn(Npoints,Nstates) + im*randn(Npoints,Nstates)
    psi = ortho_gram_schmidt(psi)
    #
    #evals, evecs = diag_lobpcg( pw, Vpot, psi, verbose=true, tol_avg=1e-10 )

    psi, Etot = Sch_solve_Emin_sd( pw, Vpot, psi, NiterMax=10 )
    psi, Etot = Sch_solve_Emin_cg( pw, Vpot, psi, NiterMax=1000 )
    #
    Y = ortho_gram_schmidt(psi)
    mu = Y' * op_H( pw, Vpot, Y )
    evals, evecs = eigen(mu)
    psi = Y*evecs

    for ist = 1:Nstates
        @printf("State # %d, Energy = %18.10f\n", ist, real(evals[ist]))
    end

    @printf("E_nn    = %18.10f\n", E_nn)
    @printf("E total = %18.10f\n", E_nn + Etot)

end

@time test_main( 64,64,64 )
@time test_main( 64,64,64 )


text:
#=
Model: Motor (8 variables, 2 inputs)
=#
using Reachability, SparseArrays

motor(o::Pair{Symbol, <:Any}...) = motor(Options(Dict{Symbol, Any}(o)))

function motor(input_options::Options)
    # =====================
    # Problem specification
    # =====================
    I = [1, 2, 2, 3, 3, 3, 3, 4, 5, 6, 6, 7, 7, 7, 7, 8]
    J = [2, 3, 2, 1, 2, 3, 4, 1, 6, 7, 6, 5, 6, 7, 8, 5]
    vals = [1, 8487.2, -1.0865, -2592.1, -21.119, -698.91, -141399.0, 1.0, 1.0,
           8487.2, -1.0865, -2592.1, -21.119, -698.91, -141399.0, 1.0]
    A = sparse(I, J, vals)

    # initial set
    X0 = Hyperrectangle([0.00225, 0.0, 0.0, 0.0, 0.00125, 0.0, 0.0, 0.0],
                        [0.00025, 0.0, 0.0, 0.0, 0.00025, 0.0, 0.0, 0.0])

    # input set
    B = sparse([4, 8], [1, 2], [-1.0, -1.0])
    U = Hyperrectangle([0.23, 0.3], [0.07, 0.1])

    # instantiate continuous LTI system
    S = InitialValueProblem(
        ConstrainedLinearControlContinuousSystem(A, B, nothing, U), X0)

    # property: x1 < 0.35 || x5 < 0.45
    property = Disjunction([
        SafeStatesProperty(HalfSpace([1.; zeros(7)], 0.35)),
        SafeStatesProperty(HalfSpace([zeros(4); 1.; zeros(3)], 0.45))])

    # =======================
    # Problem default options
    # =======================
    partition = [(2*i-1:2*i) for i in 1:4] # 2D blocks

    if input_options[:mode] == "reach"
        problem_options = Options(:vars => [5],
                                  :partition => partition,
                                  :plot_vars => [0, 5])
    elseif input_options[:mode] == "check"
        problem_options = Options(:vars => [1, 5],
                                  :partition => partition,
                                  :property => property)
    end

    return (S, merge(problem_options, input_options))
end

ext: jl

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
# ------------------------------------------------------------------
# Copyright (c) 2017, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

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

A polyalgorithm Kriging estimation solver.

Each pair `var=>param` specifies the `KrigingParam` `param`
for the Kriging variable `var`. In order to avoid boilerplate
code, the constructor expects pairs of `Symbol` and `NamedTuple`
instead.

## Parameters

* `variogram` - Variogram model (default to `GaussianVariogram()`)
* `mean`      - Simple Kriging mean
* `degree`    - Universal Kriging degree
* `drifts`    - External Drift Kriging drift functions

Latter options override former options. For example, by specifying
`drifts`, the user is telling the algorithm to ignore `degree` and
`mean`. If no option is specified, Ordinary Kriging is used by
default with the `variogram` only.

## Examples

Solve the variable `:var₁` with Simple Kriging by specifying
the `mean`, and the variable `:var₂` with Universal Kriging
by specifying the `degree` and the `variogram` model.

```julia
julia> Kriging(
  :var₁ => (mean=1.,),
  :var₂ => (degree=1, variogram=SphericalVariogram(range=20.))
)
```

Solve all variables of the problem with the default parameters
(i.e. Ordinary Kriging with unit Gaussian variogram):

```julia
julia> Kriging()
```
"""
@estimsolver Kriging begin
  @param variogram = GaussianVariogram()
  @param mean = nothing
  @param degree = nothing
  @param drifts = nothing
end

function solve(problem::EstimationProblem, solver::Kriging)
  # sanity checks
  @assert keys(solver.params) ⊆ keys(variables(problem)) "invalid variable names in solver parameters"

  # determine problem coordinate type
  T = coordtype(data(problem))

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

  # loop over target variables
  for (var,V) in variables(problem)
    # get user parameters
    if var ∈ keys(solver.params)
      varparams = solver.params[var]
    else
      varparams = KrigingParam()
    end

    # determine which Kriging variant to use
    if varparams.drifts ≠ nothing
      estimator = ExternalDriftKriging{T,V}(varparams.variogram, varparams.drifts)
    elseif varparams.degree ≠ nothing
      estimator = UniversalKriging{T,V}(varparams.variogram, varparams.degree)
    elseif varparams.mean ≠ nothing
      estimator = SimpleKriging{T,V}(varparams.variogram, varparams.mean)
    else
      estimator = OrdinaryKriging{T,V}(varparams.variogram)
    end

    # perform estimation
    varμ, varσ = solve(problem, var, estimator)

    # save result for variable
    push!(μs, var => varμ)
    push!(σs, var => varσ)
  end

  EstimationSolution(domain(problem), Dict(μs), Dict(σs))
end

function solve(problem::EstimationProblem, var::Symbol, estimator::E) where {E<:KrigingEstimator}
  # retrieve problem info
  pdata = data(problem)
  pdomain = domain(problem)

  # find valid data for variable
  X, z = valid(pdata, var)

  # fit estimator to data
  fit!(estimator, X, z)

  # pre-allocate memory for result
  varμ = Vector{eltype(z)}(undef, npoints(pdomain))
  varσ = Vector{eltype(z)}(undef, npoints(pdomain))

  # pre-allocate memory for coordinates
  coords = MVector{ndims(pdomain),coordtype(pdomain)}(undef)

  # estimation loop
  for location in SimplePath(pdomain)
    coordinates!(coords, pdomain, location)

    μ, σ² = estimate(estimator, coords)

    varμ[location] = μ
    varσ[location] = σ²
  end

  # return mean and variance
  varμ, varσ
end


text:
using FractionalSystems
using Documenter

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

makedocs(;
    modules=[FractionalSystems],
    authors="Qingyu Qu",
    repo="https://github.com/SciFracX/FractionalSystems.jl/blob/{commit}{path}#{line}",
    sitename="FractionalSystems.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://SciFracX.github.io/FractionalSystems.jl",
        assets = ["assets/favicon.ico"],
    ),
    pages = [
        "FractionalSystems.jl" => "index.md",
        "Get Started" => "get_started.md",
        "Root Locus" => "rlocus.md",
        "Frequency Analysis" => "frequency.md",
        "Lyapunov Exponents Spectrum" => "flyapunovexp.md",
        "Conversation" => "conversation.md"
    ],
)

deploydocs(;
    repo="github.com/SciFracX/FractionalSystems.jl",
    devbranch="master",
)


text:
module MixedTests

using Interpolations, Base.Test

A = rand(6,5)
knots = (collect(linspace(1,size(A,1),size(A,1))),collect(linspace(1,size(A,2),size(A,2))))
itp = @inferred(interpolate(knots, A, (Gridded(Linear()),NoInterp())))
@inferred(getindex(itp, 2, 2))
@inferred(getindex(itp, CartesianIndex((2,2))))
for j = 2:size(A,2)-1, i = 2:size(A,1)-1
    @test_approx_eq itp[i,j] A[i,j]
    @test_approx_eq itp[CartesianIndex((i,j))] A[i,j]
end
@test_approx_eq itp[knots...] A
@inferred(getindex(itp, knots...))

A = rand(8,20)
knots = ([x^2 for x = 1:8], [0.2y for y = 1:20])
itp = interpolate(knots, A, Gridded(Linear()))
@test_approx_eq itp[4,1.2] A[2,6]

@test_throws  ErrorException interpolate(knots, A, (Gridded(Linear()),NoInterp()))
end


text:
function _gradient_type(
  ::Type{A},::Type{<:SVector{D}}) where A<:StaticArray{S,T,N} where {S,T,N,D}
  SG = _gradient_size(Size(A),Val(D))
  TG = T
  NG = N+1
  LG = _gradient_length(Size(A),Val(D))
  SArray{SG,TG,NG,LG}
end

function _gradient_type(
  ::Type{A},::Type{<:VectorValue{D}}) where A<:MultiValue{S,T,N,L} where {S,T,N,L,D}
  G = _gradient_type(SArray{S,T,N,L},SVector{D,T})
  _to_multivalue_type(G)
end

function _gradient_type(::Type{T},::Type{<:VectorValue{D}}) where {T<:Real,D}
  VectorValue{D,T}
end

function _gradient_type(::Type{T},::Type{<:SVector{D}}) where {T<:Real,D}
  SVector{D,T}
end

function _to_multivalue_type(::Type{SArray{S,T,N,L}}) where {S,T,N,L}
  MultiValue{S,T,N,L}
end

@generated function _gradient_size(::Size{B},::Val{D}) where {B,D}
  str = join(["$b," for b in B])
  Meta.parse("Tuple{$D,$str}")
end

function _gradient_length(::Size{B},::Val{D}) where {B,D}
  prod((D,B...))
end

_mutable(::Type{SArray{S,T,N,L}}) where {S,T,N,L} = MArray{S,T,N,L}

_mutable(::Type{MultiValue{S,T,N,L}}) where {S,T,N,L} = MArray{S,T,N,L}

_length(x) = length(x)

_length(::Type{<:Real}) = 1

macro abstractmethod()
  quote
    error("This function belongs to an interface definition and cannot be used.")
  end
end

macro notimplemented()
  quote
    error("This function in not yet implemented")
  end
end

macro notimplementedif(condition)
  quote
    if $(esc(condition))
      @notimplemented
    end
  end
end

macro unreachable()
  quote
    error("This line of code cannot be reached")
  end
end


text:
"""ARMA process"""
function ARMA(X⃗::AbstractVector, ξ⃗::AbstractVector, ϕ⃗::AbstractVector, θ⃗::AbstractVector) # No constant, too boring
    # Just one update of an ARMA process
    if !all(length(X⃗)-1 .== (length(ξ⃗)-1, length(ϕ⃗), length(θ⃗)))
        p = max(length(X⃗)-1, length(ξ⃗)-1, length(ϕ⃗), length(θ⃗))
        X⃗ = vcat(X⃗, zeros(p+1-length(X⃗)))
        ξ⃗ = vcat(ξ⃗, randn(p+1-length(ξ⃗)))
        ϕ⃗ = vcat(ϕ⃗, zeros(p-length(ϕ⃗)))
        θ⃗ = vcat(θ⃗, zeros(p-length(θ⃗)))
    end
    X⃗[2:end] = X⃗[1:end-1]
    ξ⃗[2:end] = ξ⃗[1:end-1]
    ξ⃗[1] = randn()
    X⃗[1] = ξ⃗[1] + sum(ϕ⃗.*X⃗[2:end]) + sum(θ⃗.*ξ⃗[2:end])
    return (X⃗, ξ⃗)
end

function ARMA(X⃗::AbstractVector, ξ⃗::AbstractVector, ϕ₁::Float64=0.0, θ₁::Float64=0.0, args...)
    # Args are higher order parameters, alternating between ϕᵢ and θᵢ
    # E.g. ARMA(X⃗, ξ⃗, ϕ₁, θ₁, ϕ₂, θ₂)
    ϕ⃗ = vcat(ϕ₁, args[1:2:end])
    θ⃗ = vcat(θ₁, args[2:2:end])
    ARMA(X⃗, ξ⃗, ϕ⃗, θ⃗)
end

function ARMA(X0::AbstractVector, ξ0::AbstractVector, ϕ::Vector{Float64}, θ::Vector{Float64}, T::Int)
    X = zeros(T, length(X0))
    X[1, :] = X0
    for t ∈ 2:T
        X[t, :], ξ0 = ARMA(X[t-1, :], ξ0, ϕ, θ)
    end
    return X
end

function ARMA(X0::AbstractVector, ϕ::Vector{Float64}, θ::Vector{Float64}, T::Int)
    X = zeros(T, length(X0))
    X[1, :] = X0
    ξ = randn(length(X0))
    ARMA(X0, ξ, ϕ, θ, T)
end

"""AR process"""
function AR(X⃗::AbstractVector, ξ⃗::AbstractVector, ϕ⃗::AbstractVector)
    θ⃗ = zeros(length(ϕ⃗))
    ARMA(X⃗, ξ⃗, ϕ⃗, θ⃗)
end

function AR(X⃗::AbstractVector, ξ⃗::AbstractVector, ϕ₁::Float64=0.0, args...)
    ϕ⃗ = vcat(ϕ₁, args)
    AR(X⃗, ξ⃗, ϕ⃗)
end

function AR(X0::AbstractVector, ξ0::AbstractVector, ϕ::Vector{Float64}, T::Int)
    θ = zeros(length(ϕ))
    ARMA(X0, ξ0, ϕ, θ, T)
end

function AR(X0::AbstractVector, ϕ::Vector{Float64}, T::Int)
    θ = zeros(length(ϕ))
    ARMA(X0, ξ, ϕ, θ, T)
end

"""MA process"""
function MA(X⃗::AbstractVector, ξ⃗::AbstractVector, θ⃗::AbstractVector)
    ϕ⃗ = zeros(length(θ⃗))
    ARMA(X⃗, ξ⃗, ϕ⃗, θ⃗)
end

function MA(X⃗::AbstractVector, ξ⃗::AbstractVector, θ₁::Float64=0.0, args...)
    θ⃗ = vcat(θ₁, args)
    MA(X⃗, ξ⃗, θ⃗)
end

function MA(X0::AbstractVector, ξ0::AbstractVector, ϕ::Vector{Float64}, T::Int)
    ϕ = zeros(length(θ))
    ARMA(X0, ξ0, ϕ, θ, T)
end

function MA(X0::AbstractVector, ϕ::Vector{Float64}, T::Int)
    ϕ = zeros(length(θ))s
    ARMA(X0, ξ, ϕ, θ, T)
end

arma(args...) = ARMA(args...)
ar(args...) = AR(args...)
ma(args...) = MA(args...)


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

struct HyperbolicMesher
    topo::Symbol
    step::StepSize
    bc0::BoundaryCondition
    bc1::BoundaryCondition
    f_dissi::Float64
    f_blend::Float64
    n_smooth::Float64
    f_smooth::Float64
    wake::WakeInfo
    function HyperbolicMesher(
        topo::Symbol,
        step,
        bc0::BoundaryCondition=FreeBoundary(),
        bc1::BoundaryCondition=FreeBoundary();
        f_dissi::Real=1,
        f_blend::Real=0,
        n_smooth::Int=0,
        f_smooth::Real=1,
        wake_len::Real=0,
        wake_aoa::Real=0,
        wake_ratio::Real=1.1,
        wake_maxl::Real=Inf,
    )
        @assert topo in (:C, :O, :H)
        @assert 0 <= f_dissi <= 1
        @assert 0 <= f_blend <= 1
        @assert 0 <= f_smooth <= 1
        if topo == :C && wake_len > 0
            @assert wake_ratio >= 1.01
            @assert wake_maxl > 0
        else
            @assert wake_len <= 0 "'wake_*' is only supported in :C type grids"
        end

        return new(
            topo,
            StepSize(step),
            bc0,
            bc1,
            f_dissi,
            f_blend,
            n_smooth,
            f_smooth,
            WakeInfo(wake_len, wake_ratio, wake_maxl, wake_aoa),
        )
    end
end
include("mesh.jl")


text:
# Idea: Use LASSO for screening the active factors in an optimization
# problem.

using GLMNet

# Say that we want to optimize a function which has many factors but where
# only very few of them are active in calculating the outcome.
N = 1000
K = 3
active_factors = shuffle(collect(1:N))[1:K]
beta_intercept = 10.0
beta_slope = 1.0
betas = beta_intercept + beta_slope * randn(K, 1)
ofunc(x) = sum(x[active_factors] .* betas)

# Now lets sample and calc the function values for samples
S = int(10*K*log10(N));
xs = randn(N, S);
y = zeros(S)
[(y[i] = ofunc(xs[:,i])) for i in 1:S]

# Now do lasso regression:
#path = glmnet(xs', y; standardize = false)
# cv = glmnetcv(xs', y; standardize = false)

function active_factors_from_lasso(xs, y; alpha = 1.0)
  n = size(xs, 2)
  # Find the best lambda using cross-validation. We limit the size of the predictors.
  #cv = glmnetcv(xs, y) # ; constraints = hcat(-1000.0*ones(n), 1000.0*ones(n)))
  cv = glmnetcv(xs, y; standardize = false, alpha = alpha)
  num_steps = size(cv.path.betas, 2)
  res = Any[]
  for(i in 1:num_steps)
    pbetas = cv.path.betas[:, i]
    pbwi = [(pbetas[i], i) for i in 1:length(pbetas)]
    push!(res, filter((p) -> abs(p[1]) > 0.001, pbwi))
  end
  return res
end

# Count how many times each factor was found by the lasso.
function most_active_factors_from_lasso(xs, y; alpha = 1.0)
  lasso_factors = active_factors_from_lasso(xs, y; alpha = alpha)
  factor_counts = Dict{Int64, Int64}()
  for(i in 1:length(lasso_factors))
    factors = map((t) -> t[2], lasso_factors[i])
    for(f in factors)
      if haskey(factor_counts, f)
        factor_counts[f] += 1
      else
        factor_counts[f] = 1
      end
    end
  end
  return factor_counts, lasso_factors
end

function eval_if_lasso_finds_active_factors(N, K, S, reps = 100;
  beta_intercept = 0.0, beta_slope = 1.0, alpha = 1.0)

  active_factors = shuffle(collect(1:N))[1:K]
  print("Active factors:"); show(active_factors)
  betas = beta_intercept + beta_slope * randn(K, 1)
  print("\nBetas:"); show(betas)

  objective_func(x) = begin
    active = x[active_factors]
    sum(active .* betas)
  end

  # We want to count how many times the predicted factors (ranked based on how
  # often they were in the lasso factor set) are among the actually active factors.
  factor_found_counts = Dict{Int64, Int64}()

  for(rep in 1:reps)
    # Now lets sample
    xs = randn(N, S)

    # And calc the function values for them:
    y = zeros(S)
    for(i in 1:S)
      y[i] = objective_func(xs[:,i])
    end

    # Get the active factors from lasso reg
    factor_counts, factor_path = most_active_factors_from_lasso(xs', y; alpha = alpha)
    print("\nFactor counts:"); show(factor_counts)

    # Rank factors based on how often they were predicted as active
    factors_ranked = map((t) -> t[1], sort(collect(factor_counts), by = (t) -> t[2], rev = true))

    for(i in 1:length(factors_ranked))
      if in(factors_ranked[i], active_factors)
        if haskey(factor_found_counts, i)
          factor_found_counts[i] += 1
        else
          factor_found_counts[i] = 1
        end
      end
    end
  end

  return factor_found_counts
end

# It seems that 5*K*log10(N) samples is often always enough:
N = 10000
K = 3
@time eval_if_lasso_finds_active_factors(N, K, int(5*K*log10(N)), 10)
@time eval_if_lasso_finds_active_factors(N, K, int(5*K*log10(N)), 10; beta_intercept = 0.0)

# To make this even more systematic we need to be able to generate new objective functions
# involving a variety of factors, interactions and so forth.
type RandObjectiveFunction
  factors
  betas

  # Constructor creates a random function where we can control number of factors
  # per order, max order, total num vars as well as if effect heritage is active:
  #  N     = total number of variables
  #  maxOrder = maximum number of interaction terms in one factor (a product) of the sum that is the objective func
  #  sparsitiesPerOrder = array that gives the number of factors for that order.
  #  effect_heritage = true iff the factors of higher order can only involve vars used for order 1 (main effects)
  RandObjectiveFunction(N, maxOrder = 1, factorsPerOrder = [2];
    effect_heritage = true) = begin
    all_vars = 1:N

    # One factor is represented by the indices of vars that should be multiplied to create it.

    # Create main vars first since we might need them for heritance.
    main_vars = shuffle(collect(all_vars))[1:factorsPerOrder[1]]
    factors = map((v) -> [v], main_vars)

    # Add factors for higher orders
    for(order in 2:length(factorsPerOrder))
      for(i in 1:factorsPerOrder[order])
        factor = []
        for(j in 1:order)
          if effect_heritage
            vars_to_sample = main_vars
          else
            vars_to_sample = all_vars
          end
          push!(factor, shuffle(vars_to_sample)[1])
        end
        push!(factors, factor)
      end
    end

    # Create random betas for each factor, i.e. their weight.
    betas = randn(length(factors), 1)

    new(factors, betas)
  end
end

# Apply function f to an array a.
function eval(f::RandObjectiveFunction, a)
  sum = 0.0
  index = 1
  for(factor in f.factors)
    sum += f.betas[index] * a
    index += 1
  end
  s = size(xs, 2)
  y = zeros(s)
  for(i in 1:s)
    y[i] =
  end
end


text:
using Makie
using PortAudio
using DSP

"""
Slide the values in the given matrix to the right by 1.
The rightmosts column is discarded and the leftmost column is
left alone.
"""
function shift1!(buf::AbstractMatrix)
    for col in size(buf, 2):-1:2
        @. buf[:, col] = buf[:, col - 1]
    end
end

"""
takes a block of audio, FFT it, and write it to the beginning of the buffer
"""
function processbuf!(readbuf, win, dispbuf, fftbuf, fftplan)
    readbuf .*= win
    A_mul_B!(fftbuf, fftplan, readbuf)
    shift1!(dispbuf)
    @. dispbuf[end:-1:1, 1] = log(clamp(abs(fftbuf[1:D]), 0.0001, Inf))
end

function processblock!(src, buf, win, dispbufs, fftbuf, fftplan)
    read!(src, buf)
    for dispbuf in dispbufs
        processbuf!(buf, win, dispbuf, fftbuf, fftplan)
    end
end

N = 1024 # size of audio read
N2 = N ÷ 2 + 1 # size of rfft output
D = 200 # number of bins to display
M = 200 # amount of history to keep
src = PortAudioStream(1, 2)
buf = Array{Float32}(N) # buffer for reading
fftplan = plan_rfft(buf; flags = FFTW.EXHAUSTIVE)
fftbuf = Array{Complex{Float32}}(N2) # destination buf for FFT
dispbufs = [zeros(Float32, D, M) for i in 1:5, j in 1:5] # STFT bufs
win = gaussian(N, 0.125)

scene = Scene(resolution = (1000, 1000))

#pre-fill the display buffer so we can do a reasonable colormap
for _ in 1:M
    processblock!(src, buf, win, dispbufs, fftbuf, fftplan)
end

heatmaps = map(enumerate(IndexCartesian(), dispbufs)) do ibuf
    i = ibuf[1]
    buf = ibuf[2]

    # some function of the 2D index and the value
    heatmap(buf, offset = (i[2] * size(buf, 2), i[1] * size(buf, 1)))
end

center!(scene)

while isopen(scene[:screen])
    processblock!(src, buf, dispbufs, fftbuf, fftplan)
    for (hm, db) in zip(heatmaps, dispbufs)
        hm[:heatmap] = db
    end
    render_frame(scene)
end


text:
module ManifoldLearning

    using MultivariateStats

    import Base: show, dump, eigvals
    import MultivariateStats: outdim, projection, transform

    export

    ## common
    outdim,             # the output dimension of the transformation
    projection,         # the projection matrix
    eigvals,            # eigenvalues from the spectral analysis
    neighbors,          # the number of nearest neighbors used for building local structure
    ccomponent,         # point indexes of the largest connected component

    transform,          # perform the manifold learning

    # Transformation types
    Isomap,             # Type: Isomap model
    HLLE,               # Type: Hessian Eigenmaps model
    LLE,                # Type: Locally Linear Embedding model
    LTSA,               # Type: Local Tangent Space Alignment model
    LEM,                # Type: Laplacian Eigenmaps model
    DiffMap,            # Type: Diffusion maps model

    # example dataset
    swiss_roll          # swiss roll dataset generator

    include("types.jl")
    include("utils.jl")
    include("transformations.jl")
    include("isomap.jl")
    include("hlle.jl")
    include("lle.jl")
    include("ltsa.jl")
    include("lem.jl")
    include("diffmaps.jl")
end


text:
using LinearAlgebra
using MixedModels
using StableRNGs
using SparseArrays
using Test

using MixedModels: isconstant, average, densify, dataset
using StatsModels: FormulaTerm

@isdefined(io) || const global io = IOBuffer()
include("modelcache.jl")

@testset "average" begin
	@test average(1.1, 1.2) == 1.15
end

@testset "densify" begin
	@test densify(sparse(1:5, 1:5, ones(5))) == Diagonal(ones(5))
	rsparsev = SparseVector(float.(rand(StableRNG(123454321), Bool, 20)))
	@test densify(rsparsev) == Vector(rsparsev)
	@test densify(Diagonal(rsparsev)) == Diagonal(Vector(rsparsev))
end

@testset "isconstant" begin
	@test isconstant((true, true, true))
	@test isconstant([true, true, true])
	@test !isconstant((true, false, true))
	@test !isconstant([true, false, true])
	@test !isconstant(collect(1:4))
	@test isconstant((false, false, false))
	@test isconstant([false, false, false])
	@test isconstant(ones(3))
	@test isconstant(1, 1, 1)
	# equality of arrays with broadcasting
	@test isconstant(["(Intercept)", "days"], ["(Intercept)", "days"])
	# arrays or tuples with missing values
	@test !isconstant([missing, 1])
	@test isconstant(Int[])
	@test isconstant(Union{Int,Missing}[missing, missing, missing])
end

@testset "threaded_replicate" begin
	rng = StableRNG(42);
	single_thread = replicate(10;use_threads=false) do; only(randn(rng, 1)) ; end
	rng = StableRNG(42);
	multi_thread = replicate(10;use_threads=true) do
		if Threads.threadid() % 2 == 0
			sleep(0.001)
		end
		r = only(randn(rng, 1));
	end

	@test all(sort!(single_thread) .≈ sort!(multi_thread))
end

@testset "datasets" begin
	@test isa(MixedModels.datasets(), Vector{String})
	@test length(MixedModels.dataset(:dyestuff)) == 2
	@test length(MixedModels.dataset("dyestuff")) == 2
	dyestuff = MixedModels.dataset(:dyestuff);
	@test keys(dyestuff) == [:batch, :yield]
	@test length(dyestuff.batch) == 30
	@test_throws ArgumentError MixedModels.dataset(:foo)
end

@testset "PCA" begin
	io = IOBuffer()
	pca = models(:kb07)[3].PCA.item

	show(io, pca, covcor=true, loadings=false)
	str = String(take!(io))
	@test !isempty(findall("load: yes", str))

	show(io, pca, covcor=false, loadings=true)
	str = String(take!(io))
	@test !isempty(findall("PC1", str))
	@test !isempty(findall("load: yes", str))
end

@testset "formula" begin
	@test formula(first(models(:sleepstudy))) isa FormulaTerm
	@test formula(first(models(:contra))) isa FormulaTerm
end


text:
# Cairo backend for compose
import Cairo

import Cairo.CairoContext, Cairo.CairoSurface, Cairo.CairoARGBSurface,
       Cairo.CairoEPSSurface, Cairo.CairoPDFSurface, Cairo.CairoSVGSurface

abstract ImageBackend
abstract PNGBackend <: ImageBackend

abstract VectorImageBackend <: ImageBackend
abstract SVGBackend <: VectorImageBackend
abstract PDFBackend <: VectorImageBackend
abstract PSBackend  <: VectorImageBackend
abstract CairoBackend <: VectorImageBackend

type ImagePropertyState
    stroke::Maybe(ColorValue)
    fill::Maybe(ColorValue)
    fill_opacity::Float64
    stroke_opacity::Float64
    stroke_dash::Array{Float64,1}
    stroke_linecap::LineCap
    stroke_linejoin::LineJoin
    visible::Bool
end

type ImagePropertyFrame
    # Vector properties in this frame.
    vector_properties::Dict{Type, Property}

    # True if this property frame has scalar properties. Scalar properties are
    # emitted as a group (<g> tag) that must be closed when the frame is popped.
    has_scalar_properties::Bool

    function ImagePropertyFrame()
        return new(Dict{Type, Property}(), false)
    end
end

type Image{B <: ImageBackend} <: Backend
    out::IO
    surface::CairoSurface
    ctx::CairoContext

    width::Float64
    height::Float64

    # Current state
    stroke::Maybe(ColorValue)
    fill::Maybe(ColorValue)
    fill_opacity::Float64 # in [0,1]
    stroke_opacity::Float64 # in [0,1]
    stroke_dash::Array{Float64,1}
    stroke_linecap::LineCap
    stroke_linejoin::LineJoin
    visible::Bool

    # Keep track of property
    state_stack::Vector{ImagePropertyState}
    property_stack::Vector{ImagePropertyFrame}
    vector_properties::Dict{Type, Union(Nothing, Property)}

    # Close the surface when finished
    owns_surface::Bool

    # Backend is responsible for opening/closing the file
    ownedfile::Bool

    # Filename when ownedfile is true
    filename::Union(String, Nothing)

    # True when finish has been called and no more drawing should occur
    finished::Bool

    # Emit on finish
    emit_on_finish::Bool

    function Image(surface::CairoSurface, ctx::CairoContext, out::IO)
        img = new()
        img.out = out
        img.width = 0
        img.height = 0
        img.surface = surface
        img.ctx = ctx
        img.stroke = default_stroke_color
        img.fill   = default_fill_color
        img.fill_opacity = 1.0
        img.stroke_opacity = 1.0
        img.stroke_dash = []
        img.stroke_linecap = LineCapButt()
        img.stroke_linejoin = LineJoinMiter()
        img.visible = true
        img.state_stack = Array(ImagePropertyState, 0)
        img.property_stack = Array(ImagePropertyFrame, 0)
        img.vector_properties = Dict{Type, Union(Nothing, Property)}()
        img.owns_surface = false
        img.ownedfile = false
        img.filename = nothing
        img.finished = false
        img.emit_on_finish = false
        img
    end

    function Image(surface::CairoSurface, ctx::CairoContext)
        Image{B}(surface, ctx, IOBuffer())
    end

    Image(surface::CairoSurface) = Image{B}(surface, CairoContext(surface))

    function Image(out::IO,
                   width::MeasureOrNumber,
                   height::MeasureOrNumber,
                   emit_on_finish::Bool=true)

        width = size_measure(width)
        height = size_measure(height)

        if !isabsolute(width) || !isabsolute(height)
            error("Image size must be specificed in absolute units.")
        end

        width = absolute_native_units(B, width.abs)
        height = absolute_native_units(B, height.abs)

        surface = newsurface(B, out, width, height)
        img = Image{B}(surface, CairoContext(surface), out)
        img.width = width
        img.height = height
        img.owns_surface = true
        img.emit_on_finish = emit_on_finish
        img
    end

    function Image(filename::String,
                   width::MeasureOrNumber,
                   height::MeasureOrNumber)
        img = Image{B}(open(filename, "w"), width, height)
        img.ownedfile = true
        img.filename = filename
        img
    end

    function Image(width::MeasureOrNumber,
                   height::MeasureOrNumber,
                   emit_on_finish::Bool=true)
        img = Image{B}(IOBuffer(), width, height, emit_on_finish)
        img
    end

    function Image(c::CairoSurface)
        img = Image{B}(c,CairoContext(c))
        img
    end
end

typealias PNG Image{PNGBackend}
typealias PDF Image{PDFBackend}
typealias PS  Image{PSBackend}
typealias CAIROSURFACE  Image{CairoBackend}

# convert compose absolute units (millimeters) to the absolute units used by the
# cairo surface (pixels for PNG, points for all others)
function absolute_native_units(::Type{PNGBackend}, u::Float64)
    assumed_ppmm * u
end

function absolute_native_units{B <: ImageBackend}(::Type{B},
                                                  u::Float64)
    72 * u / 25.4
end

function absolute_native_units{B}(img::Image{B}, u::Float64)
    absolute_native_units(B, u)
end

# Width and height of a backend in absolute units
function width(img::PNG)
    Cairo.width(img.surface) / assumed_ppmm
end

function width(img::Image)
    25.4 * Cairo.width(img.surface) / 72
end

function height(img::PNG)
    Cairo.height(img.surface) / assumed_ppmm
end

function height(img::Image)
    25.4 * Cairo.height(img.surface) / 72
end

function iswithjs(img::Image)
    return false
end

function iswithousjs(img::Image)
    return true
end

finish{B <: ImageBackend}(::Type{B}, img::Image) = nothing

function finish(::Type{PNGBackend}, img::Image)
    Cairo.write_to_png(img.surface, img.out)
end

function finish{B <: ImageBackend}(img::Image{B})
    if img.finished
        return
    end

    if(img.owns_surface)
        Cairo.destroy(img.ctx)
    end
    finish(B, img)
    if(img.owns_surface)
        Cairo.destroy(img.surface)
    end

    img.finished = true

    if img.emit_on_finish && typeof(img.out) == IOBuffer
        display(img)
    end

    if method_exists(flush, (typeof(img.out),))
        flush(img.out)
    end

    if img.ownedfile
        close(img.out)
    end
end

function newsurface{B}(::Type{B}, out, width, height)
        local surface::CairoSurface
        if B == SVGBackend
            surface = CairoSVGSurface(out, width, height)
        elseif B == PNGBackend
            surface = CairoARGBSurface(iround(width), iround(height))
        elseif B == PDFBackend
            surface = CairoPDFSurface(out, width, height)
        elseif B == PSBackend
            surface = CairoEPSSurface(out, width, height)
        elseif B == CairoBackend
            surface = out
        else
            error("Unkown Cairo backend.")
        end

        if Cairo.status(surface) != Cairo.STATUS_SUCCESS
            error("Unable to create cairo surface.")
        end

        surface
end

function reset{B}(img::Image{B})
    if !img.owns_surface
        error("Backend can't be reused since an external cairo surface is being used.")
    end

    if img.ownedfile
        img.out = open(img.filename, "w")
    else
        try
            seekstart(img.out)
        catch
            error("Backend can't be reused, since the output stream is not seekable.")
        end
    end

    img.surface = newsurface(B, img.out, img.width, img.height)
    img.ctx = CairoContext(img.surface)
    img.finished = false
end

function isfinished(img::Image)
    img.finished
end

function root_box(img::Image)
    AbsoluteBoundingBox(0.0, 0.0, width(img), height(img))
end

function writemime(io::IO, ::MIME"image/png", img::PNG)
    write(io, takebuf_string(img.out))
end

function writemime(io::IO, ::MIME"application/pdf", img::PDF)
    write(io, takebuf_string(img.out))
end

function writemime(io::IO, ::MIME"application/postscript", img::PS)
    write(io, takebuf_string(img.out))
end

# Applying Properties
# -------------------

function push_property_frame(img::Image, properties::Vector{Property})
    if isempty(properties)
        return
    end

    frame = ImagePropertyFrame()
    applied_properties = Set{Type}()
    scalar_properties = Array(Property, 0)
    for property in properties
        if isscalar(property) && !(typeof(property) in applied_properties)
            push!(scalar_properties, property)
            push!(applied_properties, typeof(property))
            frame.has_scalar_properties = true
        elseif !isscalar(property)
            frame.vector_properties[typeof(property)] = property
            img.vector_properties[typeof(property)] = property
        end
    end
    push!(img.property_stack, frame)
    if isempty(scalar_properties)
        return
    end

    save_property_state(img)
    for property in scalar_properties
        apply_property(img, property.primitives[1])
    end
end

function pop_property_frame(img::Image)
    @assert !isempty(img.property_stack)
    frame = pop!(img.property_stack)

    if frame.has_scalar_properties
        restore_property_state(img)
    end

    for (propertytype, property) in frame.vector_properties
        img.vector_properties[propertytype] = nothing
        for i in length(img.property_stack):-1:1
            if haskey(img.property_stack[i].vector_properties, propertytype)
                img.vector_properties[propertytype] =
                    img.property_stack[i].vector_properties[propertytype]
            end
        end
    end
end

function save_property_state(img::Image)
    push!(img.state_stack,
        ImagePropertyState(
            img.stroke,
            img.fill,
            img.fill_opacity,
            img.stroke_opacity,
            img.stroke_dash,
            img.stroke_linecap,
            img.stroke_linejoin,
            img.visible))
    Cairo.save(img.ctx)
end

function restore_property_state(img::Image)
    state = pop!(img.state_stack)
    img.stroke = state.stroke
    img.fill = state.fill
    img.fill_opacity = state.fill_opacity
    img.stroke_opacity = state.stroke_opacity
    img.stroke_dash = state.stroke_dash
    img.stroke_linecap = state.stroke_linecap
    img.stroke_linejoin = state.stroke_linejoin
    img.visible = state.visible
    Cairo.restore(img.ctx)
end

function push_vector_properties(img::Image, idx::Int)
    save_property_state(img)
    for (propertytype, property) in img.vector_properties
        if property === nothing
            continue
        end
        if idx > length(property.primitives)
            error("Vector form and vector property differ in length. Can't distribute.")
        end
        apply_property(img, property.primitives[idx])
    end
end

function pop_vector_properties(img::Image)
    restore_property_state(img)
end

function apply_property(img::Image, p::StrokePrimitive)
    if isa(p.color, AlphaColorValue)
        img.stroke = p.color.c
        img.stroke_opacity = p.color.alpha
    else
        img.stroke = p.color
    end
end

function apply_property(img::Image, p::FillPrimitive)
    if isa(p.color, AlphaColorValue)
        img.fill = p.color.c
        img.fill_opacity = p.color.alpha
    else
        img.fill = p.color
    end
end

function apply_property(img::Image, p::FillOpacityPrimitive)
    img.fill_opacity = p.value
end

function apply_property(img::Image, p::StrokeOpacityPrimitive)
    img.stroke_opacity = p.value
end

function apply_property(img::Image, p::StrokeDashPrimitive)
    img.stroke_dash = map(v -> absolute_native_units(img, v.abs), p.value)
end

function apply_property(img::Image, p::StrokeLineCapPrimitive)
    img.stroke_linecap = p.value
end

function apply_property(img::Image, p::StrokeLineJoinPrimitive)
    img.stroke_linejoin = p.value
end

function apply_property(img::Image, p::VisiblePrimitive)
    img.visible = p.value
end

function apply_property(img::Image, property::LineWidthPrimitive)
    Cairo.set_line_width(
        img.ctx,
        absolute_native_units(img, property.value.abs))
end

function apply_property(img::Image, property::FontPrimitive)
    font_desc = ccall((:pango_layout_get_font_description, Cairo._jl_libpango),
                      Ptr{Void}, (Ptr{Void},), img.ctx.layout)

    if font_desc == C_NULL
        size = absolute_native_units(img, default_font_size.abs)
    else
        size = ccall((:pango_font_description_get_size, Cairo._jl_libpango),
                     Cint, (Ptr{Void},), font_desc)
    end

    Cairo.set_font_face(img.ctx,
        @sprintf("%s %0.2fpx", property.family, size / PANGO_SCALE))
end

function apply_property(img::Image, property::FontSizePrimitive)
    font_desc = ccall((:pango_layout_get_font_description, Cairo._jl_libpango),
                      Ptr{Void}, (Ptr{Void},), img.ctx.layout)

    if font_desc == C_NULL
        family = "sans"
    else
        family = ccall((:pango_font_description_get_family, Cairo._jl_libpango),
                       Ptr{Uint8}, (Ptr{Void},), font_desc)
        family = bytestring(family)
    end

    Cairo.set_font_face(img.ctx,
        @sprintf("%s %.2fpx",
            family, absolute_native_units(img, property.value.abs)))
end

function apply_property(img::Image, property::ClipPrimitive)
    if isempty(property.points); return; end
    move_to(img, property.points[1])
    for point in property.points[2:end]
        line_to(img, point)
    end
    close_path(img)
    Cairo.clip(img.ctx)
end

# No-op SVG+JS only properties
function apply_property(img::Image, property::JSIncludePrimitive)
end

function apply_property(img::Image, property::JSCallPrimitive)
end

function apply_property(img::Image, property::SVGClassPrimitive)
end

function apply_property(img::Image, property::SVGAttributePrimitive)
end

# Cairo Wrappers
# --------------

function move_to(img::Image, point::Point)
    Cairo.move_to(
        img.ctx,
        absolute_native_units(img, point.x.abs),
        absolute_native_units(img, point.y.abs))
end

function line_to(img::Image, point::Point)
    Cairo.line_to(
        img.ctx,
        absolute_native_units(img, point.x.abs),
        absolute_native_units(img, point.y.abs))
end

function rectangle(img::Image, corner::Point, width::Measure, height::Measure)
    Cairo.rectangle(img.ctx,
                    absolute_native_units(img, corner.x.abs),
                    absolute_native_units(img, corner.y.abs),
                    absolute_native_units(img, width.abs),
                    absolute_native_units(img, height.abs))
end

function circle(img::Image, center::Point, radius::Measure)
    Cairo.circle(img.ctx,
                 absolute_native_units(img, center.x.abs),
                 absolute_native_units(img, center.y.abs),
                 absolute_native_units(img, radius.abs))
end

function curve_to(img::Image, ctrl1::Point, ctrl2::Point, anchor::Point)
    Cairo.curve_to(
        img.ctx,
        absolute_native_units(img, ctrl1.x.abs),
        absolute_native_units(img, ctrl1.y.abs),
        absolute_native_units(img, ctrl2.x.abs),
        absolute_native_units(img, ctrl2.y.abs),
        absolute_native_units(img, anchor.x.abs),
        absolute_native_units(img, anchor.y.abs))
end

function close_path(img::Image)
    Cairo.close_path(img.ctx)
end

function arc(img::Image, x::Float64, y::Float64, radius::Float64,
                angle1::Float64, angle2::Float64)
    Cairo.arc(
        img.ctx,
        absolute_native_units(img, x),
        absolute_native_units(img, y),
        absolute_native_units(img, radius),
        absolute_native_units(img, angle1),
        absolute_native_units(img, angle2))
end

function translate(img::Image, tx::Float64, ty::Float64)
    Cairo.translate(
        img.ctx,
        absolute_native_units(img, tx),
        absolute_native_units(img, ty))
end

function scale(img::Image, sx::Float64, sy::Float64)
    Cairo.scale(img.ctx, sx, sy)
end

function rotate(img::Image, theta::Float64)
    Cairo.rotate(img.ctx, theta)
end

function rotate(img::Image, theta::Float64, x::Float64, y::Float64)
    ct = cos(theta)
    st = sin(theta)

    x′ = x - (ct * x - st * y)
    y′ = y - (st * x + ct * y)

    translate(img, x′, y′)
    rotate(img, theta)
end

# Convert native linecap/linejoin enums to the Cairo values.
cairo_linecap(::LineCapButt) = Cairo.CAIRO_LINE_CAP_BUTT
cairo_linecap(::LineCapRound) = Cairo.CAIRO_LINE_CAP_ROUND
cairo_linecap(::LineCapSquare) = Cairo.CAIRO_LINE_CAP_SQUARE
cairo_linejoin(::LineJoinMiter) = Cairo.CAIRO_LINE_JOIN_MITER
cairo_linejoin(::LineJoinBevel) = Cairo.CAIRO_LINE_JOIN_BEVEL
cairo_linejoin(::LineJoinRound) = Cairo.CAIRO_LINE_JOIN_ROUND

function fillstroke(img::Image)
    if img.fill != nothing && img.fill_opacity > 0.0 && img.visible
        rgb = convert(RGB, img.fill)
        Cairo.set_source_rgba(img.ctx, rgb.r, rgb.g, rgb.b, img.fill_opacity)

        if img.stroke != nothing
            Cairo.fill_preserve(img.ctx)
        else
            Cairo.fill(img.ctx)
        end
    end

    if img.stroke != nothing && img.stroke_opacity > 0.0
        rgb = convert(RGB, img.stroke)
        Cairo.set_source_rgba(img.ctx, rgb.r, rgb.g, rgb.b, img.stroke_opacity)
        Cairo.set_dash(img.ctx, img.stroke_dash)
        Cairo.set_line_cap(img.ctx, cairo_linecap(img.stroke_linecap))
        Cairo.set_line_join(img.ctx, cairo_linejoin(img.stroke_linejoin))

        Cairo.stroke(img.ctx)
    end

    # if the path wasn't stroked or filled, we should still clear it
    Cairo.new_path(img.ctx)
end

# Form Drawing
# ------------

function draw(img::Image, form::Form)
    for (idx, primitive) in enumerate(form.primitives)
        push_vector_properties(img, idx)
        draw(img, primitive)
        pop_vector_properties(img)
    end
end

function draw(img::Image, prim::RectanglePrimitive)
    rectangle(img, prim.corner, prim.width, prim.height)
    fillstroke(img)
end

function draw(img::Image, prim::PolygonPrimitive)
    if isempty(prim.points); return; end

    paths = make_paths(prim.points)
    for path in paths
        move_to(img, path[1])
        for point in path[2:end]
            line_to(img, point)
        end
        close_path(img)
        fillstroke(img)
    end
end

function draw(img::Image, prim::CirclePrimitive)
    circle(img, prim.center, prim.radius)
    fillstroke(img)
end

function draw(img::Image, prim::EllipsePrimitive)
    cx = prim.center.x.abs
    cy = prim.center.y.abs
    rx = sqrt((prim.x_point.x.abs - cx)^2 +
              (prim.x_point.y.abs - cy)^2)
    ry = sqrt((prim.y_point.x.abs - cx)^2 +
              (prim.y_point.y.abs - cy)^2)
    theta = atan2(prim.x_point.y.abs - cy,
                  prim.x_point.x.abs - cx)

    if !all(isfinite([cx, cy, rx, ry, theta]))
        return
    end

    save_property_state(img)
    translate(img, cx, cy)
    rotate(img, theta)
    translate(img, -rx, -ry)
    scale(img, 2rx, 2ry)
    arc(img, 0.5, 0.5, 0.5, 0.0, 2pi)
    restore_property_state(img)
    fillstroke(img)
end

function draw(img::Image, prim::LinePrimitive)
    if length(prim.points) <= 1; return; end

    paths = make_paths(prim.points)
    for path in paths
        move_to(img, path[1])
        for point in path[2:end]
            line_to(img, point)
        end
        fillstroke(img)
    end
end

function get_layout_size(img::PNG)
    width, height = Cairo.get_layout_size(img.ctx)
    width / assumed_ppmm, height / assumed_ppmm
end

function get_layout_size(img::Image)
    width, height = Cairo.get_layout_size(img.ctx)
    25.4 * width / 72, 25.4 * height / 72
end

function draw(img::Image, prim::TextPrimitive)
    if !img.visible || ((img.fill_opacity == 0.0 || img.fill === nothing) &&
            (img.stroke_opacity == 0.0 || img.stroke === nothing))
        return
    end

    pos = copy(prim.position)
    Cairo.set_text(img.ctx, prim.value, true)
    width, height = get_layout_size(img)
    pos = Point(pos.x, Measure(abs=pos.y.abs - height))

    if prim.halign != hleft || prim.valign != vbottom
        if prim.halign == hcenter
            pos = Point(Measure(abs=pos.x.abs - width/2), pos.y)
        elseif prim.halign == hright
            pos = Point(Measure(abs=pos.x.abs - width), pos.y)
        end

        if prim.valign == vcenter
            pos = Point(pos.x, Measure(abs=pos.y.abs + height/2))
        elseif prim.valign == vtop
            pos = Point(pos.x, Measure(abs=pos.y.abs + height))
        end
    end

    rgb = convert(RGB, img.fill)
    Cairo.set_source_rgba(img.ctx, rgb.r, rgb.g, rgb.b, img.fill_opacity)

    if prim.rot.theta != 0.0
        save_property_state(img)
        rotate(img, prim.rot.theta,
               prim.rot.offset.x.abs, prim.rot.offset.y.abs)
    end

    move_to(img, pos)
    Cairo.show_layout(img.ctx)

    if prim.rot.theta != 0.0
        restore_property_state(img)
    end
end

function draw(img::Image, prim::CurvePrimitive)
    move_to(img, prim.anchor0)
    curve_to(img, prim.ctrl0, prim.ctrl1, prim.anchor1)
end

function draw(img::Image, prim::BitmapPrimitive)
    error("Embedding bitmaps in Cairo backends (i.e. PNG, PDF, PS) is not supported.")
end


text:
"""
	openp : example from Bazinga.jl

	Parametric constrained minimization of the Rosenbrock function, from [SFP20].

	Original formulations:

	minimize      sum( p2 ( x[i+1] - x[i]^2 ) + ( p1 - x[i] )^2 ) for i=1:4
	subject to    || x || <= p4
	              p3 sin( x1 ) = cos( x2 + x3 )
	              x3 + x4 <= p5
    with          p = [1, 50, 1.5, 0.73, 0.2]

	Reformulation as a constrained structured problem in the form

	minimize     f(x) + g(x)
	subject to   c(x) in S

	where
	    f(x) = sum( p2 ( x[i+1] - x[i]^2 )^2 + ( p1 - x[i] )^2 ) for i=1:4
	    g(x) = IndBall2[radius=p4](x)
	    c(x) = [ p3 sin( x1 ) - cos( x2 + x3 ) ]
	           [ x3 + x4 - p5                  ]
	    S    = { c | c1 = 0, c2 <= 0 }

    References:
    [SFP20]     Sopasakis, Fresk, Patrinos, "OpEn: Code Generation for Embedded
                Nonconvex Optimization" (2020).
                21st IFAC World Congress: Proceedings
                arXiv:2003.00292v1 [math.OC] 29 Feb 2020
"""

using Bazinga, OptiMo
using Random, LinearAlgebra
using Printf

###################################################################################
# problem definition
###################################################################################
mutable struct OPENP <: AbstractOptiModel
  meta::OptiModelMeta
  p::Vector
end

function OPENP(; p::Vector=Vector([1.0, 50.0, 1.5, 0.73, 0.2]))
    @assert length(p) == 5
    @assert 0 < p[4]
    name = "OpEnParamExample"
    meta = OptiModelMeta( 5, 2, x0=zeros(Float64, 5), name=name)
    return OPENP( meta, p )
end

# necessary methods:
# obj, grad!: cons!, jprod!, jtprod!, proj!, prox!, objprox!
function OptiMo.obj( prob::OPENP, x::AbstractVector )
    OptiMo.@lencheck prob.meta.nvar x
    fx = 0
    for i=1:4
        fx += prob.p[2] * (x[i+1] - x[i]^2)^2 + (prob.p[1] - x[i])^2
    end
    return fx
end

function OptiMo.grad!( prob::OPENP, x::AbstractVector, dfx::AbstractVector )
    OptiMo.@lencheck prob.meta.nvar x
    dfx[1] = - 4.0 * prob.p[2] * x[1] * (x[2] - x[1]^2) - 2.0 * (prob.p[1] - x[1])
    dfx[2] = 2.0 * prob.p[2] * (x[2] - x[1]^2) - 4.0 * prob.p[2] * x[2] * (x[3] - x[2]^2) - 2.0 * (prob.p[1] - x[2])
    dfx[3] = 2.0 * prob.p[2] * (x[3] - x[2]^2) - 4.0 * prob.p[2] * x[3] * (x[4] - x[3]^2) - 2.0 * (prob.p[1] - x[3])
    dfx[4] = 2.0 * prob.p[2] * (x[4] - x[3]^2) - 4.0 * prob.p[2] * x[4] * (x[5] - x[4]^2) - 2.0 * (prob.p[1] - x[4])
    dfx[5] = 2.0 * prob.p[2] * (x[5] - x[4]^2)
    return nothing
end

function OptiMo.cons!( prob::OPENP, x::AbstractVector, cx::AbstractVector )
    OptiMo.@lencheck prob.meta.nvar x
    OptiMo.@lencheck prob.meta.ncon cx
    cx.= [ prob.p[3] * sin( x[1] ) - cos( x[2] + x[3] ); x[3] + x[4] - prob.p[5] ]
    return nothing
end

function OptiMo.jprod!( prob::OPENP, x::AbstractVector, v::AbstractVector, Jv::AbstractVector )
    OptiMo.@lencheck prob.meta.nvar x v
    OptiMo.@lencheck prob.meta.ncon Jv
    Jv .= [ prob.p[3] * cos( x[1] ) * v[1] + sin( x[2] + x[3] ) * ( v[2] + v[3] ); v[3] + v[4] ]
    return nothing
end

function OptiMo.jtprod!( prob::OPENP, x::AbstractVector, v::AbstractVector, Jtv::AbstractVector )
    OptiMo.@lencheck prob.meta.nvar x Jtv
    OptiMo.@lencheck prob.meta.ncon v
    sx23v1 = sin( x[2] + x[3] ) * v[1]
    Jtv .= [ prob.p[3] * cos( x[1] ) * v[1]; sx23v1; sx23v1 + v[2]; v[2]; 0.0]
    return nothing
end

function OptiMo.prox!( prob::OPENP, x::AbstractVector, a::Real, z::AbstractVector )
    OptiMo.@lencheck prob.meta.nvar x z
    z .= x .* min( 1.0, prob.p[4] / norm( x, 2 ) )
    return nothing
end

function OptiMo.objprox!( prob::OPENP, x::AbstractVector, a::Real, z::AbstractVector )
    OptiMo.@lencheck prob.meta.nvar x z
    z .= x .* min( 1.0, prob.p[4] / norm( x, 2 ) )
    return 0.0
end

function OptiMo.proj!( prob::OPENP, cx::AbstractVector, px::AbstractVector )
    OptiMo.@lencheck prob.meta.ncon cx px
    px .= [ 0.0; min( 0.0, cx[2] ) ]
    return nothing
end

# problem build
p0 = [1.0; 50.0; 1.5; 0.73; 0.2]
problem = OPENP( p=p0 )

# solver build
VERBOSE = false

# warm up
solver = Bazinga.ALPX( tol_optim=1e-5, # OpEn 1e-5
                       tol_cviol=1e-4,
                       verbose=VERBOSE )
out = solver( problem )
print( out )
@printf "\n\n\n"

# OpEn setup
solver = Bazinga.ALPX( tol_optim=1e-5, # OpEn 1e-5
                       tol_cviol=1e-4, # OpEn 1e-4
                       eps_init=1e-4,  # OpEn 1e-4
                       mu_down=0.2,    # OpEn 0.2
                       verbose=false )
mu0 = 1e-3 * ones(Float64, problem.meta.ncon) # OpEn 1e-3
out = solver( problem, mu=mu0 )
print( out )

# ALPX setup, OpEn initial penalty
solver = Bazinga.ALPX( tol_optim=1e-5, # OpEn 1e-5
                       tol_cviol=1e-4, # OpEn 1e-4
                       verbose=VERBOSE )
mu0 = 1e-3 * ones(Float64, problem.meta.ncon) # OpEn 1e-3
out = solver( problem, mu=mu0 )
print( out )

# ALPX setup, OpEn tolerances
solver = Bazinga.ALPX( tol_optim=1e-5, # OpEn 1e-5
                       tol_cviol=1e-4, # OpEn 1e-4
                       verbose=VERBOSE )
out = solver( problem )
print( out )

# ALPX setup
solver = Bazinga.ALPX( verbose=VERBOSE )
out = solver( problem )
print( out )


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Mortar2D.jl/blob/master/LICENSE

using Test
using Mortar2D: calculate_mortar_matrices

# calculate mortar matrices
elements = Dict(
    1 => [1, 2],
    2 => [3, 4])
element_types = Dict(
    1 => :Seg2,
    2 => :Seg2)
coords = Dict(
    1 => [1.0, 2.0],
    2 => [3.0, 2.0],
    3 => [2.0, 2.0],
    4 => [0.0, 2.0])
normals = Dict(
    1 => [0.0, -1.0],
    2 => [0.0, -1.0])
segmentation = Dict(1 => [(2, [-1.0, 0.0])])
De, Me = calculate_mortar_matrices(1, elements, element_types,
                                   coords, normals, segmentation)
De_expected = 1/24*[14 4; 4 2]
Me_expected = 1/24*[13 5; 5 1]
@test isapprox(De, De_expected)
@test isapprox(Me[2], Me_expected)


text:
baremodule ShowGraphviz

macro derive end
macro deriveall end

struct DOT
    source::String
    DOT(source::AbstractString) = new(source)
end

macro dot_str end

function show end

function getoption end
function setoption end
function addoption end

module Internal

import ..ShowGraphviz: @derive, @deriveall, @dot_str
using ..ShowGraphviz: ShowGraphviz

# using Graphviz_jll: dot

include("utils.jl")
include("config.jl")
include("derive.jl")

# Use README as the docstring of the module:
function define_docstring()
    path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    doc = replace(read(path, String), r"^```julia"m => "```jldoctest README")
    @eval ShowGraphviz $Base.@doc $doc ShowGraphviz
end

end  # module Internal

const CONFIG = Internal.CONFIG

Internal.define_docstring()

end  # baremodule ShowGraphviz


text: AbstractCTM = Union{CTM, gpuCTM}


text:
##### Beginning of file

module Utils # Begin submodule MirrorUpdater.Utils

__precompile__(true)

import ..delayederror

include("git.jl")

include("misc.jl")

end # End submodule MirrorUpdater.Utils

##### End of file


text:
"""
    AbstractBoundaryConditionClassification

Abstract supertype for boundary condition types.
"""
abstract type AbstractBoundaryConditionClassification end

"""
    struct Periodic <: AbstractBoundaryConditionClassification

A classification specifying a periodic boundary condition.

A condition may not be specified with a `Periodic` boundary condition.
"""
struct Periodic <: AbstractBoundaryConditionClassification end

"""
    struct Flux <: AbstractBoundaryConditionClassification

A classification specifying a boundary condition on the flux of a field.

The sign convention is such that a positive flux represents the flux of a quantity in the
positive direction. For example, a positive vertical flux implies a quantity is fluxed
upwards, in the ``+z`` direction.

Due to this convention, a positive flux applied to the top boundary specifies that a quantity
is fluxed upwards across the top boundary and thus out of the domain. As a result, a positive
flux applied to a top boundary leads to a reduction of that quantity in the interior of the
domain; for example, a positive, upwards flux of heat at the top of the domain acts to cool
the interior of the domain. Conversely, a positive flux applied to the bottom boundary leads
to an increase of the quantity in the interior of the domain. The same logic holds for east,
west, north, and south boundaries.
"""
struct Flux <: AbstractBoundaryConditionClassification end

"""
    struct Gradient <: AbstractBoundaryConditionClassification

A classification specifying a boundary condition on the derivative or gradient of a field. Also
called a Neumann boundary condition.
"""
struct Gradient <: AbstractBoundaryConditionClassification end

"""
    struct Value <: AbstractBoundaryConditionClassification

A classification specifying a boundary condition on the value of a field. Also called a Dirchlet
boundary condition.
"""
struct Value <: AbstractBoundaryConditionClassification end

"""
    struct Open <: AbstractBoundaryConditionClassification

A classification that specifies the halo regions of a field directly.

For fields located at Faces, Open also specifies field value _on_ the boundary.

Open boundary conditions are used to specify the component of a velocity field normal to a boundary
and can also be used to describe nested or linked simulation domains.
"""
struct Open <: AbstractBoundaryConditionClassification end


text:
function advance_state_internal!(agent::Agents.AbstractAgent, model::Agents.ABM)
    if agent.quarantined
        agent.ticks_quarantined += 1
    end
    if (agent.state == :Iu) | (agent.state == :Id)
        agent.ticks_infected += 1
    end
    if (agent.state == :E)
        if (agent.ticks_exposed >= agent.incubation_period)
            agent.state = :Iu
            agent.ticks_infected = 1
            return agent
        end
        agent.ticks_exposed += 1
    end
    if ((agent.state == :Iu) | (agent.state == :Id)) & (agent.ticks_infected >= agent.infection_duration)
        agent.state = :R
        return agent
    end
    if (agent.state == :Iu) & (agent.ticks_infected % 10 == 0) & (agent.ticks_infected < agent.infection_duration)
        if Random.rand() < model.test_prob
            agent.state = :Id
            model.new_cases_detected += 1
            agent.quarantined = true
            agent.ticks_quarantined = 1
        end
        return agent
    end
    return agent
end

function simulate_node_encounters!(model::Agents.ABM, node::Int64)
    encounters = Agents.get_node_agents(node, model)
    for candidate in encounters
        for interaction_partner in encounters
            adequate_contact!(candidate, interaction_partner, model)
        end
    end
    return model
end

function adequate_contact!(agent::Agents.AbstractAgent, encounter::Agents.AbstractAgent, model::Agents.ABM)
    if (agent.state == :S) & ((encounter.state == :Iu) | (encounter.state == :Id))
        get_exposed!(agent, model)
    end
    return agent
end

function get_exposed!(agent::Agents.AbstractAgent, model::Agents.ABM)
    if agent.behavior
        if Random.rand() < (model.transmission_prob * model.behavior_reduction_factor)
            agent.state = :E
            model.new_cases_real += 1
            agent.ticks_exposed = 1
        end
    else
        if Random.rand() < model.transmission_prob
            agent.state = :E
            model.new_cases_real += 1
            agent.ticks_exposed = 1
        end
    end
    return agent
end


text:
# This program reads a day of data in CSV form and computes H3 hash
# for each data point. The final result for the data is written to a
# feather file for consumption by Arrow (probably from Python)

using H3
using H3.API
using CSV
using DataFrames
using Formatting
using Arrow

# read one day of data by accumulating tables, one per CSV file
t = []
for i in 0:23
    tx = DataFrame(CSV.File("grib/20210806" * fmt("02d", i) * ".grib2.csv"))
    push!(t, tx)
end
# now glue those all together
allDay = vcat(t...)

# Add H3 tags to each data point
tiles = Dict{Int, Vector{H3Index}}()
for row in eachrow(allDay)
    for resolution in 3:9
        hex = geoToH3(GeoCoord(deg2rad(row[:latitude]), deg2rad(row[:longitude])), resolution)
        push!(get!(tiles, resolution, []), hex)
    end
end
for resolution in keys(tiles)
    allDay[!,"tile_$resolution"] = tiles[resolution]
end

sort!(allDay, [:tile, :latitude, :longitude, :datetime])

# A tiled frame is grouped on some column for writing out Feather or Parquet files
# this lets us use Table.partitioner
struct TiledFrame
    v::DataFrame
    col::Any
end

import Base.iterate, Base.getindex, Base.size
function Base.iterate(S::TiledFrame, i=0)
    if i > size(S.v)[1]
        return nothing
    elseif i == 0
        return Base.iterate(S, 1)
    else
        r0 = S[i, S.col]
        i0 = i
        while i <= size(S)[1] && S[i, S.col] == r0
            i += 1
        end
        return (S.v[i0:(i-1),:], i)
    end
end

function Base.getindex(S::TiledFrame, i::Any)
    return getindex(S.v, i)
end

function Base.getindex(S::TiledFrame, i::Any, j::Any)
    return getindex(S.v, i, j)
end

function size(S::TiledFrame)
    return size(S.v)
end

Arrow.write("x20210806.parquet", Tables.partitioner(TiledFrame(allDay, :tile)))

max_stars_repo_name: gunnarnewell/weather-server

max_issues_repo_name: gunnarnewell/weather-server

max_forks_repo_name: gunnarnewell/weather-server


text:
abstract type AbstractProcessingStep end

jsrender(session::Session, step::AbstractProcessingStep) = jsrender(session, step.card)

columns_in(step::AbstractProcessingStep) = columns_in(step.card)
columns_out(step::AbstractProcessingStep) = columns_out(step.card)

@enum State inactive scheduled computing done errored

shouldrun(state::State) = state ∉ (inactive, done)

struct StateTracker
    state::Observable{State}
    edited::Observable{Bool}
end

function jsrender(session::Session, tracker::StateTracker)
    class = map(session, tracker.state, tracker.edited, result=Observable{String}()) do state, edited
        baseclass = "float-right text-2xl pr-4 inline-block"
        edited && return "$(baseclass) text-yellow-600"
        state == inactive && return "$(baseclass) text-transparent"
        state in (scheduled, computing) && return "$(baseclass) text-blue-600 animate-pulse"
        state == done && return "$(baseclass) text-blue-800"
        state == errored && return "$(baseclass) text-red-800"
        throw(ArgumentError("Invalid state $state"))
    end
    ui = DOM.span("⬤", class=class[])
    onjs(session, class, js"""
        function (className) {
            $(ui).className = className;
        }
    """)
    return jsrender(session, ui)
end

struct ErrorContainer
    hidden::Observable{Bool}
    error::Observable{String}
end

function jsrender(session::Session, ec::ErrorContainer)
    p = DOM.p(ec.error[])
    ui = DOM.div(p; ec.hidden, class="p-8 bg-white border-2 border-red-800")
    onjs(session, ec.error, js"""
        function (value) {
            $(p).innerText = value;
        }
    """)
    return jsrender(session, ui)
end

struct ProcessingCard
    name::Symbol
    inputs::RichTextField
    target::Union{RichTextField, Nothing}
    outputs::RichTextField
    method::Union{RichTextField, RichEditor}
    process_button::Button
    clear_button::Button
    state::Observable{State}
    edited::Observable{Bool}
    error::Observable{String}
    run::Observable{Bool}
    destroy::Observable{Bool}
end

function autocompletes(card::ProcessingCard)
    textfields = (card.inputs, card.target, card.outputs, card.method)
    return filter(!isnothing, textfields)
end

function process!(card::ProcessingCard)
    foreach(parse!, autocompletes(card))
    card.state[] = scheduled
    card.run[] = true
end

function clear!(card::ProcessingCard)
    foreach(reset!, autocompletes(card))
    card.state[] = scheduled
    card.run[] = true
end

function ProcessingCard(name;
                        inputs,
                        target=nothing,
                        outputs,
                        method,
                        process_button=Button("Process", class=buttonclass(true)),
                        clear_button=Button("Clear", class=buttonclass(false)),
                        state=Observable(inactive),
                        edited=Observable(false),
                        error=Observable(""),
                        run=Observable(false),
                        destroy = Observable(false))

    card = ProcessingCard(
        name,
        inputs,
        target,
        outputs,
        method,
        process_button,
        clear_button,
        state,
        edited,
        error,
        run,
        destroy
    )

    on(_ ->  process!(card), process_button.value)
    on(_ ->  clear!(card), clear_button.value)
    confirmed = map(autocompletes(card)) do textfield
        return lift(==, textfield.widget.value, textfield.confirmedvalue)
    end
    map!(!all∘tuple, card.edited, confirmed...)

    return card
end

function used_columns(args::AbstractVector{Call}...)
    colnames = Symbol[]
    for calls in args
        for call in calls
            append!(colnames, Symbol.(call.positional), Symbol.(last.(call.named)))
        end
    end
    return colnames
end

function columns_in(card::ProcessingCard)
    args = [t.parsed for t in [card.inputs, card.target] if !isnothing(t)]
    return used_columns(args...)
end

function columns_out(card::ProcessingCard)
    return isempty(columns_in(card)) ? Symbol[] : used_columns(card.outputs.parsed)
end

function jsrender(session::Session, card::ProcessingCard)
    statetracker = StateTracker(card.state, card.edited)
    hide_error = map(!=(errored), session, card.state, result=Observable{Bool}())
    errorcontainer = ErrorContainer(hide_error, card.error)

    card_ui = DOM.div(
        DOM.span(string(card.name), class="text-blue-800 text-2xl font-semibold"),
        DOM.span(
            "✕",
            class="text-red-800 hover:text-red-900 text-2xl font-semibold float-right cursor-pointer",
            onclick=js"JSServe.update_obs($(card.destroy), true)"
        ),
        statetracker,
        autocompletes(card)...,
        DOM.div(class="mt-12", card.process_button, card.clear_button),
        class="select-none p-8 shadow bg-white border-2 border-transparent",
        dataId=string(objectid(card)),
        dataType="card",
        dataSelected="false",
    )
    ui = DOM.div(
        card_ui,
        errorcontainer
    )

    return jsrender(session, ui)
end


text:
using Pkg

# Activate test environment on older Julia versions
@static if VERSION < v"1.2"
    Pkg.activate(@__DIR__)
    Pkg.develop(PackageSpec(; path=dirname(@__DIR__)))
    Pkg.instantiate()
end

using ConsistencyResampling
using Distributions
using Documenter
using LinearAlgebra
using Random
using StatsBase
using Test

Random.seed!(1234)

@testset "ConsistencyResampling" begin
    @testset "Alias sampler" begin
        include("alias.jl")
    end
    @testset "Consistent" begin
        include("consistent.jl")
    end
    @testset "Sampler (Inf)" begin
        include("many.jl")
    end
    @testset "Sampler (1)" begin
        include("single.jl")
    end

    @testset "doctests" begin
        DocMeta.setdocmeta!(
            ConsistencyResampling,
            :DocTestSetup,
            :(using ConsistencyResampling);
            recursive=true,
        )
        doctest(ConsistencyResampling)
    end
end


text:
# Autogenerated wrapper script for PDFHighlights_jll for x86_64-unknown-freebsd
export PDFHighlightsWrapper

using Poppler_jll
JLLWrappers.@generate_wrapper_header("PDFHighlights")
JLLWrappers.@declare_library_product(PDFHighlightsWrapper, "PDFHighlightsWrapper.so")
function __init__()
    JLLWrappers.@generate_init_header(Poppler_jll)
    JLLWrappers.@init_library_product(
        PDFHighlightsWrapper,
        "lib/PDFHighlightsWrapper.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "math" begin
    x = collect(1:100)
    y = collect(6:9)

    @test ATools.rxy(1,1) ≈ sqrt(2.)
    @test ATools.phixy(1,1) ≈ π/4

    @test ATools.dxyz([1,1,1], [2,2,2]) ≈ sqrt(3.)

    m, s = ATools.mean_std(x, 5, 10)
    @test m ≈ mean(y)
    @test s ≈ std(y)

    norm_weights = Int.(round.(1000 * exp.(-(x .- 50).^2 / 18)))
    m, s = ATools.mean_std(x, 1, 90, w=norm_weights)
    @test m ≈ 50.0
    @test isapprox(s, 3.0, atol=0.003)

    fxy = ATools.gline2p(0, 10., 1., 2.)
    @test fxy(0) == 10.0
    @test fxy(1) == 2.0
    @test fxy(0.5) == 6.0
end

max_stars_repo_name: andLaing/ATools

max_issues_repo_name: andLaing/ATools

max_forks_repo_name: andLaing/ATools


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
Base.peek(io::Base.AbstractPipe) = Base.peek(Base.pipe_reader(io))

using BinaryBuilder

name = "TermboxBuilder"
version = v"1.1.0"

# Collection of sources required to build TermboxBuilder
sources = [
    "https://github.com/nsf/termbox/archive/v1.1.0.tar.gz" =>
    "2743ee4aeb0ff39fadbaf945b76c43e6f6bba544156f2576282b775a3067d748",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/termbox-1.1.0/
./waf configure --prefix=$prefix
./waf
./waf install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, :glibc),
    Linux(:x86_64, :glibc),
    Linux(:aarch64, :glibc),
    Linux(:armv7l, :glibc, :eabihf),
    Linux(:powerpc64le, :glibc),
    Linux(:i686, :musl),
    Linux(:x86_64, :musl),
    Linux(:aarch64, :musl),
    Linux(:armv7l, :musl, :eabihf),
    MacOS(:x86_64),
    FreeBSD(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libtermbox", :libtermbox)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: jgoldfar/TermboxBuilder

max_issues_repo_name: jgoldfar/TermboxBuilder

max_forks_repo_name: jgoldfar/TermboxBuilder


text:
@safetestset "Routing edge cases" begin

  @safetestset "Emoji routing" begin
    using Genie
    using HTTP

    port = rand(8500:8900)

    route("/✔/🧞/♥/❤") do
      "/✔/🧞/♥/❤"
    end

    server = up(port)

    response = try
      HTTP.request("GET", "http://127.0.0.1:$port/✔/🧞/♥/❤")
    catch ex
      ex.response
    end

    @test response.status == 200
    @test String(response.body) == "/✔/🧞/♥/❤"

    down()
    sleep(1)
    server = nothing
  end;

  @safetestset "Emoji routing ✔" begin
    using Genie
    using HTTP

    port = rand(8500:8900)

    route("/✔") do
      "All good"
    end

    server = up(port)

    response = try
      HTTP.request("GET", "http://127.0.0.1:$port/✔")
    catch ex
      ex.response
    end

    @test response.status == 200
    @test String(response.body) == "All good"

    down()
    sleep(1)
    server = nothing
  end;

  @safetestset "Emoji routing with params" begin
    using Genie, Genie.Requests
    using HTTP

    port = rand(8500:8900)

    route("/:check/:genie/:smallheart/:bigheart") do
      "/$(payload(:check))/$(payload(:genie))/$(payload(:smallheart))/$(payload(:bigheart))"
    end

    server = up(port)

    response = try
      HTTP.request("GET", "http://127.0.0.1:$port/✔/🧞/♥/❤")
    catch ex
      ex.response
    end

    @test response.status == 200
    @test String(response.body) == "/✔/🧞/♥/❤"

    down()
    sleep(1)
    server = nothing
  end;

end


text:
function approxextrema(I)
    arr = sample(I)
    return (minimum(arr), maximum(arr))
end

function estimatequantile(array, p)
    try
        return quantile(sample(array; n=1e5), p)
    catch
        @warn "quantile could not be estimated! (maybe only NaNs)"
        return 0
    end
end

function sample(I; n=10000)
    n = min(n, length(I))
    len = ceil(Int, √n) # take len blocks of len elements
    startindices = round.(Int, range(firstindex(I) - 1, lastindex(I) - len; length=len))
    indices = vcat((i .+ (1:len) for i in startindices)...)
    ret = filter(isfinite, I[indices])
    if isempty(ret)
        ret = filter(isfinite, I)
    end
    return ret
end

function get_corner_indices(I; max_length=10)
    d = size(I)
    n = min.(max_length, ceil.(Int, d ./ 3)) # n voxels for each dim
    getrange(num, len) = [1:num, (len-num+1):len] # first and last voxels
    return collect(Iterators.product(getrange.(n, d)...))
end

# estimate noise parameters from corner without signal
function estimatenoise(image)
    corners = get_corner_indices(image)
    (lowestmean, ind) = findmin(mean.(filter(isfinite, image[I...]) for I in corners))
    sigma = std(filter(isfinite, image[corners[ind]...]))
    if isnan(sigma) # no corner available
        # estimation that is correct if half the image is signal and half noise
        sigma = 2estimatesigma_from_quantile(image, 1/4)
        lowestmean = sigma / 2
    end
    return lowestmean, sigma
end

# sigma is only calculated for quantile (biased)
function estimatesigma_from_quantile(image, quantile)
    q = estimatequantile(image, quantile)
    samples = filter(x -> x < q, sample(image))
    return std(samples)
end

function robustmask!(image; maskedvalue=if eltype(image) <: AbstractFloat NaN else 0 end)
    image[.!robustmask(image)] .= maskedvalue
    image
end
function robustmask(weight)
    μ, σ = estimatenoise(weight)
    m = mean(filter(isfinite, weight[weight .> 5σ]))
    maximum((5σ, m/5))
    mask = weight .> maximum((5σ, m/5))
    # remove small holes and minimally grow
    boxsizes=[[3,3] for i in 1:ndims(weight)]
    return gaussiansmooth3d(mask; nbox=2, boxsizes) .> 0.55
end

getcomplex(fnmag::AbstractString, fnphase::AbstractString) = getcomplex(niread(fnmag), niread(fnphase))

function getcomplex(mag, phase)
    minp, maxp = approxextrema(phase)
    mag .* exp.((2im * pi / (maxp - minp)) .* phase)
end

function readfromtextheader(filename, searchstring)
    for line in readlines(open(filename, "r"))
        if occursin(searchstring, line)
            # regex to search for "= " or ": " and return the following non-whitespace characters
            return match(r"(?<=(= |: ))(\S+)", line).match
        end
    end
end

# root sum of squares combination
RSS(mag; dims=ndims(mag)) = dropdims(.√sum(mag.^Float32(2); dims); dims)

function getscaledimage(array, div::Number, offset = 0, type::Symbol = :trans)
    array = reshape(array, size(array)[1:2]) # drops trailing singleton dimensions
    scaled = if offset != 0
        (array .- offset) .* (1 / div) .+ 0.5
    else
        array .* (1 / div)
    end
    scaled[isnan.(scaled) .| (scaled .< 0)] .= 0
    scaled[scaled .> 1] .= 1
    if type == :trans
        scaled = reverse(permutedims(scaled, [2 1]); dims = 1)
    else
    end
    scaled
end

function getscaledimage(array, type::Symbol = :trans)
    scaled = robustrescale(array, 0, 1, threshold=true)
    getscaledimage(scaled, 1, 0, type)
end

robustrescale(array, newmin, newmax; threshold=false, mask=trues(size(array)), datatype=Float64) =
    robustrescale!(datatype.(array), newmin, newmax; threshold, mask)

function robustrescale!(array, newmin, newmax; threshold=false, mask=trues(size(array)))
    array[isnan.(array)] .= minimum(filter(isfinite, array))
    q = [0.01, 0.99] # quantiles
    oldq = estimatequantile(array[mask], q)
    oldrange = (oldq[2] - oldq[1]) / (q[2] - q[1])
    oldmin = oldq[1] - q[1] * oldrange
    newrange = newmax - newmin

    array .= (array .- oldmin) .* (newrange / oldrange) .+ newmin

    if threshold
        array[array .< newmin] .= newmin
        array[array .> newmax] .= newmax
    end
    array
end

function rescale(array, newmin, newmax; datatype=eltype(array))
    rescale!(datatype.(array), newmin, newmax)
end

function rescale!(array, newmin, newmax)
    oldmin, oldmax = approxextrema(array)
    factor = (newmax - newmin) / (oldmax - oldmin)
    array .= (array .- oldmin) .* factor .+ newmin
end

to_dim(a::Real, dim::Int) = to_dim([a], dim)
to_dim(V::AbstractVector, dim::Int) = reshape(V, ones(Int, dim-1)..., :)

"""
    getHIP(mag, phase; echoes=[1,2])
return the hermitian inner product between the specified echoes.
"""
function getHIP(mag, phase; echoes=[1,2])
    e1, e2 = echoes
    compl = zeros(ComplexF64, size(mag)[1:3])
    for iCha in 1:size(mag, 5)
        compl .+= exp.(1.0im .* (phase[:,:,:,e2,iCha] .- phase[:,:,:,e1,iCha])) .* mag[:,:,:,e1,iCha] .* mag[:,:,:,e2,iCha]
    end
    compl
end

function getHIP(compl; echoes=[1,2])
    e1, e2 = echoes
    c = zeros(eltype(compl), size(compl)[1:3])
    for iCha in 1:size(compl, 5)
        c .+=  compl[:,:,:,e2,iCha] .* conj.(compl[:,:,:,e1,iCha])
    end
    return c
end


text:
using Test
using LinearAlgebra
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCell
include("useful_routines.jl")

polyorder = 1
numqp = required_quadrature_order(polyorder)
basis = TensorProductBasis(2,polyorder)
levelset = InterpolatingPolynomial(1,basis)
levelsetcoeffs = [-1.,-1.,1.,1.]
update!(levelset,levelsetcoeffs)

xL,xR = [-1.,-1.],[1.,1.]
pquad = area_quadrature(levelset,+1,xL,xR,numqp)
nquad = area_quadrature(levelset,-1,xL,xR,numqp)
squad = surface_quadrature(levelset,xL,xR,numqp)

lambda,mu = 1.,2.
penalty = 1.0
dx = 0.1
e11 = dx / 2.0
e22 = -lambda / (lambda + 2mu) * e11
dy = e22
penalty = 1.
stiffness = CutCell.plane_strain_voigt_hooke_matrix(lambda,mu)
cellmap = CutCell.CellMap([0.,0.],[2.,1.])

pbf = CutCell.bilinear_form(basis,pquad,stiffness,cellmap)
nbf = CutCell.bilinear_form(basis,nquad,stiffness,cellmap)

normals = repeat([1.,0.],1,numqp)
top = CutCell.coherent_traction_operator(basis,squad,squad,normals,stiffness,cellmap)

facescale = CutCell.scale_area(cellmap,normals)
mm = penalty*CutCell.interface_mass_matrix(basis,squad,squad,facescale)

eta = -1.0
BF = [pbf         zeros(8,8)
      zeros(8,8)  nbf       ]
T1 = -0.5*[-top  -top
           +top  +top]
T2 = eta*T1'

MM = [+mm  -mm
      -mm  +mm]

K = BF+T1+T2+MM
@test norm(K-K') > 1.0

rhs = zeros(16)

CutCell.apply_dirichlet_bc!(K,rhs,[5,6],1,0.,2)
CutCell.apply_dirichlet_bc!(K,rhs,[5],2,0.,2)
CutCell.apply_dirichlet_bc!(K,rhs,[3,4],1,dx,2)

sol = K\rhs
disp = reshape(sol,2,:)

testdisp = [0.0 0.0 dx dx 0.0 0.0 dx dx
            0.0 dy  0. dy 0.  dy  0. dy]

@test allapprox(disp,testdisp,1e2*eps())

eta = 1.0
BF = [pbf         zeros(8,8)
      zeros(8,8)  nbf       ]
T1 = -0.5*[-top  -top
           +top  +top]
T2 = eta*T1'

MM = [+mm  -mm
      -mm  +mm]

K = BF+T1+T2+MM
@test isapprox(norm(K-K'),0.0,atol=1e1eps())

rhs = zeros(16)

CutCell.apply_dirichlet_bc!(K,rhs,[5,6],1,0.,2)
CutCell.apply_dirichlet_bc!(K,rhs,[5],2,0.,2)
CutCell.apply_dirichlet_bc!(K,rhs,[3,4],1,dx,2)

sol = K\rhs
disp = reshape(sol,2,:)

testdisp = [0.0 0.0 dx dx 0.0 0.0 dx dx
            0.0 dy  0. dy 0.  dy  0. dy]

@test allapprox(disp,testdisp,1e2*eps())

eta = 0.0
BF = [pbf         zeros(8,8)
      zeros(8,8)  nbf       ]
T1 = -0.5*[-top  -top
           +top  +top]
T2 = eta*T1'

MM = [+mm  -mm
      -mm  +mm]

K = BF+T1+T2+MM
@test norm(K-K') > 1.0

rhs = zeros(16)

CutCell.apply_dirichlet_bc!(K,rhs,[5,6],1,0.,2)
CutCell.apply_dirichlet_bc!(K,rhs,[5],2,0.,2)
CutCell.apply_dirichlet_bc!(K,rhs,[3,4],1,dx,2)

sol = K\rhs
disp = reshape(sol,2,:)

testdisp = [0.0 0.0 dx dx 0.0 0.0 dx dx
            0.0 dy  0. dy 0.  dy  0. dy]

@test allapprox(disp,testdisp,1e2*eps())


text:
for PT in (:Particles, :StaticParticles)
    @eval begin
        Base.promote_rule(::Type{Complex{S}}, ::Type{$PT{T,N}}) where {S<:Real,T<:Real,N} = Complex{$PT{promote_type(S,T),N}}

    end

end

"""
    ℂ2ℂ_function(f::Function, z::Complex{<:AbstractParticles})
applies `f : ℂ → ℂ ` to `z::Complex{<:AbstractParticles}`.
"""
function ℂ2ℂ_function(f::F, z::Complex{T}) where {F,T<:AbstractParticles}
    rz,iz = z.re,z.im
    s = map(1:length(rz.particles)) do i
        @inbounds f(complex(rz[i], iz[i]))
    end
    complex(T(real.(s)), T(imag.(s)))
end

function ℂ2ℂ_function!(f::F, s, z::Complex{T}) where {F,T<:AbstractParticles}
    rz,iz = z.re,z.im
    map!(s, 1:length(rz.particles)) do i
        @inbounds f(complex(rz[i], iz[i]))
    end
    complex(T(real.(s)), T(imag.(s)))
end

for ff in (sqrt, exp, sin, cos)
    f = nameof(ff)
    @eval Base.$f(z::Complex{<: AbstractParticles}) = ℂ2ℂ_function($f, z)
    @eval $(Symbol(f,:!))(s, z::Complex{<: AbstractParticles}) = ℂ2ℂ_function!($f, s, z)
end

function Base.:(/)(a::Complex{T}, b::Complex{T}) where T<:AbstractParticles
    are = real(a); aim = imag(a); bre = real(b); bim = imag(b)
    if mean(abs(bre)) <= mean(abs(bim))
        if isinf(bre) && isinf(bim)
            r = sign(bre)/sign(bim)
        else
            r = bre / bim
        end
        den = bim + r*bre
        Complex((are*r + aim)/den, (aim*r - are)/den)
    else
        if isinf(bre) && isinf(bim)
            r = sign(bim)/sign(bre)
        else
            r = bim / bre
        end
        den = bre + r*bim
        Complex((are + aim*r)/den, (aim - are*r)/den)
    end
end


text:
using Base.Test
using CGP
CGP.Config.init("cfg/test.yaml")

@testset "Creation tests" begin
    for ct in CGP.chromosomes
        println(ct)
        nin = rand(1:100); nout = rand(1:100);
        c = ct(nin, nout)
        cgenes = deepcopy(c.genes)
        @testset "Simple $ct" begin
            @test all(c.genes .<= 1.0)
            @test all(c.genes .>= -1.0)
            @test c.nin == nin
            @test c.nout == nout
        end
        @testset "Genes $ct" begin
            newgenes = rand(size(c.genes))
            d = ct(newgenes, nin, nout)
            @test c != d
            @test c.genes != d.genes
        end
        @testset "Gene equality $ct" begin
            for i=1:10
                d = ct(cgenes, nin, nout)
                @test c.genes == d.genes
                @test c.outputs == d.outputs
                @test length(c.nodes) == length(d.nodes)
                @test [n.connections for n in c.nodes] == [n.connections for n in d.nodes]
                @test [n.f for n in c.nodes] == [n.f for n in d.nodes]
                @test [n.output for n in c.nodes] == [n.output for n in d.nodes]
                @test [n.active for n in c.nodes] == [n.active for n in d.nodes]
            end
        end
        @testset "Gene reconstruction $ct" begin
            genes = get_genes(c, collect((c.nin+1):length(c.nodes)))
            @test length(genes) == (length(cgenes) - c.nin - c.nout)
            @test genes == cgenes[(c.nin+c.nout+1):end]
        end
    end
end

@testset "Functional tests" begin
    for ct in CGP.chromosomes
        println(ct)
        nin = rand(1:100); nout = rand(1:100)
        @testset "Process $ct" begin
            c = ct(nin, nout)
            genecopy = deepcopy(c.genes)
            funccopy = deepcopy([n.f for n in c.nodes])
            conncopy = deepcopy([n.connections for n in c.nodes])
            activecopy = deepcopy([n.active for n in c.nodes])
            for i in 1:10
                out = process(c, rand(nin))
                @test all(genecopy .== c.genes)
                @test length(out) == nout
                @test all(out .<= 1.0)
                @test all(out .>= -1.0)
                @test [n.f for n in c.nodes] == funccopy
                @test [n.connections for n in c.nodes] == conncopy
                @test [n.active for n in c.nodes] == activecopy
            end
        end
        @testset "Process equality $ct" begin
            c = ct(nin, nout)
            d = ct(deepcopy(c.genes), nin, nout)
            for i in 1:10
                inp = rand(nin)
                cout = process(c, inp)
                dout = process(d, inp)
                @test [n.connections for n in c.nodes] == [n.connections for n in d.nodes]
                @test [n.f for n in c.nodes] == [n.f for n in d.nodes]
                @test [n.active for n in c.nodes] == [n.active for n in d.nodes]
                @test [n.output for n in c.nodes] == [n.output for n in d.nodes]
                @test all(cout == dout)
            end
        end
    end
end

nothing


text:
module QFT

export gen_qft, gen_iqft
using ..GateOps

"""
    gen_qft(qubit_indices::Vector{Int})

Generates an OpenQASM representation of the quantum Fourier transform
for the given number of qubits.

"""
function gen_qft(qubit_indices::Vector{Int})
    circuit = Vector{GateOps.GateCall}()

    angle(k) = 2*pi / 2^k;
    for (iidx,ival) = Iterators.reverse(enumerate(qubit_indices))
        push!(circuit, GateOps.hadamard(qubit_indices[iidx]))
        for (jidx,jval) = Iterators.reverse(enumerate(qubit_indices[1:iidx-1]))
            theta = angle(iidx - jidx);
            push!(circuit, GateOps.c_r_phase(ival, jval, theta))
        end
    end
    for idx = 1:div(length(qubit_indices),2)
        push!(circuit, GateOps.swap(qubit_indices[idx], qubit_indices[length(qubit_indices)-idx+1]))
    end
    return circuit
end

"""
    gen_iqft(qubit_indices::Vector{Int})

Generates an OpenQASM representation of the inverse quantum Fourier transform
for the given number of qubits.

"""
function gen_iqft(qubit_indices::Vector{Int})
    circuit = Vector{GateOps.GateCall}()

    angle(k) = -2*pi / 2^k;
    for (iidx,ival) = enumerate(qubit_indices)
        for (jidx,jval) = enumerate(qubit_indices[1:iidx-1])
            theta = angle(iidx - jidx);
            push!(circuit, GateOps.c_r_phase(ival, jval, theta))
        end
        circuit *= GateOps.apply_gate_h(q_reg, qubit_indices[iidx])

    end
    for idx = 1:div(length(qubit_indices),2)
        push!(circuit, GateOps.swap(qubit_indices[idx], qubit_indices[length(qubit_indices)-idx+1]))
    end
    return circuit
end

end


text:
# general method given a reachability algorithm for the linear system
function reach_CARLIN_alg(X0, F1, F2; alg, resets, N, T, Δt0, bloat, compress)
    if resets == 0
        reach_CARLIN(X0, F1, F2; alg, N, T, Δt=Δt0, bloat, compress)
    else
        reach_CARLIN_resets(X0, F1, F2, resets; alg, N, T, Δt=Δt0, bloat, compress)
    end
end

function reach_CARLIN(X0, F1, F2; alg, N, T, Δt, bloat, compress, A=nothing)
    error_bound_func = error_bound_specabs

    # lift initial states
    if compress
        ŷ0 = lift_vector(X0, N) # see CarlemanLinearization/linearization.jl
    else
        ŷ0 = kron_pow_stack(X0, N) |> box_approximation
    end

    # solve continuous ODE
    if isnothing(A)
        A = build_matrix(F1, F2, N; compress)
    end
    prob = @ivp(ŷ' = Aŷ, ŷ(0) ∈ ŷ0)
    sol = solve(prob, T=T, alg=alg, Δt0=Δt)

    # projection onto the first n variables
    n = dim(X0)
    πsol_1n = _project(sol, 1:n)

    if !bloat
        return πsol_1n
    end

    # compute errors
    errfunc = error_bound_func(X0, Matrix(F1), Matrix(F2), N=N)

    # evaluate error bounds for each reach-set in the solution
    E = [errfunc.(tspan(R)) for R in sol]

    # if the interval is always > 0 then we can just take max(Ei)

    # symmetrize intervals
    E_rad = [symmetric_interval_hull(Interval(ei)) for ei in E]
    E_ball = [BallInf(zeros(n), max(ei)) for ei in E_rad]

    # sum the solution with the error
    fp_bloated = [ReachSet(set(Ri) ⊕ Ei, tspan(Ri)) for (Ri, Ei) in zip(πsol_1n, E_ball)] |> Flowpipe

    return fp_bloated
end

function _compute_resets(resets::Int, T)
    return mince(0 .. T, resets+1)
end

function _compute_resets(resets::Vector{Float64}, T)
    # assumes initial time is 0
    aux = vcat(0.0, resets, T)
    return [interval(aux[i], aux[i+1]) for i in 1:length(aux)-1]
end

function reach_CARLIN_resets(X0, F1, F2, resets; alg, N, T, Δt, bloat, compress)

    # build state matrix (remains unchanged upon resets)
    A = build_matrix(F1, F2, N; compress)

    # time intervals to compute
    time_intervals = _compute_resets(resets, T)

    # compute until first chunk
    T1 = sup(first(time_intervals))
    sol_1 = reach_CARLIN(X0, F1, F2; alg, N, T=T1, Δt=interval(0)+Δt, bloat, compress, A=A)

    # preallocate output flowpipe
    fp_1 = flowpipe(sol_1)
    out = Vector{typeof(fp_1)}()
    push!(out, fp_1)

    # approximate final reach-set
    Rlast = sol_1[end]
    X0 = box_approximation(set(Rlast))

    # compute remaining chunks
    for i in 2:length(time_intervals)
        T0 = T1
        Ti = sup(time_intervals[i])
        sol_i = reach_CARLIN(X0, F1, F2; alg, N, T=Ti-T0, Δt=interval(T0)+Δt, bloat, compress, A=A)
        push!(out, flowpipe(sol_i))
        X0 = box_approximation(set(sol_i[end]))
        T1 = Ti
    end

    return MixedFlowpipe(out)
end


text:
# Tests for Profit Profit Efficiency Enhanced Russell Graph Slack Based Measure DEA Models
@testset "ProfitERGDEAModel" begin

    # Test using Book data
    X = [2; 4; 8; 12; 6; 14; 14; 9.412];
    Y = [1; 5; 8; 9; 3; 7; 9; 2.353];
    W = [1; 1; 1; 1; 1; 1; 1; 1];
    P = [2; 2; 2; 2; 2; 2; 2; 2];

    # Profit Enhanced Russell Graph Slack Based Measure
    profiterg = deaprofiterg(X, Y, W, P)

    @test typeof(profiterg) == ProfitERGDEAModel

    @test nobs(profiterg) == 8
    @test ninputs(profiterg) == 1
    @test noutputs(profiterg) == 1
    @test ismonetary(profiterg) == false

    @test efficiency(profiterg) == efficiency(profiterg, :Economic)
    @test efficiency(profiterg, :Economic)   ≈ [4.0; 0.5; 0.0; 1/6; 0.8; 4/7; 2/7; 1.2706] atol = 1e-5
    @test efficiency(profiterg, :Technical)  ≈ [0.0; 0.0; 0.0; 0.0; 0.6; 11/21; 1/7; 0.8] atol = 1e-5
    @test efficiency(profiterg, :Allocative) ≈ [4.0; 0.5; 0.0; 1/6; 0.2; 1/21; 1/7; 0.4706] atol = 1e-5
    @test normfactor(profiterg) ≈ [2.0; 4.0; 8.0; 12.0; 10.0; 14.0; 14.0; 10.0]

    # With monetary = true
    profitergmon = deaprofiterg(X, Y, W, P, monetary = true)

    @test ismonetary(profitergmon) == true

    @test efficiency(profitergmon) ≈ efficiency(profiterg) .* normfactor(profiterg)
    @test efficiency(profitergmon, :Technical) ≈ efficiency(profiterg, :Technical) .* normfactor(profiterg)
    @test efficiency(profitergmon, :Allocative) ≈ efficiency(profiterg, :Allocative) .* normfactor(profiterg)

    # Print
    show(IOBuffer(), profiterg)

    # Test errors
    @test_throws DimensionMismatch deaprofiterg([1; 2 ; 3], [4; 5], [1; 1; 1], [1; 1]) #  Different number of observations
    @test_throws DimensionMismatch deaprofiterg([1; 2; 3], [4; 5; 6], [1; 2; 3; 4], [1; 1; 1]) # Different number of observation in input prices
    @test_throws DimensionMismatch deaprofiterg([1; 2; 3], [4; 5; 6], [1; 1; 1], [1; 2; 3; 4]) # Different number of observation in output prices
    @test_throws DimensionMismatch deaprofiterg([1 1; 2 2; 3 3 ], [4; 5; 6], [1 1 1; 2 2 2; 3 3 3], [1; 1; 1]) # Different number of input prices and inputs
    @test_throws DimensionMismatch deaprofiterg([1; 2; 3 ], [4; 5; 6], [1; 1; 1], [1 1 1; 2 2 2; 3 3 3]) # Different number of output prices and outputs

end


text:
function hpipm_dims()
    py_qp_dims = hp.hpipm_dims()
    ju_qp_dims = py2ju_obj(py_qp_dims)
    return ju_qp_dims
end


text:
# routes.jl
using Genie.Router
using MyLib

route("/") do
  MyLib.example()
end


text:
#===============================================================================
 Robust NonMonotone Line search
 From DF-Sane paper https://www.ams.org/journals/mcom/2006-75-255/S0025-5718-06-01840-0/S0025-5718-06-01840-0.pdf
 This line search doesn't fit the rest of the API so I need to come up with an
 interface for nonmonotone line searches (it has to accept an fbar (or Q in Hager
 Zhang's notion) and a forcing term)
===============================================================================#
struct RNMS end

function find_steplength(rnms::RNMS, φ, φ0, fbar, ηk::T, τmin, τmax) where T
  α₊ = T(1)
  α₋ = α₊
  γ = T(1)/10000
  for k = 1:100
    φα₊ = φ(α₊)
    if φα₊ ≤ fbar + ηk - γ*α₊^2*φ0
      return α₊, φα₊
    end
    φα₋ = φ(-α₋)
    if φα₋ ≤ fbar + ηk - γ*α₋^2*φ0
      return -α₋, φα₋
    end

    # update alpha+ and alpha-
    αt = α₊^2*φ0/(φα₊ + (2*α₊-1)*φ0)
    α₊ = clamp(αt, τmin*α₊, τmax*α₊)

    αt = α₋^2*φ0/(φα₋ + (2*α₋-1)*φ0)
    α₋ = clamp(αt, τmin*α₋, τmax*α₋)
  end
  throw(ErrorException("Line search failed."))
end

function safeguard_σ(σ::T, σmin, σmax, F) where T
  if abs(σ) < σmin || abs(σ) > σmax
    normF = norm(F)
    if normF > T(1)
      σ = T(1)
    elseif T(1)/10^5 ≤ normF ≤ T(1)
      σ = inv(normF)
    elseif normF < T(1)/10^5
      σ = T(10)^5
    end
  end
  return σ
end

struct DFSANE end
function solve(prob::NEqProblem, x0, method::DFSANE, options::NEqOptions)
  if !(mstyle(prob) === InPlace())
    throw(ErrorException("solve() not defined for OutOfPlace() with DFSANE"))
  end

  t0 = time()
  F = prob.R.F
  T = eltype(x0)

  σmin, σmax = 1e-5, 1e5
  τmin, τmax = T(1)/10, T(1)/2
  M = 4
  nexp = 2

  x = copy(x0)
  ρs = norm(x)
  Fx = copy(x0)
  Fx = F(Fx, x)
  y = copy(Fx)
  ρ2F0 = norm(Fx, 2)
  ρF0 = norm(Fx, Inf)
  fx = ρ2F0^nexp
  fvals = [fx]

  abstol = 1e-5
  reltol = 1e-8

  σ₀ = T(1)
  σ = safeguard_σ(σ₀, σmin, σmax, Fx)
  iter = 0
  while iter < options.maxiter
    iter += 1
    z = copy(x) # FIXME
    push!(fvals, fx)
    if length(fvals) > M
      popfirst!(fvals)
    end
    fbar = maximum(fvals)
    d = -σ*Fx
    ηk = ρ2F0/(1+iter)^2
    φ(α) = norm(F(Fx, x.+α.*d))^nexp
    φ0 = fx
    α, φα = find_steplength(RNMS(), φ, φ0, fbar, ηk, τmin, τmax)
    s = α*d
    ρs = norm(s)
    x .= x.+s
    y .= -Fx
    Fx = F(Fx, x)
    y .+= Fx
    ρ2fx = norm(Fx, 2)
    ρfx = norm(Fx ,Inf)
    fx = ρ2fx^nexp
    σ = norm(s)^2/dot(s, y)
    σ = safeguard_σ(σ, σmin, σmax, Fx)
    # use sqrt(length(x))*abs or abstol?
    if ρfx < abstol + reltol*ρF0
      break
    end
  end
  ConvergenceInfo(method, (solution=x, best_residual=Fx, ρF0=ρF0, ρ2F0=ρ2F0, ρs=ρs, iter=iter, time=time()-t0), options)
end


text:
module HWconstrained

greet() = print("Hello World!")

	using JuMP, NLopt, DataFrames, Ipopt
	using LinearAlgebra

	export data, table_NLopt, table_JuMP

	function data(a=0.5)
        na = 3
		nc = 1
		ns = 4
		nss = 16
		e = [2,0,0]
		p = [1,1,1]
		z_2 = [0.72,0.92,1.12,1.32]
		z_3 = [0.86,0.96,1.06,1.16]
		z = [[1.0,z_2[j],z_3[k]] for j in 1:4 for k in 1:4]
		pi = 1/16*ones(1,16)
		return Dict("a"=>a,"na"=>na,"nc"=>nc,"ns"=>ns,"nss"=>nss,"e"=>e,"p"=>p,"z"=>z,"pi"=>pi)
	end

	function maxJuMP(a=0.5)
	        d = data(a)
			n = Model(with_optimizer(Ipopt.Optimizer))
			@variable(n, omega[1:3])
			@variable(n, c)
			@NLobjective(n, Min, exp(-a*c)+sum(d["pi"][i]*exp((-a)*sum(d["z"][i][j]'*omega[j] for j in 1:d["na"])) for i in 1:d["nss"]))
			@constraint(n, c + d["p"]'*(omega - d["e"]) <= 0)
			JuMP.optimize!(n)
	  	return Dict("obj"=>objective_value(n),"c"=>value(c),"omegas"=>[value(omega[i]) for i in 1:length(omega)])
	end

	function table_JuMP()
		d = DataFrame(a=[0.5;1.0;5.0],c = zeros(3),omega1=zeros(3),omega2=zeros(3),omega3=zeros(3),fval=zeros(3))
		for i in 1:nrow(d)
			xx = maxJuMP(d[i,:a])
			d[i,:c] = xx["c"]
			d[i,:omega1] = xx["omegas"][1]
			d[i,:omega2] = xx["omegas"][2]
			d[i,:omega3] = xx["omegas"][3]
			d[i,:fval] = xx["obj"]
		end
		return d
	end

	function obj(x::Vector,d::Dict)
		a = d["a"]
		z = d["z"]
		pi = d["pi"]
		if length(grad)>0
			  grad[1] = sum(d["pi"][i]*(-a)*exp((-a)*z[i][1]*sum(z[i]'*x[1:3])) for i in 1:d["nss"])
			  grad[2] = sum(d["pi"][i]*(-a)*exp((-a)*z[i][2]*sum(z[i]'*x[1:3])) for i in 1:d["nss"])
			  grad[3] = sum(d["pi"][i]*(-a)*exp((-a)*z[i][3]*sum(z[i]'*x[1:3])) for i in 1:d["nss"])
			  grad[4] = (-a)*exp.(-a*x[4])
		  end
		return (exp(-a*x[4])+sum(d["pi"][i]*exp((-a)*sum(z[i]'*x[1:3])) for i in 1:d["nss"]))
	  end

	function constr(x::Vector,grad::Vector,d::Dict)
		if length(grad)>0
			   grad[1] = d["p"][1]
			   grad[2] = d["p"][2]
			   grad[3] = d["p"][3]
			   grad[4] = 1
	   end
	   return x[4] + d["p"]'*(x[1:3] - d["e"])
	end

	function max_NLopt(a=0.5)
		count = 0 # keep track of # function evaluations
		opt = Opt(:LD_MMA, 4)
		lower_bounds!(opt, [-Inf, -Inf, -Inf, 0.])
		xtol_rel!(opt, 1e-7)
		grad = zeros(4)
		println(length(grad))
		d = data(a)
		min_objective!(opt,(x,grad) -> obj(x,grad,d))
		inequality_constraint!(opt, (x,grad) -> constr(x,grad,d))
		ftol_rel!(opt,1e-9)
		NLopt.optimize(opt, [0,0,0,0])
end

	function table_NLopt()
		d = DataFrame(a=[0.5;1.0;5.0],c = zeros(3),omega1=zeros(3),omega2=zeros(3),omega3=zeros(3),fval=zeros(3))
		for i in 1:nrow(d)
			xx = max_NLopt(d[i,:a])
			for j in 2:ncol(d)-1
				d[i,j] = xx[2][j-1]
			end
			d[i,end] = xx[1]
		end
		return d
	end

end # module


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Julia 1.1.0
#     language: julia
#     name: julia-1.1
# ---

# # This notebook
#
# _Notes:_ For documentation see <https://gaelforget.github.io/MeshArrays.jl/stable/>, <https://docs.juliadiffeq.org/latest/solvers/ode_solve.html> and <https://en.wikipedia.org/wiki/Displacement_(vector)>

# ## 1. import software

using IndividualDisplacements, MeshArrays, DifferentialEquations, Plots, Statistics
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"PlotIndDisp.jl"))

# ## 2. Read gridded variables as `MeshArray`s

# Put grid variables in a dictionary.

# +
mygrid=gcmgrid("llc90_latlon/","ll",1,[(360,178)], [360 178], Float32, read, write)

GridVariables=Dict("XC" => read(mygrid.path*"XC.latlon.data",MeshArray(mygrid,Float32)),
"YC" => read(mygrid.path*"YC.latlon.data",MeshArray(mygrid,Float32)),
"XG" => read(mygrid.path*"XG.data",MeshArray(mygrid,Float32)),
"YG" => read(mygrid.path*"YG.data",MeshArray(mygrid,Float32)),
"DXC" => read(mygrid.path*"DXC.latlon.data",MeshArray(mygrid,Float32)),
"DYC" => read(mygrid.path*"DYC.latlon.data",MeshArray(mygrid,Float32)) );
# -

# Read velocity fields as `MeshArray`s.

# +
using MAT
file = matopen(mygrid.path*"uv_lonlat.mat")
u=read(file, "u")
v=read(file, "v")
close(file)

u=dropdims(mean(u,dims=3),dims=3)
v=dropdims(mean(v,dims=3),dims=3)

u=read(u,MeshArray(mygrid,Float32))
v=read(v,MeshArray(mygrid,Float32));

u[findall(isnan.(u))]=0.0
v[findall(isnan.(v))]=0.0

u0=u; u1=u;
v0=v; v1=v;
# -

# Put velocity fields and time range in a dictionary.

# +
t0=0.0; t1=86400*366*2.0; dt=3600;

u0=u0./GridVariables["DXC"]
u1=u1./GridVariables["DXC"]
v0=v0./GridVariables["DYC"]
v1=v1./GridVariables["DYC"]

uvt = Dict("u0" => u0, "u1" => u1, "v0" => v0, "v1" => v1, "t0" => t0, "t1" => t1, "dt" => dt) ;
# -

# Merge the two dictionaries and add masks

# +
uvetc=merge(uvt,GridVariables);

nr=50; kk=1;

mskW=read(mygrid.path*"hFacW.latlon.data",MeshArray(mygrid,Float32,nr))
mskW=1.0 .+ 0.0 * mask(mskW[:,kk],NaN,0.0)
mskS=read(mygrid.path*"hFacS.latlon.data",MeshArray(mygrid,Float32,nr))
mskS=1.0 .+ 0.0 * mask(mskS[:,kk],NaN,0.0)

msk=Dict("mskW" => mskW, "mskS" => mskS)

uvetc=merge(uvetc,msk);
# -

# ## 3. Compute trajectories from gridded flow fields

comp_vel=IndividualDisplacements.VelComp
get_vel=IndividualDisplacements.VelCopy

# +
ii1=0:5:360; ii2=20:2:150;

n1=length(ii1); n2=length(ii2);
uInitS=Array{Float64,2}(undef,(2,n1*n2))
for i1 in eachindex(ii1); for i2 in eachindex(ii2);
        i=i1+(i2-1)*n1
        uInitS[1,i]=ii1[i1]
        uInitS[2,i]=ii2[i2]
end; end;
du=fill(0.0,size(uInitS));
comp_vel(du,uInitS,uvetc,0.0)
du
# -

tspan = (0.0,uvt["t1"]-uvt["t0"])
prob = ODEProblem(comp_vel,uInitS,tspan,uvetc)
sol = solve(prob,Tsit5(),reltol=1e-4,abstol=1e-4)
size(sol)

# ## 4. Store trajectories in a `DataFrame` and plot

ID=collect(1:size(sol,2))*ones(1,size(sol,3))
lon=mod.(sol[1,:,:],360)
lat=mod.(sol[2,:,:],180)
df = DataFrame(ID=Int.(ID[:]), lon=lon[:], lat=lat[:])
size(df)

PyPlot.figure(); PlotBasic(df,size(sol,2),90.0)