text:
using ImplicitDomainQuadrature
import Base.==, Base.≈

function allapprox(v1, v2)
    return all(v1 .≈ v2)
end

function allapprox(v1, v2, atol)
    err = maximum(abs.(v1 - v2))
    # if 10err > atol
    #     @warn "Approximation is very tight!"
    # end
    return length(v1) == length(v2) &&
           all([isapprox(v1[i], v2[i], atol = atol) for i = 1:length(v1)])
end

function allequal(v1, v2)
    return all(v1 .== v2)
end

function Base.isequal(c1::CutCellDG.CellMap, c2::CutCellDG.CellMap)
    return allequal(c1.yL, c2.yL) && allequal(c1.yR, c2.yR)
end

function ==(c1::CutCellDG.CellMap, c2::CutCellDG.CellMap)
    return isequal(c1, c2)
end

function plane_distance_function(coords, normal, x0)
    return (coords .- x0)' * normal
end

function closest_point_on_plane(querypoints,normal,x0)
    normalcomp = vec(normal'*(querypoints .- x0))
    disp = hcat([c*normal for c in normalcomp]...)
    return querypoints - disp
end

function circle_distance_function(coords, center, radius)
    difference = (coords .- center) .^ 2
    distance = radius .- sqrt.(vec(mapslices(sum, difference, dims = 1)))
    return distance
end

function closest_point_on_arc(querypoints,center,radius)
    relpos = querypoints .- center
    cqpoints = relpos[1,:] + im*relpos[2,:]
    theta = angle.(cqpoints)
    v = vcat(cos.(theta)',sin.(theta)')
    return center .+ radius*v
end

function corner_distance_function(x::V, xc) where {V<:AbstractVector}
    v = xc - x
    if all(x .<= xc)
        minimum(v)
    elseif all(x .> xc)
        return -sqrt(v' * v)
    elseif x[2] > xc[2]
        return v[2]
    else
        return v[1]
    end
end

function corner_distance_function(points::M, xc) where {M<:AbstractMatrix}
    return vec(
        mapslices(x -> corner_distance_function(x, xc), points, dims = 1),
    )
end

function ≈(q1::QuadratureRule, q2::QuadratureRule)
    flag = allapprox(q1.points, q2.points) && allapprox(q1.weights, q2.weights)
end

function required_quadrature_order(polyorder)
    ceil(Int, 0.5 * (2polyorder + 1))
end

function add_cell_error_squared!(
    err,
    interpolater,
    exactsolution,
    cellmap,
    quad,
)
    detjac = CutCellDG.determinant_jacobian(cellmap)
    for (p, w) in quad
        numsol = interpolater(p)
        exsol = exactsolution(cellmap(p))
        err .+= (numsol - exsol) .^ 2 * detjac * w
    end
end

function mesh_L2_error(nodalsolutions, exactsolution, basis, cellquads, mesh)
    ndofs = size(nodalsolutions)[1]
    err = zeros(ndofs)
    interpolater = InterpolatingPolynomial(ndofs, basis)
    ncells = CutCellDG.number_of_cells(mesh)
    for cellid = 1:ncells
        cellsign = CutCellDG.cell_sign(mesh, cellid)
        CutCellDG.check_cellsign(cellsign)
        if cellsign == +1 || cellsign == 0
            cellmap = CutCellDG.cell_map(mesh, +1, cellid)
            nodeids = CutCellDG.nodal_connectivity(mesh, +1, cellid)
            elementsolution = nodalsolutions[:, nodeids]
            update!(interpolater, elementsolution)
            quad = cellquads[+1, cellid]
            add_cell_error_squared!(
                err,
                interpolater,
                exactsolution,
                cellmap,
                quad,
            )
        end
        if cellsign == -1 || cellsign == 0
            cellmap = CutCellDG.cell_map(mesh, -1, cellid)
            nodeids = CutCellDG.nodal_connectivity(mesh, -1, cellid)
            elementsolution = nodalsolutions[:, nodeids]
            update!(interpolater, elementsolution)
            quad = cellquads[-1, cellid]
            add_cell_error_squared!(
                err,
                interpolater,
                exactsolution,
                cellmap,
                quad,
            )
        end
    end
    return sqrt.(err)
end

function cellwise_L2_error(
    nodalsolutions,
    exactsolution,
    basis,
    cellquads,
    mesh,
)

    ndofs = size(nodalsolutions)[1]
    interpolater = InterpolatingPolynomial(ndofs, basis)
    ncells = CutCellDG.number_of_cells(mesh)
    err = zeros(ndofs, 2, ncells)

    for cellid = 1:ncells
        cellsign = CutCellDG.cell_sign(mesh, cellid)
        CutCellDG.check_cellsign(cellsign)
        if cellsign == +1 || cellsign == 0
            cellmap = CutCellDG.cell_map(mesh, +1, cellid)
            nodeids = CutCellDG.nodal_connectivity(mesh, +1, cellid)
            elementsolution = nodalsolutions[:, nodeids]
            update!(interpolater, elementsolution)
            quad = cellquads[+1, cellid]
            cellerr = zeros(ndofs)
            add_cell_error_squared!(
                cellerr,
                interpolater,
                exactsolution,
                cellmap,
                quad,
            )
            err[:, 1, cellid] = cellerr
        end
        if cellsign == -1 || cellsign == 0
            cellmap = CutCellDG.cell_map(mesh, -1, cellid)
            nodeids = CutCellDG.nodal_connectivity(mesh, -1, cellid)
            elementsolution = nodalsolutions[:, nodeids]
            update!(interpolater, elementsolution)
            quad = cellquads[-1, cellid]
            cellerr = zeros(ndofs)
            add_cell_error_squared!(
                cellerr,
                interpolater,
                exactsolution,
                cellmap,
                quad,
            )
            err[:, 2, cellid] = cellerr
        end
    end
    return sqrt.(err)
end

function average(v)
    return sum(v)/length(v)
end

function convergence_rate(dx, err)
    return diff(log.(err)) ./ diff(log.(dx))
end

function add_cell_norm_squared!(vals, func, cellmap, quad)
    detjac = CutCellDG.determinant_jacobian(cellmap)
    for (p, w) in quad
        v = func(cellmap(p))
        vals .+= v .^ 2 * detjac * w
    end
end

function integral_norm_on_mesh(func, cellquads, mesh, ndofs)
    vals = zeros(ndofs)
    ncells = CutCellDG.number_of_cells(mesh)
    for cellid = 1:ncells
        s = CutCellDG.cell_sign(mesh, cellid)
        @assert s == -1 || s == 0 || s == 1
        if s == 1 || s == 0
            cellmap = CutCellDG.cell_map(mesh, +1, cellid)
            pquad = cellquads[+1, cellid]
            add_cell_norm_squared!(vals, func, cellmap, pquad)
        end
        if s == -1 || s == 0
            cellmap = CutCellDG.cell_map(mesh, -1, cellid)
            nquad = cellquads[-1, cellid]
            add_cell_norm_squared!(vals, func, cellmap, nquad)
        end
    end
    return sqrt.(vals)
end

function add_interface_error_squared!(
    err,
    interpolater,
    exactsolution,
    cellmap,
    quad,
    facescale,
)
    @assert length(quad) == length(facescale)
    for (i, (p, w)) in enumerate(quad)
        numsol = interpolater(p)
        exsol = exactsolution(cellmap(p))
        err .+= (numsol - exsol) .^ 2 * facescale[i] * w
    end
end

function interface_L2_error(
    nodalsolutions,
    exactsolution,
    levelsetsign,
    basis,
    interfacequads,
    mesh,
)
    ndofs = size(nodalsolutions)[1]
    err = zeros(ndofs)
    interpolater = InterpolatingPolynomial(ndofs, basis)

    ncells = CutCellDG.number_of_cells(mesh)

    for cellid = 1:ncells
        if CutCellDG.cell_sign(mesh, cellid) == 0
            quad = interfacequads[levelsetsign, cellid]
            cellmap = CutCellDG.cell_map(mesh, levelsetsign, cellid)
            scalearea = CutCellDG.interface_scale_areas(interfacequads, cellid)

            nodeids = CutCellDG.nodal_connectivity(mesh, levelsetsign, cellid)
            elementsolution = nodalsolutions[:, nodeids]
            update!(interpolater, elementsolution)

            add_interface_error_squared!(
                err,
                interpolater,
                exactsolution,
                cellmap,
                quad,
                scalearea,
            )
        end
    end
    return sqrt.(err)
end

function add_interface_norm_squared!(vals, func, cellmap, quad, facescale)
    @assert length(quad) == length(facescale)
    for (i, (p, w)) in enumerate(quad)
        v = func(cellmap(p))
        vals .+= v .^ 2 * facescale[i] * w
    end
end

function integral_norm_on_interface(
    func,
    interfacequads,
    levelsetsign,
    mesh,
    ndofs,
)
    vals = zeros(ndofs)
    ncells = CutCellDG.number_of_cells(mesh)

    for cellid = 1:ncells
        if CutCellDG.cell_sign(mesh, cellid) == 0
            cellmap = CutCellDG.cell_map(mesh, levelsetsign, cellid)
            quad = interfacequads[levelsetsign, cellid]
            facescale = CutCellDG.interface_scale_areas(interfacequads, cellid)

            add_interface_norm_squared!(vals, func, cellmap, quad, facescale)
        end
    end
    return sqrt.(vals)
end

function update_maxnorm_error!(
    globalerror,
    interpolater,
    exactsolution,
    cellmap,
    quad,
)
    ndofs = length(globalerror)
    for (p, w) in quad
        numsol = interpolater(p)
        exsol = exactsolution(cellmap(p))
        pointerror = abs.(numsol - exsol)
        for i = 1:ndofs
            globalerror[i] = max(globalerror[i], pointerror[i])
        end
    end
end

function interface_maxnorm_error(
    nodalsolutions,
    exactsolution,
    levelsetsign,
    basis,
    interfacequads,
    mesh,
)
    ndofs = size(nodalsolutions)[1]
    err = zeros(ndofs)
    interpolater = InterpolatingPolynomial(ndofs, basis)

    ncells = CutCellDG.number_of_cells(mesh)

    for cellid = 1:ncells
        if CutCellDG.cell_sign(mesh, cellid) == 0
            quad = interfacequads[levelsetsign, cellid]
            cellmap = CutCellDG.cell_map(mesh, levelsetsign, cellid)

            nodeids = CutCellDG.nodal_connectivity(mesh, levelsetsign, cellid)
            elementsolution = nodalsolutions[:, nodeids]
            update!(interpolater, elementsolution)

            update_maxnorm_error!(
                err,
                interpolater,
                exactsolution,
                cellmap,
                quad,
            )
        end
    end
    return err
end

function update_maxnorm!(vals, func, cellmap, quad)
    ndofs = length(vals)
    for (p, w) in quad
        v = abs.(func(cellmap(p)))
        for i = 1:ndofs
            vals[i] = max(vals[i], v[i])
        end
    end
end

function maxnorm_on_interface(func, interfacequads, levelsetsign, mesh, ndofs)
    vals = zeros(ndofs)
    ncells = CutCellDG.number_of_cells(mesh)

    for cellid = 1:ncells
        if CutCellDG.cell_sign(mesh, cellid) == 0
            cellmap = CutCellDG.cell_map(mesh, levelsetsign, cellid)
            quad = interfacequads[levelsetsign, cellid]

            update_maxnorm!(vals, func, cellmap, quad)
        end
    end
    return vals
end

function uniform_mesh_L2_error(nodalsolutions, exactsolution, basis, quad, mesh)
    ndofs, nnodes = size(nodalsolutions)
    err = zeros(ndofs)
    interpolater = InterpolatingPolynomial(ndofs, basis)
    ncells = CutCellDG.number_of_cells(mesh)

    for cellid = 1:ncells
        cellmap = CutCellDG.cell_map(mesh, cellid)
        nodeids = CutCellDG.nodal_connectivity(mesh, cellid)
        elementsolution = nodalsolutions[:, nodeids]
        update!(interpolater, elementsolution)
        add_cell_error_squared!(err, interpolater, exactsolution, cellmap, quad)
    end
    return sqrt.(err)
end

function integral_norm_on_uniform_mesh(func, quad, mesh, ndofs)
    numcells = CutCellDG.number_of_cells(mesh)
    vals = zeros(ndofs)
    for cellid in 1:numcells
        cellmap = CutCellDG.cell_map(mesh,cellid)
        add_cell_norm_squared!(vals, func, cellmap, quad)
    end
    return sqrt.(vals)
end

function lame_lambda(k, m)
    return k - 2m / 3
end


text:
# Tests insertion and deletion of gates.
using Test

function test_indel( n::Int64, p::Parameters )
  @testset "testing insertion and deletion of gates" begin
    for i = 1:n
      c = random_chromosome(p)
      sav_c = deepcopy(c)
      output_values(sav_c)  # needed to set active gates to active
      insert_gate!(c)
      @test output_values(c) == output_values(sav_c)
      delete_gate!(c)
      @test output_values(c) == output_values(sav_c)
    end
  end
end

function test_delete_gate( n::Int64, p::Parameters )
  for i = 1:n
    c = random_chromosome(p)
    sav_c = deepcopy(c)
    output_values(sav_c)  # needed to set active gates to active
    #dg = rand(1:(c.params.numinteriors-p.numoutputs))
    #println("dg: ",dg,"  active: ",sav_c.interiors[dg].active)
    #delete_gate!(c,dg)
    #if sav_c.interiors[dg].active == false
    delete_gate!(c)
    if number_active_gates(sav_c) < sav_c.params.numinteriors
      println("test for i = ",i)
      if output_values(c) != output_values(sav_c)
        println((output_values(c),output_values(sav_c)))
        print_circuit(c)
        print_circuit(sav_c)
        return (sav_c,c)
      end
    end
  end
end

function test_insert_gate( n::Int64, p::Parameters )
  for i = 1:n
    c = random_chromosome(p)
    sav_c = deepcopy(c)
    insert_gate!(c)
    if output_values(c) != output_values(sav_c)
      println((output_values(c),output_values(sav_c)))
      print_circuit(c)
      print_circuit(sav_c)
      return (sav_c,c)
    end
  end
end

ext: jl


text:
# # ChainRules integration demo: Relaxed Unit Commitment

#md # [![](https://img.shields.io/badge/show-github-579ACA.svg)](@__REPO_ROOT_URL__/docs/src/examples/chainrules_unit.jl)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/chainrules_unit.ipynb)

# In this example, we will demonstrate the integration of DiffOpt with
# [ChainRulesCore.jl](https://juliadiff.org/ChainRulesCore.jl/stable/),
# the library allowing the definition of derivatives for functions
# that can then be used by automatic differentiation systems.

using JuMP
import DiffOpt
import Plots
import LinearAlgebra: ⋅
import GLPK
import ChainRulesCore

# ## Unit commitment problem

# We will consider a unit commitment problem, finding the cost-minimizing activation
# of generation units in a power network over multiple time periods.
# The considered constraints include:
# - Demand satisfaction of several loads
# - Ramping constraints
# - Generation limits.

# The decisions are:
# - ``u_{it} \in \{0,1\}``: activation of the ``i``-th unit at time ``t``
# - ``p_{it}``: power output of the ``i``-th unit at time ``t``.

# DiffOpt handles convex optimization problems only, we therefore
# relax the domain of the ``u_{it}`` variables to ``\left[0,1\right]``.

# ## Primal UC problem

# ChainRules defines the differentiation of functions.
# The actual function that is differentiated in the context of DiffOpt is the
# solution map taking in input the problem parameters and returning the solution.

function unit_commitment(
        load1_demand, load2_demand, gen_costs, noload_costs;
        model = Model(GLPK.Optimizer), silent=false)
    MOI.set(model, MOI.Silent(), silent)

    ## Problem data
    units = [1, 2] # Generator identifiers
    load_names = ["Load1", "Load2"] # Load identifiers
    n_periods = 4 # Number of time periods
    Pmin = Dict(1 => fill(0.5, n_periods), 2 => fill(0.5, n_periods)) # Minimum power output (pu)
    Pmax = Dict(1 => fill(3.0, n_periods), 2 => fill(3.0, n_periods)) # Maximum power output (pu)
    RR = Dict(1 => 0.25, 2 => 0.25) # Ramp rates (pu/min)
    P0 = Dict(1 => 0.0, 2 => 0.0) # Initial power output (pu)
    D = Dict("Load1" => load1_demand, "Load2" => load2_demand) # Demand (pu)
    Cp = Dict(1 => gen_costs[1], 2 => gen_costs[2]) # Generation cost coefficient ($/pu)
    Cnl = Dict(1 => noload_costs[1], 2 => noload_costs[2]) # No-load cost ($)

    ## Variables
    ## Note: u represents the activation of generation units.
    ## Would be binary in the typical UC problem, relaxed here to u ∈ [0,1]
    ## for a linear relaxation.
    @variable(model, 0 <= u[g in units, t in 1:n_periods] <= 1) # Commitment
    @variable(model, p[g in units, t in 1:n_periods] >= 0) # Power output (pu)

    ## Constraints

    ## Energy balance
    @constraint(
        model,
        energy_balance_cons[t in 1:n_periods],
        sum(p[g, t] for g in units) == sum(D[l][t] for l in load_names),
    )

    ## Generation limits
    @constraint(model, [g in units, t in 1:n_periods], Pmin[g][t] * u[g, t] <= p[g, t])
    @constraint(model, [g in units, t in 1:n_periods], p[g, t] <= Pmax[g][t] * u[g, t])

    ## Ramp rates
    @constraint(model, [g in units, t in 2:n_periods], p[g, t] - p[g, t - 1] <= 60 * RR[g])
    @constraint(model, [g in units], p[g, 1] - P0[g] <= 60 * RR[g])
    @constraint(model, [g in units, t in 2:n_periods], p[g, t - 1] - p[g, t] <= 60 * RR[g])
    @constraint(model, [g in units], P0[g] - p[g, 1] <= 60 * RR[g])

    ## Objective
    @objective(
        model,
        Min,
        sum((Cp[g] * p[g, t]) + (Cnl[g] * u[g, t]) for g in units, t in 1:n_periods),
    )

    optimize!(model)
    ## asserting finite optimal value
    @assert termination_status(model) == MOI.OPTIMAL
    ## converting to dense matrix
    return JuMP.value.(p.data)
end

m = Model(GLPK.Optimizer)
@show unit_commitment(
    [1.0, 1.2, 1.4, 1.6], [1.0, 1.2, 1.4, 1.6],
    [1000.0, 1500.0], [500.0, 1000.0],
    model=m, silent=true
)

# ## Perturbation of a single input parameter

# Let us vary the demand at the second time frame on both loads:

demand_values = 0.05:0.05:3.0
pvalues = map(demand_values) do di
    unit_commitment(
        [1.0, di, 1.4, 1.6], [1.0, di, 1.4, 1.6],
        [1000.0, 1500.0], [500.0, 1000.0],
        silent=true,
    )
end
pflat = [getindex.(pvalues, i) for i in eachindex(pvalues[1])];

# The influence of this variation of the demand is piecewise linear on the
# generation at different time frames:

Plots.scatter(demand_values, pflat, xaxis = ("Demand"), yaxis = ("Generation"))
Plots.title!("Different time frames and generators")
Plots.xlims!(0.0, 3.5)

# ## Forward Differentiation

# Forward differentiation rule for the solution map of the unit commitment problem.
# It takes as arguments:
# 1. the perturbations on the input parameters
# 2. the differentiated function
# 3. the primal values of the input parameters,

# and returns a tuple `(primal_output, perturbations)`, the main primal result
# and the perturbation propagated to this result:

function ChainRulesCore.frule(
        (_, Δload1_demand, Δload2_demand, Δgen_costs, Δnoload_costs),
        ::typeof(unit_commitment),
        load1_demand, load2_demand, gen_costs, noload_costs;
        optimizer=GLPK.Optimizer,
        )
    ## creating the UC model with a DiffOpt optimizer wrapper around GLPK
    model = Model(() -> DiffOpt.diff_optimizer(optimizer))
    ## building and solving the main model
    pv = unit_commitment(
        load1_demand, load2_demand, gen_costs, noload_costs, model=model)
    energy_balance_cons = model[:energy_balance_cons]

    ## Setting some perturbation of the energy balance constraints
    ## Perturbations are set as MOI functions
    Δenergy_balance = [
        convert(MOI.ScalarAffineFunction{Float64}, d1 + d2)
        for (d1, d2) in zip(Δload1_demand, Δload2_demand)
    ]
    MOI.set.(
        model,
        DiffOpt.ForwardInConstraint(), energy_balance_cons,
        Δenergy_balance,
    )

    p = model[:p]
    u = model[:u]

    ## setting the perturbation of the linear objective
    Δobj = sum(Δgen_costs ⋅ p[:,t] + Δnoload_costs ⋅ u[:,t] for t in size(p, 2))
    MOI.set(model, DiffOpt.ForwardInObjective(), Δobj)
    DiffOpt.forward(JuMP.backend(model))
    ## querying the corresponding perturbation of the decision
    Δp = MOI.get.(model, DiffOpt.ForwardOutVariablePrimal(), p)
    return (pv, Δp.data)
end

# We can now compute the perturbation of the output powers `Δpv`
# for a perturbation of the first load demand at time 2:

load1_demand = [1.0, 1.0, 1.4, 1.6]
load2_demand = [1.0, 1.0, 1.4, 1.6]
gen_costs = [1000.0, 1500.0]
noload_costs = [500.0, 1000.0];

# all input perturbations are 0
# except first load at time 2
Δload1_demand = 0 * load1_demand
Δload1_demand[2] = 1.0
Δload2_demand = 0 * load2_demand
Δgen_costs = 0 * gen_costs
Δnoload_costs = 0 * noload_costs
(pv, Δpv) = ChainRulesCore.frule(
    (nothing, Δload1_demand, Δload2_demand, Δgen_costs, Δnoload_costs),
    unit_commitment,
    load1_demand, load2_demand, gen_costs, noload_costs,
)

Δpv

# The result matches what we observe in the previous figure:
# the generation of the first generator at the second time frame (third element on the plot).

# # Reverse-mode differentiation of the solution map

# The `rrule` returns the primal and a pullback.
# The pullback takes a seed for the optimal solution `̄p` and returns
# derivatives with respect to each input parameter of the function.

function ChainRulesCore.rrule(
        ::typeof(unit_commitment),
        load1_demand, load2_demand, gen_costs, noload_costs;
        optimizer=GLPK.Optimizer,
        silent=false)
    model = Model(() -> DiffOpt.diff_optimizer(optimizer))
    ## solve the forward UC problem
    pv = unit_commitment(
        load1_demand, load2_demand, gen_costs, noload_costs,
        model=model, silent=silent)
    function pullback_unit_commitment(pb)
        p = model[:p]
        u = model[:u]
        energy_balance_cons = model[:energy_balance_cons]

        MOI.set.(model, DiffOpt.BackwardInVariablePrimal(), p, pb)
        DiffOpt.backward(JuMP.backend(model))

        obj = MOI.get(model, DiffOpt.BackwardOutObjective())

        ## computing derivative wrt linear objective costs
        dgen_costs = similar(gen_costs)
        dgen_costs[1] = sum(JuMP.coefficient.(obj, p[1,:]))
        dgen_costs[2] = sum(JuMP.coefficient.(obj, p[2,:]))

        dnoload_costs = similar(noload_costs)
        dnoload_costs[1] = sum(JuMP.coefficient.(obj, u[1,:]))
        dnoload_costs[2] = sum(JuMP.coefficient.(obj, u[2,:]))

        ## computing derivative wrt constraint constant
        dload1_demand = JuMP.constant.(
            MOI.get.(model, DiffOpt.BackwardOutConstraint(), energy_balance_cons))
        dload2_demand = copy(dload1_demand)
        return (dload1_demand, dload2_demand, dgen_costs, dnoload_costs)
    end
    return (pv, pullback_unit_commitment)
end

# We can set a seed of one on the power of the first generator at the second time frame and zero for all other
# parts of the solution:

(pv, pullback_unit_commitment) = ChainRulesCore.rrule(
    unit_commitment,
    load1_demand, load2_demand, gen_costs, noload_costs,
    optimizer=GLPK.Optimizer,
    silent=true,
)
dpv = 0 * pv
dpv[1,2] = 1
dargs = pullback_unit_commitment(dpv)
(dload1_demand, dload2_demand, dgen_costs, dnoload_costs) = dargs;

# The sensitivities with respect to the load demands are:
dload1_demand

# and:
dload2_demand

# The sensitivity of the generation is propagated to the sensitivity of both
# loads at the second time frame.

# This example integrating ChainRules was designed with support
# from [Invenia Technical Computing](https://www.invenia.ca/).


text:
using Juqst

potentialSingles = [
                    [[0,0],[0,1]], # IX
                    [[0,0],[1,0]], # IY
                    [[0,0],[1,1]], # IZ
                    ]

all2QlMuBs =  [  [[0,0,0,0],[1,1,0,1],[1,0,1,1],[0,1,1,0]], #II ZX YZ XY
                 [[0,0,0,0],[1,1,1,0],[0,1,1,1],[1,0,0,1]], #II ZY XZ YX
                 [[0,0,0,0],[0,0,0,1],[0,1,0,0],[0,1,0,1]], #II IX XI XX
                 [[0,0,0,0],[0,0,1,0],[1,0,0,0],[1,0,1,0]], #II IY YI YY
                 [[0,0,0,0],[0,0,1,1],[1,1,0,0],[1,1,1,1]]] #II IZ ZI ZZ

# We only want to choose the first two types for the initial runs
potentialMuBs = [all2QlMuBs[1],all2QlMuBs[2]]

pX = [0 1;1 0]
pY = [0 -im;im 0]
pZ = [1 0;0 -1]
pI = [1 0;0 1]
paulis = [pI,pX,pY,pZ]
superPaulis = [makeSuper(x) for x in paulis];

"""
generateSensibleSubsamples(experiments)
## Arguments
	experiments: A list of tuples, showing how to locally divide the qubits.
Each experiment needs to be for the same number of qubits.
## Example
If you pass in [(2,2,2),(1,2,2,1)]
Then this would generage the needed data for two expeiments, each for 6 qubits.
- The first would use two qubit gates on qubits 0&1, 2&3 and 4&5.
- The second would use single qubit gates on qubits 0 and 5, and two qubit gates on 1&2 and 3&4.

With three qubits, we might call it thus:
```
(experiments,paulisAll,ds) = generateSensibleSubsamples([(2,1),(1,2)])
print("\$experiments\\n")
print("\$paulisAll\\n")
```

```
Any[Any[(2, 1), (1, 1)], Any[(1, 3), (2, 2)]]
Any[Any[[[0, 0, 0, 0], [1, 1, 0, 1], [1, 0, 1, 1], [0, 1, 1, 0]], [[0, 0], [0, 1]]], Any[[[0, 0], [1, 1]], [[0, 0, 0, 0], [1, 1, 1, 0], [0, 1, 1, 1], [1, 0, 0, 1]]]]
```

The Experiment numbers are defined like this:

```
potentialSingles = [
                    [[0,0],[0,1]], # IX
                    [[0,0],[1,0]], # IY
                    [[0,0],[1,1]], # IZ
                    ]

all2QlMuBs =  [  [[0,0,0,0],[1,1,0,1],[1,0,1,1],[0,1,1,0]], #II ZX YZ XY
                 [[0,0,0,0],[1,1,1,0],[0,1,1,1],[1,0,0,1]], #II ZY XZ YX
                 [[0,0,0,0],[0,0,0,1],[0,1,0,0],[0,1,0,1]], #II IX XI XX
                 [[0,0,0,0],[0,0,1,0],[1,0,0,0],[1,0,1,0]], #II IY YI YY
                 [[0,0,0,0],[0,0,1,1],[1,1,0,0],[1,1,1,1]]] #II IZ ZI ZZ
```

## Returns
The following tuple (chosenExperiments, chosenPaulis, ds), where
- chosenExperiments is a list of experiment list, where each experiment list contains a tuple of qubit size and experiment number.
- chosenPaulis are the Paulis that will be generated by each tuple in an experiment list (using the 01->X 10->Y 11->Z mapping)
- ds is just a convenient bit-array of offsets of 2*noOfQubits, used in working out offsets for the decoder.
"""
function generateSensibleSubsamples(experiments)
	qubits = sum(experiments[1])
    for i in experiments
        @assert(qubits == sum(i),"All experiments have to use the same number of qubits")
    end
    chosenExperiments = []
    chosenPaulis = []
    for e in experiments
        thisExperiment =[]
        thisPaulis = []
        for q in e
            if q ==2
                chosen = rand(1:2)
                push!(thisExperiment,(2,chosen))
                push!(thisPaulis,all2QlMuBs[chosen])
            elseif q == 1
                chosen = rand(1:3)
                push!(thisExperiment,(1,chosen))
                push!(thisPaulis,potentialSingles[chosen])
            else
                @assert(false,"We only deal with one and two qubit splis just now\n")
            end
        end
        push!(chosenExperiments,thisExperiment)
        push!(chosenPaulis,thisPaulis)
    end
    ds = vcat(
        [[0 for _ = 1:2*qubits]],
        [map(x->parse(Int,x),collect(reverse(string(i,base=2,pad=2*qubits)))) for i in [2^b for b=0:2*qubits-1]]);
    return (chosenExperiments,chosenPaulis,ds)

end

"""
addCircuit(circuit,experiment_type,qubits_id)

## Arguments
   circuit - the qiskit circuit we will add gates to.
   experiment-type: (Integer,Integer) - the stabilizer group we wont to generate (see below)
   qubits_id: Array of Int - the 'number' of the qubits the circuit is to be on.

The stabilizer circuits vary depending on the experiment-type. This is a tuple, the first number indicating the number of
qubits (currently 1 or 2), the second which experiment currently 1:3 for single qubits and 1:5 for two qubits.

The qubits_id will be used to work out which qubits in the circuit to use. They are "JULIA INDEXED".

The qubits are identified via the circuit ```qubits = circuit.qubits```. But because of the way PyCall works, they will be
in a 1 indexed array. e.g. qiskitits qubit_0 will be qubits[1].

Therefore if you want the circuit to be on qubits 0 and 1, you pass in [1,2] into qubits_id.

The circuits that will be added will generate the following for two qubits its:

1.    [[0,0,0,0],[1,1,0,1],[1,0,1,1],[0,1,1,0]], #II ZX YZ XY
2.    [[0,0,0,0],[1,1,1,0],[0,1,1,1],[1,0,0,1]], #II ZY XZ YX
3.    [[0,0,0,0],[0,0,0,1],[0,1,0,0],[0,1,0,1]], #II IX XI XX
4.    [[0,0,0,0],[0,0,1,0],[1,0,0,0],[1,0,1,0]], #II IY YI YY
5.    [[0,0,0,0],[0,0,1,1],[1,1,0,0],[1,1,1,1]]] #II IZ ZI ZZ

For one qubit its

1.    [[0,0],[0,1]], # IX
2.    [[0,0],[1,0]], # IY
3.    [[0,0],[1,1]], # IZ

## Returns nothing

The qiskit circuit will have had the gates added.
"""
function addCircuit(circuit,experiment_type,qubits_id)
    qubits = circuit.qubits
    (q,no) = experiment_type
    if q == 2 # Two qubit MUBS
        q1 = qubits[qubits_id[1]]
        q2 = qubits[qubits_id[2]]
        if no == 1
            circuit.h(q1)
            circuit.h(q2)
            circuit.s(q2)
            circuit.cx(q2,q1)
            circuit.s(q1)
            circuit.cx(q2,q1)
            return
        elseif no == 2
            circuit.h(q1)
            circuit.h(q2)
            circuit.s(q1)
            circuit.cx(q2,q1)
            circuit.s(q1)
            circuit.cx(q2,q1)
            return
        elseif no ==3
           circuit.h(q1)
           circuit.h(q2)
           return
        elseif no == 4
            circuit.h(q1)
            circuit.h(q2)
            circuit.s(q1)
            circuit.s(q2)
            return
        else
            return # no circuit needed for Z
        end
    elseif q == 1 # One Qubits MUBS
        q1 = qubits[qubits_id[1]]
        if no == 1
            circuit.h(q1)
            return
        elseif no ==2
            circuit.h(q1)
            circuit.s(q1)
            return
        else
            return # Nothing needed for z
        end
    else
        @warn("passed in an invalid experiment, qubit number $q")
    end

end

"""
addReverseCircuit(circuit,experiment_type,qubits_id)

## Arguments
   circuit - the qiskit circuit we will add gates to.
   experiment-type: (Integer,Integer) - the stabilizer group we wont to generate (see below)
   qubits_id: Array of Int - the 'number' of the qubits the circuit is to be on.

Effectively this 'reverses' the circuit added by addCircuit.

e.g.

```
addCircuit(my_circuit,(2,1),[1,2])
addReverseCircuit(my_circuit,(2,1),[1,2,])
```

Will create a 'nothing' circuit i.e. self inverting.

The stabilizer circuits vary depending on the experiment-type. This is a tuple, the first number indicating the number of
qubits (currently 1 or 2), the second which experiment currently 1:3 for single qubits and 1:5 for two qubits.

The qubits_id will be used to work out which qubits in the circuit to use. They are "JULIA INDEXED".

The qubits are identified via the circuit ```qubits = circuit.qubits```. But because of the way PyCall works, they will be
in a 1 indexed array. e.g. qiskitits qubit_0 will be qubits[1].

Therefore if you want the circuit to be on qubits 0 and 1, you pass in [1,2] into qubits_id.

The circuits that will be added will generate the following for two qubits its:

1.    [[0,0,0,0],[1,1,0,1],[1,0,1,1],[0,1,1,0]], #II ZX YZ XY
2.    [[0,0,0,0],[1,1,1,0],[0,1,1,1],[1,0,0,1]], #II ZY XZ YX
3.    [[0,0,0,0],[0,0,0,1],[0,1,0,0],[0,1,0,1]], #II IX XI XX
4.    [[0,0,0,0],[0,0,1,0],[1,0,0,0],[1,0,1,0]], #II IY YI YY
5.    [[0,0,0,0],[0,0,1,1],[1,1,0,0],[1,1,1,1]]] #II IZ ZI ZZ

For one qubit its

1.    [[0,0],[0,1]], # IX
2.    [[0,0],[1,0]], # IY
3.    [[0,0],[1,1]], # IZ

## Returns nothing

The qiskit circuit will have had the gates added.
"""
function addReverseCircuit(circuit,experiment_type,qubits_id)
# Should probably actually reverse the circuit - this is
# possibly simpler to understand but will be harder to maintain.
    qubits = circuit.qubits
    (q,no) = experiment_type
    if q == 2 # Two qubit MUBS
#              [[0,0,0,0],[1,1,0,1],[1,0,1,1],[0,1,1,0]], #II ZX YZ XY
#              [[0,0,0,0],[1,1,1,0],[0,1,1,1],[1,0,0,1]], #II ZY XZ YX
#              [[0,0,0,0],[0,0,0,1],[0,1,0,0],[0,1,0,1]], #II IX XI XX
#              [[0,0,0,0],[0,0,1,0],[1,0,0,0],[1,0,1,0]], #II IY YI YY
#              [[0,0,0,0],[0,0,1,1],[1,1,0,0],[1,1,1,1]]] #II IZ ZI ZZ

        q1 = qubits[qubits_id[1]]
        q2 = qubits[qubits_id[2]]
        if no == 1
            circuit.cx(q2,q1)
            circuit.s(q1)
            circuit.s(q1)
            circuit.s(q1)
            circuit.cx(q2,q1)
            circuit.s(q2)
            circuit.s(q2)
            circuit.s(q2)
            circuit.h(q1)
            circuit.h(q2)
            return
        elseif no == 2
            circuit.cx(q2,q1)
            circuit.s(q1)
            circuit.s(q1)
            circuit.s(q1)
            circuit.cx(q2,q1)
            circuit.s(q1)
            circuit.s(q1)
            circuit.s(q1)
            circuit.h(q1)
            circuit.h(q2)
            return
        elseif no ==3
           circuit.h(q1)
           circuit.h(q2)
           return
        elseif no == 4
            circuit.s(q1)
            circuit.s(q2)
            circuit.s(q1)
            circuit.s(q2)
            circuit.s(q1)
            circuit.s(q2)
            circuit.h(q1)
            circuit.h(q2)
            return
        else
            return # no circuit needed for Z
        end
    elseif q == 1 # One Qubits MUBS
        q1 = qubits[qubits_id[1]]
#                     [[0,0],[0,1]], # IX
#                     [[0,0],[1,0]], # IY
#                     [[0,0],[1,1]], # IZ
        if no == 1
            circuit.h(q1)
            return
        elseif no ==2
            circuit.s(q1)
            circuit.s(q1)
            circuit.s(q1)
            circuit.h(q1)
            return
        else
            return # Nothing needed for z
        end
    else
        @warn("passed in an invalid experiment, qubit number $q")
    end

end

"""
setUpExperiment(experiment,circuit;reverse=false)

## Arguments
    experiment: List of tuples - the stabilizers to be extracted
    circuit: circuit to be done
    reverse: apply or removing? Boolean.

Cycles through the experiment applying the correct gates in the circuit.
The experiment is set up as a list of tuples (q_no,exp_no), representing
the number of qubits (currently 1 or 2) and the type of experiment.

For example [(1,3),(2,1),(2,3)] would be a 5 qubit experiment.
-  The 1 qubit experiment 3 would be applied to the first qubit.
-  The 2 qubit experiment 1 would be applied to the second and third qubits
-  The 2 qubit experiment 3 would be applied to the fourth and fifth qubits.

## Returns

Nothing, but amends the supplied circuit.
"""
function setUpExperiment(experiment,circuit;reverse=false)
    numberOfQubits = foldl((a,b)->a[1]+b[1],experiment)
    start = 1
    for (q,no) in experiment
        if q ==2
            ids = [start,start+1]
            start +=2
            if reverse
                addReverseCircuit(circuit,(q,no),ids)
            else
                addCircuit(circuit,(q,no),ids)
            end
        elseif q==1
            ids = [start]
            start +=1
            if reverse
                addReverseCircuit(circuit,(q,no),ids)
            else
                addCircuit(circuit,(q,no),ids)
            end
        else
            @warn("Currently we only deal with local mubs, not on $q qubits")
        end
    end
end

"""
genPTwirl(circuit,len,qubit_id)

generates a Pauli twirl on the circuit consisting of len gates + inversion gate. On qubit_id

## Arguments
	circuit - the qiskit circuit to alter
	len - an integer representing the number of gates in the twirl
	qubit_id - an integer representing the qubit. Index from 1 (Julia), the qubit register is retrieved from circuit and is not needed.
## Returns
	Nothing - the circuit is altered.
"""
function genPTwirl(circuit,len,qubit_id)
    qubits = circuit.qubits
    list = rand(1:4,len)
    total = foldr(*,map(x->superPaulis[x],list))
    push!(list,findfirst(x->x==total,superPaulis))
    q1 = qubits[qubit_id]
    for i in list
        if i == 2
            circuit.x(q1)
        elseif i == 3
            circuit.y(q1)
        elseif i == 4
            circuit.z(q1)
        else
            circuit.id(q1)
        end
        circuit.barrier(q1)
    end

end

"""
getCircuit(experiment,noOfGates)

## Arguments
	experiment: a list of tuples (such as that obtained from generateSensibleSubsamples)
	noOfGates: the number of Paulis we want in the middle (to twirl)

Generates a qiskit circuit that implements the experiment passed in (e.g [(2,1) (2,2)] will generate
a cricuit that uses two 2 qubit MUBS at the beginning and end (1 and 2 respectively) and then performs
a Pauli Twirl in the middle). The circuits are self-inverting.

Appropriate barriers are inserted.

## Example
```
(experiments,paulisAll,ds) = generateSensibleSubsamples([(2,1),(1,2)]) # Get some experiments
circuit = getCircuit(experiments[1],12) # twelve pauli twirl circuit.
qiskit.execute(circuit,qiskit.Aer.get_backend("qasm_simulator"),shots=100).result().get_counts() # noisless simulator.
circuit.draw() # Will draw the circuit
```

## Returns
 	the qiskit circuit.
 """
function getCircuit(experiment,noOfGates)
	numberOfQubits = -1
	if length(experiment) == 1
		numberOfQubits = experiment[1][1]
	else
    	numberOfQubits = foldl((a,b)->a[1]+b[1],experiment)
    end
    circuit =qiskit.QuantumCircuit(numberOfQubits,numberOfQubits)
    setUpExperiment(experiment,circuit)
    circuit.barrier()
    for i in 1:numberOfQubits
        genPTwirl(circuit,noOfGates,i)
    end
    setUpExperiment(experiment,circuit,reverse=true)
    circuit.barrier()
    circuit.measure(0:numberOfQubits-1,0:numberOfQubits-1)
    return circuit
end

"""
doASeries(experiment,lengths,sequences,shots,noise_mode)
## Arguments
    experiment - the experiment to perform (see getCircuit for detailed explanation)
    lengths - a list of the lengths to perform the experiment on.
	sequences - the number of randomised sequences to perform at each length.
	shots - the number of measurements of each sequences
	noise-model - the noise model you want passed to the 'Aer' simulator.

A convenience function to simulate the basic information gathering algorithm.
We assume a an Aer simulator, and you supply the noise model.
For edifferent simulators or to run on a real device you will have to re-implement this function, probably
just change the execute line.
## Returns
	a list of the a list of the results of each sequence.(ordered by the same way as lengths).
Note on a real device you probably want to randomise the order that different lengths are called i.e. don't do all the short ones, slightly longer and then longer ones in a row.
"""
function doASeries(experiment,lengths,sequences,shots,noise_model)
    results = []
    for l in lengths
        perSeqResults = []
        for seq = 1:sequences
            circuit = getCircuit(experiment,l)
            counts = qiskit.execute(circuit,qiskit.Aer.get_backend("qasm_simulator"),noise_model=noise_model,shots=shots).result().get_counts()
            push!(perSeqResults,counts)
        end
        push!(results,perSeqResults)
    end
    return results
end

"""
getCounts(results,lengths,noOfQubtis)

Convenience funciton that takes the results from doASeries and consolidates the counts accross different sequences
i.e. if we have 10 sequences at lengths 10,20 and 40, groups all the sequences at length 10 together, the seqeunces at length 20 ... etc.
It then transforms into a probability vector (i.e. the vector for each length by the total number of counts).

NOTE: we assume that the machine is able to handle a vector of size 2^noOfQubits

##  Returns
	A list of probability vectors, one for each length.
"""
function getCounts(results,lengths,noOfQubits)
    consolidatedCounts = [zeros(2^noOfQubits) for _ in 1:length(lengths)];
    for (ix,seqs) in enumerate(results)
        for s in seqs
            for k in keys(s)
                consolidatedCounts[ix][parse(Int,k,base=2)+1] += s[k]
            end
        end
    end
    #Normalise to a percentage.
    for i  in 1:length(consolidatedCounts)
        consolidatedCounts[i]  = consolidatedCounts[i]./sum(consolidatedCounts[i])
    end
    return consolidatedCounts
end

"""
extractEigenvalues(counts,lengths)

Uses the functionality of Juqst - quantumNoise algorithms to extract the eigenvalues from a vector of counts.
For documentation on fitTheFidelities, see rharper2/Juqst.jl on gitHub.
To see how the counts should look see getCounts or the example workbooks.

## Returns
	the Eigenvalues for the supplied probability vector.
"""
function extractEigenvalues(counts,lengths)
    (params,l, failed) = fitTheFidelities(lengths,counts)
    return vcat(1,[p[2] for p in params]) # We don't fit the first one, it is always 1 for CPTP maps
end

function binaryArrayToNumber(x)
    return foldl((y1,y2)->y1<<1+y2,x)
end

function generateEigenvalueRecur(stabs,current)
    if length(stabs) != 1
        current = generateEigenvalueRecur(stabs[2:end],current)
        q = stabs[1]
        multiplier = 2^length(q)
        newCurrent = []
        for c in current
            for i in [binaryArrayToNumber(x) for x in stabs[1]]
                push!(newCurrent,multiplier*c+i)
            end
        end
        return newCurrent
    else
        return [binaryArrayToNumber(x) for x in stabs[1]]
    end
end

"""
generateEigenvalues(stabs)

## Arguments
	stabs - a list of the Paulis used in the experiment, this can be obtained from the generateSensibleSubsamples function.

Uses a recursive algorithm to work out which global eigenvalues will be learnt through the supplied Paulis/experiment. You can use
PEEL.fidelityLabels(ix-1) to see the Pauli representation of these eigenvalues. Note the need to subtract 1, to mvoe from Julia 1 index
to a 0 index label (i.e. in Julia II...I is 1, whereas fidelityLabels expects this to be 0).
## Returns
	a list of eigenvalues (Julia 1 indexed).
"""
function generateEigenvalues(stabs)
    x = generateEigenvalueRecur(stabs,[])
    return x.+1 # Julia indexing.
end

"""
getStabiliserForExperiment(experiment)

For a given experiment (see generateSensibleSubsamples)  returns the Paulis (in bit form) that are accessed via the experiment.
"""
function getStabilizerForExperiment(experiment)
    paulis =[]
    for (q,s) in experiment
        if q == 2
            push!(paulis,all2QlMuBs[s])
        else
            push!(paulis,potentialSingles[s])
        end
    end
    return paulis
end


text:
#!/usr/bin/env julia
doc""" Script to test KullbackLeibler.jl

- FIXME write it!
- XXX Or do it automatically with jldoctest (https://juliadocs.github.io/Documenter.jl/latest/man/doctests/#Doctests-1).

- MIT Licensed, (C) 2018 Lilian Besson (Naereen)
  https://GitHub.com/Naereen/KullbackLeibler.jl
"""

__author__ = "Lilian Besson (Naereen)"
__version__ = "0.1"

include("../julia-src/KullbackLeibler.jl")
# include("KullbackLeibler.jl")
using KullbackLeibler
using Base.Test

# utility for the tests

tolerance = 1e-6

function isclose(x, y, tol=tolerance)
    return abs(x - y) <= tol
end

# write your own tests here

@test isclose(KullbackLeibler.klBern(0.5, 0.5), 0.0)
@test isclose(KullbackLeibler.klBern(0.1, 0.9), 1.757779)
# And this KL is symmetric
@test isclose(KullbackLeibler.klBern(0.9, 0.1), 1.757779)
@test isclose(KullbackLeibler.klBern(0.4, 0.5), 0.020135)
@test isclose(KullbackLeibler.klBern(0.01, 0.99), 4.503217)


text:
# ----------------------------------------------------------
# ------------------------  Show  --------------------------
# ----------------------------------------------------------
# I think it is pretty much self-explenatory

Base.show(io::IO, FO::FunctionOperator) =
    print(io, "FunctionOperator{$(eltype(FO))}($(FO.name), $(FO.inDims), $(FO.outDims))")

Base.show(io::IO, ::MIME"text/plain", FO::FunctionOperator) =
    print(io, """
FunctionOperator with eltype $(eltype(FO))
    Name: $(FO.name)
    Input dimensions: $(FO.inDims)
    Output dimensions: $(FO.outDims)""")

Base.show(io::IO, FO::FunctionOperatorComposite) =
    print(io, "FunctionOperatorComposite{$(eltype(FO))}($(FO.name), $(FO.inDims), $(FO.outDims), $(FO.plan_string))")

Base.show(io::IO, ::MIME"text/plain", FO::FunctionOperatorComposite) =
    print(io, """
FunctionOperatorComposite with eltype $(eltype(FO))
    Name: $(FO.name)
    Input dimensions: $(FO.inDims)
    Output dimensions: $(FO.outDims)
    Plan: $(FO.plan_string)""")


text:
using QuantEcon
using NLsolve
using NLopt

import QuantEcon.simulate

mutable struct Model{TF <: AbstractFloat,
                     TM <: AbstractMatrix{TF},
                     TV <: AbstractVector{TF}}
    β::TF
    Π::TM
    G::TV
    Θ::TV
    transfers::Bool
    U::Function
    Uc::Function
    Ucc::Function
    Un::Function
    Unn::Function
    n_less_than_one::Bool
end

"""
{{ class_word }} returns planner's allocation as a function of the multiplier
on the implementability constraint μ
"""
struct SequentialAllocation{TP <: Model,
                            TI <: Integer,
                            TV <: AbstractVector}
    model::TP
    mc::MarkovChain
    S::TI
    cFB::TV
    nFB::TV
    ΞFB::TV
    zFB::TV
end

"""
Initializes the {{ class_word }} from the calibration model
"""
function SequentialAllocation(model::Model)
    β, Π, G, Θ = model.β, model.Π, model.G, model.Θ
    mc = MarkovChain(Π)
    S = size(Π, 1)   # Number of states
    # Now find the first best allocation
    cFB, nFB, ΞFB, zFB = find_first_best(model, S, 1)

    return SequentialAllocation(model, mc, S, cFB, nFB, ΞFB, zFB)
end

"""
Find the first best allocation
"""
function find_first_best(model::Model, S::Integer, version::Integer)
    if version != 1 && version != 2
        throw(ArgumentError("version must be 1 or 2"))
    end
    β, Θ, Uc, Un, G, Π =
        model.β, model.Θ, model.Uc, model.Un, model.G, model.Π
    function res!(z, out)
        c = z[1:S]
        n = z[S+1:end]
        out[1:S] = Θ .* Uc(c, n) + Un(c, n)
        out[S+1:end] = Θ .* n - c - G
    end
    res = nlsolve(res!, 0.5 * ones(2 * S))

    if converged(res) == false
        error("Could not find first best")
    end

    if version == 1
        cFB = res.zero[1:S]
        nFB = res.zero[S+1:end]
        ΞFB = Uc(cFB, nFB)         # Multiplier on the resource constraint
        zFB = vcat(cFB, nFB, ΞFB)
        return cFB, nFB, ΞFB, zFB
    elseif version == 2
        cFB = res.zero[1:S]
        nFB = res.zero[S+1:end]
        IFB = Uc(cFB, nFB) .* cFB + Un(cFB, nFB) .* nFB
        xFB = \(eye(S) - β * Π, IFB)
        zFB = [vcat(cFB[s], xFB[s], xFB) for s in 1:S]
        return cFB, nFB, IFB, xFB, zFB
    end
end

"""
Computes optimal allocation for time ``t\geq 1`` for a given ``\mu``
"""
function time1_allocation(pas::SequentialAllocation, μ::Real)
    model, S = pas.model, pas.S
    Θ, β, Π, G, Uc, Ucc, Un, Unn =
        model.Θ, model.β, model.Π, model.G,
        model.Uc, model.Ucc, model.Un, model.Unn
    function FOC!(z::Vector, out)
        c = z[1:S]
        n = z[S+1:2S]
        Ξ = z[2S+1:end]
        out[1:S] = Uc(c, n) - μ * (Ucc(c, n) .* c + Uc(c, n)) - Ξ            # FOC c
        out[S+1:2S] = Un(c, n) - μ * (Unn(c, n) .* n + Un(c, n)) + Θ .* Ξ    # FOC n
        out[2S+1:end] = Θ .* n - c - G                                       # Resource constraint
        return out
    end
    # Find the root of the FOC
    res = nlsolve(FOC!, pas.zFB)
    if res.f_converged == false
        error("Could not find LS allocation.")
    end
    z = res.zero
    c, n, Ξ = z[1:S], z[S+1:2S], z[2S+1:end]
    # Now compute x
    I  = Uc(c, n) .* c +  Un(c, n) .* n
    x = \(eye(S) - β * model.Π, I)
    return c, n, x, Ξ
end

"""
Finds the optimal allocation given initial government debt `B_` and state `s_0`
"""
function time0_allocation(pas::SequentialAllocation,
                          B_::AbstractFloat, s_0::Integer)
    model = pas.model
    Π, Θ, G, β =
        model.Π, model.Θ, model.G, model.β
    Uc, Ucc, Un, Unn =
        model.Uc, model.Ucc, model.Un, model.Unn
    # First order conditions of planner's problem
    function FOC!(z, out)
        μ, c, n, Ξ = z[1], z[2], z[3], z[4]
        xprime = time1_allocation(pas, μ)[3]
        out .= vcat(
        Uc(c, n) .* (c - B_) + Un(c, n) .* n + β * dot(Π[s_0, :], xprime),
        Uc(c, n) - μ * (Ucc(c, n) .* (c - B_) + Uc(c, n)) - Ξ,
        Un(c, n) - μ * (Unn(c, n) .* n + Un(c, n)) + Θ[s_0] .* Ξ,
        (Θ .* n - c - G)[s_0]
        )
    end
    # Find root
    res = nlsolve(FOC!, [0.0, pas.cFB[s_0], pas.nFB[s_0], pas.ΞFB[s_0]])
    if res.f_converged == false
        error("Could not find time 0 LS allocation.")
    end
    return (res.zero...)
end

"""
Find the value associated with multiplier `μ`
"""
function time1_value(pas::SequentialAllocation, μ::Real)
    model = pas.model
    c, n, x, Ξ = time1_allocation(pas, μ)
    U_val = model.U.(c, n)
    V = \(eye(pas.S) - model.β*model.Π, U_val)
    return c, n, x, V
end

"""
Computes Τ given `c`, `n`
"""
function Τ(model::Model, c::Union{Real,Vector}, n::Union{Real,Vector})
    Uc, Un = model.Uc.(c, n), model.Un.(c, n)
    return 1+Un./(model.Θ .* Uc)
end

"""
Simulates planners policies for `T` periods
"""
function simulate(pas::SequentialAllocation,
                  B_::AbstractFloat, s_0::Integer, T::Integer,
                  sHist::Union{Vector, Void}=nothing)
    model = pas.model
    Π, β, Uc = model.Π, model.β, model.Uc
    if sHist == nothing
        sHist = QuantEcon.simulate(pas.mc, T, init=s_0)
    end
    cHist = zeros(T)
    nHist = zeros(T)
    Bhist = zeros(T)
    ΤHist = zeros(T)
    μHist = zeros(T)
    RHist = zeros(T-1)
    # time 0
    μ, cHist[1], nHist[1], _  = time0_allocation(pas, B_, s_0)
    ΤHist[1] = Τ(pas.model, cHist[1], nHist[1])[s_0]
    Bhist[1] = B_
    μHist[1] = μ
    # time 1 onward
    for t in 2:T
        c, n, x, Ξ = time1_allocation(pas,μ)
        u_c = Uc(c,n)
        s = sHist[t]
        ΤHist[t] = Τ(pas.model, c, n)[s]
        Eu_c = dot(Π[sHist[t-1],:], u_c)
        cHist[t], nHist[t], Bhist[t] = c[s], n[s], x[s] / u_c[s]
        RHist[t-1] = Uc(cHist[t-1], nHist[t-1]) / (β * Eu_c)
        μHist[t] = μ
    end
    return cHist, nHist, Bhist, ΤHist, sHist, μHist, RHist
end

"""
Bellman equation for the continuation of the Lucas-Stokey Problem
"""
mutable struct BellmanEquation{TP <: Model,
                               TI <: Integer,
                               TV <: AbstractVector,
                               TM <: AbstractMatrix{TV},
                               TVV <: AbstractVector{TV}}
    model::TP
    S::TI
    xbar::TV
    time_0::Bool
    z0::TM
    cFB::TV
    nFB::TV
    xFB::TV
    zFB::TVV
end

"""
Initializes the {{ class_word }} from the calibration `model`
"""
function BellmanEquation(model::Model, xgrid::AbstractVector, policies0::Vector)
    S = size(model.Π, 1) # Number of states
    xbar = [minimum(xgrid), maximum(xgrid)]
    time_0 = false
    cf, nf, xprimef = policies0
    z0 = [vcat(cf[s](x), nf[s](x), [xprimef[s, sprime](x) for sprime in 1:S])
                        for x in xgrid, s in 1:S]
    cFB, nFB, IFB, xFB, zFB = find_first_best(model, S, 2)
    return BellmanEquation(model, S, xbar, time_0, z0, cFB, nFB, xFB, zFB)
end

"""
Finds the optimal policies
"""
function get_policies_time1(T::BellmanEquation,
                        i_x::Integer, x::AbstractFloat,
                        s::Integer, Vf::AbstractArray)
    model, S = T.model, T.S
    β, Θ, G, Π = model.β, model.Θ, model.G, model.Π
    U, Uc, Un = model.U, model.Uc, model.Un

    function objf(z::Vector, grad)
        c, xprime = z[1], z[2:end]
        n=c+G[s]
        Vprime = [Vf[sprime](xprime[sprime]) for sprime in 1:S]
        return -(U(c, n) + β * dot(Π[s, :], Vprime))
    end
    function cons(z::Vector, grad)
        c, xprime = z[1], z[2:end]
        n=c+G[s]
        return x - Uc(c, n) * c - Un(c, n) * n - β * dot(Π[s, :], xprime)
    end
    lb = vcat(0, T.xbar[1] * ones(S))
    ub = vcat(1 - G[s], T.xbar[2] * ones(S))
    opt = Opt(:LN_COBYLA, length(T.z0[i_x, s])-1)
    min_objective!(opt, objf)
    equality_constraint!(opt, cons)
    lower_bounds!(opt, lb)
    upper_bounds!(opt, ub)
    maxeval!(opt, 300)
    maxtime!(opt, 10)
    init = vcat(T.z0[i_x, s][1], T.z0[i_x, s][3:end])
    for (i, val) in enumerate(init)
        if val > ub[i]
            init[i] = ub[i]
        elseif val < lb[i]
            init[i] = lb[i]
        end
    end
    (minf, minx, ret) = optimize(opt, init)
    T.z0[i_x, s] = vcat(minx[1], minx[1] + G[s], minx[2:end])
    return vcat(-minf, T.z0[i_x, s])
end
"""
Finds the optimal policies
"""
function get_policies_time0(T::BellmanEquation,
                        B_::AbstractFloat, s0::Integer, Vf::Array)
    model, S = T.model, T.S
    β, Θ, G, Π = model.β, model.Θ, model.G, model.Π
    U, Uc, Un = model.U, model.Uc, model.Un
    function objf(z, grad)
        c, xprime = z[1], z[2:end]
        n = c+G[s0]
        Vprime = [Vf[sprime](xprime[sprime]) for sprime in 1:S]
        return -(U(c, n) + β * dot(Π[s0, :], Vprime))
    end
    function cons(z::Vector, grad)
        c, xprime = z[1], z[2:end]
        n = c + G[s0]
        return -Uc(c, n) * (c - B_) - Un(c, n) * n - β * dot(Π[s0, :], xprime)
    end
    lb = vcat(0, T.xbar[1] * ones(S))
    ub = vcat(1-G[s0], T.xbar[2] * ones(S))
    opt = Opt(:LN_COBYLA, length(T.zFB[s0])-1)
    min_objective!(opt, objf)
    equality_constraint!(opt, cons)
    lower_bounds!(opt, lb)
    upper_bounds!(opt, ub)
    maxeval!(opt, 300)
    maxtime!(opt, 10)
    init = vcat(T.zFB[s0][1], T.zFB[s0][3:end])
    for (i, val) in enumerate(init)
        if val > ub[i]
            init[i] = ub[i]
        elseif val < lb[i]
            init[i] = lb[i]
        end
    end
    (minf, minx, ret) = optimize(opt, init)
    return vcat(-minf, vcat(minx[1], minx[1]+G[s0], minx[2:end]))
end


text:
using ADMPS
using Test

@testset "ADMPS.jl" begin
    @testset "cuda_patch.jl" begin
        println("cuda_patch tests running...")
        include("cuda_patch.jl")
    end

    @testset "hamiltonianmodels" begin
        println("hamiltonianmodels tests running...")
        include("hamiltonianmodels.jl")
    end

    @testset "exampletensors" begin
        println("exampletensors tests running...")
        include("exampletensors.jl")
    end

    @testset "environment" begin
        println("environment tests running...")
        include("environment.jl")
    end

    @testset "autodiff" begin
        println("autodiff tests running...")
        include("autodiff.jl")
    end

    @testset "variationalmps" begin
        println("variationalmps tests running...")
        include("variationalmps.jl")
    end
end


text:
using MAT
using LinearAlgebra
using StatsBase
using Plots
include("../src/SparseCard.jl")
include("image_graphs_utils.jl")

## Load optimal solution to the problem, for comparison
numsuperpix = 200

if numsuperpix == 200
    lambda = 0.002
    K = 13
else
    lambda = 0.008
    K = 16
end

dataset = "001_oct"
short = "oct"
G = matread("GraphPotentials/Graph_noclique_$(dataset)_$numsuperpix.mat")
An = G["A"]
tvec = G["tvec"]
svec = G["svec"]
n = size(An,1)
In,Jn,Wn = findnz(An)
sourceside = sum(svec)
Edges, EdgesW = load_superpixel_cliques(dataset,lambda,numsuperpix)

Opt = matread("OptimalSolutions/Optimal_$(dataset)_$numsuperpix.mat")
S_opt = Opt["S"]
opttime = Opt["opttime"]
total_edges = Opt["total_edges"]
opt = eval_fS_image(Edges,EdgesW,svec,tvec,An,S_opt)

## Load and plot results

pl = plot()

c1 = 10     # ACDM default
c2 = 25     # ACDM best
T = 600
Col = [1 0 0; 0 .75 0; 0 0 1; 1 .5 0]
timelimit = 62
methods = ["APcompact_T_$T", "RCDM_$K", "ACDM_$(K)_$c1","ACDM_$(K)_$c2"]
methodname = ["IAP","RCDM", "ACDM-default","ACDM-best",]

for j = 1:4
col = RGBA(Col[j,1],Col[j,2],Col[j,3],1)
ms = 3.5
method = methods[j]
if j == 1
    jump = 6
else
    jump = 20
end

method = methods[j]
C = matread("ContOutput/ContOutput_$(short)$(numsuperpix)/$(method)_1.mat")
levelcost = vec(C["levelcost"])
L = length(levelcost)
Runs = zeros(5,L)
Objs = zeros(5,L)

for jj = 1:5
    C = matread("ContOutput/ContOutput_$(short)$(numsuperpix)/$(method)_$jj.mat")
    levelcost = vec(C["levelcost"])
    timevec = vec(C["timevec"])/1000
    Runs[jj,:] = timevec

    # objective is levelcost + weight of source-side edges
    Objs[jj,:] = levelcost .+ sourceside
end

timevec = vec(mean(Runs,dims = 1))
Fvals = vec(mean(Objs,dims = 1))
Upper = vec(maximum(Objs,dims = 1))./opt .-1
Lower = vec(minimum(Objs,dims = 1))./opt .-1
Approx = Fvals./opt .- 1
xrange = 1:jump:length(Approx)

# Truncate if over time limit
if maximum(timevec)>timelimit
    lst = findfirst(x->x>timelimit,timevec)
    st = rem(lst,jump)
    timevec = timevec[1:lst]
    Approx = Approx[1:lst]
    Lower = Lower[1:lst]
    Upper = Upper[1:lst]
    xrange = 1:jump:lst
end

# Truncate and adjust if exact convergence is reached
if minimum(Approx) <= 0
    converged = findfirst(x->x<= 0,Approx)
    # avoids log(0) errors without changing plot
    Approx[converged] = 1e-12
    Lower[converged] = 1e-12
    Upper[converged] = 1e-12
    st = rem(converged,jump)
    xrange = st:jump:converged
end

label = methodname[j]
alp = 1.0
plot!(pl,timevec[xrange],Approx[xrange],
    label = label,alpha = alp,
    markershape = :circle,markersize = ms,markerstrokewidth = 0.0,
    markercolor = col, markerstrokecolor = col, color = col,yscale = :log10)

plot!(pl,timevec[xrange],Lower[xrange],alpha = 0,
    fillrange=Upper[xrange], fillalpha=0.3,
    markerstrokecolor = col, color = col,linewidth = 0.0,label = "")

end
pl

## Plot sparsecard results
ms = 3.5
M = matread("Output/SparseCard_$(dataset)_$(numsuperpix)_4.mat")
runtimes = M["runtimes"]
objvals = M["objvals"]
epsilons = M["epsilons"]
flow_approx = (objvals./opt) .- 1
J = findall(x->x==0.0,flow_approx)
flow_approx[J] .= 1e-12

p = sortperm(runtimes)
plot!(pl, runtimes[p[1:end-2]], flow_approx[p[1:end-2]], label = "SparseCard",
    xlabel = "Runtime (seconds)",ylabel = "Approx Ratio - 1",yscale = :log10,
    markershape = :circle,markersize = ms,markerstrokewidth = 0.0,
    markercolor = :black, color = :black)

## Adjust size and save
s1 = 300
s2 = 250
gfs = 12
tfs = 12
plot!(pl,xlimit = [0,timelimit+1],ylimit = [5e-7,.5],yscale = :log10,
        grid = false, size = (s1,s2),fontfamily = "helvetica",
        title = "\n oct, $numsuperpix superpixels",
        legend = false,tickfontsize=tfs,guidefontsize = gfs)

# savefig("Figures/RunApprox_$(dataset)_$(numsuperpix).pdf")

max_stars_repo_name: nveldt/SparseCardDSFM

max_issues_repo_name: nveldt/SparseCardDSFM

max_forks_repo_name: nveldt/SparseCardDSFM


text:
using Test
using OMEinsum
using OMEinsum: IndexSize, get_size_dict

@testset "index size" begin
    is = IndexSize(('i', 'j'), (2, 3))
    @test is['i'] == 2
    is = is + is
    @test is == IndexSize(('i', 'j', 'i', 'j'), (2, 3, 2, 3))
    @test is['i'] == 2
    a, b, c = randn(2,3), randn(3,4), randn(2,4)
    xs = (a, b, c)
    ixs = (('i', 'j'), ('j', 'k'), ('i', 'k'))
    is = get_size_dict(ixs, xs)
    @test is == IndexSize(('i', 'j', 'j', 'k', 'i', 'k'), (2,3,3,4,2,4))
    @test is['i'] == 2
    @test is['j'] == 3
    @test is['k'] == 4

    @test IndexSize('i'=>2) == IndexSize(('i',), (2,))
    noindex = IndexSize((), ())
    @test noindex isa IndexSize
    @test noindex + is == is
end

@testset "unspecified index sizes" begin
    v = randn(5)
    a = randn(6,6)
    xs = (v, a)
    ixs = (('j',), ('k', 'k'))
    @test einsum(ein"j,kk->j", xs, get_size_dict(ixs, xs)) ≈ ein"j,kk->j"(xs...)
    res = ein"j,kk->iij"(xs...; size_info=IndexSize('i'=>9))
    @test size(res) == (9,9,5)
end


text:
#data set has names at the beginning of it. Since the parts of the names have spaces in between them and a space is the seperator for the
#rest of the file, it gets confused how many colums exist

dataset_path = joinpath(dirname(@__FILE__), "immigrant-salaries.data.orig")

f = open("$dataset_path.custom", "w")
for line in readlines(dataset_path)
  line = line[30:end]
  line = lstrip(line)
  write(f, "$line\n")
end
close(f)

max_stars_repo_name: CoventryResearch/uci-data

max_issues_repo_name: CoventryResearch/uci-data

max_forks_repo_name: CoventryResearch/uci-data


text:
using DashBootstrapComponents, DashHtmlComponents

row = html_div([
    dbc_row(dbc_col(html_div("A single, half-width column"), width=6)),
    dbc_row(dbc_col(html_div("An automatically sized column"), width="auto")),
    dbc_row([
        dbc_col(html_div("One of three columns"), width=3),
        dbc_col(html_div("One of three columns")),
        dbc_col(html_div("One of three columns"), width=3),
    ]),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
using Test
using KnetMetrics
using Random

@testset "pairwise-metrics" begin
    Random.seed!(42)
    u = rand(100)
    Random.seed!(123)
    v = rand(100)

    @testset "minkowski-distance" begin
        @test isapprox(minkowski_distance(u,v), 19.192858878774214)
    end

    @testset "minkowski-distance-p3" begin
        @test isapprox(minkowski_distance(u,v, p=3), 11.841710005717449)
    end

    @testset "euclidian-distance" begin
        @test isapprox(euclidian_distance(u,v), 19.192858878774214)
    end

    @testset "manhattan-distance" begin
        @test isapprox(manhattan_distance(u,v), 36.666029519718066)
    end

    @testset "chebyshev-distance" begin
        @test isapprox(chebyshev_distance(u,v), 0.8624100780860133)
    end

    @testset "braycurtis-distance" begin
        @test isapprox(braycurtis_distance(u,v), 0.3653169890723984)
    end

    @testset "canberra-distance" begin
        @test isapprox(canberra_distance(u,v), 32.09459846314555)
    end

    @testset "cityblock-distance" begin
        @test isapprox(cityblock_distance(u,v), 36.666029519718066)
    end

    @testset "correlation" begin
        @test isapprox(correlation(u,v), 1.1508533815000055)
    end

    @testset "cosine-similarity" begin
        @test isapprox(cosine_similarity(u,v), -0.1508533815000055)
    end

end


text:
mutable struct StyleAndLanguageManager
    languageDefinitions::Dict{AbstractString, GtkSourceWidget.GtkSourceLanguage}
    main_style::GtkSourceWidget.GtkSourceStyleScheme
    fontsize::Int
    fontCss::String
    style_provider::GtkCssProvider
    styles::Dict{String, GtkSourceWidget.GtkSourceStyle}

    function StyleAndLanguageManager()

        #FIXME this should be in GtkSourceWidget
        sourceStyleManager = GtkSourceStyleSchemeManager()
        GtkSourceWidget.set_search_path(sourceStyleManager,
            Any[joinpath(GtkIDE.HOMEDIR, "../share/styles/"), C_NULL])

        languageDefinitions = Dict{AbstractString, GtkSourceWidget.GtkSourceLanguage}()
        sourceLanguageManager = GtkSourceWidget.sourceLanguageManager
        GtkSourceWidget.set_search_path(sourceLanguageManager,
            Any[joinpath(GtkIDE.HOMEDIR, "../share/language-specs/"), C_NULL])

        languageDefinitions[".jl"] = GtkSourceWidget.language(sourceLanguageManager, "julia")
        languageDefinitions[".md"] = GtkSourceWidget.language(sourceLanguageManager, "markdown")

        if Sys.iswindows()
            main_style = style_scheme(sourceStyleManager, "visualcode")
            fontsize = opt("fontsize")
            fontCss =  """button, entry, window, sourceview, textview  {
                font-family: Consolas, Courier, monospace;
                font-size: $(fontsize)px;
            }"""
        end
        if Sys.isapple()
            main_style = style_scheme(sourceStyleManager, "visualcode")
            fontsize = opt("fontsize")
            fontCss =  "button, entry, window, sourceview, textview {
                font-family: Monaco, Consolas, Courier, monospace;
                font-size: $(fontsize)px;
            }"
        end

        if Sys.islinux()
            main_style = style_scheme(sourceStyleManager, "visualcode")
            fontsize = opt("fontsize")-1
            fontCss =  """GtkButton, GtkEntry, GtkWindow, GtkSourceView, GtkTextView {
                font-family: Consolas, Courier, monospace;
                font-size: $(fontsize)px;
            }"""
        end

        #provider = GtkExtensions.default_css_provider
        #provider = GtkCssProviderFromData!(provider, data=fontCss)
        provider = GtkCssProvider(data = fontCss)
        Gtk.icon_theme_append_search_path(Gtk.icon_theme_get_default(), joinpath(HOMEDIR, "../icons/"))

        # I'm getting duplicated GObject when calling style several times,
        # so let's call it only once.
        styles = Dict(
            "text" => GtkSourceWidget.style(main_style, "text"),
            "def:note" => GtkSourceWidget.style(main_style, "def:note")
        )

        new(
            languageDefinitions,
            main_style,
            fontsize,
            fontCss,
            provider,
            styles
        )
    end
end

get_style(mng::StyleAndLanguageManager, style_id::String) =
    haskey(mng.styles, style_id) ? mng.styles[style_id] : GtkSourceWidget.style(mng.main_style, style_id)


text:
@mkCppMethod IVRSystem GetRecommendedRenderTargetSize
@mkCppMethod IVRSystem GetProjectionMatrix
@mkCppMethod IVRSystem GetProjectionRaw
@mkCppMethod IVRSystem ComputeDistortion
@mkCppMethod IVRSystem GetEyeToHeadTransform
@mkCppMethod IVRSystem GetTimeSinceLastVsync
@mkCppMethod IVRSystem GetD3D9AdapterIndex
@mkCppMethod IVRSystem GetDXGIOutputInfo
@mkCppMethod IVRSystem GetOutputDevice
@mkCppMethod IVRSystem IsDisplayOnDesktop
@mkCppMethod IVRSystem SetDisplayVisibility
@mkCppMethod IVRSystem GetDeviceToAbsoluteTrackingPose
@mkCppMethod IVRSystem ResetSeatedZeroPose
@mkCppMethod IVRSystem GetSeatedZeroPoseToStandingAbsoluteTrackingPose
@mkCppMethod IVRSystem GetRawZeroPoseToStandingAbsoluteTrackingPose
@mkCppMethod IVRSystem GetSortedTrackedDeviceIndicesOfClass
@mkCppMethod IVRSystem GetTrackedDeviceActivityLevel
@mkCppMethod IVRSystem ApplyTransform
@mkCppMethod IVRSystem GetTrackedDeviceIndexForControllerRole
@mkCppMethod IVRSystem GetControllerRoleForTrackedDeviceIndex
@mkCppMethod IVRSystem GetTrackedDeviceClass
@mkCppMethod IVRSystem IsTrackedDeviceConnected
@mkCppMethod IVRSystem GetBoolTrackedDeviceProperty
@mkCppMethod IVRSystem GetFloatTrackedDeviceProperty
@mkCppMethod IVRSystem GetInt32TrackedDeviceProperty
@mkCppMethod IVRSystem GetUint64TrackedDeviceProperty
@mkCppMethod IVRSystem GetMatrix34TrackedDeviceProperty
@mkCppMethod IVRSystem GetArrayTrackedDeviceProperty
@mkCppMethod IVRSystem GetStringTrackedDeviceProperty
@mkCppMethod IVRSystem GetPropErrorNameFromEnum
@mkCppMethod IVRSystem PollNextEvent
@mkCppMethod IVRSystem PollNextEventWithPose
@mkCppMethod IVRSystem GetEventTypeNameFromEnum
@mkCppMethod IVRSystem GetHiddenAreaMesh
@mkCppMethod IVRSystem GetControllerState
@mkCppMethod IVRSystem GetControllerStateWithPose
@mkCppMethod IVRSystem TriggerHapticPulse
@mkCppMethod IVRSystem GetButtonIdNameFromEnum
@mkCppMethod IVRSystem GetControllerAxisTypeNameFromEnum
@mkCppMethod IVRSystem IsInputAvailable
@mkCppMethod IVRSystem IsSteamVRDrawingControllers
@mkCppMethod IVRSystem ShouldApplicationPause
@mkCppMethod IVRSystem ShouldApplicationReduceRenderingWork
# @mkCppMethod IVRSystem DriverDebugRequest
@mkCppMethod IVRSystem PerformFirmwareUpdate
@mkCppMethod IVRSystem AcknowledgeQuit_Exiting
# @mkCppMethod IVRSystem AcknowledgeQuit_UserPrompt

@mkCppMethod IVRApplications AddApplicationManifest
@mkCppMethod IVRApplications RemoveApplicationManifest
@mkCppMethod IVRApplications IsApplicationInstalled
@mkCppMethod IVRApplications GetApplicationCount
@mkCppMethod IVRApplications GetApplicationKeyByIndex
@mkCppMethod IVRApplications GetApplicationKeyByProcessId
@mkCppMethod IVRApplications LaunchApplication
@mkCppMethod IVRApplications LaunchTemplateApplication
@mkCppMethod IVRApplications LaunchApplicationFromMimeType
@mkCppMethod IVRApplications LaunchDashboardOverlay
@mkCppMethod IVRApplications CancelApplicationLaunch
@mkCppMethod IVRApplications IdentifyApplication
@mkCppMethod IVRApplications GetApplicationProcessId
@mkCppMethod IVRApplications GetApplicationsErrorNameFromEnum
@mkCppMethod IVRApplications GetApplicationPropertyString
@mkCppMethod IVRApplications GetApplicationPropertyBool
@mkCppMethod IVRApplications GetApplicationPropertyUint64
@mkCppMethod IVRApplications SetApplicationAutoLaunch
@mkCppMethod IVRApplications GetApplicationAutoLaunch
@mkCppMethod IVRApplications SetDefaultApplicationForMimeType
@mkCppMethod IVRApplications GetDefaultApplicationForMimeType
@mkCppMethod IVRApplications GetApplicationSupportedMimeTypes
@mkCppMethod IVRApplications GetApplicationsThatSupportMimeType
@mkCppMethod IVRApplications GetApplicationLaunchArguments
@mkCppMethod IVRApplications GetStartingApplication
# @mkCppMethod IVRApplications GetTransitionState
@mkCppMethod IVRApplications PerformApplicationPrelaunchCheck
# @mkCppMethod IVRApplications GetApplicationsTransitionStateNameFromEnum
# @mkCppMethod IVRApplications IsQuitUserPromptRequested
@mkCppMethod IVRApplications LaunchInternalProcess
@mkCppMethod IVRApplications GetCurrentSceneProcessId

@mkCppMethod IVRSettings GetSettingsErrorNameFromEnum
# @mkCppMethod IVRSettings Sync
@mkCppMethod IVRSettings SetBool
@mkCppMethod IVRSettings SetInt32
@mkCppMethod IVRSettings SetFloat
@mkCppMethod IVRSettings SetString
@mkCppMethod IVRSettings GetBool
@mkCppMethod IVRSettings GetInt32
@mkCppMethod IVRSettings GetFloat
@mkCppMethod IVRSettings GetString
@mkCppMethod IVRSettings RemoveSection
@mkCppMethod IVRSettings RemoveKeyInSection

@mkCppMethod CVRSettingHelper GetSettingsErrorNameFromEnum
# @mkCppMethod CVRSettingHelper Sync
@mkCppMethod CVRSettingHelper SetBool
@mkCppMethod CVRSettingHelper SetInt32
@mkCppMethod CVRSettingHelper SetFloat
@mkCppMethodCast CVRSettingHelper SetString "void" "(const char *, const char *, const char *, vr::EVRSettingsError *)"
@mkCppMethodCast CVRSettingHelper SetString "void" "(const std::string &, const std::string &, const std::string &, vr::EVRSettingsError *)"
@mkCppMethod CVRSettingHelper GetBool
@mkCppMethod CVRSettingHelper GetInt32
@mkCppMethod CVRSettingHelper GetFloat
@mkCppMethodCast CVRSettingHelper GetString "void"        "(const char *, const char *, char *, uint32_t, vr::EVRSettingsError *)"
@mkCppMethodCast CVRSettingHelper GetString "std::string" "(const std::string &, const std::string &, vr::EVRSettingsError *)"
@mkCppMethod CVRSettingHelper RemoveSection
@mkCppMethod CVRSettingHelper RemoveKeyInSection

@mkCppMethod IVRChaperone GetCalibrationState
@mkCppMethod IVRChaperone GetPlayAreaSize
@mkCppMethod IVRChaperone GetPlayAreaRect
@mkCppMethod IVRChaperone ReloadInfo
@mkCppMethod IVRChaperone SetSceneColor
@mkCppMethod IVRChaperone GetBoundsColor
@mkCppMethod IVRChaperone AreBoundsVisible
@mkCppMethod IVRChaperone ForceBoundsVisible

@mkCppMethod IVRChaperoneSetup CommitWorkingCopy
@mkCppMethod IVRChaperoneSetup RevertWorkingCopy
@mkCppMethod IVRChaperoneSetup GetWorkingPlayAreaSize
@mkCppMethod IVRChaperoneSetup GetWorkingPlayAreaRect
@mkCppMethod IVRChaperoneSetup GetWorkingCollisionBoundsInfo
@mkCppMethod IVRChaperoneSetup GetLiveCollisionBoundsInfo
@mkCppMethod IVRChaperoneSetup GetWorkingSeatedZeroPoseToRawTrackingPose
@mkCppMethod IVRChaperoneSetup GetWorkingStandingZeroPoseToRawTrackingPose
@mkCppMethod IVRChaperoneSetup SetWorkingPlayAreaSize
@mkCppMethod IVRChaperoneSetup SetWorkingCollisionBoundsInfo
@mkCppMethod IVRChaperoneSetup SetWorkingPerimeter
@mkCppMethod IVRChaperoneSetup SetWorkingSeatedZeroPoseToRawTrackingPose
@mkCppMethod IVRChaperoneSetup SetWorkingStandingZeroPoseToRawTrackingPose
@mkCppMethod IVRChaperoneSetup ReloadFromDisk
@mkCppMethod IVRChaperoneSetup GetLiveSeatedZeroPoseToRawTrackingPose
@mkCppMethod IVRChaperoneSetup ExportLiveToBuffer
@mkCppMethod IVRChaperoneSetup ImportFromBufferToWorking
@mkCppMethod IVRChaperoneSetup ShowWorkingSetPreview
@mkCppMethod IVRChaperoneSetup HideWorkingSetPreview

@mkCppMethod IVRCompositor SetTrackingSpace
@mkCppMethod IVRCompositor GetTrackingSpace
@mkCppMethod IVRCompositor WaitGetPoses
@mkCppMethod IVRCompositor GetLastPoses
@mkCppMethod IVRCompositor GetLastPoseForTrackedDeviceIndex
@mkCppMethod IVRCompositor Submit
@mkCppMethod IVRCompositor ClearLastSubmittedFrame
@mkCppMethod IVRCompositor PostPresentHandoff
@mkCppMethod IVRCompositor GetFrameTiming
@mkCppMethod IVRCompositor GetFrameTimings
@mkCppMethod IVRCompositor GetFrameTimeRemaining
@mkCppMethod IVRCompositor GetCumulativeStats
@mkCppMethod IVRCompositor FadeToColor
@mkCppMethod IVRCompositor GetCurrentFadeColor
@mkCppMethod IVRCompositor FadeGrid
@mkCppMethod IVRCompositor GetCurrentGridAlpha
@mkCppMethod IVRCompositor SetSkyboxOverride
@mkCppMethod IVRCompositor ClearSkyboxOverride
@mkCppMethod IVRCompositor CompositorBringToFront
@mkCppMethod IVRCompositor CompositorGoToBack
@mkCppMethod IVRCompositor CompositorQuit
@mkCppMethod IVRCompositor IsFullscreen
@mkCppMethod IVRCompositor GetCurrentSceneFocusProcess
@mkCppMethod IVRCompositor GetLastFrameRenderer
@mkCppMethod IVRCompositor CanRenderScene
@mkCppMethod IVRCompositor ShowMirrorWindow
@mkCppMethod IVRCompositor HideMirrorWindow
@mkCppMethod IVRCompositor IsMirrorWindowVisible
@mkCppMethod IVRCompositor CompositorDumpImages
@mkCppMethod IVRCompositor ShouldAppRenderWithLowResources
@mkCppMethod IVRCompositor ForceInterleavedReprojectionOn
@mkCppMethod IVRCompositor ForceReconnectProcess
@mkCppMethod IVRCompositor SuspendRendering
@mkCppMethod IVRCompositor GetMirrorTextureD3D11
@mkCppMethod IVRCompositor ReleaseMirrorTextureD3D11
@mkCppMethod IVRCompositor GetMirrorTextureGL
@mkCppMethod IVRCompositor ReleaseSharedGLTexture
@mkCppMethod IVRCompositor LockGLSharedTextureForAccess
@mkCppMethod IVRCompositor UnlockGLSharedTextureForAccess
@mkCppMethod IVRCompositor GetVulkanInstanceExtensionsRequired
#   // .method("GetVulkanDeviceExtensionsRequired", &vr::IVRCompositor::GetVulkanDeviceExtensionsRequired) // VkPhysicalDevice_T ↯
@mkCppMethod IVRCompositor SetExplicitTimingMode
@mkCppMethod IVRCompositor SubmitExplicitTimingData
@mkCppMethod IVRCompositor IsMotionSmoothingEnabled

@mkCppMethod IVRNotifications CreateNotification
@mkCppMethod IVRNotifications RemoveNotification

@mkCppMethod IVROverlay FindOverlay
@mkCppMethod IVROverlay CreateOverlay
@mkCppMethod IVROverlay DestroyOverlay
# @mkCppMethod IVROverlay SetHighQualityOverlay
# @mkCppMethod IVROverlay GetHighQualityOverlay
@mkCppMethod IVROverlay GetOverlayKey
@mkCppMethod IVROverlay GetOverlayName
@mkCppMethod IVROverlay SetOverlayName
@mkCppMethod IVROverlay GetOverlayImageData
@mkCppMethod IVROverlay GetOverlayErrorNameFromEnum
@mkCppMethod IVROverlay SetOverlayRenderingPid
@mkCppMethod IVROverlay GetOverlayRenderingPid
@mkCppMethod IVROverlay SetOverlayFlag
@mkCppMethod IVROverlay GetOverlayFlag
@mkCppMethod IVROverlay SetOverlayColor
@mkCppMethod IVROverlay GetOverlayColor
@mkCppMethod IVROverlay SetOverlayAlpha
@mkCppMethod IVROverlay GetOverlayAlpha
@mkCppMethod IVROverlay SetOverlayTexelAspect
@mkCppMethod IVROverlay GetOverlayTexelAspect
@mkCppMethod IVROverlay SetOverlaySortOrder
@mkCppMethod IVROverlay GetOverlaySortOrder
@mkCppMethod IVROverlay SetOverlayWidthInMeters
@mkCppMethod IVROverlay GetOverlayWidthInMeters
# @mkCppMethod IVROverlay SetOverlayAutoCurveDistanceRangeInMeters
# @mkCppMethod IVROverlay GetOverlayAutoCurveDistanceRangeInMeters
@mkCppMethod IVROverlay SetOverlayTextureColorSpace
@mkCppMethod IVROverlay GetOverlayTextureColorSpace
@mkCppMethod IVROverlay SetOverlayTextureBounds
@mkCppMethod IVROverlay GetOverlayTextureBounds
@mkCppMethod IVROverlay GetOverlayRenderModel
@mkCppMethod IVROverlay SetOverlayRenderModel
@mkCppMethod IVROverlay GetOverlayTransformType
@mkCppMethod IVROverlay SetOverlayTransformAbsolute
@mkCppMethod IVROverlay GetOverlayTransformAbsolute
@mkCppMethod IVROverlay SetOverlayTransformTrackedDeviceRelative
@mkCppMethod IVROverlay GetOverlayTransformTrackedDeviceRelative
@mkCppMethod IVROverlay SetOverlayTransformTrackedDeviceComponent
@mkCppMethod IVROverlay GetOverlayTransformTrackedDeviceComponent
@mkCppMethod IVROverlay GetOverlayTransformOverlayRelative
@mkCppMethod IVROverlay SetOverlayTransformOverlayRelative
@mkCppMethod IVROverlay ShowOverlay
@mkCppMethod IVROverlay HideOverlay
@mkCppMethod IVROverlay IsOverlayVisible
@mkCppMethod IVROverlay GetTransformForOverlayCoordinates
@mkCppMethod IVROverlay PollNextOverlayEvent
@mkCppMethod IVROverlay GetOverlayInputMethod
@mkCppMethod IVROverlay SetOverlayInputMethod
@mkCppMethod IVROverlay GetOverlayMouseScale
@mkCppMethod IVROverlay SetOverlayMouseScale
@mkCppMethod IVROverlay ComputeOverlayIntersection
@mkCppMethod IVROverlay IsHoverTargetOverlay
# @mkCppMethod IVROverlay GetGamepadFocusOverlay
# @mkCppMethod IVROverlay SetGamepadFocusOverlay
# @mkCppMethod IVROverlay SetOverlayNeighbor
# @mkCppMethod IVROverlay MoveGamepadFocusToNeighbor
@mkCppMethod IVROverlay SetOverlayDualAnalogTransform
@mkCppMethod IVROverlay GetOverlayDualAnalogTransform
@mkCppMethod IVROverlay SetOverlayTexture
@mkCppMethod IVROverlay ClearOverlayTexture
@mkCppMethod IVROverlay SetOverlayRaw
@mkCppMethod IVROverlay SetOverlayFromFile
@mkCppMethod IVROverlay GetOverlayTexture
@mkCppMethod IVROverlay ReleaseNativeOverlayHandle
@mkCppMethod IVROverlay GetOverlayTextureSize
@mkCppMethod IVROverlay CreateDashboardOverlay
@mkCppMethod IVROverlay IsDashboardVisible
@mkCppMethod IVROverlay IsActiveDashboardOverlay
@mkCppMethod IVROverlay SetDashboardOverlaySceneProcess
@mkCppMethod IVROverlay GetDashboardOverlaySceneProcess
@mkCppMethod IVROverlay ShowDashboard
@mkCppMethod IVROverlay GetPrimaryDashboardDevice
@mkCppMethod IVROverlay ShowKeyboard
@mkCppMethod IVROverlay ShowKeyboardForOverlay
@mkCppMethod IVROverlay GetKeyboardText
@mkCppMethod IVROverlay HideKeyboard
@mkCppMethod IVROverlay SetKeyboardTransformAbsolute
@mkCppMethod IVROverlay SetKeyboardPositionForOverlay
@mkCppMethod IVROverlay SetOverlayIntersectionMask
@mkCppMethod IVROverlay GetOverlayFlags
@mkCppMethod IVROverlay ShowMessageOverlay
@mkCppMethod IVROverlay CloseMessageOverlay

@mkCppMethod IVRRenderModels LoadRenderModel_Async
@mkCppMethod IVRRenderModels FreeRenderModel
@mkCppMethod IVRRenderModels LoadTexture_Async
@mkCppMethod IVRRenderModels FreeTexture
@mkCppMethod IVRRenderModels LoadTextureD3D11_Async
@mkCppMethod IVRRenderModels LoadIntoTextureD3D11_Async
@mkCppMethod IVRRenderModels FreeTextureD3D11
@mkCppMethod IVRRenderModels GetRenderModelName
@mkCppMethod IVRRenderModels GetRenderModelCount
@mkCppMethod IVRRenderModels GetComponentCount
@mkCppMethod IVRRenderModels GetComponentName
@mkCppMethod IVRRenderModels GetComponentButtonMask
@mkCppMethod IVRRenderModels GetComponentRenderModelName
@mkCppMethod IVRRenderModels GetComponentStateForDevicePath
@mkCppMethod IVRRenderModels GetComponentState
@mkCppMethod IVRRenderModels RenderModelHasComponent
@mkCppMethod IVRRenderModels GetRenderModelThumbnailURL
@mkCppMethod IVRRenderModels GetRenderModelOriginalPath
@mkCppMethod IVRRenderModels GetRenderModelErrorNameFromEnum

@mkCppMethod IVRExtendedDisplay GetWindowBounds
@mkCppMethod IVRExtendedDisplay GetEyeOutputViewport
@mkCppMethod IVRExtendedDisplay GetDXGIOutputInfo

@mkCppMethod IVRTrackedCamera GetCameraErrorNameFromEnum
@mkCppMethod IVRTrackedCamera HasCamera
@mkCppMethod IVRTrackedCamera GetCameraFrameSize
@mkCppMethod IVRTrackedCamera GetCameraIntrinsics
@mkCppMethod IVRTrackedCamera GetCameraProjection
@mkCppMethod IVRTrackedCamera AcquireVideoStreamingService
@mkCppMethod IVRTrackedCamera ReleaseVideoStreamingService
@mkCppMethod IVRTrackedCamera GetVideoStreamFrameBuffer
@mkCppMethod IVRTrackedCamera GetVideoStreamTextureSize
@mkCppMethod IVRTrackedCamera GetVideoStreamTextureD3D11
@mkCppMethod IVRTrackedCamera GetVideoStreamTextureGL
@mkCppMethod IVRTrackedCamera ReleaseVideoStreamTextureGL

@mkCppMethod IVRScreenshots RequestScreenshot
@mkCppMethod IVRScreenshots HookScreenshot
@mkCppMethod IVRScreenshots GetScreenshotPropertyType
@mkCppMethod IVRScreenshots GetScreenshotPropertyFilename
@mkCppMethod IVRScreenshots UpdateScreenshotProgress
@mkCppMethod IVRScreenshots TakeStereoScreenshot
@mkCppMethod IVRScreenshots SubmitScreenshot

@mkCppMethod IVRResources LoadSharedResource
@mkCppMethod IVRResources GetResourceFullPath

@mkCppMethod IVRDriverManager GetDriverCount
@mkCppMethod IVRDriverManager GetDriverName
@mkCppMethod IVRDriverManager GetDriverHandle

@mkCppMethod IVRInput SetActionManifestPath
@mkCppMethod IVRInput GetActionSetHandle
@mkCppMethod IVRInput GetActionHandle
@mkCppMethod IVRInput GetInputSourceHandle
@mkCppMethod IVRInput UpdateActionState
@mkCppMethod IVRInput GetDigitalActionData
@mkCppMethod IVRInput GetAnalogActionData
# @mkCppMethod IVRInput GetPoseActionData
@mkCppMethod IVRInput GetSkeletalActionData
@mkCppMethod IVRInput GetBoneCount
@mkCppMethod IVRInput GetBoneHierarchy
@mkCppMethod IVRInput GetBoneName
@mkCppMethod IVRInput GetSkeletalReferenceTransforms
@mkCppMethod IVRInput GetSkeletalTrackingLevel
@mkCppMethod IVRInput GetSkeletalBoneData
@mkCppMethod IVRInput GetSkeletalSummaryData
@mkCppMethod IVRInput GetSkeletalBoneDataCompressed
@mkCppMethod IVRInput DecompressSkeletalBoneData
@mkCppMethod IVRInput TriggerHapticVibrationAction
@mkCppMethod IVRInput GetActionOrigins
@mkCppMethod IVRInput GetOriginLocalizedName
@mkCppMethod IVRInput GetOriginTrackedDeviceInfo
@mkCppMethod IVRInput ShowActionOrigins
@mkCppMethod IVRInput ShowBindingsForActionSet

@mkCppMethod IVRIOBuffer Open
@mkCppMethod IVRIOBuffer Close
@mkCppMethod IVRIOBuffer Read
@mkCppMethod IVRIOBuffer Write
@mkCppMethod IVRIOBuffer PropertyContainer
# @mkCppMethod IVRIOBuffer HasReaders

@mkCppMethod IVRSpatialAnchors CreateSpatialAnchorFromDescriptor
@mkCppMethod IVRSpatialAnchors CreateSpatialAnchorFromPose
@mkCppMethod IVRSpatialAnchors GetSpatialAnchorPose
@mkCppMethod IVRSpatialAnchors GetSpatialAnchorDescriptor

@mkCppMethodInline COpenVRContext Clear
@mkCppMethod COpenVRContext CheckClear
@mkCppMethod COpenVRContext VRSystem
@mkCppMethod COpenVRContext VRChaperone
@mkCppMethod COpenVRContext VRChaperoneSetup
@mkCppMethod COpenVRContext VRCompositor
@mkCppMethod COpenVRContext VROverlay
@mkCppMethod COpenVRContext VRResources
@mkCppMethod COpenVRContext VRScreenshots
@mkCppMethod COpenVRContext VRRenderModels
@mkCppMethod COpenVRContext VRExtendedDisplay
@mkCppMethod COpenVRContext VRSettings
@mkCppMethod COpenVRContext VRApplications
@mkCppMethod COpenVRContext VRTrackedCamera
@mkCppMethod COpenVRContext VRDriverManager
@mkCppMethod COpenVRContext VRInput
@mkCppMethod COpenVRContext VRIOBuffer
@mkCppMethod COpenVRContext VRSpatialAnchors
# @mkCppMethod COpenVRContext VRNotifications


text:
"""
    Pump(pump_efficiency=0.7, cost_per_kW=0.28, derating=0.75)

An on-farm pump to access water.
"""
@with_kw mutable struct Pump <: Infrastructure

    @infrastructure_fields

    pump_efficiency::Union{Float64, AgParameter} = 0.7  # Efficiency of pump. Defaults to 0.7 (70%)
    cost_per_kW::Union{Float64, AgParameter} = 0.28  # cost in dollars/kW. Defaults to 0.28/kW

    # Accounts for efficiency losses between the energy required at the pump
    # shaft and the total energy required. Defaults to 0.75
    derating::Union{Float64, AgParameter} = 0.75
end

function cost_per_ha(p::Pump, year_step::Int64, area::Float64)::Float64
    return maintenance_cost(p, year_step) / area
end

function subtotal_costs(p::Pump, year_step::Int64)::Float64
    return cost_per_ha(p, year_step, 1.0)
end

"""
    pumping_costs_per_ML(p::Pump, flow_rate_Lps::Float64, head_pressure::Float64)::Float64

Calculate pumping cost per ML for a given flow rate and head pressure.

``P(Kw) = (H * Q) / ((102 * Ep) * D)``

where:
* `H` is head pressure in metres (m)
* `Q` is Flow in Litres per Second
* `Ep` is Pump efficiency (defaults to 0.7)
* `D` is the derating factor
* `102` is a constant, as given in Velloti & Kalogernis (2013)

See:
    * [Robinson, D. W. (2002)](http://www.clw.csiro.au/publications/technical2002/tr20-02.pdf)
    * [Vic. Dept. Agriculture (2006)](http://agriculture.vic.gov.au/agriculture/farm-management/soil-and-water/irrigation/border-check-irrigation-design>)
    * [Vellotti & Kalogernis (2013)](http://irrigation.org.au/wp-content/uploads/2013/06/Gennaro-Velloti-and-Kosi-Kalogernis-presentation.pdf>)

Parameters
----------
`flow_rate_Lps` : required flow rate in Litres per second over the irrigation duration

`head_pressure` : Head pressure of pumping system in metres.

Parameters taken from `Pump`
----------------------------
`pump_efficiency` : Efficiency of pump. Defaults to 0.7 (70%)

`derating` : Accounts for efficiency losses between the energy required at the pump
                    shaft and the total energy required. Defaults to 0.75

Returns
-------
float, cost_per_ML

- link to [`Pump(pump_efficiency, cost_per_kW, derating)`](@ref)
"""
function pumping_costs_per_ML(p::Pump, flow_rate_Lps::Float64,
                              head_pressure::Float64)::Float64
    if flow_rate_Lps <= 0.0
        return 0.0
    end

    c::Float64 = 102.0
    pe::Float64 = p.pump_efficiency
    dr::Float64 = p.derating
    energy_required_kW::Float64 = (head_pressure * flow_rate_Lps) / ((c * pe) * dr)

    # Litres / minutes in hour / seconds in minute
    hours_per_ML::Float64 = (1e6 / flow_rate_Lps) / 60.0^2

    cost_per_hour::Float64 = p.cost_per_kW * energy_required_kW
    cost_per_ML::Float64 = (cost_per_hour * hours_per_ML)

    @assert cost_per_ML >= 0.0 "Pumping costs cannot be negative.
    $cost_per_ML
    $flow_rate_Lps
    $head_pressure"

    return cost_per_ML
end

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl


text:
function abs(obj::SparseGradient)
    return jcall(obj, "abs", SparseGradient, ())
end

function acos(obj::SparseGradient)
    return jcall(obj, "acos", SparseGradient, ())
end

function acosh(obj::SparseGradient)
    return jcall(obj, "acosh", SparseGradient, ())
end

function add(obj::SparseGradient, arg0::SparseGradient)
    return jcall(obj, "add", SparseGradient, (SparseGradient,), arg0)
end

function add(obj::SparseGradient, arg0::jdouble)
    return jcall(obj, "add", SparseGradient, (jdouble,), arg0)
end

function add_in_place(obj::SparseGradient, arg0::SparseGradient)
    return jcall(obj, "addInPlace", void, (SparseGradient,), arg0)
end

function asin(obj::SparseGradient)
    return jcall(obj, "asin", SparseGradient, ())
end

function asinh(obj::SparseGradient)
    return jcall(obj, "asinh", SparseGradient, ())
end

function atan(obj::SparseGradient)
    return jcall(obj, "atan", SparseGradient, ())
end

function atan2(::Type{SparseGradient}, arg0::SparseGradient, arg1::SparseGradient)
    return jcall(SparseGradient, "atan2", SparseGradient, (SparseGradient, SparseGradient), arg0, arg1)
end

function atan2(obj::SparseGradient, arg0::SparseGradient)
    return jcall(obj, "atan2", SparseGradient, (SparseGradient,), arg0)
end

function atanh(obj::SparseGradient)
    return jcall(obj, "atanh", SparseGradient, ())
end

function cbrt(obj::SparseGradient)
    return jcall(obj, "cbrt", SparseGradient, ())
end

function ceil(obj::SparseGradient)
    return jcall(obj, "ceil", SparseGradient, ())
end

function compose(obj::SparseGradient, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "compose", SparseGradient, (jdouble, jdouble), arg0, arg1)
end

function copy_sign(obj::SparseGradient, arg0::SparseGradient)
    return jcall(obj, "copySign", SparseGradient, (SparseGradient,), arg0)
end

function copy_sign(obj::SparseGradient, arg0::jdouble)
    return jcall(obj, "copySign", SparseGradient, (jdouble,), arg0)
end

function cos(obj::SparseGradient)
    return jcall(obj, "cos", SparseGradient, ())
end

function cosh(obj::SparseGradient)
    return jcall(obj, "cosh", SparseGradient, ())
end

function create_constant(::Type{SparseGradient}, arg0::jdouble)
    return jcall(SparseGradient, "createConstant", SparseGradient, (jdouble,), arg0)
end

function create_variable(::Type{SparseGradient}, arg0::jint, arg1::jdouble)
    return jcall(SparseGradient, "createVariable", SparseGradient, (jint, jdouble), arg0, arg1)
end

function divide(obj::SparseGradient, arg0::SparseGradient)
    return jcall(obj, "divide", SparseGradient, (SparseGradient,), arg0)
end

function divide(obj::SparseGradient, arg0::jdouble)
    return jcall(obj, "divide", SparseGradient, (jdouble,), arg0)
end

function equals(obj::SparseGradient, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function exp(obj::SparseGradient)
    return jcall(obj, "exp", SparseGradient, ())
end

function expm1(obj::SparseGradient)
    return jcall(obj, "expm1", SparseGradient, ())
end

function floor(obj::SparseGradient)
    return jcall(obj, "floor", SparseGradient, ())
end

function get_derivative(obj::SparseGradient, arg0::jint)
    return jcall(obj, "getDerivative", jdouble, (jint,), arg0)
end

function get_field(obj::SparseGradient)
    return jcall(obj, "getField", Field, ())
end

function get_real(obj::SparseGradient)
    return jcall(obj, "getReal", jdouble, ())
end

function get_value(obj::SparseGradient)
    return jcall(obj, "getValue", jdouble, ())
end

function hash_code(obj::SparseGradient)
    return jcall(obj, "hashCode", jint, ())
end

function hypot(::Type{SparseGradient}, arg0::SparseGradient, arg1::SparseGradient)
    return jcall(SparseGradient, "hypot", SparseGradient, (SparseGradient, SparseGradient), arg0, arg1)
end

function hypot(obj::SparseGradient, arg0::SparseGradient)
    return jcall(obj, "hypot", SparseGradient, (SparseGradient,), arg0)
end

function linear_combination(obj::SparseGradient, arg0::SparseGradient, arg1::SparseGradient, arg2::SparseGradient, arg3::SparseGradient)
    return jcall(obj, "linearCombination", SparseGradient, (SparseGradient, SparseGradient, SparseGradient, SparseGradient), arg0, arg1, arg2, arg3)
end

function linear_combination(obj::SparseGradient, arg0::SparseGradient, arg1::SparseGradient, arg2::SparseGradient, arg3::SparseGradient, arg4::SparseGradient, arg5::SparseGradient)
    return jcall(obj, "linearCombination", SparseGradient, (SparseGradient, SparseGradient, SparseGradient, SparseGradient, SparseGradient, SparseGradient), arg0, arg1, arg2, arg3, arg4, arg5)
end

function linear_combination(obj::SparseGradient, arg0::SparseGradient, arg1::SparseGradient, arg2::SparseGradient, arg3::SparseGradient, arg4::SparseGradient, arg5::SparseGradient, arg6::SparseGradient, arg7::SparseGradient)
    return jcall(obj, "linearCombination", SparseGradient, (SparseGradient, SparseGradient, SparseGradient, SparseGradient, SparseGradient, SparseGradient, SparseGradient, SparseGradient), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function linear_combination(obj::SparseGradient, arg0::Vector{SparseGradient}, arg1::Vector{SparseGradient})
    return jcall(obj, "linearCombination", SparseGradient, (Vector{SparseGradient}, Vector{SparseGradient}), arg0, arg1)
end

function linear_combination(obj::SparseGradient, arg0::Vector{jdouble}, arg1::Vector{SparseGradient})
    return jcall(obj, "linearCombination", SparseGradient, (Vector{jdouble}, Vector{SparseGradient}), arg0, arg1)
end

function linear_combination(obj::SparseGradient, arg0::jdouble, arg1::SparseGradient, arg2::jdouble, arg3::SparseGradient)
    return jcall(obj, "linearCombination", SparseGradient, (jdouble, SparseGradient, jdouble, SparseGradient), arg0, arg1, arg2, arg3)
end

function linear_combination(obj::SparseGradient, arg0::jdouble, arg1::SparseGradient, arg2::jdouble, arg3::SparseGradient, arg4::jdouble, arg5::SparseGradient)
    return jcall(obj, "linearCombination", SparseGradient, (jdouble, SparseGradient, jdouble, SparseGradient, jdouble, SparseGradient), arg0, arg1, arg2, arg3, arg4, arg5)
end

function linear_combination(obj::SparseGradient, arg0::jdouble, arg1::SparseGradient, arg2::jdouble, arg3::SparseGradient, arg4::jdouble, arg5::SparseGradient, arg6::jdouble, arg7::SparseGradient)
    return jcall(obj, "linearCombination", SparseGradient, (jdouble, SparseGradient, jdouble, SparseGradient, jdouble, SparseGradient, jdouble, SparseGradient), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function log(obj::SparseGradient)
    return jcall(obj, "log", SparseGradient, ())
end

function log10(obj::SparseGradient)
    return jcall(obj, "log10", SparseGradient, ())
end

function log1p(obj::SparseGradient)
    return jcall(obj, "log1p", SparseGradient, ())
end

function multiply(obj::SparseGradient, arg0::SparseGradient)
    return jcall(obj, "multiply", SparseGradient, (SparseGradient,), arg0)
end

function multiply(obj::SparseGradient, arg0::jdouble)
    return jcall(obj, "multiply", SparseGradient, (jdouble,), arg0)
end

function multiply(obj::SparseGradient, arg0::jint)
    return jcall(obj, "multiply", SparseGradient, (jint,), arg0)
end

function multiply_in_place(obj::SparseGradient, arg0::SparseGradient)
    return jcall(obj, "multiplyInPlace", void, (SparseGradient,), arg0)
end

function negate(obj::SparseGradient)
    return jcall(obj, "negate", SparseGradient, ())
end

function new_instance(obj::SparseGradient, arg0::jdouble)
    return jcall(obj, "newInstance", SparseGradient, (jdouble,), arg0)
end

function num_vars(obj::SparseGradient)
    return jcall(obj, "numVars", jint, ())
end

function pow(::Type{SparseGradient}, arg0::jdouble, arg1::SparseGradient)
    return jcall(SparseGradient, "pow", SparseGradient, (jdouble, SparseGradient), arg0, arg1)
end

function pow(obj::SparseGradient, arg0::SparseGradient)
    return jcall(obj, "pow", SparseGradient, (SparseGradient,), arg0)
end

function pow(obj::SparseGradient, arg0::jdouble)
    return jcall(obj, "pow", SparseGradient, (jdouble,), arg0)
end

function pow(obj::SparseGradient, arg0::jint)
    return jcall(obj, "pow", SparseGradient, (jint,), arg0)
end

function reciprocal(obj::SparseGradient)
    return jcall(obj, "reciprocal", SparseGradient, ())
end

function remainder(obj::SparseGradient, arg0::SparseGradient)
    return jcall(obj, "remainder", SparseGradient, (SparseGradient,), arg0)
end

function remainder(obj::SparseGradient, arg0::jdouble)
    return jcall(obj, "remainder", SparseGradient, (jdouble,), arg0)
end

function rint(obj::SparseGradient)
    return jcall(obj, "rint", SparseGradient, ())
end

function root_n(obj::SparseGradient, arg0::jint)
    return jcall(obj, "rootN", SparseGradient, (jint,), arg0)
end

function scalb(obj::SparseGradient, arg0::jint)
    return jcall(obj, "scalb", SparseGradient, (jint,), arg0)
end

function signum(obj::SparseGradient)
    return jcall(obj, "signum", SparseGradient, ())
end

function sin(obj::SparseGradient)
    return jcall(obj, "sin", SparseGradient, ())
end

function sin_cos(obj::SparseGradient)
    return jcall(obj, "sinCos", FieldSinCos, ())
end

function sinh(obj::SparseGradient)
    return jcall(obj, "sinh", SparseGradient, ())
end

function sqrt(obj::SparseGradient)
    return jcall(obj, "sqrt", SparseGradient, ())
end

function subtract(obj::SparseGradient, arg0::SparseGradient)
    return jcall(obj, "subtract", SparseGradient, (SparseGradient,), arg0)
end

function subtract(obj::SparseGradient, arg0::jdouble)
    return jcall(obj, "subtract", SparseGradient, (jdouble,), arg0)
end

function tan(obj::SparseGradient)
    return jcall(obj, "tan", SparseGradient, ())
end

function tanh(obj::SparseGradient)
    return jcall(obj, "tanh", SparseGradient, ())
end

function taylor(obj::SparseGradient, arg0::Vector{jdouble})
    return jcall(obj, "taylor", jdouble, (Vector{jdouble},), arg0)
end

function to_degrees(obj::SparseGradient)
    return jcall(obj, "toDegrees", SparseGradient, ())
end

function to_radians(obj::SparseGradient)
    return jcall(obj, "toRadians", SparseGradient, ())
end


text:
#! /usr/bin/julia

# Rosetta Code, Zig-zag matrix

immutable ZigZag
    m::Int
    n::Int
    diag::Array{Int,1}
    cmax::Int
    numd::Int
    lohi::(Int,Int)
end

function zigzag(m::Int, n::Int)
    0<m && 0<n || error("The matrix dimensions must be positive.")
    ZigZag(m, n, [-1,1], m*n, m+n-1, extrema([m,n]))
end
zigzag(n::Int) = zigzag(n, n)

type ZZState
    cnt::Int
    cell::Array{Int,1}
    dir::Int
    dnum::Int
    dlen::Int
    dcnt::Int
end

Base.length(zz::ZigZag) = zz.cmax
Base.start(zz::ZigZag) = ZZState(1, [1,1], 1, 1, 1, 1)
Base.done(zz::ZigZag, zzs::ZZState) = zzs.cnt > zz.cmax

function Base.next(zz::ZigZag, zzs::ZZState)
    s = sub2ind((zz.m, zz.n), zzs.cell[1], zzs.cell[2])
    if zzs.dcnt == zzs.dlen
        if isodd(zzs.dnum)
            if zzs.cell[2] < zz.n
                zzs.cell[2] += 1
            else
                zzs.cell[1] += 1
            end
        else
            if zzs.cell[1] < zz.m
                zzs.cell[1] += 1
            else
                zzs.cell[2] += 1
            end
        end
        zzs.dcnt = 1
        zzs.dnum += 1
        zzs.dir = -zzs.dir
        if zzs.dnum <= zz.lohi[1]
            zzs.dlen += 1
        elseif zz.lohi[2] < zzs.dnum
            zzs.dlen -= 1
        end
    else
        zzs.cell += zzs.dir*zz.diag
        zzs.dcnt += 1
    end
    zzs.cnt += 1
    return (s, zzs)
end

using Formatting

function width{T<:Integer}(n::T)
    w = ndigits(n)
    n < 0 || return w
    return w + 1
end

function pretty{T<:Integer}(a::Array{T,2}, indent::Int=4)
    lo, hi = extrema(a)
    w = max(width(lo), width(hi))
    id = " "^indent
    fe = FormatExpr(@sprintf(" {:%dd}", w))
    s = id
    nrow = size(a)[1]
    for i in 1:nrow
        for j in a[i,:]
            s *= format(fe, j)
        end
        i != nrow || continue
        s *= "\n"*id
    end
    return s
end

n = 5
println("The n = ", n, " zig-zag matrix:")
a = zeros(Int, (n, n))
for (i, s) in enumerate(zigzag(n))
    a[s] = i-1
end
println(pretty(a))

m = 3
println()
println("Generalize to a non-square matrix (", m, "x", n, "):")
a = zeros(Int, (m, n))
for (i, s) in enumerate(zigzag(m, n))
    a[s] = i-1
end
println(pretty(a))

p = primes(10^3)
n = 7
println()
println("An n = ", n, " prime spiral matrix:")
a = zeros(Int, (n, n))
for (i, s) in enumerate(zigzag(n))
    a[s] = p[i]
end
println(pretty(a))

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code