text:
export lsmr, lsmr!

using LinearAlgebra

"""
    lsmr(A, b; kwrags...) -> x, [history]

Same as [`lsmr!`](@ref), but allocates a solution vector `x` initialized with zeros.
"""
lsmr(A, b; kwargs...) = lsmr!(zerox(A, b), A, b; kwargs...)

"""
    lsmr!(x, A, b; kwargs...) -> x, [history]

Minimizes ``\\|Ax - b\\|^2 + \\|λx\\|^2`` in the Euclidean norm. If multiple solutions
exists the minimum norm solution is returned.

The method is based on the Golub-Kahan bidiagonalization process. It is
algebraically equivalent to applying MINRES to the normal equations
``(A^*A + λ^2I)x = A^*b``, but has better numerical properties,
especially if ``A`` is ill-conditioned.

# Arguments
- `x`: Initial guess, will be updated in-place;
- `A`: linear operator;
- `b`: right-hand side.

## Keywords

- `λ::Number = 0`: lambda.
- `atol::Number = 1e-6`, `btol::Number = 1e-6`: stopping tolerances. If both are
  1.0e-9 (say), the final residual norm should be accurate to about 9 digits.
  (The final `x` will usually have fewer correct digits,
  depending on `cond(A)` and the size of damp).
- `conlim::Number = 1e8`: stopping tolerance. `lsmr` terminates if an estimate
  of `cond(A)` exceeds conlim.  For compatible systems Ax = b,
  conlim could be as large as 1.0e+12 (say).  For least-squares
  problems, conlim should be less than 1.0e+8.
  Maximum precision can be obtained by setting
- `atol` = `btol` = `conlim` = zero, but the number of iterations
  may then be excessive.
- `maxiter::Int = maximum(size(A))`: maximum number of iterations.
- `log::Bool`: keep track of the residual norm in each iteration;
- `verbose::Bool`: print convergence information during the iterations.

# Return values

**if `log` is `false`**

- `x`: approximated solution.

**if `log` is `true`**

- `x`: approximated solution.
- `ch`: convergence history.

**ConvergenceHistory keys**

- `:atol` => `::Real`: atol stopping tolerance.
- `:btol` => `::Real`: btol stopping tolerance.
- `:ctol` => `::Real`: ctol stopping tolerance.
- `:anorm` => `::Real`: anorm.
- `:rnorm` => `::Real`: rnorm.
- `:cnorm` => `::Real`: cnorm.
- `:resnom` => `::Vector`: residual norm at each iteration.
"""
function lsmr!(x, A, b;
    maxiter::Int = maximum(size(A)),
    log::Bool=false, kwargs...
    )
    history = ConvergenceHistory(partial=!log)
    reserve!(history,[:anorm,:rnorm,:cnorm],maxiter)

    T = Adivtype(A, b)
    m, n = size(A, 1), size(A, 2)
    btmp = similar(b, T)
    copyto!(btmp, b)
    v, h, hbar = similar(x, T), similar(x, T), similar(x, T)
    lsmr_method!(history, x, A, btmp, v, h, hbar; maxiter=maxiter, kwargs...)
    log && shrink!(history)
    log ? (x, history) : x
end

#########################
# Method Implementation #
#########################

function lsmr_method!(log::ConvergenceHistory, x, A, b, v, h, hbar;
    atol::Number = 1e-6, btol::Number = 1e-6, conlim::Number = 1e8,
    maxiter::Int = maximum(size(A)), λ::Number = 0,
    verbose::Bool=false
    )
    verbose && @printf("=== lsmr ===\n%4s\t%7s\t\t%7s\t\t%7s\n","iter","anorm","cnorm","rnorm")

    # Sanity-checking
    m = size(A, 1)
    n = size(A, 2)
    length(x) == n || error("x has length $(length(x)) but should have length $n")
    length(v) == n || error("v has length $(length(v)) but should have length $n")
    length(h) == n || error("h has length $(length(h)) but should have length $n")
    length(hbar) == n || error("hbar has length $(length(hbar)) but should have length $n")
    length(b) == m || error("b has length $(length(b)) but should have length $m")

    T = Adivtype(A, b)
    Tr = real(T)
    normrs = Tr[]
    normArs = Tr[]
    conlim > 0 ? ctol = convert(Tr, inv(conlim)) : ctol = zero(Tr)
    # form the first vectors u and v (satisfy  β*u = b,  α*v = A'u)
    tmp_u = similar(b)
    tmp_v = similar(v)
    mul!(tmp_u, A, x)
    b .-= tmp_u
    u = b
    β = norm(u)
    u .*= inv(β)
    mul!(v, adjoint(A), u)
    α = norm(v)
    v .*= inv(α)

    log[:atol] = atol
    log[:btol] = btol
    log[:ctol] = ctol

    # Initialize variables for 1st iteration.
    ζbar = α * β
    αbar = α
    ρ = one(Tr)
    ρbar = one(Tr)
    cbar = one(Tr)
    sbar = zero(Tr)

    copyto!(h, v)
    fill!(hbar, zero(Tr))

    # Initialize variables for estimation of ||r||.
    βdd = β
    βd = zero(Tr)
    ρdold = one(Tr)
    τtildeold = zero(Tr)
    θtilde  = zero(Tr)
    ζ = zero(Tr)
    d = zero(Tr)

    # Initialize variables for estimation of ||A|| and cond(A).
    normA, condA, normx = -one(Tr), -one(Tr), -one(Tr)
    normA2 = abs2(α)
    maxrbar = zero(Tr)
    minrbar = 1e100

    # Items for use in stopping rules.
    normb = β
    istop = 0
    normr = β
    normAr = α * β
    iter = 0
    # Exit if b = 0 or A'b = 0.

    log.mvps=1
    log.mtvps=1
    if normAr != 0
        while iter < maxiter
            nextiter!(log,mvps=1)
            iter += 1
            mul!(tmp_u, A, v)
            u .= tmp_u .+ u .* -α
            β = norm(u)
            if β > 0
                log.mtvps+=1
                u .*= inv(β)
                mul!(tmp_v, adjoint(A), u)
                v .= tmp_v .+ v .* -β
                α = norm(v)
                v .*= inv(α)
            end

            # Construct rotation Qhat_{k,2k+1}.
            αhat = hypot(αbar, λ)
            chat = αbar / αhat
            shat = λ / αhat

            # Use a plane rotation (Q_i) to turn B_i to R_i.
            ρold = ρ
            ρ = hypot(αhat, β)
            c = αhat / ρ
            s = β / ρ
            θnew = s * α
            αbar = c * α

            # Use a plane rotation (Qbar_i) to turn R_i^T to R_i^bar.
            ρbarold = ρbar
            ζold = ζ
            θbar = sbar * ρ
            ρtemp = cbar * ρ
            ρbar = hypot(cbar * ρ, θnew)
            cbar = cbar * ρ / ρbar
            sbar = θnew / ρbar
            ζ = cbar * ζbar
            ζbar = - sbar * ζbar

            # Update h, h_hat, x.
            hbar .= hbar .* (-θbar * ρ / (ρold * ρbarold)) .+ h
            x .+= (ζ / (ρ * ρbar)) * hbar
            h .= h .* (-θnew / ρ) .+ v

            ##############################################################################
            ##
            ## Estimate of ||r||
            ##
            ##############################################################################

            # Apply rotation Qhat_{k,2k+1}.
            βacute = chat * βdd
            βcheck = - shat * βdd

            # Apply rotation Q_{k,k+1}.
            βhat = c * βacute
            βdd = - s * βacute

            # Apply rotation Qtilde_{k-1}.
            θtildeold = θtilde
            ρtildeold = hypot(ρdold, θbar)
            ctildeold = ρdold / ρtildeold
            stildeold = θbar / ρtildeold
            θtilde = stildeold * ρbar
            ρdold = ctildeold * ρbar
            βd = - stildeold * βd + ctildeold * βhat

            τtildeold = (ζold - θtildeold * τtildeold) / ρtildeold
            τd = (ζ - θtilde * τtildeold) / ρdold
            d += abs2(βcheck)
            normr = sqrt(d + abs2(βd - τd) + abs2(βdd))

            # Estimate ||A||.
            normA2 += abs2(β)
            normA  = sqrt(normA2)
            normA2 += abs2(α)

            # Estimate cond(A).
            maxrbar = max(maxrbar, ρbarold)
            if iter > 1
                minrbar = min(minrbar, ρbarold)
            end
            condA = max(maxrbar, ρtemp) / min(minrbar, ρtemp)

            ##############################################################################
            ##
            ## Test for convergence
            ##
            ##############################################################################

            # Compute norms for convergence testing.
            normAr  = abs(ζbar)
            normx = norm(x)

            # Now use these norms to estimate certain other quantities,
            # some of which will be small near a solution.
            test1 = normr / normb
            test2 = normAr / (normA * normr)
            test3 = inv(condA)
            push!(log, :cnorm, test3)
            push!(log, :anorm, test2)
            push!(log, :rnorm, test1)
            verbose && @printf("%3d\t%1.2e\t%1.2e\t%1.2e\n",iter,test2,test3,test1)

            t1 = test1 / (one(Tr) + normA * normx / normb)
            rtol = btol + atol * normA * normx / normb
            # The following tests guard against extremely small values of
            # atol, btol or ctol.  (The user may have set any or all of
            # the parameters atol, btol, conlim  to 0.)
            # The effect is equivalent to the normAl tests using
            # atol = eps,  btol = eps,  conlim = 1/eps.
            if iter >= maxiter istop = 7; break end
            if 1 + test3 <= 1 istop = 6; break end
            if 1 + test2 <= 1 istop = 5; break end
            if 1 + t1 <= 1 istop = 4; break end
            # Allow for tolerances set by the user.
            if test3 <= ctol istop = 3; break end
            if test2 <= atol istop = 2; break end
            if test1 <= rtol  istop = 1; break end
        end
    end
    verbose && @printf("\n")
    setconv(log, istop ∉ (3, 6, 7))
    x
end


text:
plt = Observable{Any}(plot_M0(seq;darkmode))
ui = dom"div"(plt)
map!(p->plot_M0(p;darkmode), plt, seq_obs)
content!(w, "div#content", ui)

ext: jl


text:
# Benchmark

using StaticArrays
using Libdl: dlopen, dlsym, dlext
find_grad_fn(path_basename,extension) = dlsym(
               dlopen(joinpath(@__DIR__,
                      "$path_basename.$extension")),
               :density_gradient)

function density_gradient_rust(wfn; nsamples=10)

   point = ones(3)
   ngtoH = size(wfn.nuexp, 3)
   ∇ρ = zeros(3)

   for i in 1:nsamples
      @ccall $RUST_gradient(
                 point      :: Ptr{Cdouble},
                 wfn.nmo    :: Clong,
                 wfn.natm   :: Clong,
                 wfn.nprims :: Clong,
                 wfn.maxgrp :: Clong,
                 ngtoH      :: Clong,
                 wfn.ngroup :: Ptr{Clong},
                 wfn.ityp   :: Ptr{Clong},
                 wfn.nzexp  :: Ptr{Clong},
                 wfn.nlm    :: Ptr{Clong},
                 wfn.nuexp  :: Ptr{Clong},
                 wfn.occ    :: Ptr{Cdouble},
                 wfn.oexp   :: Ptr{Cdouble},
                 wfn.coords :: Ptr{Cdouble},
                 wfn.rcutte :: Ptr{Cdouble},
                 wfn.coef   :: Ptr{Cdouble},
                 ∇ρ         :: Ptr{Cdouble}, # inout
             )::Cvoid
   end
end

function density_gradient_cpp(wfn; nsamples=10)

   point = ones(3)
   ngtoH = size(wfn.nuexp, 3)
   ∇ρ = zeros(3)

   for i in 1:nsamples
      @ccall $CPP_gradient(
                 point      :: Ptr{Cdouble},
                 wfn.nmo    :: Clong,
                 wfn.natm   :: Clong,
                 wfn.nprims :: Clong,
                 wfn.maxgrp :: Clong,
                 ngtoH      :: Clong,
                 wfn.ngroup :: Ptr{Clong},
                 wfn.ityp   :: Ptr{Clong},
                 wfn.nzexp  :: Ptr{Clong},
                 wfn.nlm    :: Ptr{Clong},
                 wfn.nuexp  :: Ptr{Clong},
                 wfn.occ    :: Ptr{Cdouble},
                 wfn.oexp   :: Ptr{Cdouble},
                 wfn.coords :: Ptr{Cdouble},
                 wfn.rcutte :: Ptr{Cdouble},
                 wfn.coef   :: Ptr{Cdouble},
                 ∇ρ         :: Ptr{Cdouble}, # inout
             )::Cvoid
   end
end

function density_gradient_fortran(wfn; nsamples=10)

   point = ones(3)
   ngtoH = size(wfn.nuexp, 3)
   ∇ρ = zeros(3)

   for i in 1:nsamples
      @ccall $FORTRAN_gradient(
                 point      :: Ptr{Cdouble},
                 wfn.nmo    :: Ref{Clong},
                 wfn.natm   :: Ref{Clong},
                 wfn.nprims :: Ref{Clong},
                 wfn.maxgrp :: Ref{Clong},
                 ngtoH      :: Ref{Clong},
                 wfn.ngroup :: Ptr{Clong},
                 wfn.ityp   :: Ptr{Clong},
                 wfn.nzexp  :: Ptr{Clong},
                 wfn.nlm    :: Ptr{Clong},
                 wfn.nuexp  :: Ptr{Clong},
                 wfn.occ    :: Ptr{Cdouble},
                 wfn.oexp   :: Ptr{Cdouble},
                 wfn.coords :: Ptr{Cdouble},
                 wfn.rcutte :: Ptr{Cdouble},
                 wfn.coef   :: Ptr{Cdouble},
                 ∇ρ         :: Ptr{Cdouble}, # inout
             )::Cvoid
   end
end

# Pure Julia
function density_gradient_julia(wfn; nsamples=10)

   point = ones(3)
   ngtoH = size(wfn.nuexp, 3)
   ∇ρ = zeros(3)

   for i in 1:nsamples
      density_gradient!(∇ρ, point, wfn)
   end
   nothing
end

# ------------------------------------------------
#
# Prepare

cd("fortran")
run(`make`)
cd("../cpp")
run(`make`)
cd("../rust")
run(`make`)
cd("..")

const RUST_gradient    = find_grad_fn("rust/target/release/libexp_chemint_rust", dlext)
const CPP_gradient     = find_grad_fn("cpp/lib", "so")
const FORTRAN_gradient = find_grad_fn("fortran/libeval_gradrho", "so")

struct AIMWFN_Mol_Gamess_HF
   nmo       :: Int64
   nocc      :: Int64
   nprims    :: Int64
   natm      :: Int64
   atnam     :: Vector{String}
   coords    :: Matrix{Float64}
   charges   :: Vector{Float64}
   icen      :: Vector{Int64}
   ityp      :: Vector{Int64}
   oexp      :: Vector{Float64}
   occ       :: Vector{Float64}
   eorb      :: Vector{Float64}
   coef      :: Matrix{Float64}
   tote      :: Float64
   virial    :: Float64
   npc       :: Vector{Int64}
   icenat    :: Matrix{Int64}
   nlm       :: Matrix{Int64}
   ngroup    :: Vector{Int64}
   nzexp     :: Matrix{Int64}
   nuexp     :: Array{Int64,3}
   npcant    :: Int64
   maxgrp    :: Int64
   numshells :: Int64
   cuttz     :: Float64
   rcutte    :: Matrix{Float64}
   chkfile   :: String
end
include("julia/eval_gradrho.jl")

using Serialization
using BenchmarkTools
#using ChemInt

@doc """
    bench_library(lang=:f)

Benchmark C++ (:cpp), Rust (:rs), Fortran(:f), or Julia (:jl)
"""
function bench_library(lang=:f,system=:ch4)

   wfnpath = joinpath(@__DIR__, "$(string(system)).wfn")
   wfn = if :ChemInt in names(Main, imported = true)
            parse_wavefunction(wfnpath)
         else
            deserialize(wfnpath * ".ser")
         end

   formula = if system == :ch4
            "CH₄"
         elseif system == :c2h4
            "C₂H₄"
         elseif system == :imidazol
            "imidazol"
         end

   println()
   println("Calculating the density gradient of $formula at a point (x10) [$(string(lang))]:")
   if lang == :cpp
      return @btime density_gradient_cpp($wfn)
   elseif lang == :rs
      return @btime density_gradient_rust($wfn)
   elseif lang == :f
      return @btime density_gradient_fortran($wfn)
   elseif lang == :jl
      return @btime density_gradient_julia($wfn)
   else
      @warn "lang must be one of :cpp, :rs, :f, or :jl"
   end
end

for system in [:ch4, :c2h4, :imidazol]
   for lang in [:f,:rs,:cpp,:jl]
      bench_library(lang, system)
   end
   println()
end


text:
using CautiousEngine

@testset "dfa function tests" begin
    # Simple reach-avoid Specification
    dfa_states = collect(1:3)
    dfa_props = ["a", "b"]
    dfa_transitions = [(1, "!a∧!b", 1),
                        (1, "a∧!b", 2),
                        (2, "true", 2),
                        (1, "!a∧b", 3),
                        (3, "true", 3)]
    dfa_accepting_state = 2
    dfa_sink_state = 3
    dfa_initial_state = 1

    dfa = CautiousEngine.DFA(dfa_states, dfa_props, dfa_transitions, dfa_accepting_state, dfa_sink_state, dfa_initial_state)

    @test CautiousEngine.get_next_dfa_labels(dfa, 1) == ["a∧!b"]
    @test CautiousEngine.get_next_dfa_labels(dfa, 2) == []
    @test CautiousEngine.get_next_dfa_labels(dfa, 3) == []

end


text:
using EmpiricalRiskMinimization
srand(123)

# data
n = 2000; k = 30;
d = k + 1;
X = randn(n, k); theta = randn(d);
y = [ones(n) X] * theta + 0.5 * randn(n);

# model
M = Model(X, y, embedall=true)
train(M);
status(M);
theta = thetaopt(M);


text:
const ParamCacheKey = NamedTuple{(:problem, :type, :name), NTuple{3, Symbol}}

make_cache_key(problem, type, name) = (problem = problem, type = type, name = name)

# Priority for keeping data in cache to serve reads. Currently unused.
IS.@scoped_enum(CachePriority, LOW = 1, MEDIUM = 2, HIGH = 3,)

struct CacheFlushRule
    keep_in_cache::Bool
    priority::CachePriority
end

CacheFlushRule() = CacheFlushRule(false, CachePriority.LOW)

"""
Informs the flusher on what data to keep in cache.
"""
struct CacheFlushRules
    data::Dict{ParamCacheKey, CacheFlushRule}
    min_flush_size::Int
    max_size::Int
end

const MIN_CACHE_FLUSH_SIZE_MiB = MiB

function CacheFlushRules(; max_size = GiB, min_flush_size = MIN_CACHE_FLUSH_SIZE_MiB)
    return CacheFlushRules(Dict{ParamCacheKey, CacheFlushRule}(), min_flush_size, max_size)
end

function add_rule!(rules::CacheFlushRules, problem, type, name, keep_in_cache, priority)
    key = make_cache_key(problem, type, name)
    rules.data[key] = CacheFlushRule(keep_in_cache, priority)
end

get_rule(x::CacheFlushRules, problem, type, name) =
    get_rule(x, make_cache_key(problem, type, name))
get_rule(x::CacheFlushRules, key) = x.data[key]

mutable struct CacheStats
    hits::Int
    misses::Int
end

CacheStats() = CacheStats(0, 0)

function get_cache_hit_percentage(x::CacheStats)
    total = x.hits + x.misses
    total == 0 && return 0.0

    return x.hits / (x.hits + x.misses) * 100
end


text:
using OOPMacro
import OOPMacro: findFnCall, findFnSelfArgNameSymbol, setFnSelf!, deleteFnSelf!, setFnName!, getFnName, setFnSelfArgType!

testFiles = ("fnUtil.jl",
"basic.jl",
"constructor.jl",
"inheritence.jl",
"multi-module.jl")

@testset "Testing: $filename" for filename in testFiles
    include(filename)
end


text:
using DeformableRegistration: Distance, ImageProcessing, Transformation, Interpolation, regOptions
 using Base.Test
 using MicroLogging
 using Interpolations

include("helpers/checkDerivative.jl")

@testset "Distance: distance derivatives - SSD" begin
    # create reference image
    refImg = createImage(100*rand(50,60))

    # measure distance and check derivative (nonparametric)
    centeredGrid = getCellCenteredGrid(refImg)
     centeredGrid.data[:] = centeredGrid.data[:] + 0.3 *rand(size(centeredGrid.data))
     D,dD,d2D = ssdDistance(refImg,refImg,centeredGrid.data,doDerivative=true,doHessian=true)
     Dfunc(x) = ssdDistance(refImg,refImg,x)[1]
     errlin,errquad = checkDerivative(Dfunc,dD',centeredGrid.data,doPlot=false)
     @test checkErrorDecay(errquad)

    options = regOptions()
     options.matrixFree = true;
     D,dD,d2D = ssdDistance(refImg,refImg,centeredGrid.data,doDerivative=true,doHessian=true,options=options)
     Dfunc(x) = ssdDistance(refImg,refImg,x)[1]
     errlin,errquad = checkDerivative(Dfunc,dD',centeredGrid.data)
     @test checkErrorDecay(errquad)

    # measure distance and check derivative (parametic)
    centeredGrid = getCellCenteredGrid(refImg)
     options = DeformableRegistration.regOptions()
     options.parametricOnly=true
     evaluationPoint = [1,0,0,0,1,0]+0.1*rand(6)
     D,dD,d2D = ssdDistance(refImg,refImg,transformGridAffine(centeredGrid,evaluationPoint).data,doDerivative=true,doHessian=true,options=options)
     Dfunc(p) = ssdDistance(refImg,refImg,transformGridAffine(centeredGrid,p).data,options=options)[1]
     errlin,errquad = checkDerivative(Dfunc,dD',evaluationPoint)
     @test checkErrorDecay(errquad)

    options.matrixFree = true;
     D,dD,d2D = ssdDistance(refImg,refImg,transformGridAffine(centeredGrid,evaluationPoint).data,doDerivative=true,doHessian=true,options=options)
     Dfunc(p) = ssdDistance(refImg,refImg,transformGridAffine(centeredGrid,p).data,options=options)[1]
     errlin,errquad = checkDerivative(Dfunc,dD',evaluationPoint)
     @test checkErrorDecay(errquad)
end

@testset "NGF" begin
    @testset "Distance: distance derivatives - NGF - parametric" begin
        refImg = createImage(100*rand(50,60))

        # measure distance and check derivative (nonparametric)
        centeredGrid = getCellCenteredGrid(refImg)
        centeredGrid.data[:] = centeredGrid.data[:] + 0.3 *rand(size(centeredGrid.data))
        D,dD,d2D = ngfDistance(refImg,refImg,centeredGrid.data,doDerivative=true,doHessian=true)
        Dfunc(x) = ngfDistance(refImg,refImg,x)[1]
        errlin,errquad = checkDerivative(Dfunc,dD',centeredGrid.data,doPlot=false)
        @test checkErrorDecay(errquad)

        evaluationPoint = [1,0,0,0,1,0.0]+0.1*rand(6)
        centeredGrid = getCellCenteredGrid(refImg)
        options = DeformableRegistration.regOptions()
        options.edgeParameterR = 0.1
        options.edgeParameterT = 0.1
        options.matrixFree = false
        options.parametricOnly = true
        D,dD,d2D = ngfDistance(refImg,refImg,transformGridAffine(centeredGrid,evaluationPoint).data,doDerivative=true,doHessian=true,options=options)
        Dfunc(p) = ngfDistance(refImg,refImg,transformGridAffine(centeredGrid,p).data,options=options)[1]
        errlin,errquad = checkDerivative(Dfunc,dD',evaluationPoint,doPlot=false)
        @test checkErrorDecay(errquad)

        evaluationPoint = [1,0,0,0,1,0.0]+0.1*rand(6)
        centeredGrid = getCellCenteredGrid(refImg)
        options = DeformableRegistration.regOptions()
        options.edgeParameterR = 0.1
        options.edgeParameterT = 0.1
        options.matrixFree = true
        options.parametricOnly = true
        D,dD,d2D = ngfDistance(refImg,refImg,transformGridAffine(centeredGrid,evaluationPoint).data,doDerivative=true,doHessian=true,options=options)
        Dfunc(p) = ngfDistance(refImg,refImg,transformGridAffine(centeredGrid,p).data,options=options)[1]
        errlin,errquad = checkDerivative(Dfunc,dD',evaluationPoint,doPlot=false)
        @test checkErrorDecay(errquad)
    end

    @testset "Distance: distance derivatives - NGF - nonparametric" begin
        refImg = createImage(100*rand(50,60))

        ##
        centeredGrid = getCellCenteredGrid(refImg)
        centeredGrid.data[:] = centeredGrid.data[:] + 0.3 *rand(size(centeredGrid.data))
        options = DeformableRegistration.regOptions()
        options.edgeParameterR = 0.1
        options.edgeParameterT = 0.1
        options.matrixFree = false
        options.parametricOnly = false
        D,dD,d2D = ngfDistance(refImg,refImg,centeredGrid.data,doDerivative=true,doHessian=true,options=options)
        Dfunc(grid) = ngfDistance(refImg,refImg,grid,options=options)[1]
        errlin,errquad = checkDerivative(Dfunc,dD',centeredGrid.data,doPlot=false)
        @test checkErrorDecay(errquad)

        centeredGrid = getCellCenteredGrid(refImg)
        centeredGrid.data[:] = centeredGrid.data[:] + 0.3 *rand(size(centeredGrid.data))
        options = DeformableRegistration.regOptions()
        options.edgeParameterR = 0.1
        options.edgeParameterT = 0.1
        options.matrixFree = true
        options.parametricOnly = false
        D,dD,d2D = ngfDistance(refImg,refImg,centeredGrid.data,doDerivative=true,doHessian=true,options=options)
        Dfunc(grid) = ngfDistance(refImg,refImg,grid,options=options)[1]
        errlin,errquad = checkDerivative(Dfunc,dD',centeredGrid.data,doPlot=false)
        @test checkErrorDecay(errquad)

    end

    @testset "Distance: ngf epsilon estimation" begin
        testimage = dirname(Base.source_path()) * "/testdata/luebeck.jpg"
        refImg = loadImage(testimage)
        for r = refImg.data
             r = r + 0.01*rand()
        end
        warn("Skipping NGF ϵ estimation test.")
        @test_skip estimateNGFEpsilon(refImg,cutoffPercent=80)[1] >= 0 # seems to run forever
    end
end

@testset "SSD arbitrary grid" begin
    @testset "fixedGrid_IdenticalImagesZeroDeformation" begin
        n = 10

        imageData = ones(n,n)
        referenceImage = createImage(imageData, voxelsize = [1.0,1.0], shift = [0.0, 0.0])
        templateImage = createImage(imageData, voxelsize = [1.0,1.0], shift = [0.0, 0.0])
        transformedGrid = getCellCenteredGrid(referenceImage)
        # voxelsize::Array{Float64,1},shift::Array{Float64,1}, gridSize::Tuple{Int64,Int64}
        functionValue,dFunctionValue,d2FunctionValue,d_transformedImage =
            ssdDistanceArbitraryGrid(referenceImage, templateImage, transformedGrid)
                             # doDerivative::Bool=false,doHessian::Bool=false,
                             # options::regOptions=regOptions(),centeredGrid::scaledArray=scaledArray(zeros(1),(0,),[0,],[0,])))
        @test functionValue ≈ 0 atol=0.0001

        functionValue,dFunctionValue,d2FunctionValue,d_transformedImage =
            ssdDistanceArbitraryGrid(referenceImage, templateImage, transformedGrid, doDerivative=true)
        @test functionValue ≈ 0 atol=0.0001

        imageData = rand(n,n)
        referenceImage = createImage(imageData, voxelsize = [1.0,1.0], shift = [0.0, 0.0])
        templateImage = createImage(imageData, voxelsize = [1.0,1.0], shift = [0.0, 0.0])
        transformedGrid = getCellCenteredGrid(referenceImage)
        # voxelsize::Array{Float64,1},shift::Array{Float64,1}, gridSize::Tuple{Int64,Int64}
        functionValue,dFunctionValue,d2FunctionValue,d_transformedImage =
            ssdDistanceArbitraryGrid(referenceImage, templateImage, transformedGrid)
        @test functionValue ≈ 0 atol=0.0001

        functionValue,dFunctionValue,d2FunctionValue,d_transformedImage =
            ssdDistanceArbitraryGrid(referenceImage, templateImage, transformedGrid, doDerivative=true)
        @test functionValue ≈ 0 atol=0.0001

        imageData = rand(n,n)
        referenceImage = createImage(imageData, voxelsize = [1.0,1.0], shift = [0.0, 0.0])
        templateImage = createImage(imageData, voxelsize = [1.0,1.0], shift = [0.0, 0.0])
        transformedGrid = getCellCenteredGrid(restrictResolutionToLevel(referenceImage,2))
        @debug "-------- transformedGrid"
        @debug transformedGrid
        # voxelsize::Array{Float64,1},shift::Array{Float64,1}, gridSize::Tuple{Int64,Int64}
        functionValue,dFunctionValue,d2FunctionValue,d_transformedImage =
            ssdDistanceArbitraryGrid(referenceImage, templateImage, transformedGrid)
        @test functionValue ≈ 0 atol=0.0001

    end

    @testset "fixedGrid_IdenticalImagesFormats" begin
        n = 10

        imageData = ones(n,n)
        referenceImage = createImage(imageData, voxelsize = [1.0,1.0], shift = [0.0, 0.0])
        templateImage = createImage(imageData, voxelsize = [1.0,1.0], shift = [0.0, 0.0])

        gridSize = (5,5)
        voxelsize = n/gridSize[1]
        transformedGrid = getCellCenteredGrid([voxelsize,voxelsize],[0.0,0.0],gridSize)
        # voxelsize::Array{Float64,1},shift::Array{Float64,1}, gridSize::Tuple{Int64,Int64}
        functionValue,dFunctionValue,d2FunctionValue,d_transformedImage =
            ssdDistanceArbitraryGrid(referenceImage, templateImage, transformedGrid, doDerivative=true)
                             # doDerivative::Bool=false,doHessian::Bool=false,
                             # options::regOptions=regOptions(),centeredGrid::scaledArray=scaledArray(zeros(1),(0,),[0,],[0,])))
        @test functionValue ≈ 0 atol=0.0001
        @test size(dFunctionValue)... == prod(gridSize)*2
    ##
        imageData = rand(n,n)
        referenceImage = createImage(imageData, voxelsize = [1.0,1.0], shift = [0.0, 0.0])
        templateImage = createImage(imageData, voxelsize = [1.0,1.0], shift = [0.0, 0.0])

    ##
        gridSize = (10,20)
        voxelsize = 1.0
        transformedGrid = getCellCenteredGrid([voxelsize,voxelsize],[0.0,0.0],gridSize)
        transformedGrid.data[1:200] = transformedGrid.data[1:200] + 4.0
        # voxelsize::Array{Float64,1},shift::Array{Float64,1}, gridSize::Tuple{Int64,Int64}
        functionValue,dFunctionValue,d2FunctionValue,d_transformedImage =
            ssdDistanceArbitraryGrid(referenceImage, templateImage, transformedGrid, doDerivative=true)
    ##
        @test norm(dFunctionValue[105:195]) ≈ 0
        @test norm(dFunctionValue[305:395]) ≈ 0
        #@test norm(dFunctionValue[5:95]) > 1
        @test size(dFunctionValue)... == prod(gridSize)*2

    end

    @testset "ssdResampling_derivatives" begin

    include("../test/helpers/checkDerivative.jl")
    #include("../src/distances/ssdDistance.jl")

    smoothrand(x::Array{Float64,1}) = smoothArray(rand(size(x)),3, 1.0)

    refImg = createImage( smoothArray(100.0*rand(99,111),3, 3.0) )
        #refImg.data.data = smoothArray(refImg.data.data,11, 1.0)
        centeredGrid = getCellCenteredGrid(refImg)
        centeredGrid.data[:] = centeredGrid.data[:] + 5.0 *smoothrand(centeredGrid.data)
        options = regOptions()
        options.matrixFree = true;
        D,dD,d2D = Distance.ssdDistanceArbitraryGrid(refImg,refImg,centeredGrid, doDerivative=true,
            doHessian=true,options=options, interpolationScheme=BSpline(Linear()))
        Dfunc(x) = Distance.ssdDistanceArbitraryGrid(refImg,refImg,x, doDerivative=false, interpolationScheme=BSpline(Linear()))[1]
        errlin,errquad = checkDerivative(Dfunc,dD',centeredGrid, doPlot=false)
        @test_skip checkErrorDecay(errquad)

    # The transposed interpolation operator is not implemented correctly which is
    # why in the case where the grid resolution changes, the derivative checks are failing.

    gridSize = (25,30)
        voxelsize = size(refImg.data)./gridSize
        vs = voxelsize
        centeredGrid = getCellCenteredGrid([vs[1]; vs[2]],[0.0,0.0],gridSize)
        centeredGrid.data[:] = centeredGrid.data[:] + 5 *smoothrand(centeredGrid.data)
        D,dD,d2D = ssdDistanceArbitraryGrid(refImg,refImg,centeredGrid,doDerivative=true,doHessian=true,options=options, interpolationScheme=BSpline(Linear()))
        Dfunc(x) = ssdDistanceArbitraryGrid(refImg,refImg,x, doDerivative=false, interpolationScheme=BSpline(Cubic(Line())))[1]
        errlin,errquad = checkDerivative(Dfunc,dD',centeredGrid, doPlot=false)
        @test_skip checkErrorDecay(errquad)

    gridSize = (10,10)
        voxelsize = size(refImg.data)./gridSize
        vs = voxelsize
        centeredGrid = getCellCenteredGrid([vs[1]; vs[2]],[0.0,0.0],gridSize)
        centeredGrid.data[:] = centeredGrid.data[:] + 10 * smoothrand(centeredGrid.data)
        D,dD,d2D = ssdDistanceArbitraryGrid(refImg,refImg,centeredGrid,doDerivative=true,doHessian=true,options=options)
        Dfunc(x) = ssdDistanceArbitraryGrid(refImg,refImg,x)[1]
        errlin,errquad = checkDerivative(Dfunc,dD',centeredGrid, doPlot=false)
        @test_skip checkErrorDecay(errquad)

    gridSize = (25,30)
        vs = size(refImg.data)./gridSize
        centeredGrid = getCellCenteredGrid([vs[1]; vs[2]],[0.0,0.0],gridSize)
        disp = 10 * ones(gridSize)
        disp[4:6,4:6]=0.0
        testGrid = 0.0 * centeredGrid + vcat(disp[:], disp[:])
        testGridInterp = interpolateDeformationField(testGrid, getCellCenteredGrid(refImg), interpolationScheme=BSpline(Linear()))
        testGridInterp2 = interpolateDeformationField(testGridInterp, getCellCenteredGrid([vs[1]; vs[2]],[0.0,0.0],gridSize), interpolationScheme=BSpline(Linear()))
         #
         # using PyPlot
         # figure()
         # subplot(1,3,1)
         # imshow(reshape(testGrid.data[1:Int(end/2)],testGrid.dimensions))
         # clim([0,10])
         # subplot(1,3,2)
         # imshow(reshape(testGridInterp.data[1:Int(end/2)],testGridInterp.dimensions))
         # clim([0,10])
         # subplot(1,3,3)
         # imshow(reshape(testGridInterp2.data[1:Int(end/2)],testGridInterp2.dimensions))
         # clim([0,10])

        @test_skip sum(testGrid.data[:]) ≈ sum(testGridInterp2.data[:]) atol=0.1

    end

end


text: get_core_id() = Int(@ccall sched_getcpu()::Cint)

max_stars_repo_name: carstenbauer/BandwidthBenchmark

max_issues_repo_name: carstenbauer/BandwidthBenchmark

max_forks_repo_name: carstenbauer/BandwidthBenchmark


text:
# Dispatch types to call different problem builders
abstract type AbstractExperiment end
abstract type AbstractDynamicExperiment <: AbstractExperiment end
abstract type AbstractStochasticExperiment <: AbstractExperiment end
abstract type AbstractSteadyStateExperiment <: AbstractExperiment end
struct DynamicLL <: AbstractDynamicExperiment end
struct DynamicNLSE <: AbstractDynamicExperiment end
struct DynamicIkeda <: AbstractDynamicExperiment end

struct StochasticNLSE <: AbstractStochasticExperiment end

struct SteadyStateLL <: AbstractSteadyStateExperiment end
struct SteadyStateNLSE <: AbstractSteadyStateExperiment end
struct SteadyStateIkeda <: AbstractSteadyStateExperiment end

# Abstraction layers over DiffEq solution/problem types
abstract type AbstractNLOProblem end
mutable struct DynamicNLSEProblem{fType, f0Type, fftType, ifftType, meshType, opType} <: AbstractNLOProblem
    prob::ODEProblem
    model::AbstractModel
    ω::fType
    ω0::f0Type
    tmesh::meshType
    planned_fft::fftType
    planned_ifft::ifftType
    Doperator::opType
end
mutable struct StochasticNLSEProblem{fType, f0Type, fftType, ifftType, meshType, opType} <: AbstractNLOProblem
    prob::SDEProblem
    model::AbstractModel
    ω::fType
    ω0::f0Type
    tmesh::meshType
    planned_fft::fftType
    planned_ifft::ifftType
    Doperator::opType
end
mutable struct DynamicLLProblem{fType, f0Type, fftType, ifftType, meshType, opType} <: AbstractNLOProblem
    prob::ODEProblem
    model::AbstractModel
    ω::fType
    ω0::f0Type
    tmesh::meshType
    planned_fft::fftType
    planned_ifft::ifftType
    Doperator::opType
end
mutable struct DynamicIkedaProblem{fType, f0Type, fftType, ifftType, meshType, opType, cbType} <: AbstractNLOProblem
    prob::ODEProblem
    model::AbstractModel
    ω::fType
    ω0::f0Type
    tmesh::meshType
    planned_fft::fftType
    planned_ifft::ifftType
    Doperator::opType
    ikeda_callback::cbType
end

abstract type AbstractNLOSolution end
mutable struct DynamicNLOSolution <: AbstractNLOSolution
    sol::DESolution
    prob::AbstractNLOProblem
end

# calling solution returns time-domain solution
function (sol::DynamicNLOSolution)(z, mode_idx::Int=1)
    structure_idx, ls = currentstructure(sol.prob.model, z)
    return sol.prob.planned_fft * (sol.sol(z/m)[:, mode_idx, structure_idx].*exp.(view(sol.prob.Doperator, :, mode_idx, structure_idx) * z))
end
function FT(sol::DynamicNLOSolution, z, mode_idx::Int=1)
    dt = sol.prob.tmesh[2]-sol.prob.tmesh[1]
    structure_idx, ls = currentstructure(sol.prob.model, z)
    return fftshift(view(sol.sol(z/m), :, mode_idx, structure_idx) .* exp.(view(sol.prob.Doperator, :, mode_idx, structure_idx) * z)) / dt
end

##Abstract interfacei forwards to DiffEq sol field
Base.getindex(A::AbstractNLOSolution,i::Int) = getindex(A.sol, i)
Base.getindex(A::AbstractNLOSolution,I::Vararg{Int, N}) where {N} = getindeX(A.sol, I)
Base.setindex!(A::AbstractNLOSolution, v, i::Int) = setindex!(A.sol, v, i)
Base.setindex!(A::AbstractNLOSolution, v, I::Vararg{Int, N}) where {N} = setindex!(A.sol, v, I)
size(A::AbstractNLOSolution) = size(A.sol)
Base.summary(A::AbstractNLOSolution) = summary(A.sol)
Base.show(io::IO, A::AbstractNLOSolution) = show(io, A.sol)
Base.show(io::IO, m::MIME"text/plain", A::AbstractNLOSolution) = show(io, m, A.sol)
tuples(sol::AbstractNLOSolution) = tuples(sol.sol)
start(sol::AbstractNLOSolution) = start(sol.sol)
next(sol::AbstractNLOSolution,state) = next(sol.sol, state)
done(sol::AbstractNLOSolution,state) = done(sol.sol, state)


text:
# this is probably a bad way of doing things but I want to test some utils
# that are not otherwise exposed to the user
include("../src/utils/utils.jl")

@testset "Monin-Obukhov length tests" begin
    @test monin_obukhov("A", 1.2) ≈ -11.4*1.2^0.10
    @test monin_obukhov("D", 1.2) ≈ Inf
    @test_throws ErrorException monin_obukhov("Q", 1.2)
end

@testset "Pasquill-Gifford dispersion tests" begin
    σy_plume = crosswind_dispersion("A")
    σy_puff = crosswind_dispersion("A"; plume=false)
    @test σy_plume(1.2) ≈ 0.423*((600/600)^0.2)*(1.2^0.9)
    @test σy_puff(1.2) ≈ 0.18*(1.2^0.92)
    @test_throws ErrorException crosswind_dispersion("Q")

    σz_plume = vertical_dispersion("A")
    σz_puff = vertical_dispersion("A"; plume=false)
    @test σz_plume(1.2) ≈ 107.7*(1.2^-1.7172)*exp(0.2770*log(1.2)^2)
    @test σz_puff(1.2) ≈ 0.60*(1.2^0.75)
    @test_throws ErrorException vertical_dispersion("Q")
end

@testset "Windspeed by powerlaw" begin
    u0, z0, p = 3.0, 1.0, 0.108
    u = windspeed(u0, z0, "A")
    @test u(10) ≈ u0*(10/z0)^p
    @test_throws ErrorException windspeed(u0, z0, "Q")
end

@testset "Windspeed by Monin-Obukhov length" begin
    ustar, zR, k = 3.0, 1.0, 0.35
    λ = monin_obukhov("A", zR)
    a(z) = (1-15*(z/λ))^0.25
    Ψ(a) = 2*log((1+a)/2) + log((1+a^2)/2) - 2*atan(a) + π/2

    uA = windspeed(ustar, zR, λ, "A")
    uD = windspeed(ustar, zR, λ, "D")
    uE = windspeed(ustar, zR, λ, "E")

    @test uA(10) ≈ (ustar/k)*(log((10+zR)/zR) - Ψ(a(10)))
    @test uD(10) ≈ (ustar/k)*(log((10+zR)/zR))
    @test uE(10) ≈ (ustar/k)*(log((10+zR)/zR) - 4.7*(10/λ))
    @test_throws ErrorException windspeed(ustar, zR, λ, "Q")
end


text:
# Julia wrapper for header: /usr/include/libavcodec/mediacodec.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_mediacodec_alloc_context,
    av_mediacodec_default_init,
    av_mediacodec_default_free,
    av_mediacodec_release_buffer,
    av_mediacodec_render_buffer_at_time

function av_mediacodec_alloc_context()
    ccall((:av_mediacodec_alloc_context, libavcodec), Ptr{AVMediaCodecContext}, ())
end

function av_mediacodec_default_init(avctx, ctx, surface)
    ccall((:av_mediacodec_default_init, libavcodec), Cint, (Ptr{AVCodecContext}, Ptr{AVMediaCodecContext}, Ptr{Cvoid}), avctx, ctx, surface)
end

function av_mediacodec_default_free(avctx)
    ccall((:av_mediacodec_default_free, libavcodec), Cvoid, (Ptr{AVCodecContext},), avctx)
end

function av_mediacodec_release_buffer(buffer, render::Integer)
    ccall((:av_mediacodec_release_buffer, libavcodec), Cint, (Ptr{AVMediaCodecBuffer}, Cint), buffer, render)
end

function av_mediacodec_render_buffer_at_time(buffer, time::Int64)
    ccall((:av_mediacodec_render_buffer_at_time, libavcodec), Cint, (Ptr{AVMediaCodecBuffer}, Int64), buffer, time)
end


text:
@timedtestset "sparse gp" begin
    x = 0:0.1:10
    xu = 0:10
    σ = 1.0
    σu = 1e-3
    f = wrap(GP(Matern32Kernel()), GPC())
    covariance_error = "The covariance matrix of a sparse GP can often be dense and " *
        "can cause the computer to run out of memory. If you are sure you have enough " *
        "memory, you can use `cov(f.fobs)`."

    @timedtestset "SparseFiniteGP Constructors" begin
        f = wrap(GP(Matern32Kernel()), GPC())
        @test SparseFiniteGP(f(x), f(xu)) == SparseFiniteGP(f(x, 1e-18), f(xu, 1e-18))
    end

    @timedtestset "SparseFiniteGP methods" begin
        f = wrap(GP(Matern32Kernel()), GPC())
        fx = f(x)
        fxu = SparseFiniteGP(f(x), f(xu))
        @test mean(fxu) == mean(fx)
        @test marginals(fxu) == marginals(fx)
        @test rand(MersenneTwister(12345), fxu) == rand(MersenneTwister(12345), fx)
        @test_throws ErrorException(covariance_error) cov(fxu)
        @test cov(fxu.fobs) == cov(fx)
    end

    @timedtestset "SparseFiniteGP inference" begin
        f = wrap(GP(Matern32Kernel()), GPC())
        fx = f(x, σ)
        fxu = SparseFiniteGP(f(x, σ), f(xu, σu))
        y = rand(MersenneTwister(12345), fxu)

        fpost1 = approx_posterior(VFE(), fxu.fobs, y, fxu.finducing)
        fpost2 = posterior(fxu, y)

        @test marginals(fpost1(x)) == marginals(fpost2(x))
        @test elbo(fxu, y) == logpdf(fxu, y)
        @test logpdf(fxu, y) == elbo(fxu.fobs, y, fxu.finducing)
        yy = rand(fxu, 10)
        @test all(logpdf(fx, yy) .> logpdf(fxu, yy))
    end
end


text:
# include("../src/LarSurf.jl")

using LarSurf
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm, SparseArrays
using Pio3d

threshold = 4000
pth = Pio3d.datasets_join_path("medical/orig/sample-data/nrn4.pklz")
datap = Pio3d.read3d(pth)

data3d = datap["data3d"]
# segmentation = convert(Array{Int8, 2}, data3d .> threshold)
segmentation = data3d .> threshold

coord_vect_of_segment = [Int8(segmentation[k]) for k=1:length(segmentation)]

# coord_vect_of_segment = [x>0.25 ? 1 : 0  for x in rand(size(∂_2,1)) ]

# V,(_,_,FV,CV) = Lar.cuboidGrid([32,32,16], true)
sz = size(data3d)
V,(_,_,FV,CV) = Lar.cuboidGrid([sz[1], sz[2], sz[3]], true)
∂_2 = Lar.u_coboundary_2( CV, FV)
coord_vect_of_all_3D_cells  = ones(size(∂_2,1),1)
coord_vect_of_boundary_2D_cells = ∂_2' * coord_vect_of_all_3D_cells .% 2
out = coord_vect_of_boundary_2D_cells
boundary_2cells = [ FV[f] for f in findnz(SparseArrays.sparse(out))[1] ]
hpc = Plasm.lar2exploded_hpc(V, boundary_2cells)(1.3,1.3,1.3)
Plasm.view(hpc)
out = ∂_2' * coord_vect_of_segment .% 2
boundary_2cells = [ FV[f] for f in findnz(SparseArrays.sparse(out))[1] ]
hpc = Plasm.lar2exploded_hpc(V, boundary_2cells)(1.,1.,1.)
Plasm.view(hpc)


text:
#=

  Set partition problem in Julia ConstraintSolver.jl

  Problem formulation from
    http://www.koalog.com/resources/samples/PartitionProblem.java.html
  """
   This is a partition problem.
   Given the set S = {1, 2, ..., n},
   it consists in finding two sets A and B such that:

    *  A U B = S
    *  |A| = |B|
    * sum(A) = sum(B)
    * sum_squares(A) = sum_squares(B)
  """

  This model was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function set_partition(n=8,num_sets=2,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            # "branch_strategy" => :IMPS, # default
                                                            "branch_strategy" => :ABS, # Activity Based Search
                                                            "activity.decay" => 0.999, # default 0.999
                                                            "activity.max_probes" => 100, # default, 10
                                                            "activity.max_confidence_deviation" => 1, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>16,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    println("n:$n num_sets:$num_sets")
    n2 = round(Int, n / num_sets)
    if n % 4 != 0
        println("n ($n) must be a multiple of 4")
        return
    end

    # Which set (1..num_sets) is  1..n assigned to
    @variable(model, 1 <= a[1:n] <= num_sets, Int)

    # number of integers in each set
    @variable(model, n2 <= counts[1:num_sets] <= n2, Int)

    @variable(model, 0 <= sums[1:num_sets] <= n*n, Int)
    @variable(model, 0 <= sum_squared[1:num_sets] <= n*n*n, Int)

    for s in 1:num_sets
        count_ctr(model, a, :(==), s, counts[s])
        b = @variable(model, [1:n],Bin)
        for i in 1:n
            @constraint(model, b[i] := {a[i] == s})
        end
        @constraint(model,sums[s] == sum(i*b[i] for i in 1:n))
        @constraint(model,sum_squared[s] == sum(i*i*b[i] for i in 1:n))
    end

    for s in 1:num_sets-1
        @constraint(model,sums[s] == sums[s+1])
        @constraint(model,sum_squared[s] == sum_squared[s+1])
    end

    # Symmetry breaking
    @constraint(model,a[1] == 1)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                a_val = convert.(Integer,JuMP.value.(a; result=sol))
                sums_val = convert.(Integer,JuMP.value.(sums; result=sol))
                sum_squared_val = convert.(Integer,JuMP.value.(sum_squared; result=sol))
                counts_val = convert.(Integer,JuMP.value.(counts; result=sol))
                println("a:$a_val")
                for s in 1:num_sets
                    println("set $s: ", [i for i in 1:n if a_val[i] == s])
                end
                println("sums:$sums_val")
                println("sum_squared:$sum_squared_val")
                println("counts:$counts_val")
                println()
            end
        end
    else
        println("status:$status")
    end

    return status
end

# @time set_partition(16,2,true,false)
function testit()
    for i in 2:10
        n = 4*i
        println("\nn:$n")
        @time status = set_partition(n,2,true,false)
        if status != MOI.OPTIMAL
            break
        end
    end
end

testit()

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using Documenter, CovarianceMatrices

makedocs(modules = [CovarianceMatrices],
         format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
         clean = true,
         sitename = "CovarianceMatrices.jl",
         authors = "Giuseppe Ragusa",
         checkdocs = :exports,
         # strict = true,
         pages = ["Introduction" => "introduction.md",
                  #"Long run covariances" => "lrcov.md",
                  #"GLM covariances" => "glmcov.md"])a
                 ])

deploydocs(repo = "github.com/gragusa/CovarianceMatrices.jl.git")


text:
export getFirstArival
function getFirstArival(dataTraces::Array{Array{Float32,2},1},mask,x_rcv_loc::Array{Float64},x_src_loc::Array{Float64},vmax::AbstractFloat,vmin::AbstractFloat,window::Int64,dt::Float64,offset)

dobsEikPicked = zeros(Float64,size(mask));

for k = 1:length(dataTraces)
	Dk = dataTraces[k];
	sigma = maximum(abs.(Dk[:]));
	# print("Mean Signal: ")
	# println(sigma)

	dobsEik_k = vec(dobsEikPicked[:,k]);
	Dk = copy(Dk);

	Dk[:,3:end-2] = (Dk[:,1:end-4] + 1.5*Dk[:,2:end-3] + 2*Dk[:,3:end-2] + 1.5*Dk[:,4:end-1] + Dk[:,5:end] )./7;
	Dk[:,3:end-2] = (Dk[:,1:end-4] + 1.5*Dk[:,2:end-3] + 2*Dk[:,3:end-2] + 1.5*Dk[:,4:end-1] + Dk[:,5:end] )./7;
	dobsEikPicked[:,k] = pickTime(Dk,mask[:,k],x_rcv_loc,x_src_loc[k],vmax,vmin,window,dt,offset); # Try to compute travel time
end

return dobsEikPicked;
end

function pickTime(dataTrace::Array,mask::Array,x_rcv_loc::Array{Float64},x_src_loc::Float64,vmax::AbstractFloat,vmin::AbstractFloat,window::Int64,dt::Float64,zeroTimeOffset::Int64)
dpredEik = zeros(Int64,size(dataTrace,1));
println("window = ",window);
println("zeroTimeOffset = ",zeroTimeOffset);
println("vmax = ",vmax);
println("vmin = ",vmin);
window = window;
halfWindow = div(window,2);

for ii = 1:length(dpredEik)

	if mask[ii]==0 || x_rcv_loc[ii] < x_src_loc
		continue;
	end
	windowIdx = ii;
	lookAtPrev = ii > 1 && mask[ii-1]!=0;
	lookAtNext = ii < length(dpredEik) &&  mask[ii+1]!=0;
	if lookAtPrev && lookAtNext
		windowIdx = (ii-1):(ii+1);
	elseif lookAtPrev
		windowIdx = (ii-1):ii;
	elseif lookAtNext
		windowIdx = ii:(ii+1);
	end

	jjmin = 3*window;
	jjmax = size(dataTrace,2)-3*window;
	determineZeroTimeOffset = 0;
	if lookAtPrev && x_rcv_loc[ii-1] >= x_src_loc
		# jjmin = dpredEik[ii-1] + convert(Int64,round((x_rcv_loc[ii] - x_rcv_loc[ii-1])/(vmax*dt)));
		# jjmax = dpredEik[ii-1] + convert(Int64,round((x_rcv_loc[ii] - x_rcv_loc[ii-1])/(vmin*dt)));
	else
		# jjmin = convert(Int64,round((x_rcv_loc[ii] - x_src_loc)/(vmax*dt)))+zeroTimeOffset;
		determineZeroTimeOffset=1;
		# we know that this point is very close to the source and we know the velocity in the sea.
		# jjmax = convert(Int64,round((x_rcv_loc[ii] - x_src_loc)/(vmin*dt)))+zeroTimeOffset;
		# jjmin = jjmax;
	end

	vals = zeros(max(jjmax - jjmin +1,0));
	valsDeriv = zeros(max(jjmax - jjmin +1,0));
	# println(size(s_tag))

	for jj = jjmin:jjmax
		nearWindowSignal = dataTrace[windowIdx,(jj-window+1):jj];
		farWindowSignal = dataTrace[windowIdx,(jj-2*window):(jj-window)];
		vals[jj - jjmin + 1 ] = norm(nearWindowSignal)/norm(farWindowSignal);

		forwardWindowSignal = sum(vec(dataTrace[windowIdx,(jj+1):jj+halfWindow+1]));
		backwardWindowSignal = sum(vec(dataTrace[windowIdx,(jj-halfWindow):jj]));
		valsDeriv[jj - jjmin + 1] = forwardWindowSignal - backwardWindowSignal;

		# println("placing vals: ")
		if jj > jjmin
			if vals[jj - jjmin + 1 ] > 2.0 && valsDeriv[jj - jjmin + 1]*valsDeriv[jj - jjmin] < 0.0
				dpredEik[ii] = jj;
				break;
			end
		end

	end

	if determineZeroTimeOffset==1
		zeroTimeOffset = dpredEik[ii] - convert(Int64,round((x_rcv_loc[ii] - x_src_loc)/(vmin*dt)));
	end
	if dpredEik[ii] == 0
		warn("Could not find travel time - maybe recorded time is too short..., Aborting");
		break;
	end

end

for ii = length(dpredEik):-1:1
	if mask[ii]==0 || x_rcv_loc[ii] >= x_src_loc
		continue;
	end
	windowIdx = ii;
	lookAtPrev = ii > 1 && mask[ii-1]!=0;
	lookAtNext = ii < length(dpredEik) &&  mask[ii+1]!=0;
	if lookAtPrev && lookAtNext
		windowIdx = (ii-1):(ii+1);
	elseif lookAtPrev
		windowIdx = (ii-1):ii;
	elseif lookAtNext
		windowIdx = ii:(ii+1);
	end

	jjmin = 3*window;
	jjmax = size(dataTrace,2)-3*window;
	determineZeroTimeOffset = 0;
	if lookAtNext && x_rcv_loc[ii+1] < x_src_loc
	else
		determineZeroTimeOffset=1;
	end

	vals = zeros(max(jjmax - jjmin +1,0));
	valsDeriv = zeros(max(jjmax - jjmin +1,0));
	for jj = jjmin:jjmax
		nearWindowSignal = dataTrace[windowIdx,(jj-window+1):jj];
		farWindowSignal = dataTrace[windowIdx,(jj-2*window):(jj-window)];
		vals[jj - jjmin + 1 ] = norm(nearWindowSignal)/norm(farWindowSignal);

		forwardWindowSignal = sum(vec(dataTrace[windowIdx,(jj+1):jj+halfWindow+1]));
		backwardWindowSignal = sum(vec(dataTrace[windowIdx,(jj-halfWindow):jj]));
		valsDeriv[jj - jjmin + 1] = forwardWindowSignal - backwardWindowSignal;

		# println("placing vals: ")
		if jj > jjmin
			if vals[jj - jjmin + 1 ] > 2.0 && valsDeriv[jj - jjmin + 1]*valsDeriv[jj - jjmin] < 0.0
				dpredEik[ii] = jj;
				break;
			end
		end
	end

	if determineZeroTimeOffset==1
		zeroTimeOffset = dpredEik[ii] - convert(Int64,round(-(x_rcv_loc[ii] - x_src_loc)/(vmin*dt)));
	end
	if dpredEik[ii] == 0
		warn("Could not find travel time - maybe recorded time is too short..., Aborting");
		break;
	end

end

dpredEik = convert(Array{Float64},dpredEik);
for ii=1:length(dpredEik)
	dpredEik[ii] = max(dpredEik[ii] - zeroTimeOffset,0)*dt;
end

return dpredEik
end

function pickTimeSimple(dataTrace::Array,mask::Array,x_rcv_loc::Array{Float64},x_src_loc::Float64,vmax::AbstractFloat,vmin::AbstractFloat,window::Int64,dt::Float64,offset::Int64)
dpredEik = zeros(Int64,size(dataTrace,1));
window = 2*window;
for ii = 1:length(dpredEik)
	if mask[ii]==0
		continue;
	end
	windowIdx = ii;
	lookAtPrev = ii > 1 && mask[ii-1]!=0;
	lookAtNext = ii < length(dpredEik) &&  mask[ii+1]!=0;
	if lookAtPrev && lookAtNext
		windowIdx = (ii-1):(ii+1);
	elseif lookAtPrev
		windowIdx = (ii-1):ii;
	elseif lookAtNext
		windowIdx = ii:(ii+1);
	end
	jjmin = 2*window+1;
	jjmax = size(dataTrace,2) - 4*window

	vals = zeros(jjmax - jjmin + 1);
	for jj = jjmin : jjmax
		windowSignal = vec(dataTrace[windowIdx,jj-2*window:(jj-window)]);
		wideWindowSignal = vec(dataTrace[windowIdx,(jj-window+1):jj]);
		stdShort = norm(windowSignal);
		stdLong = norm(wideWindowSignal);
		vals[jj - jjmin + 1 ] = abs(stdLong)./abs(stdShort);
		if vals[jj - jjmin + 1 ] > 10.0
			dpredEik[ii] = jj;
			found = 1;
			break;
		end
	end
	dpredEik[ii] = jjmin + indmax(vals) + 1;
end
dpredEik = convert(Array{Float64},dpredEik);
for ii=1:length(dpredEik)
	dpredEik[ii] = max(dpredEik[ii] - offset,0)*dt;
end

return dpredEik
end

############ OLD IDEA FOR BACKUP

# for ii = length(dpredEik):-1:1
	# if mask[ii]==0 || x_rcv_loc[ii] >= x_src_loc
		# continue;
	# end
	# windowIdx = ii;
	# lookAtPrev = ii > 1 && mask[ii-1]!=0;
	# lookAtNext = ii < length(dpredEik) &&  mask[ii+1]!=0;
	# if lookAtPrev && lookAtNext
		# windowIdx = (ii-1):(ii+1);
	# elseif lookAtPrev
		# windowIdx = (ii-1):ii;
	# elseif lookAtNext
		# windowIdx = ii:(ii+1);
	# end

	# jjmin = 1;
	# jjmax = size(dataTrace,2);
	# if lookAtNext && x_rcv_loc[ii+1] < x_src_loc
		# jjmin = dpredEik[ii+1] + convert(Int64,round((x_rcv_loc[ii+1] - x_rcv_loc[ii])/(vmax*dt)));
		# jjmax = dpredEik[ii+1] + convert(Int64,round((x_rcv_loc[ii+1] - x_rcv_loc[ii])/(vmin*dt)));
	# else
		### jjmin = convert(Int64,round(-(x_rcv_loc[ii] - x_src_loc)/(vmax*dt)))+zeroTimeOffset-10;
		# jjmax = convert(Int64,round(-(x_rcv_loc[ii] - x_src_loc)/(vmin*dt)))+zeroTimeOffset;
		# jjmin = jjmax;
	# end
	# found = 0;
	# jjmax = min(jjmax,size(dataTrace,2));
	# found = 0;
	# vals = zeros(max(jjmax - jjmin +1,0));
	# for jj = jjmin:jjmax
		# farWindowSignal = vec(dataTrace[windowIdx,(jj-2*window):(jj-window)]);
		# nearWindowSignal = vec(dataTrace[windowIdx,(jj-window+1):jj]);
		# stdNear = norm(nearWindowSignal);
		# stdFar = norm(farWindowSignal);
		# vals[jj - jjmin + 1 ] = (abs(stdNear)./abs(stdFar));
		# if vals[jj - jjmin + 1 ] > 4
			# dpredEik[ii] = jj;
			# found = 1;
			# break;
		# end
	# end
	### println(vals);
	# if found == 0 && jjmax >= jjmin
		# dpredEik[ii] = jjmin + indmax(vals) + 1;
	# end

	# if dpredEik[ii] == 0
		# warn("Could not find travel time - maybe recorded time is too short..., Aborting");
		# break;
	# end

# end


text:
const TOKEN_STORAGE_MAGIC = b"HDTS"
const TOKEN_STORAGE_VERSION = UInt8(0)

# token protocol: https://github.com/hanborq/hadoop/blob/master/src/core/org/apache/hadoop/security/Credentials.java
struct Credentials
    tokens::Dict{AbstractString,TokenProto}
    secrets::Dict{AbstractString,Vector{UInt8}}

    function Credentials()
        new(Dict{AbstractString,TokenProto}(), Dict{AbstractString,Vector{UInt8}}())
    end
end

function read_credentials!(token_file::String; credentials::Credentials=Credentials())
    open(token_file, "r") do io
        read_credentials!(io; credentials=credentials)
    end
    credentials
end

function read_credentials!(io::IO; credentials::Credentials=Credentials())
    hdts = Vector{UInt8}(undef, 4)
    readbytes!(io, hdts)
    (hdts == TOKEN_STORAGE_MAGIC) || error("Bad header found in token storage")
    version = read(io, UInt8)
    (version == TOKEN_STORAGE_VERSION) || error("Unknown version $version in token storage (expected $TOKEN_STORAGE_VERSION)")

    # read tokens
    size = ProtoBuf.read_varint(io, Int32)
    for idx in 1:size
        alias = ProtoBuf.read_string(io)

        identifier = ProtoBuf.read_bytes(io)
        password = ProtoBuf.read_bytes(io)
        kind = ProtoBuf.read_string(io)
        service = ProtoBuf.read_string(io)
        token = TokenProto(; identifier=identifier, password=password, kind=kind, service=service)

        if isempty(alias)
            alias = string(token.service, token.kind)
        end

        credentials.tokens[alias] = token
    end

    # read secrets
    size = ProtoBuf.read_varint(io, Int32)
    for idx in 1:size
        alias = ProtoBuf.read_string(io)
        secret = ProtoBuf.read_bytes(io)
        credentials.secrets[alias] = secret
    end
    credentials
end

add_token!(credentials::Credentials, token::TokenProto) = add_token!(credentials, string(token.service, token.kind), token)
add_token!(credentials::Credentials, alias::AbstractString, token::TokenProto) = (credentials.tokens[alias] = token; nothing)

_match_alias(alias, pattern) = isempty(pattern) || (pattern == alias)
_match_kind(token, pattern) = isempty(pattern) || (pattern == token.kind)
_match_token(dict_entry, alias_pattern, kind_pattern) = _match_alias(first(dict_entry), alias_pattern) && _match_kind(last(dict_entry), kind_pattern)
find_tokens(credentials::Credentials; alias::AbstractString="", kind::AbstractString="") = collect(values(filter(x->_match_token(x, alias, kind), credentials.tokens)))


text:
export
    OptimalRHS,
    optimize_rhs,
    DPSolution,
    value_iteration,
    value_residual,
    solve_residual,
    nonlinear_solve_value,
    bellman_residual_plot

"""
Object characterizing the optimal solution (a value and a policy) to
the right hand side of a Bellman equation at a given state.
"""
immutable OptimalRHS{TV,TP}
    value::TV
    policy::TP
end

"""
Accessor function for the value.
"""
value(o::OptimalRHS) = o.value

"""
Accessor function for the optimal policy.
"""
policy(o::OptimalRHS) = o.policy

"""
`optimize_rhs(model, value, state)` is the core function for solving
dynamic programming problems.

The user should implement this function for each type of `model` that
will be solved by the library. For a given value function `value` and
state `state`, it should return an OptimalRHS object, containing the
optimal value and the optimal policy.
"""
function optimize_rhs end

"""
Solution to a dynamic programming problem.
"""
@with_kw immutable DPSolution
    model
    value
    policy
    iterations
    converged
end

value(dp::DPSolution) = dp.value
policy(dp::DPSolution) = dp.policy

"""
Solve a dynamic programming problem by value iteration.
"""
function value_iteration(model, V, options = IterationOptions())
    @unpack max_iter = options
    basis = V.basis
    states = collocation_points(basis)
    for iteration in 1:max_iter
        VP = [optimize_rhs(model, V, s) for s in states]
        V′points = value.(VP)
        Δ = V′points - collocation_values(V)
        V′ = basis \ V′points
        converged = _converged(Δ, options)
        if converged || iteration == max_iter
            return DPSolution(model = model,
                              value = V′,
                              policy = basis \ policy.(VP),
                              iterations = iteration,
                              converged = converged)
        else
            V = V′
        end
    end
end

function value_iteration(solution::DPSolution,
                         options = IterationOptions())
    @unpack model, value = solution
    value_iteration(model, value, options)
end

"""
Calculate the residual between the right- and left-hand sides of a
Bellman equation at the given state.
"""
function value_residual(model, value, state)
    optimize_rhs(model, value, state).value - value(state)
end

function value_residual(solution::DPSolution, state)
    @unpack model, value = solution
    value_residual(model, value, state)
end

"""
Solve for `residual(f) = 0` in a linearly interpolated space of
functions.

The collocation points are determined from `initial_f`, which is also
used as a starting point.

`residual_function(model, f, s)` will be called with all `s`
collocation points, and the result is used as the residual. A
nonlinear solver is used to find `f`.

`f` and the result of `nlsolve` are returned as two values.
"""
function solve_residual(model,
                        residual_function,
                        initial_f;
                        nlsolve_options...)
    @unpack basis, α = initial_f
    states = collocation_points(basis)
    function f!(α, residual)
        f = basis * α
        residual .= [residual_function(model, f, state) for state in states]
    end
    result = nlsolve(f!, α; nlsolve_options...)
    basis * result.zero, result
end

function nonlinear_solve_value(solution::DPSolution; nlsolve_options...)
    value, root = solve_residual(solution.model, value_residual,
                                 solution.value; nlsolve_options...)
    policy(state) = optimize_rhs(solution.model, value, state).policy
    DPSolution(model = solution.model,
               value = value,
               policy = value.basis \ policy,
               iterations = root.iterations,
               converged = NLsolve.converged(root))
end

"""
Plot the residual of the Bellman equation of `solution`, using `N`
points, for a univariate solution. `statename` labels the x axis.
"""
function bellman_residual_plot(solution; statename="state", N=100)
    basis = solution.value.basis
    dense_z = linspace(basis, N)
    residual = [value_residual(solution, z) for z in dense_z]
    plot(dense_z, residual, xlab=statename, ylab="residual of Bellman equation",
         legend = false)
    scatter!(collocation_points(basis), zeros(degf(basis)))
end


text:
include( "uses.jl")
include( "types.jl")
include( "logger.jl")
## FIXME either use or replace
include( "examples.jl")
include( "display_constants.jl")
include( "static_texts.jl")
## FIXME shouldn't be needed
include( "text_html_libs.jl")
include( "runner_libs.jl")
include( "table_libs.jl")
include( "base_results.jl")

include( "dash_libs.jl")

const PREAMBLE = """

Every year the Scottish Government [sets its annual budget](https://www.gov.scot/publications/scottish-budget-2021-22/).

This page lets you use a microsimulation tax-benefit model to experiment with some of the most important things that can be changed in the budget,
and with some that currently can't because they are reserved to Westminster (reserved items are shown with a grey background).

You can experiment with the difficult choices involved in balancing the need for fairness and equality against the need not to discourage people from working and saving.

For simplicity, only a few key items can be changed on this page. The [full model](https://github.com/grahamstark/ScottishTaxBenefitModel.jl) allows changing
practically all aspects of the Scottish fiscal system that directly affect individuals. Full instructions on installing and using the full model on your own computer will follow presently.

"""

"""
Create the block of sliders and radios on the LHS
"""
function get_input_block()
	return dbc_form([
		dbc_row([
			dbc_col( it_fieldset() ),
			dbc_col( ni_fieldset() ),
			dbc_col( ben_fieldset() )
		])
		dbc_row([
			dbc_col([
				submit_button()
			]) # col
		]) # row
	]) # form
end

app = dash(
	external_stylesheets=[dbc_themes.UNITED],
	url_base_pathname="/scotbudg/" )
# BOOTSTRAP|SIMPLEX|MINTY|COSMO|SANDSTONE|UNITED|SLATE|SOLAR|UNITED|

app.layout = dbc_container(fluid=true, className="p-5") do
	html_title( "A Budget for Scotland")
	html_h1("A Budget for Scotland"),
	dcc_markdown( PREAMBLE ),
	get_input_block(),
	dcc_loading(
		id="model_running",
		type="default",
		children = [
			html_div(
				id="loading-output-1"
			) # div
		] # children
	), # loading
	html_div( id="output-block" ),
	html_div( dcc_markdown( ENDNOTES ))
end # layout

function do_output( br, hr, tr, pa, ni_prim, ni_sec, cb, pen, uct, ucs, wtcb, scp, scp_age )
	results = nothing
	cache_key = "$br, $hr, $tr, $pa, $ni_prim, $ni_sec, $cb, $pen, $uct, $ucs, $wtcb, $scp, $scp_age"
	sys = deepcopy( BASE_PARAMS )
	settings = deepcopy( BASE_SETTINGS )
	settings.uuid = UUIDs.uuid4()
	# settings.dump_frames = true
	if (br != 20) || (hr !=41)||(tr !=46)||(uct != 55 )||(cb != 21.15)||(pen!= 179.60)||(scp!=10)||(pa!=12_570)||(ni_prim!=12)||(ni_sec!=13.8)||(ucs!=324.84)||(wtcb!=2_005) || (scp_age != 5)
		br /= 100.0
		hr /= 100.0
		tr /= 100.0
		uct /= 100.0
		pa /= WEEKS_PER_YEAR
		wtcb /= WEEKS_PER_YEAR
		ucs /= WEEKS_PER_MONTH

		ni_prim /= 100.0
		ni_sec /= 100.0

		if br == 0
			sys.it.non_savings_rates[1:3] .= 0.0
		else
			bincr = br-sys.it.non_savings_rates[2]
			sys.it.non_savings_rates[1:3] .+= bincr
			sys.it.non_savings_rates[1] = max(0, sys.it.non_savings_rates[1])
		end
		sys.it.non_savings_rates[4] = hr
		sys.it.non_savings_rates[5] = tr
		sys.it.personal_allowance = pa
		sys.uc.taper = uct
		sys.lmt.working_tax_credit.basic = wtcb

		if ucs == 0
			sys.uc.age_25_and_over = 0.0
			sys.uc.age_18_24 = 0.0
			sys.uc.couple_both_under_25 = 0.0
			sys.uc.couple_oldest_25_plus = 0.0
		else
			ucsd = ucs - sys.uc.age_25_and_over
			# move main uc allows equally, as in covid uplift
			sys.uc.age_25_and_over = max(0.0, ucsd+sys.uc.age_25_and_over)
			sys.uc.age_18_24 = max(0.0, ucsd+sys.uc.age_18_24)
			sys.uc.couple_both_under_25 = max(0.0, sys.uc.couple_both_under_25+ucsd)
			sys.uc.couple_oldest_25_plus = max(0.0,sys.uc.couple_oldest_25_plus+ucsd)
		end

		@debug "sys.nmt_bens.child_benefit.first_child was $(sys.nmt_bens.child_benefit.first_child) now $cb"
		sys.nmt_bens.child_benefit.first_child = cb
		@debug "sys.nmt_bens.pensions.new_state_pension was $(sys.nmt_bens.pensions.new_state_pension) now $pen"
		sys.nmt_bens.pensions.new_state_pension = pen
		sys.scottish_child_payment.amount = scp
		sys.scottish_child_payment.maximum_age = scp_age
		sys.ni.primary_class_1_rates[3] = ni_prim
		sys.ni.secondary_class_1_rates[2:3] .= ni_sec

		results = do_run_a( cache_key, sys, settings )
	else
		@debug "returning base results"
		results = BASE_RESULTS
	end
	println("sys.it.non_savings_rates $(sys.it.non_savings_rates)")
	println("BASE_PARAMS.it.non_savings_rates $(BASE_PARAMS.it.non_savings_rates)")
	return make_output_table( results, sys )
end

function no_nothings( things ...)::Bool
	for thing in things
		if isnothing( thing )
			return false
		end
	end
	return true
end

callback!(
    app,
    Output("model_running",  "children"),
	Output("output-block", "children"),
	Input("submit-button", "n_clicks"),

	State( "basic_rate", "value"),
	State( "higher_rate", "value"),
	State( "top_rate", "value"),
	State( "pa", "value"),

	State( "ni_prim", "value"),
	State( "ni_sec", "value"),

	State( "cb", "value"),
	State( "pen", "value"),
	State( "uctaper", "value"),
	State( "ucs", "value"),
	State( "wtcb", "value"),
	State( "scp", "value"),
	State( "scp_age", "value")

	) do n_clicks, basic_rate, higher_rate, top_rate, pa,
	     ni_prim, ni_sec,
		 cb, pen, uctaper, ucs, wtcb, scp, scp_age

	println( "n_clicks = $n_clicks")
	# will return 'nothing' if something is out-of-range or not a number, or if no clicks on submit
	if no_nothings( n_clicks, basic_rate, higher_rate, top_rate, pa, ni_prim, ni_sec, cb, pen, uctaper, ucs, wtcb, scp, scp_age )
		println( "running the live calc version")
		return [nothing, do_output( basic_rate, higher_rate, top_rate, pa, ni_prim, ni_sec, cb, pen, uctaper, ucs, wtcb, scp, scp_age )]
	end
	println( "doing the do-nothing version")
	[nothing, do_output( 20, 41, 46, 12_570, 12, 13.8, 21.15, 179.60, 55, 324.84, 2_005, 10.0, 5 )]
end

run_server( app, "0.0.0.0", 8052; debug=true )


text:
# Find a way to run these example in isolation
using Femshop

init_femshop("poisson1d")

useLog("poisson1dlog")

domain(1)

functionSpace(order=3)
mesh(LINEMESH, elsperdim=20)

u = variable("u")
testSymbol("v")
boundary(u, 1, DIRICHLET, 0)

coefficient("f", "-100*pi*pi*sin(10*pi*x)*sin(pi*x) -
             pi*pi*sin(10*pi*x)*sin(pi*x) + 20*pi*pi*cos(10*pi*x)*cos(pi*x)")
weakForm(u, "-grad(u)*grad(v) - f*v")

import Femshop: grid_data, mesh_data, refel
dofs_per_node = 1

N1 = size(grid_data.allnodes,2);
Nn = dofs_per_node * N1;
Np = refel.Np;
nel = mesh_data.nel;

rhsvec = zeros(Nn);
lhsmatI = zeros(Int, nel*dofs_per_node*Np*dofs_per_node*Np);
lhsmatJ = zeros(Int, nel*dofs_per_node*Np*dofs_per_node*Np);
lhsmatV = zeros(nel*dofs_per_node*Np*dofs_per_node*Np);
allocated_vecs = [rhsvec, lhsmatI, lhsmatJ, lhsmatV];

Femshop.CGSolver.assemble(u, Femshop.bilinears[u.index], Femshop.linears[u.index], allocated_vecs)

print("END")

#SUITE["poisson1d-assemble"] = @benchmarkable Femshop.CGSolver.assemble($u,
                         #$Femshop.bilinears[u.index], $Femshop.linears[u.index])