text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

relaxed_tols = (default_tol_relax = 100,)
insts = OrderedDict()
insts["minimal"] = [
    ((2, 3, 4, 2, true, false, false), nothing, relaxed_tols),
    ((2, 3, 4, 2, true, false, false), :SOCExpPSD),
    ((2, 3, 4, 2, false, true, false),),
    ((2, 3, 4, 2, false, true, false), :ExpPSD),
    ((2, 3, 4, 2, false, false, true),),
    ((2, 3, 4, 2, false, false, true), :ExpPSD),
    ]
insts["fast"] = [
    ((3, 5, 7, 2, true, false, false),),
    ((3, 5, 7, 2, true, false, false), :SOCExpPSD),
    ((3, 5, 7, 2, false, true, false),),
    ((3, 5, 7, 2, false, true, false), :ExpPSD),
    ((3, 5, 7, 2, false, false, true),),
    ((3, 5, 7, 2, false, false, true), :ExpPSD),
    ((5, 15, 25, 5, true, false, false),),
    ((5, 15, 25, 5, false, true, false),),
    ((5, 15, 25, 5, false, false, true),),
    ((10, 30, 50, 5, true, false, false),),
    ((10, 30, 50, 5, false, true, false),),
    ((10, 30, 50, 5, false, false, true),),
    ((25, 75, 125, 10, true, false, false),),
    ((25, 75, 125, 10, false, true, false),),
    ((25, 75, 125, 10, false, false, true),),
    ]
insts["various"] = [
    ((25, 45, 75, 10, true, false, false),),
    ((25, 45, 75, 10, true, false, false), :SOCExpPSD),
    ((25, 45, 75, 10, false, true, false),),
    ((25, 45, 75, 10, false, true, false), :ExpPSD),
    ((25, 45, 75, 10, false, true, false), :SOCExpPSD),
    ((25, 45, 75, 10, false, false, true),),
    ((25, 45, 75, 10, false, false, true), :ExpPSD),
    ((25, 45, 75, 10, false, false, true), :SOCExpPSD),
    ((50, 90, 150, 10, true, false, false),),
    ((50, 90, 150, 10, true, false, false), :SOCExpPSD),
    ((50, 90, 150, 10, false, true, false),),
    ((50, 90, 150, 10, false, true, false), :ExpPSD),
    ((50, 90, 150, 10, false, true, false), :SOCExpPSD),
    ((50, 90, 150, 10, false, false, true),),
    ((100, 180, 300, 10, true, false, false),),
    ((100, 180, 300, 10, false, true, false),),
    ]
return (DOptimalDesignJuMP, insts)


text:
# Test to see if the rejection algorithm is working !

push!(LOAD_PATH, "../models")
push!(LOAD_PATH, "../")
using SparseInverseProblems
using SuperResModels
using TestCases
using Utils

using PyCall
@pyimport numpy as np
@pyimport matplotlib.pyplot as plt

### Static parameters and container for use.
f_c = 20
x_max = 1.0
num_x = 10*f_c
filter_x =  ones(2*f_c+1);
model_static = SuperResModels.Fourier1d(f_c, x_max, filter_x, num_x)
### Dynamic parameters and container for use.
K = 2
tau = 1.0/(K*1.0)
# the maximum velocity on the domain is given by x_max/2 + K tau v_max = 1
v_max = 0.5
num_v = 10
model_dynamic = SuperResModels.DynamicFourier1d(model_static, v_max, tau, K, num_v*K)

# How particles are generated
test_case = () -> TestCases.cloud_1d_full(x_max, v_max, 1, 1, K, tau, rand(4:10))

### We generate particles and we evaluate their separation
N = 10000
separations = zeros(N)
dataFolder = "../data"
for i = 1:N
    (thetas, weights) = Utils.Rejection_sampling(test_case,dataFolder, K,tau,x_max)
    separations[i] = Utils.separation_val(thetas,K,tau,x_max)
end

plt.figure()
plt.hist(separations,100)
plt.show()

ext: jl


text:
# Analytical solution for systems with no coupling between surfaces.

"""
Exact solution for an uncoupled system.
"""
struct Analytical <: Solution
    "Partition function."
    Z::Float64
    "Energy."
    E::Float64
    "Heat capacity."
    Cv::Float64
end

"""
    Analytical(sys::DiagonalSystem{S,M}, beta::Float64)

Calculate the solution for `sys` at `beta`.
"""
function Analytical(sys::DiagonalSystem{S,M}, beta::Float64) where {S,M}
    Zs = exp.(-beta * diag(sys.energy))
    E1s = zero(Zs)
    E2s = zero(Zs)

    for s in 1:S
        # Find the normal mode frequencies for this surface.
        freq_sqrt = sqrt.(sys.freq[:, s])
        lin_p = sys.lin[:, s, s] .* freq_sqrt
        A = diagm(0 => sys.freq[:, s].^2) + sys.quad[:, :, s, s] .* (freq_sqrt * freq_sqrt')
        issymmetric(A) || @warn "Asymmetric A"
        F = eigen(Symmetric(A))
        any(F.values .< 0) && error("Imaginary normal mode frequencies")
        lambdas = sqrt.(F.values)
        T = F.vectors'
        lin_pp = T * lin_p

        E_ = sys.energy[s, s]

        for m in 1:M
            x = exp(-0.5 * beta * lambdas[m])
            E_ -= 0.5 * (lin_pp[m] / lambdas[m])^2

            Zs[s] *= exp(0.5 * beta * (lin_pp[m] / lambdas[m])^2)
            Zs[s] *= x / (1-x^2)

            E1s[s] += lambdas[m] * (x^2+1) / (2 * (1-x^2))
        end

        E1s[s] += E_
        E2s[s] += E_^2

        for m1 in 1:M
            x1 = exp(-0.5 * beta * lambdas[m1])

            E2s[s] += 2.0 * E_ * lambdas[m1] * (x1^2+1) / (2 * (1-x1^2))
            E2s[s] += lambdas[m1]^2 * (x1^4+6x1^2+1) / (4 * (1-x1^2)^2)

            for m2 in 1:M
                m1 == m2 && continue

                x2 = exp(-0.5 * beta * lambdas[m2])

                E2s[s] += lambdas[m1] * (x1^2+1) / (2 * (1-x1^2)) * lambdas[m2] * (x2^2+1) / (2 * (1-x2^2))
            end
        end
    end

    Z = sum(Zs)
    E = sum(E1s .* Zs / Z)
    Cv = (sum(E2s .* Zs / Z) - E^2) * beta^2

    Analytical(Z, E, Cv)
end


text:
using QuantumNLDiffEq
using Test
using Yao: dispatch, EasyBuild, put, Z, chain, Ry, parameters, zero_state, expect
using Zygote: gradient

≃(a, b) = round(a, digits = 2) ==  round(b, digits = 2)

@testset "Check if chainrules_patch works" begin
	N = 6
	cost = sum([put(N, i=>Z) for i in 1:N])
	theta = rand(96)
	new_var = dispatch(EasyBuild.variational_circuit(N,5), :random);
	function new_loss(theta, j)
		new_circ = dispatch(new_var, theta)
		diff_evalue = expect(cost, zero_state(N) => new_circ)
		return real.(diff_evalue)*j
	end

	p = gradient(_theta->new_loss(_theta, 1), theta)[1][3]
	d = gradient(_theta->new_loss(_theta, 3), theta)[1][3]
	@test d ≃ p*3
end

@testset "Tests for Chebyshev sparse mapping" begin
	M = range(0, stop=0.9, length=20)
        cost = sum([put(6, i=>Z) for i in 1:6])
        var = dispatch(EasyBuild.variational_circuit(6,5), :random)
        quantum_feature_map_circuit = chain(6, [put(i=>Ry(0)) for i in 1:6])
        params = parameters(var)
        u_0 = 1.0
        N = 6
        train!(quantum_feature_map_circuit, cost, var, M, N, 1.0, params)

	@testset "Check if FD works" begin
		t = 6
		evalue(x) = expect(cost, zero_state(N)=>QuantumNLDiffEq.new_circuit(quantum_feature_map_circuit, var, x, params, N, 2; mapping=QuantumNLDiffEq.ChebyshevSparse()))

		x_fd = M[t]
		dx = 0.0001
		x_fd += dx
		loss_grad_fd = (evalue(x_fd) - evalue(M[t]))/dx

		@test QuantumNLDiffEq.calculate_diff_evalue(quantum_feature_map_circuit, cost, var, N, params, M[t], 2; mapping=QuantumNLDiffEq.ChebyshevSparse()) ≃ loss_grad_fd
	end

	@test loss(quantum_feature_map_circuit, cost, var, M, N, u_0, params, 2; mapping=QuantumNLDiffEq.ChebyshevSparse()) < 1.0
end


text:
using OrdinaryDiffEq, Test

function d_alembert(du,u,p,t)
  du[1] = p[1] - p[2]*u[1] + p[3]*t
end

function d_alembert_jac(J,u,p,t)
  J[1] = -p[2]
end

function d_alembert_analytic(u0,p,t::Number)
    a,b,c = p
    ebt = exp(b*t)
    @. exp(-b*t)*(-a*b + c + ebt*(a*b + c*(b*t - 1)) + b^2 * u0)/(b^2)
end

p = (1., 2., 3.)
u0 = [1.0]
tspan = (0.0,10.0)
prob = ODEProblem(ODEFunction(d_alembert,
                              jac = d_alembert_jac,
                              analytic=d_alembert_analytic),
                              u0,tspan,p)

sol = solve(prob,Tsit5(),abstol=1e-10,reltol=1e-10)
@test sol.errors[:l2] < 1e-7
sol = solve(prob,Rosenbrock23(),abstol=1e-8,reltol=1e-8)
@test sol.errors[:l2] < 1e-7
sol = solve(prob,Rodas4(),abstol=1e-10,reltol=1e-10)
@test sol.errors[:l2] < 1e-7
sol = solve(prob,Veldd4(),abstol=1e-10,reltol=1e-10)
@test sol.errors[:l2] < 1e-7
sol = solve(prob,Rodas5(),abstol=1e-10,reltol=1e-10)
@test sol.errors[:l2] < 1e-7
sol = solve(prob,TRBDF2(),abstol=1e-10,reltol=1e-10)
@test sol.errors[:l2] < 2e-6
sol = solve(prob,Trapezoid(),abstol=1e-10,reltol=1e-10)
@test sol.errors[:l2] < 2e-6
sol = solve(prob,KenCarp3(),abstol=1e-10,reltol=1e-10)
@test sol.errors[:l2] < 8e-4
sol = solve(prob,KenCarp4(),abstol=1e-10,reltol=1e-10)
@test sol.errors[:l2] < 1e-7

using ModelingToolkit
function lotka(du,u,p,t)
  x = u[1]
  y = u[2]
  du[1] = p[1]*x - p[2]*x*y
  du[2] = -p[3]*y + p[4]*x*y
end

prob = ODEProblem(lotka,[1.0,1.0],(0.0,1.0),[1.5,1.0,3.0,1.0])
de = ModelingToolkit.modelingtoolkitize(prob)
prob2 = remake(prob, f=ODEFunction(de...; jac=true, Wfact=true))

sol = solve(prob, TRBDF2())

for Alg in [Rodas5, Rosenbrock23, TRBDF2, KenCarp4]
  @test Array( solve(prob2, Alg(), tstops=sol.t, adaptive=false) ) ≈ Array( solve(prob, Alg(), tstops=sol.t, adaptive=false) ) atol=1e-4
end


text:
module Fractal

using Distributions

export IFS, WIFS, IFSNonlinear
export itrPtsProb

include("IFS.jl")

end


text:
"""
Subroutine ihsort

Given an integer array keys of length n, this subroutine uses
a heap sort to sort the keys in increasing order.

This subroutine is a minor modification of code written by
Mark Jones and Paul Plassmann.

MINPACK-2 Project. March 1998.
Argonne National Laboratory.
Chih-Jen Lin and Jorge J. More'.
"""
function ihsort(n,keys)
    if  n <= 1
        return
    end

    # Build the heap.

    mid = floor(Int,n/2)
    for k=mid:-1:1
        x = keys[k]
        lheap = k
        rheap = n
        m = lheap*2
        while m <= rheap
            if m < rheap
                if keys[m] < keys[m+1]
                    m = m + 1
                end
            end
            if x >= keys[m]
                m = rheap + 1
            else
                keys[lheap] = keys[m]
                lheap = m
                m = 2*lheap
            end
        end
        keys[lheap] = x
    end

    # Sort the heap.

    for k=n:-1:2
        x = keys[k]
        keys[k] = keys[1]
        lheap = 1
        rheap = k-1
        m = 2
        while m <= rheap
            if m < rheap
                if keys[m] < keys[m+1]
                    m = m+1
                end
            end
            if x >= keys[m]
                m = rheap + 1
            else
                keys[lheap] = keys[m]
                lheap = m
                m = 2*lheap
            end
        end
        keys[lheap] = x
    end

    return
end


text:
using Dates

out_dir() = mkpath(normpath(joinpath(@__FILE__, "..", "..", "..", "out")))

snapshot_name(time::DateTime, step::Int) = "bert-abs-$(Dates.format(time, "yyyy-mm-dd-HH:MM"))-step-$(string(step, pad=6))"

snapshot_file(time::DateTime, step::Int, name::String) = joinpath(out_dir(), "$(snapshot_name(time, step))-$name")

function snapshot_files(time::DateTime, inc_steps::Int, name::String)::Tuple{Vector{String}, Int}
    step = inc_steps
    files = []
    push!(files, snapshot_file(time, step, name))
    while isfile(files[end])
        step += inc_steps
        push!(files, snapshot_file(time, step, name))
    end
    return files[1:end-1], step - inc_steps
end


text:
module RateLimiter

using Dates

export AbstractRateLimiter, @rate_limit, NoLimitRateLimiter, TokenBucketRateLimiter

#TODO: should the Manifest be committed?

""" Abstract base type for all rate limiters. """
abstract type AbstractRateLimiter end

"""
    !(limiter, cost::T) where {T <: Real}

Wait until the method can be called again.  This method is thread safe.

# Arguments
- `limiter`: rate limiter
- `cost`: cost of executing the method.
"""
wait!(limiter::AbstractRateLimiter, cost::T) where {T <: Real} = error("Function is not defined for AbstractRateLimiter: type=$(typeof(limiter)) function=wait!")

"""
    @rate_limit(limiter, cost, expr)

A macro to limit the rate an expression can be called.

# Arguments
    - `limiter`: rate limiter
    - `cost`: cost of executing the method.
    - `expr`: expression to be evaluated
"""
macro rate_limit(limiter, cost, expr)
    return esc(quote
        RateLimiter.wait!($limiter, $cost)
        $expr
    end)
end

#######################################################################

"""
A rate limiter which does not limit the rate.
"""
struct NoLimitRateLimiter <: AbstractRateLimiter end

"""
    !(limiter, cost::T) where {T <: Real}

Wait until the method can be called again.  This method is thread safe.

# Arguments
- `limiter`: rate limiter
- `cost`: cost of executing the method.
"""
wait!(limiter::NoLimitRateLimiter, cost::T) where {T <: Real} = return

#######################################################################

"""
Token-Bucket rate limiter.  This limiter is used to control for constraints such as bandwidth and burstiness.

See:  https://en.wikipedia.org/wiki/Token_bucket#Algorithm
"""
mutable struct TokenBucketRateLimiter <: AbstractRateLimiter
    tokens_per_second::Float64
    max_tokens::Float64
    lock::ReentrantLock

    last_update::DateTime
    tokens::Float64

    """
    Token-Bucket rate limiter.  This limiter is used to control for constraints such as bandwidth and burstiness.

    See:  https://en.wikipedia.org/wiki/Token_bucket#Algorithm

    # Arguments
    - `tokens_per_second`: number of tokens added per second.
    - `max_tokens`: maximum number of tokens.
    - `initial_tokens`: initial number of tokens.
    - `sleep_seconds`: number of seconds to sleep before checking to see if execution is possible.
    """
    function TokenBucketRateLimiter(tokens_per_second, max_tokens, initial_tokens)
        tokens_per_second <= 0 && error("tokens_per_second must be positive: tokens_per_second=$tokens_per_second")
        max_tokens <= 0 && error("max_tokens must be positive: max_tokens=$max_tokens")
        return new(tokens_per_second, max_tokens, ReentrantLock(), now(), initial_tokens)
    end
end

"""
    !(limiter, cost::T) where {T <: Real}

Wait until the method can be called again.  This method is thread safe.

# Arguments
- `limiter`: rate limiter
- `cost`: cost of executing the method.
"""
function wait!(limiter::TokenBucketRateLimiter, cost::T) where {T <: Real}
    # error check
    cost < 0 && error("Cost must be greater than zero: cost=$cost")
    cost > limiter.max_tokens && error("Cost must be less than the max number of tokens: cost=$cost max_tokens=$(limiter.max_tokens)")

    lock(limiter.lock) do
        while true
            # fill
            t = now()
            limiter.tokens = min(limiter.max_tokens, limiter.tokens + Dates.toms(t-limiter.last_update) / Dates.toms(Second(1)) * limiter.tokens_per_second)
            limiter.last_update = t

            # take
            if limiter.tokens >= cost
                limiter.tokens -= cost
                return
            end

            # retry
            sleep((cost-limiter.tokens)/limiter.tokens_per_second)
        end
    end
end

#######################################################################

end # module


text:
using Kalman

k = 1.0
p0 = 0.1*k
r = 10.0
deltaT = 0.0005
logit(p0,k,x) = (k*p0*exp(x))/(k+p0*(exp(x)-1))

function an(x,k,r,deltaT)
    x1 = zeros(x)
    x1[1] = x[1]
    x1[2] = logit(x[2],k,r*deltaT)
    return x1
end

fekf(x) = an(x,k,x[1],deltaT)
gekf(x) = eye(2)
hekf(x) = x[2]

Q = [1e-10 0.0; 0.0 1e-10]
R = [1e-2]'

x0 = State([5.0,0.1],[1e2 0.0; 0.0 1e-10])

ekff = NonlinearModel(fekf,gekf,Q)
z = NonlinearObservationModel(hekf,R)

kf0 = BasicExtendedKalmanFilter(x0,ekff,z)
kf = kf0

t = [0:deltaT:deltaT*300]
z = map(t->logit(p0,k,r*t),t)
ys = z.+R[1]*randn(301)

ps = zeros(t)
rs = zeros(t)

kf1 = predictupdate(kf,Observation([ys[1]]))
predictupdate!(kf,Observation([ys[1]]))
@test kf1.x == kf.x

for i in 2:length(t)
    kf = predictupdate(kf,Observation([ys[i]]))
    ps[i] = kf.x.x[2]
    rs[i] = kf.x.x[1]
end


text:
using TrafficAssignment
using Test
using Revise
using BenchmarkTools
using DataFrames
using CSV
using Random
using Printf
Random.seed!(1104)

testnetworks = ["Anaheim"]
algorithms = [:FW, :fukushimaFW, :conjugateFW, :TAPAS]
@testset "TrafficAssignment" begin
    for network in testnetworks
        for algorithm in algorithms
            metadata, report, solution = assigntraffic(network, method=algorithm, assignment=:UE, tol=1e-10, maxiters=50, maxruntime=300, log=:off)
            refsol = DataFrame(CSV.File(joinpath(@__DIR__, "Solution/$network.csv")))
            rg₁ = 10^report[end, :LOG₁₀RG]
            rg₂ = abs(TrafficAssignment.getrg(network, refsol))
            println("$network - $algorithm")
            @printf("   Solution relative gap : %.3e",  rg₁)
            @printf("\n   Benchmark relative gap: %.3e\n\n",  rg₂)
            @test rg₁ ≤ rg₂
        end
    end
end
# TODO: Add more test networks


text:
"""
    abstract type FEOperator <: GridapType

A `FEOperator` contains finite element problem,
that is assembled as far as possible and ready to be solved.
See also [FETerm](@ref)
"""
abstract type FEOperator <: GridapType end

"""
"""
function get_test(op::FEOperator)
  @abstractmethod
end

"""
"""
function get_trial(op::FEOperator)
  @abstractmethod
end

"""
"""
function allocate_residual(op::FEOperator,u)
  @assert is_a_fe_function(u)
  @abstractmethod
end

"""
    $(SIGNATURES)

Inplace version of [`residual`](@ref).
"""
function residual!(b::AbstractVector,op::FEOperator,u)
  @assert is_a_fe_function(u)
  @abstractmethod
end

"""
    $(SIGNATURES)

Compute the residual of `op` at `u`. See also [`residual_and_jacobian`](@ref)
"""
function residual(op::FEOperator,u)
  @assert is_a_fe_function(u)
  b = allocate_residual(op,u)
  residual!(b,op,u)
  b
end

"""
"""
function allocate_jacobian(op::FEOperator,u)
  @assert is_a_fe_function(u)
  @abstractmethod
end

"""
    $(SIGNATURES)

Inplace version of [`jacobian`](@ref).
"""
function jacobian!(A::AbstractMatrix,op::FEOperator,u)
  @assert is_a_fe_function(u)
  @abstractmethod
end

"""
    $(SIGNATURES)

Compute the jacobian of an operator `op`.
See also [`get_algebraic_operator`](@ref), [`residual_and_jacobian!`](@ref).
"""
function jacobian(op::FEOperator,u)
  @assert is_a_fe_function(u)
  A = allocate_jacobian(op,u)
  jacobian!(A,op,u)
  A
end

"""
    $(SIGNATURES)

Inplace version of [`residual_and_jacobian`](@ref).
"""
function residual_and_jacobian!(b::AbstractVector,A::AbstractMatrix,op::FEOperator,u)
  residual!(b,op,u)
  jacobian!(A,op,u)
  (b,A)
end

"""
    residual, jacobian = $(SIGNATURES)

Compute the residual and jacobian of an operator `op` at a given point `u`.
Depending on the nature of `op` the point `u` can either be a plain array or a `FEFunction`.

See also [`jacobian`](@ref), [`residual`](@ref), [`get_algebraic_operator`](@ref).
"""
function residual_and_jacobian(op::FEOperator,u)
  b = residual(op,u)
  A = jacobian(op,u)
  (b,A)
end

"""
"""
function test_fe_operator(op::FEOperator,args...;kwargs...)
  test = get_test(op)
  trial = get_trial(op)
  @test isa(test,FESpace)
  @test isa(trial,FESpace)
  u = zero(trial)
  @test is_a_fe_function(u)
  b = allocate_residual(op,u)
  @test isa(b,AbstractVector)
  residual!(b,op,u)
  b2 = residual(op,u)
  @test isa(b2,AbstractVector)
  A = allocate_jacobian(op,u)
  @test isa(A,AbstractMatrix)
  jacobian!(A,op,u)
  A2 = jacobian(op,u)
  @test isa(A2,AbstractMatrix)
  residual_and_jacobian!(b,A,op,u)
  b, A = residual_and_jacobian(op,u)
  @test isa(b,AbstractVector)
  @test isa(A,AbstractMatrix)
  _op = get_algebraic_operator(op)
  test_nonlinear_operator(_op,args...;kwargs...)
end

"""
    $(SIGNATURES)

Return an "algebraic view" of an operator. Algebraic
means, that the resulting operator acts on plain arrays,
instead of `FEFunctions`. This can be useful for solving
with external tools like `NLsolve.jl`.
See also [`FEOperator`](@ref).
"""
function get_algebraic_operator(feop::FEOperator)
  AlgebraicOpFromFEOp(feop)
end

struct AlgebraicOpFromFEOp <: NonlinearOperator
  feop::FEOperator
end

function allocate_residual(op::AlgebraicOpFromFEOp,x::AbstractVector)
  trial = get_trial(op.feop)
  u = EvaluationFunction(trial,x)
  allocate_residual(op.feop,u)
end

function residual!(b::AbstractVector,op::AlgebraicOpFromFEOp,x::AbstractVector)
  trial = get_trial(op.feop)
  u = EvaluationFunction(trial,x)
  residual!(b,op.feop,u)
end

function residual(op::AlgebraicOpFromFEOp,x::AbstractVector)
  trial = get_trial(op.feop)
  u = EvaluationFunction(trial,x)
  residual(op.feop,u)
end

function allocate_jacobian(op::AlgebraicOpFromFEOp,x::AbstractVector)
  trial = get_trial(op.feop)
  u = EvaluationFunction(trial,x)
  allocate_jacobian(op.feop,u)
end

function jacobian!(A::AbstractMatrix,op::AlgebraicOpFromFEOp,x::AbstractVector)
  trial = get_trial(op.feop)
  u = EvaluationFunction(trial,x)
  jacobian!(A,op.feop,u)
end

function jacobian(op::AlgebraicOpFromFEOp,x::AbstractVector)
  trial = get_trial(op.feop)
  u = EvaluationFunction(trial,x)
  jacobian(op.feop,u)
end

function residual_and_jacobian!(b::AbstractVector,A::AbstractMatrix,op::AlgebraicOpFromFEOp,x::AbstractVector)
  trial = get_trial(op.feop)
  u = EvaluationFunction(trial,x)
  residual_and_jacobian!(b,A,op.feop,u)
end

function residual_and_jacobian(op::AlgebraicOpFromFEOp,x::AbstractVector)
  trial = get_trial(op.feop)
  u = EvaluationFunction(trial,x)
  residual_and_jacobian(op.feop,u)
end

function zero_initial_guess(op::AlgebraicOpFromFEOp)
  trial = get_trial(op.feop)
  x = zero_free_values(trial)
end


text:
__precompile__(true)

module LAPACK

end #module


text:
# Some TODO:
# Rename:
#   LinearFilter -> AbstractLinearFilter
#   LinearFilter_Default -> LinearFilter
#   LinearFilter_DefaultCentered -> LinearFilter_Centered (or only compute the right offset)
#
export LinearFilter
export fcoef, length, offset, range

import Base: length, range, isapprox, show

"""
    abstract type LinearFilter{T<:Number} end

Abstract type defining a linear filter.
"""
abstract type LinearFilter{T<:Number} end

"""
    fcoef(c::LinearFilter)

Returns filter coefficients
"""
fcoef(c::LinearFilter) = c._fcoef

"""
    length(c::LinearFilter)::Int

Returns filter length
"""
length(c::LinearFilter)::Int = length(fcoef(c))

"""
    offset(c::LinearFilter)::Int

Returns filter offset

**Caveat:** the first position is **0** (and not **1**)
"""
offset(c::LinearFilter)::Int = c._offset

# Internal
#
# Computes [[range_filter][]] using primitive types.
# This allows reuse by =directConv!= for instance.
#
# *Caveat:* do not overload Base.range !!!
filter_range(size::Int,offset::Int)::UnitRange = UnitRange(-offset,size-offset-1)

"""
    range(c::LinearFilter)::UnitRange

Returns filter range Ω

Filter support of a filter α is defined by Ω = [ - offset(α), length(α) -  offset(α) - 1 ]
"""
range(c::LinearFilter)::UnitRange = filter_range(length(c),offset(c))

# For convenience only, used in utests
#
function isapprox(f::LinearFilter{T},v::AbstractArray{T,1}) where {T<:Number}
    return isapprox(fcoef(f),v)
end

#
# Pretty print
#
# I defined this afterward to avoid jldoctest to fail because of
# different rounding errors occuring for different architectures.
#
# Now filters coefficient are rounded before printing.
#
function Base.show(io::IO, f::LinearFilter)
    r = range(f)
    coef = fcoef(f)
    print(io,"Filter(r=$r,c=")
    print(io, round.(coef; sigdigits=4))
    println(io,")")
end

"""
    struct LinearFilter_Default{T<:Number,N}

Default linear filter.

You can create a filter as follows

```jldoctest
julia> linear_filter=LinearFilter([1,-2,1],1)
Filter(r=-1:1,c=[1.0, -2.0, 1.0])

julia> offset(linear_filter)
1

julia> range(linear_filter)
-1:1
```
"""
struct LinearFilter_Default{T<:Number,N} <: LinearFilter{T}
    _fcoef::SVector{N,T}
    _offset::Int
end

#+LinearFilter,Internal
# Creates a linear filter from a coefficient vector and its associated offset
#
# *Example:*
#!linear_filter=LinearFilter(rand(3),5)
#!offset(linear_filter)
#!range(linear_filter)
#
function LinearFilter_Default(c::AbstractArray{T,1},offset::Int)  where {T<:Number}
    v=SVector{length(c),T}(c)
    return LinearFilter_Default{T,length(c)}(v,offset)
end

#+LinearFilter,Internal
#
# Default *centered* linear filter
#
# Array length has to be odd, 2n+1. Filter offset is n by construction.
#
struct LinearFilter_DefaultCentered{T<:Number,N} <: LinearFilter{T}
    _fcoef::SVector{N,T}
end

#+LinearFilter,Internal
function LinearFilter_DefaultCentered(c::AbstractArray{T,1}) where {T<:Number}
    N = length(c)
    @assert isodd(length(c)) "Centered filters must have an odd number of coefficients, $N is even"
    return LinearFilter_DefaultCentered{T,N}(SVector{N,T}(c))
end

#+LinearFilter,Internal
offset(f::LinearFilter_DefaultCentered{T,N}) where {T<:Number,N} = (N-1)>>1

# Once that we have defined LinearFilter_Default and
# LinearFilter_DefaultCentered we can unify construction. Switch to
# the right type is decided according to arguments

#+LinearFilter
#
# Creates a linear filter from its coefficients and an offset
#
# The *offset* is the position of the filter coefficient to be aligned with zero, see [[range_filter][]].
#
# *Example:*
#!f=LinearFilter([0:5;],4);
#!hcat([range(f);],fcoef(f))
#
function LinearFilter(c::AbstractArray{T,1},offset::Int)::LinearFilter  where {T}
    return LinearFilter_Default(c,offset)
end

#+LinearFilter
#
# Creates a centered linear filter, it must have an odd number of
# coefficients, $2n+1$ and is centered by construction (offset=n)
#
# *Example:*
#!f=LinearFilter([0:4;]);
#!hcat([range(f);],fcoef(f))
#
function LinearFilter(c::AbstractArray{T,1})::LinearFilter  where {T}
    return LinearFilter_DefaultCentered(c)
end

max_stars_repo_name: vincent-picaud/DirectConvolution

max_issues_repo_name: vincent-picaud/DirectConvolution

max_forks_repo_name: vincent-picaud/DirectConvolution


text:
using Test
using Unfold

firbase = firbasis((-1, 1), 10)

# test optional call
@test Unfold.kernel(firbase) == Unfold.kernel(firbasis((-1, 1), 10))

@test typeof(Unfold.name(firbase)) <: String
# test basics of basisfunction
@test length(collect(Unfold.colnames(firbase))) == 21
@test unique(Unfold.kernel(firbase)(1)) == [1.0, 0.0]

# test length consistency
@test length(Unfold.colnames(firbase)) ==
      length(Unfold.times(firbase)) ==
      size(Unfold.kernel(firbase)(3.1))[2]

# testing the non-sampling rate samples
@test Unfold.kernel(firbase)(0.5)[1:3, 1:3] == [0.5 0.0 0.0; 0.5 0.5 0.0; 0.0 0.5 0.5]

boldbase = hrfbasis(2.0, name = "test")

@test Unfold.name(boldbase) == "test"
@test Unfold.kernel(boldbase)(0) == Unfold.kernel(boldbase)(1)
@test Unfold.kernel(boldbase)(0.1) != Unfold.kernel(boldbase)(1) # testing fractional kernel generation
@test findmax(Unfold.kernel(boldbase)(0.3))[2][1] == 4


text:
abstract type AbstractSpace end

"""
    sample(s::AbstractSpace)

Get a random sample from `s`.
"""
function sample end

"""
    occursin(s::AbstractSpace, x)

Return wheather `x` is a valid sample in space `s`.
"""
function occursin end


text:
function version()
  version = "19.151"
  println("Version: $version")
end


text:
using DataFrames
using CSV
using ArgParse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--datafile"
            help = "Path to ngsim_point.csv"
        "--vehiclenum", "-v"
            help = "Vehicle UniqueID number to focus on"
            arg_type = Int
            default = 3000
        "--radius"
            help = "radius around focus vehicle to include"
            arg_type = Float64
            default = 100.0
        "--trim"
            help = "Number of timesteps to remove from begin and end of focal dataset"
            default = 50
            arg_type = Int
        "--long"
            help = "Whether or not to tidy the dataframe"
            action= :store_true
        "vars"
            help = "Which data-containing variables to include"
            default = ["Local_X", "Local_Y"]
        "--output"
            help = "Path to write output dataset"
            default = "./output_dataset.csv"
    end

    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    println("Parsed args:")
    for (arg,val) in parsed_args
        println("  $arg  =>  $val")
    end
    println("Reading dataset...")
    dataset = DataFrame(CSV.File(parsed_args["datafile"]))
    unique!(dataset, ["timestep", "Unique_Veh_ID"])
    dataset.timestep = round.(Int, dataset.timestep) # stored as floats
    focaldf = filter(:Unique_Veh_ID => x -> x == parsed_args["vehiclenum"], dataset)
    maxtime = maximum(focaldf[!, :timestep])
    mintime = minimum(focaldf[!, :timestep])

    # Select rows
    outdf = filter!(:timestep => x -> mintime <= x <= maxtime, dataset)
    rename!(focaldf, :Global_X => :Focal_X, :Global_Y => :Focal_Y)
    outdf = leftjoin(outdf, focaldf[:, ["timestep", "Focal_X", "Focal_Y"]], on="timestep")
    outdf.dist_to_focal = sqrt.((outdf.Global_X - outdf.Focal_X).^2 + (outdf.Global_Y - outdf.Focal_Y).^2)
    filter!(:dist_to_focal => x -> x <= parsed_args["radius"], outdf)
    # Select columsn
    vars_to_include = copy(parsed_args["vars"])
    pushfirst!(vars_to_include, "timestep")
    pushfirst!(vars_to_include, "Unique_Veh_ID")
    select!(outdf, vars_to_include)
    if parsed_args["long"] == true
        outdf = stack(outdf, parsed_args["vars"], ["timestep", "Unique_Veh_ID"];
                    variable_name="channel", value_name="value")
        rename!(outdf, :timestep => :sample_id, :Unique_Veh_ID => :observation_id)
    end
    CSV.write(parsed_args["output"], outdf)
end

main()


text:
""" Solve the classical GP Regression (especially valid for augmented likelihoods) """
mutable struct Analytic{T<:Real} <: Inference{T}
    ϵ::T #Convergence criteria
    nIter::Integer #Number of steps performed
    Stochastic::Bool #Use of mini-batches
    nSamples::Int64
    nSamplesUsed::Int64 #Size of mini-batches
    MBIndices::Vector{Int64} #Indices of the minibatch
    ρ::T #Stochastic Coefficient
    HyperParametersUpdated::Bool #To know if the inverse kernel matrix must updated
    xview::SubArray{T,2,Matrix{T}}
    yview::SubArray
    function Analytic{T}(ϵ::T,nIter::Integer,Stochastic::Bool,nSamples::Integer,MBIndices::AbstractVector,ρ::T,flag::Bool) where T
        return new{T}(ϵ,nIter,Stochastic,nSamples,nSamples,MBIndices,ρ,flag)
    end
end

"""`Analytic(;ϵ::T=1e-5)`

Analytic inference structure for the classical GP regression

**Keyword arguments**

    - `ϵ::T` : convergence criteria, which can be user defined
"""
function Analytic(;ϵ::T=1e-5) where {T<:Real}
    Analytic{T}(ϵ,0,false,1,collect(1:1),1.0,true)
end

function Base.show(io::IO,inference::Analytic{T}) where T
    print(io,"Analytic Inference")
end

function init_inference(inference::Analytic{T},nLatent::Integer,nFeatures::Integer,nSamples::Integer,nSamplesUsed::Integer) where {T<:Real}
    inference.nSamples = nSamples
    inference.nSamplesUsed = nSamples
    inference.MBIndices = 1:nSamples
    inference.ρ = one(T)
    return inference
end

function analytic_updates!(model::GP{T}) where {T}
    if !isnothing(model.likelihood.opt_noise)
        model.likelihood.σ² = mean(abs2,model.y.-first(model.f).μ)
    end
    model.f[1].μ = first(model.f).K\(model.y - first(model.f).μ₀(model.X))
end


text:
using FixedSizeArrays
using PlotlyJS
include("../src/CHull2D.jl")

#
# LineSegment Intersection
#

# Simple example
p1 = Point(0.0, 0.0)
p2 = Point(0.5, 0.5)
p3 = Point(0.0, 0.5)
p4 = Point(0.5, 0.0)
l1 = CHull2D.LineSegment(p1, p2)
l2 = CHull2D.LineSegment(p3, p4)
println("Should be true", ":", intersect(l1, l2))

# TODO: Should be able to handle Inf instead of 1e12
p5 = Point(0.0, 0.25)
p6 = Point(1e12, 0.25)
l3 = CHull2D.LineSegment(p5, p6)
println("Should be true", ":", intersect(l1, l3))

#
# Convex Hull with Quadrant
#
p = [Point(randn(), randn()) for i=1:25]
ch = CHull2D.convexhull(p)
q = CHull2D.Quadrant(Point(0.0, 0.0), Inf, Inf)
chiq = intersect(ch, q)

t1 = scatter(;x=[foo[1] for foo in chiq], y=[foo[2] for foo in chiq])
t2 = scatter(;x=[foo[1] for foo in ch.extremepoints], y=[foo[2] for foo in ch.extremepoints])
Plot([t1, t2])


text:
function mandel_iter(c::Complex,max_iter::Int,thresh::Int)

    i = 0
    z = 0

    while (i < max_iter) & (abs2(z) < thresh)
        z = z^2 + c
        i += 1

    end

    return i
end

function mandelbrot_protein(codon::Array{Float64,1},thresh::Int,max_iter::Int,n::Int)

    x,y,w = codon

    protein = Array{Int,2}(undef, n, n)

    x_arr = range(x - w/2,x + w/2,length=n)

    y_arr = range(y - w/2,y + w/2, length=n)

    @inbounds for y_i in 1:n

        @inbounds for x_i in 1:n

            protein[x_i,y_i] =  floor(Int,(mandel_iter(Complex{Float32}(x_arr[x_i],y_arr[y_i]),max_iter,thresh) * 255 / max_iter))

        end

    end

    return protein

end

max_stars_repo_name: harrybooth/FractalDevelopment

max_issues_repo_name: harrybooth/FractalDevelopment

max_forks_repo_name: harrybooth/FractalDevelopment


text:
# Contains various membership function types
# ------------------------------------------

"""	 Membership function type
"""
abstract type MF end

"""	 Triangular membership function type

	 TriangularMF(l_vertex, center, r_vertex)

	 Properties
	 ----------
	 `l_vertex`, `center` and `r_vertex` are the vertices of the triangle, in order

	 `eval` function returns membership value at a point
	 `mean_at` function returns mean value at line clipped by given firing strength
"""
mutable struct TriangularMF<:MF

	l_vertex::Real
	center::Real
	r_vertex::Real

	eval::Function
	mean_at::Function

	function TriangularMF(l_vertex::Real, center::Real, r_vertex::Real)

		if l_vertex <= center <= r_vertex

			this = new()

			this.l_vertex = l_vertex
			this.center = center
			this.r_vertex = r_vertex

			this.eval = function eval(x)
				maximum([minimum([((x - this.l_vertex) / (this.center - this.l_vertex)), ((this.r_vertex - x) / (this.r_vertex - this.center))]), 0])
			end

			this.mean_at = function mean_at(firing_strength)

				if firing_strength != 1
					p1 = (this.center - this.l_vertex) * firing_strength + this.l_vertex
					p2 = (this.center - this.r_vertex) * firing_strength + this.r_vertex
					(p1 + p2) / 2
				elseif firing_strength == 1
					return this.center
				end

			end

			this

		else

			error("invalid vertices")

		end

	end

end

"""
 Gaussian membership function type

	GaussianMF(center, sigma)

	 Properties
	 ----------
	 `center` is the center of the distribution
	 `sigma` determines width of the distribution

	 `eval` function returns membership value at a point
	 `mean_at` function returns mean value at line clipped by given firing strength

"""
mutable struct GaussianMF<:MF

	center::Real
	sigma::Real

	eval::Function
	mean_at::Function

	function GaussianMF(center::Real, sigma::Real)

		this = new()

		this.center = center
		this.sigma = sigma

		this.eval = function eval(x)
			exp( - 0.5 * ((x - this.center) / this.sigma) ^ 2)
		end

		this.mean_at = function mean_at(firing_strength)
			this.center
		end

		this

	end

end

"""
	 Generalised Bell membership function type

	BellMF(a, b, c)

	 Properties
	 ----------
	 `a`, `b` and `c` the usual bell parameters with `c` being the center

	 `eval` function returns membership value at a point
	 `mean_at` function returns mean value at line clipped by given firing strength

"""
mutable struct BellMF<:MF

	a::Real
	b::Real
	c::Real

	eval::Function
	mean_at::Function

	function BellMF(a::Real, b::Real, c::Real)

		this = new()

		this.a = a
		this.b = b
		this.c = c

		this.eval = function eval(x)
			1 / (1 + abs((x - this.c) / this.a) ^ (2 * this.b))
		end

		this.mean_at = function mean_at(firing_strength)
			this.c
		end

		this

	end

end

"""
	 Trapezoidal membership function type

	TrapezoidalMF(l_bottom_vertex, l_top_vertex, r_top_vertex, r_bottom_vertex)

	 Properties
	 ----------
	 `l_bottom_vertex`, `l_top_vertex`, `r_top_vertex` and `r_bottom_vertex` are the vertices of the trapezoid, in order

	 `eval` function returns membership value at a point
	 `mean_at` function returns mean value at line clipped by given firing strength

"""
mutable struct TrapezoidalMF<:MF

	l_bottom_vertex::Real
	l_top_vertex::Real
	r_top_vertex::Real
	r_bottom_vertex::Real

	eval::Function
	mean_at::Function

	function TrapezoidalMF(l_bottom_vertex::Real, l_top_vertex::Real, r_top_vertex::Real, r_bottom_vertex::Real)

		if l_bottom_vertex <= l_top_vertex <= r_top_vertex <= r_bottom_vertex

			this = new()

			this.l_bottom_vertex = l_bottom_vertex
			this.l_top_vertex = l_top_vertex
			this.r_top_vertex = r_top_vertex
			this.r_bottom_vertex = r_bottom_vertex

			this.eval = function eval(x)
				maximum([minimum([((x - this.l_bottom_vertex) / (this.l_top_vertex - this.l_bottom_vertex)), 1, ((this.r_bottom_vertex - x) / (this.r_bottom_vertex - this.r_top_vertex))]), 0])
			end

			this.mean_at = function mean_at(firing_strength)
				p1 = (this.l_top_vertex - this.l_bottom_vertex) * firing_strength + this.l_bottom_vertex
				p2 = (this.r_top_vertex - this.r_bottom_vertex) * firing_strength + this.r_bottom_vertex
				(p1 + p2) / 2
			end

			this

		else

			error("invalid vertices")

		end

	end

end

"""
	 Sigmoid membership function type

	SigmoidMF(a, c, limit)

	 Properties
	 ----------
	 `a` controls slope
	 `c` is the crossover point
	 `limit` sets the extreme limit

	 `eval` function returns membership value at a point
	 `mean_at` function returns mean value at line clipped by given firing strength

"""
mutable struct SigmoidMF<:MF

	a::Real
	c::Real
	limit::Real

	eval::Function
	mean_at::Function

	function SigmoidMF(a::Real, c::Real, limit::Real)

		if (a > 0 && limit > c) || (a < 0 && limit < c)

			this = new()

			this.a = a
			this.c = c
			this.limit = limit

			this.eval = function eval(x)
				1 / (1 + exp(-this.a * (x - this.c)))
			end

			this.mean_at = function mean_at(firing_strength)

				if firing_strength == 1
					p_firing_strength = 0.999
				elseif firing_strength == 0
					p_firing_strength = 0.001
				else
					p_firing_strength = firing_strength
				end

				p1 = -log((1 / p_firing_strength) - 1) / this.a + this.c
				p2 = this.limit
				(p1 + p2) / 2

			end

			this

		else

			error("invalid parameters")

		end

	end

end


text:
include("../server/simulation/simulation.jl")

using .Simulation

# Focus on synapse #3 for this test.

const SYNAPSE = 3
const TOTAL_SPIKES = 25

function test_merger_stream()
    println("Testing regular pattern stream")

    pat_stream = Simulation.RegularPatternStream()

    println(@__DIR__)
    # Load pattern stream
    file = "../../src/data/stim_2.data"

    open(file, "r") do f
        patterns = readlines(f)

        duration = length(patterns[1])
        syn_lanes = length(patterns)
        stimulus = zeros(UInt8, syn_lanes, duration)

        row = 1
        for pat_samples in patterns
            col = 1
            # Format:
            # ..|..|.....| etc
            # println(pat_samples)
            for c in pat_samples
                if c == '|'
                    stimulus[row, col] = UInt8(1)
                end
                col += 1
            end
            # println(samples[row, :])
            row += 1
        end

        Simulation.config_stream!(pat_stream, stimulus, 50)
    end

    # Create poisson stream to merge with.
    poi_stream = Simulation.PoissonStream(UInt64(13163), 0.2)

    merger = Simulation.StreamMerger()

    Simulation.add_stream!(merger, pat_stream)
    Simulation.add_stream!(merger, poi_stream)

    # 0000010010000001000000000000010000010000000000000000011100000000100000110000000000010010100010100010 = poisson
    # 0000100000000000000000001000000000000000000010000000000000000000100000000000000000001000000000000000 = pattern
    # 0000110010000001000000001000010000010000000010000000011100000000100000110000000000011010100010100010 = merge

    expected = [0,0,0,0,1,1,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,0,1,0,0,0,1,0,1,0,0,0,1,0];

    for t in 1:100
        @assert Simulation.output(merger) == expected[t] "Expected output did not match at t = ($t)"
        print(Simulation.output(merger))
        Simulation.step!(merger)
    end
    println("")

    println("Passed")
end

test_merger_stream()


text:
function []=PlotFaceKiteGrid(Grid,fig)
figure(fig)
clf(fig)
hold on;
for iF=1:size(Grid.Faces,2)
  for iK=1:size(Grid.Faces(iF).Kite.Faces,2)
    X=zeros(size(Grid.Faces(iF).Kite.Faces(iK).N,2),3);
    for i=1:size(Grid.Faces(iF).Kite.Faces(iK).N,2)
      X(i,:)=Grid.Faces(iF).Kite.Nodes(Grid.Faces(iF).Kite.Faces(iK).N(i)).P;
    end
    fill3(X(:,1),X(:,2),X(:,3),0);
  end
end
hold off
end


text:
struct MINLP <: GroundStructureSolver
    solver             # to solve MINLP model
    contz::Bool
    timelimit::Float64 # seconds
    writemodels::Bool

    function MINLP(solver; contz=false, timelimit=Inf, writemodels=false)
        new(solver, contz, timelimit, writemodels)
    end
end

"""
Build MINLP model for instance on given ground structure.

To be solved with given MPB solver (capable of NLP).
Use continuous variables 0 ≤ z ≤ 1 if `contz` is true.
"""
function make_minlp(inst::Instance, topo::Topology, optimizer; contz=false)
    nodes, nnodes = topo.nodes, length(topo.nodes)
    arcs, narcs = topo.arcs, length(topo.arcs)
    terms, nterms = inst.nodes, length(inst.nodes)
    termidx = termindex(nodes, terms)
    ndiams = length(inst.diameters)

    # demand for all nodes, including junctions
    dem = fill(0.0, nnodes)
    for t in 1:nterms
        dem[termidx[t]] = inst.demand[t]
    end

    # pressure bounds for all nodes
    πlb_min = minimum([b.lb for b in inst.pressure])^2
    πub_max = maximum([b.ub for b in inst.pressure])^2
    πlb = fill(πlb_min, nnodes)
    πub = fill(πub_max, nnodes)
    πlb[termidx] = [b.lb^2 for b in inst.pressure]
    πub[termidx] = [b.ub^2 for b in inst.pressure]

    # adjacency lists
    inarcs, outarcs = [Int[] for v in 1:nnodes], [Int[] for v in 1:nnodes]
    for a in 1:narcs
        tail, head = arcs[a]
        push!(inarcs[head], a)
        push!(outarcs[tail], a)
    end
    tail = [arcs[a].tail for a in 1:narcs]
    head = [arcs[a].head for a in 1:narcs]

    # "big-M" bound for flow on arcs
    qmax = 0.5 * sum(abs.(inst.demand))

    L = pipelengths(topo)
    c = [diam.cost for diam in inst.diameters]
    Dm5 = [diam.value^(-5) for diam in inst.diameters]
    C = inst.ploss_coeff * L

    # always use SCIP directly
    model = JuMP.direct_model(MOI.instantiate(optimizer))

    # select arcs from topology
    @variable(model, y[1:narcs], Bin)

    # select diameter on arcs
    if contz # continuous
        @variable(model, 0 ≤ z[1:narcs, 1:ndiams] ≤ 1)
    else # discrete
        @variable(model, z[1:narcs, 1:ndiams], Bin)
    end

    # flow through arcs
    @variable(model, 0 <= q[1:narcs] <= qmax)

    # squared pressure on nodes
    @variable(model, πlb[v] ≤ π[v=1:nnodes] ≤ πub[v])

    # mass flow balance at nodes
    @constraint(model, flowbalance[v=1:nnodes],
                sum(q[a] for a=inarcs[v]) - sum(q[a] for a=outarcs[v]) == dem[v])

    # allow flow only for active arcs
    @constraint(model, active[a=1:narcs], q[a] <= qmax*y[a])

    # (conditional) pressure loss constraint (awkward nonlinear formulation)
    @NLconstraint(model, ploss[a=1:narcs],
                  y[a]*(π[tail[a]] - π[head[a]]) ==
                  C[a]*sum(Dm5[i]*z[a,i] for i=1:ndiams)*q[a]^2 )

    # choose diameter on active arcs
    @constraint(model, choice[a=1:narcs], sum(z[a,i] for i=1:ndiams) == y[a])

    # minimize total construction cost
    @objective(model, Min, sum(c[i] * L[a] * z[a,i] for a=1:narcs for i=1:ndiams))

    return model, y, z, q, π
end

function PipeLayout.optimize(inst::Instance, topo::Topology, solver::MINLP)
    if solver.timelimit < Inf
        MOI.set(solver.solver, MOI.TimeLimitSec, solver.timelimit)
    end
    model, y, z, q, π = make_minlp(inst, topo, solver.solver, contz=solver.contz)
    JuMP.optimize!(model)
    status = JuMP.termination_status(model)

    if solver.writemodels
        MOI.writeproblem(internalmodel(model), "minlp.orig.cip")
    end

    if status == MOI.INFEASIBLE
        return Result(status, nothing, Inf, Inf, 0)
    end

    zsol = round.(Bool, JuMP.value.(z) .>= 0.5)
    qsol = JuMP.value.(q)
    bestsol = CandSol(zsol, qsol, qsol.^2)

    primal = JuMP.objective_value(model)
    dual = JuMP.objective_bound(model)

    Result(status, bestsol, primal, dual, 0)
end


text:
using FactCheck, Base.Test, Images, Colors, FixedPointNumbers

facts("Corner") do
    A = zeros(41,41)
    A[16:26,16:26] = 1

    context("Corners API") do
    	img = zeros(20, 20)
    	img[4:17, 4:17] = 1
    	img[8:13, 8:13] = 0

    	expected_corners = zeros(20, 20)
    	ids = map(CartesianIndex{2}, [(4, 4), (4, 17), (17, 4), (17, 17), (8, 8), (8, 13), (13, 8), (13, 13)])
    	for id in ids expected_corners[id] = 1 end
    	corners = imcorner(img, method = harris)
    	for id in ids @fact corners[id] --> true end
    	@fact sum(corners .!= expected_corners) < 3 --> true
    	corners = imcorner(img, method = shi_tomasi)
    	for id in ids @fact corners[id] --> true end
    	@fact sum(corners .!= expected_corners) < 3 --> true

    	ids = map(CartesianIndex{2}, [(4, 4), (4, 17), (17, 4), (17, 17)])
    	corners = imcorner(img, method = kitchen_rosenfeld)
    	for id in ids @fact corners[id] --> true end
    	ids = map(CartesianIndex{2}, [(8, 8), (8, 13), (13, 8), (13, 13)])
    	for id in ids expected_corners[id] = 0 end
    	@fact sum(corners .!= expected_corners) < 3 --> true
    end

	context("Harris") do
		Ac = imcorner(A, 0.99, true, method = harris)
		# check corners
		@fact Ac[16,16] --> true
		@fact Ac[16,26] --> true
		@fact Ac[26,16] --> true
		@fact Ac[26,26] --> true

		@fact all(Ac[1:14, :]) --> false
		@fact all(Ac[:, 1:14]) --> false
		@fact all(Ac[28:end, :]) --> false
		@fact all(Ac[:, 28:end]) --> false
		@fact all(Ac[18:24, 18:24]) --> false

		Ac = harris(A)

		@fact Ac[16,16] --> roughly(maximum(Ac))
		@fact Ac[16,26] --> roughly(Ac[16,16])
		@fact Ac[26,16] --> roughly(Ac[16,16])
		@fact Ac[26,26] --> roughly(Ac[16,16])

		# check edge intensity
		@fact Ac[16,21] --> roughly(minimum(Ac))
		@fact Ac[21,16] --> roughly(Ac[16,21])
		@fact Ac[21,26] --> roughly(Ac[16,21])
		@fact Ac[26,21] --> roughly(Ac[16,21])

		# check background intensity
		@fact Ac[5,5] --> roughly(0)
		@fact Ac[5,36] --> roughly(0)
		@fact Ac[36,5] --> roughly(0)
		@fact Ac[36,36] --> roughly(0)

	end

	context("Shi-Tomasi") do
		Ac = imcorner(A, 0.99, true, method = shi_tomasi)
		# check corners
		@fact Ac[16,16] --> true
		@fact Ac[16,26] --> true
		@fact Ac[26,16] --> true
		@fact Ac[26,26] --> true

		@fact all(Ac[1:14, :]) --> false
		@fact all(Ac[:, 1:14]) --> false
		@fact all(Ac[28:end, :]) --> false
		@fact all(Ac[:, 28:end]) --> false
		@fact all(Ac[18:24, 18:24]) --> false

		Ac = shi_tomasi(A)
		# check corner intensity
		@fact Ac[16,16] --> roughly(maximum(Ac))
		@fact Ac[16,26] --> roughly(Ac[16,16])
		@fact Ac[26,16] --> roughly(Ac[16,16])
		@fact Ac[26,26] --> roughly(Ac[16,16])

		# check edge intensity
		@fact Ac[16,21] --> roughly(0)
		@fact Ac[21,16] --> roughly(0)
		@fact Ac[21,26] --> roughly(0)
		@fact Ac[26,21] --> roughly(0)

		# check background intensity
		@fact Ac[5,5] --> roughly(0)
		@fact Ac[5,36] --> roughly(0)
		@fact Ac[36,5] --> roughly(0)
		@fact Ac[36,36] --> roughly(0)

	end

	context("Kitchen-Rosenfeld") do
		A[10:30, 10:30] = 1
		A[15:25, 15:25] = 0
		Ac = imcorner(A, 0.99, true, method = kitchen_rosenfeld)
		@fact Ac[10, 10] --> true
		@fact Ac[10, 30] --> true
		@fact Ac[30, 10] --> true
		@fact Ac[30, 30] --> true

		@fact all(Ac[1:9, :]) --> false
		@fact all(Ac[:, 1:9]) --> false
		@fact all(Ac[31:end, :]) --> false
		@fact all(Ac[:, 31:end]) --> false
		@fact all(Ac[12:28, 12:28]) --> false
	end

	context("Fast Corners") do

		img = reshape(1:1:49, 7, 7)

		for i in 1:8
			corners = fastcorners(img, i)
			@fact all(corners) --> true
		end

		corners = fastcorners(img, 9)
		@fact all(corners[4, 1:end - 1]) --> false
		corners[4, 1:end - 1] = true
		@fact all(corners) --> true

		corners = fastcorners(img, 10)
		@fact all(corners[3:4, 1:end - 2]) --> false
		@fact corners[4, end - 1] --> false
		corners[3:4, 1:end - 2] = true
		corners[4, end - 1] = true
		@fact all(corners) --> true

		corners = fastcorners(img, 11)
		@fact all(corners[2:5, 1:end - 3]) --> false
		@fact all(corners[3:5, end - 2]) --> false
		@fact corners[4, end - 1] --> false
		corners[2:5, 1:end - 3] = true
		corners[3:5, 1:end - 2] = true
		corners[4, end - 1] = true
		@fact all(corners) --> true

		corners = fastcorners(img, 12)
		@fact all(corners[1:6, 1:end - 3]) --> false
		@fact all(corners[3:5, end - 2]) --> false
		@fact corners[4, end - 1] --> false
		corners[1:6, 1:end - 3] = true
		corners[3:5, 1:end - 2] = true
		corners[4, end - 1] = true
		@fact all(corners) --> true

		img = gaussian2d(1.4)
		img = vcat(img, img)
		img = hcat(img, img)
		corners = fastcorners(img, 12, 0.05)

		@fact corners[5, 5] --> true
		@fact corners[5, 14] --> true
		@fact corners[14, 5] --> true
		@fact corners[14, 14] --> true

		@fact all(corners[:, 1:3] .== false) --> true
		@fact all(corners[1:3, :] .== false) --> true
		@fact all(corners[:, 16:end] .== false) --> true
		@fact all(corners[16:end, :] .== false) --> true
		@fact all(corners[6:12, 6:12] .== false) --> true

		img = 1 - img

		corners = fastcorners(img, 12, 0.05)

		@fact corners[5, 5] --> true
		@fact corners[5, 14] --> true
		@fact corners[14, 5] --> true
		@fact corners[14, 14] --> true

		@fact all(corners[:, 1:3] .== false) --> true
		@fact all(corners[1:3, :] .== false) --> true
		@fact all(corners[:, 16:end] .== false) --> true
		@fact all(corners[16:end, :] .== false) --> true
		@fact all(corners[6:12, 6:12] .== false) --> true

	end

end


text:
export Sequential

"""
$(TYPEDEF)

Represents a sequential (feed-forward) neural net, with `layers` ordered from input
to output.

## Fields:
$(FIELDS)
"""
@auto_hash_equals struct Sequential <: NeuralNet
    layers::Array{Layer, 1}
    UUID::String
end

function Base.show(io::IO, p::Sequential)

    println(io, "sequential net $(p.UUID)")
    for (index, value) in enumerate(p.layers)
        println(io, "  ($index) $value")
    end
end

(p::Sequential)(x::Array{<:JuMPReal}) = x |> p.layers


text:
__precompile__()
module iteratedintegral
export init, get_ito, get_strat
#################################################################
function init(fcn, opt)
  #=
  fcn=SDE definition (for dimension m and noise type)
  opt= options (for max_p)
  =#
  S=Dict("KL_C"=>1, "m"=>fcn.m, "NoiseType"=>fcn.noise_type)
  try
    S["max_p"]=opt["max_p"]
  catch
    S["max_p"]=1000
  end
  # get tail of sum_{p=P]^infty 1/p^2
  no_terms=100
  H2_tail=Array{Float16}(undef,no_terms)
  H2_tail[1]=pi^2/6 # zeta(2)
  for i=2:no_terms
    H2_tail[i]=H2_tail[i-1]-1/(i^2)
  end
  #S["H2_tail"]=H2_tail
  # map for upper-triangular to vector
  T=Int16
  m=fcn.m
  m2=round(T,m*(m-1)/2)
  ri=Array{T}(undef,m2)
  ci=Array{T}(undef,m2)
  k=1
  for i = 1:(m-1)
      for j = (i+1):m
        ri[k]=i
        ci[k]=j
        k=k+1
      end
  end
  S=merge(S,Dict("ri"=>ri, "ci"=>ci, "m2"=>m2))
  # initialise storage
  S["tmp_m2"]=Array{Float64}(undef,(m,m))
  S["tail"]=Array{Float64}(undef,(m,m))
  S["J"]=Array{Float64}(undef,m)
  S["K"]=Array{Float64}(undef,(m,m))
  #
  return S
end
#####
function __init__()
end
######################
function get_p(dt,S)
  #=
  get number of terms for a KL expansion in terms of dt
  and KL_C
  (not actually used as tail correction is implemented here)
  =#
  thresh=2*pi^2*S["KL_C"]^2 * dt
  p=findfirst(x->x<thresh,S["H2_tail"])
  if p<S["max_p"]
      p= 1 / (2 * (pi* S["KL_C"])^2 * dt)
  end
  return round(Int,p)
end
####################################
function Gamma_1!(out,b,S)
  #=
  helper function: see notes \Gamma_1
  in: b m2-vector
      S dictionary (see init function for defs of m and m2 and ri,ci)
  out: out (m x m)-matrix
  =#
  r=S["ri"]
  c=S["ci"]
  for i in 1:length(b)
    out[r[i],c[i]]=b[i]
    out[c[i],r[i]]=-b[i]
  end
end
######################################
function Gamma_2!(out,A,dW,S)
  #=
  helper function: see notes
  in: A (m x m)-matrix
  dW m-vector
  S dictionary (see init function)
  out: out (m x m)-matrix
  =#
  for i=1:length(dW)
    for j=1:length(dW)
      out[i,j]=(A[i,:]'*dW)*dW[j]
    end
  end
end
########################################
function Gamma_3!(out,A,S)
  #=
  helper function: see notes
  in: A (m x m)-matrix
  S dictionary (see init function)
  out: m2-vector
  =#
  r=S["ri"]
  c=S["ci"]
  for i in 1:length(r)
    out[i]=A[r[i],c[i]] - A[c[i],r[i]]
  end
end
####################################
function compute_tail!(tail,dt,dW,ap,S)
  #=
  Tail approximation for KL expansion of stochastic integrals
  in: dt timestep
  dW m-vector
  ap coefficient (see notes)
  S dictionary (see init function)
  out: tail (m x m)-matrix
  =#
  scale=dt^2*ap/(2*pi^2)
  alpha=sqrt(1+sum(dW.*dW))
  noiseVec=randn(S["m2"])*sqrt(scale)
  tmp_M=noiseVec/((1+alpha)*dt)
  tmp_m2=S["tmp_m2"]
  Gamma_1!(tmp_m2,tmp_M,S)
  Gamma_2!(tail,tmp_m2,dW,S)
  Gamma_3!(tmp_M,tail,S)
  Gamma_1!(tail,tmp_M+noiseVec,S)
  for i = 1:S["m"]
    tail[i,i] = 0
  end
end
#
function get_strat!(dt,S)
  #=
  Main routine for generating integrals (Stratonvich)
  in: dt timestep
  S (see init function)
  out: J,K updated in S
  =#
  m=S["m"]
  J=S["J"]
  K=S["K"]
  if S["NoiseType"]==2 # diagonal
    J[:]=sqrt(dt)*randn(m)
    fill!(K,0)
    for i=1:m
      K[i,i]=0.5*J[i]^2
    end
  elseif S["NoiseType"]==3 # commutative
    J[:]=sqrt(dt)*randn(m)
    for i=1:m
      for j=1:m
       K[i,j]=0.5*J[i]*J[j]
    end
  end
  else # unstructured
    BBT!(dt,S)
  end
end
######################################
function get_ito!(dt,S)
  #=
  Main routine for generating integrals (Ito)
  in: dt timestep
  S (see init function)
  out: J,K updated in S
  =#
  get_strat!(dt,S)
  K=S["K"]
  for i=1:S["m"]
    K[i,i]-=0.5*dt
  end
end
########################################
function BBT!(dt,S)
  #=
  Routine for generating integrals (Stratonovich) unstructured case
  in: dt timestpe
  S dictionary (see init function)
  out: J,K updated in S
  =#
  m=S["m"]
  J=S["J"]
  K=S["K"]
  #
  p_const=sqrt(m*(m-1)/(24*dt)/(pi*S["KL_C"]))
  a_p=pi^2/6
  J[:]=sqrt(dt)*randn(m)
  for i=1:m
    K[i,i]=0.5*J[i]^2
  end
  p=round(p_const*sqrt(m*4*sum(J.*J)/dt))
  for r=1:(p+1)
    a_p=a_p-(1/r^2)
    sd_sq=dt/(2*pi*r)
    sd=sqrt(sd_sq)
    eta=randn(m)*sd
    zeta=randn(m)*sd
    K[:,1]-= zeta*(2/sqrt(pi*r))
    for i=2:m
      for j=1:i-1
        K[i,j]+=zeta[i]*eta[j]-eta[i]*zeta[j]
      end
    end
  end
  for i=2:m
    for j=1:i-1
      tmp1=K[i,j]-0.5*(K[j,1]*J[i]-K[i,1]*J[j])
      tmp2=0.5*J[i]*J[j]
      K[i,j]=tmp2+tmp1
      K[j,i]=tmp2-tmp1
    end
  end
  tail=S["tail"]
  compute_tail!(tail,dt,J,a_p,S)
  K+=tail
end
end


text:
@testset "utils/utils.jl" begin
    @test AGP.jitt isa AGP.Jittering
    @test Float64(AGP.jitt) ≈ 1e-4
    @test Float32(AGP.jitt) ≈ 1e-3
    @test Float16(AGP.jitt) ≈ 1e-2

    @test AGP.δ(0, 1) == 0.0
    @test AGP.δ(1, 1) == 1.0

    A = rand(2, 2)
    B = rand(2, 2)
    x = rand(2)
    @test AGP.hadamard(A, B) == A .* B

    X = copy(A)
    AGP.add_transpose!(X)
    @test X == A + A'

    D = A * A' + I
    C = cholesky(D)
    @test AGP.invquad(C, x) ≈ dot(x, D \ x)

    @test AGP.trace_ABt(A, B) ≈ tr(A * B')
    @test AGP.diag_ABt(A, B) ≈ diag(A * B')
    @test AGP.diagv_B(x, B) ≈ Diagonal(x) * B
    @test AGP.κdiagθκ(A, x) ≈ A' * Diagonal(x) * A
    @test AGP.ρκdiagθκ(2.0, A, x) ≈ 2.0 * A' * Diagonal(x) * A
    @test AGP.opt_add_diag_mat(x, A) ≈ A + Diagonal(x)
    @test AGP.safe_expcosh(2.0, 1.0) ≈ exp(2.0) / cosh(1.0)
    @test AGP.logcosh(2.0) ≈ log(cosh(2.0))

    #TODO test
    #=
    function symcat(S::Symmetric, v::AbstractVector, vv::Real)
        S = vcat(S, v')
        S = hcat(S, vcat(v, vv))
        return Symmetric(S)
    end

    function make_grid(range1, range2)
        return hcat(
            [i for i in range1, j in range2][:],
            [j for i in range1, j in range2][:],
        )
    end
    =#
    @test (2.0 * C).L ≈ cholesky(2.0 * D).L
    @test C * x ≈ D * x
    @test (C + 2.0 * I).L ≈ cholesky(D + 2.0 * I).L
end