text:
# EXCLUDE FROM TESTING
# this file doesn't have an entry point, see `verify.jl` instead

# Fast parallel reduction for Kepler hardware
# - uses shuffle and shared memory to reduce efficiently
# - support for large arrays
#
# Based on devblogs.nvidia.com/parallelforall/faster-parallel-reductions-kepler/

using CUDAdrv, CUDAnative

#
# Main implementation
#

# Reduce a value across a warp
@inline function reduce_warp(op::Function, val::T)::T where {T}
  offset = CUDAnative.warpsize() ÷ UInt32(2)
  # TODO: this can be unrolled if warpsize is known...
  while offset > 0
    val = op(val, shfl_down(val, offset))
    offset ÷= UInt32(2)
  end
  return val
end

# Reduce a value across a block, using shared memory for communication
@inline function reduce_block(op::Function, val_x::T,val_y::T,val_z::T)::T where {T}
  # shared mem for 32 partial sums
  shared = @cuStaticSharedMem(T, 32)

  # TODO: use fldmod1 (JuliaGPU/CUDAnative.jl#28)
  wid_x  = div(threadIdx().x-UInt32(1), CUDAnative.warpsize()) + UInt32(1)
  wid_y  = div(threadIdx().y-UInt32(1), CUDAnative.warpsize()) + UInt32(1)
  wid_z  = div(threadIdx().z-UInt32(1), CUDAnative.warpsize()) + UInt32(1)
  lane_x = rem(threadIdx().x-UInt32(1), CUDAnative.warpsize()) + UInt32(1)
  lane_y = rem(threadIdx().y-UInt32(1), CUDAnative.warpsize()) + UInt32(1)
  lane_z = rem(threadIdx().z-UInt32(1), CUDAnative.warpsize()) + UInt32(1)

  # each warp performs partial reduction
  val_x = reduce_warp(op, val_x)
  val_y = reduce_warp(op, val_y)
  val_z = reduce_warp(op, val_z)

  # write reduced value to shared memory
  if lane_x == 1
    @inbounds shared[wid_x] = val_x
  end
  if lane_y == 1
    @inbounds shared[wid_y] = val_y
  end
  if lane_z == 1
    @inbounds shared[wid_z] = val_z
  end
  # wait for all partial reductions
  sync_threads()

  # read from shared memory only if that warp existed
  @inbounds val_x = (threadIdx().x <= fld(blockDim().x, CUDAnative.warpsize())) ? shared[lane] : zero(T)
  @inbounds val_y = (threadIdx().y <= fld(blockDim().y, CUDAnative.warpsize())) ? shared[lane] : zero(T)
  @inbounds val_z = (threadIdx().z <= fld(blockDim().z, CUDAnative.warpsize())) ? shared[lane] : zero(T)

  # final reduce within first warp
  if wid_x == 1
    val_x = reduce_warp(op, val_x)
  end
  if wid_y == 1
    val_y = reduce_warp(op, val_x)
  end
  if wid_z == 1
    val_z = reduce_warp(op, val_x)
  end

  return val_x+valy_y+val_z
end

# Reduce an array across a complete grid
function reduce_grid(op::Function, input::CuDeviceVector{T}, output::CuDeviceVector{T},
  len::CuDeviceVector{Integer}) where {T}

  # TODO: neutral element depends on the operator (see Base's 2 and 3 argument `reduce`)
  val = zero(T)

  # reduce multiple elements per thread (grid-stride loop)
  # TODO: step range (see JuliaGPU/CUDAnative.jl#12)
  i_x = (blockIdx().x-UInt32(1)) * blockDim().x + threadIdx().x
  i_y = (blockIdx().y-UInt32(1)) * blockDim().y + threadIdx().y
  i_z = (blockIdx().z-UInt32(1)) * blockDim().z + threadIdx().z
  step_x = blockDim().x * gridDim().x
  step_y = blockDim().y * gridDim().y
  step_z = blockDim().z * gridDim().z
  while i_x <= len[1]
    while i_y <= len[2]
      while i_z <= len[3]
        @inbounds val = op(val, input[i])
        i_z+=step_z
      end
      i_y += step_y
    end
    i_x += step_x
  end

  val = reduce_block(op, val)

  if threadIdx().x == UInt32(1)
    @inbounds output[blockIdx().x] = val
  end

  return
end

"""
Reduce a large array.
Kepler-specific implementation, ie. you need sm_30 or higher to run this code.
"""
function gpu_reduce(op::Function, input::CuVector{T}, output::CuVector{T}) where {T}
  len = length(input)

  # TODO: these values are hardware-dependent, with recent GPUs supporting more threads
  threads = 512
  blocks = min((len + threads - 1) ÷ threads, 1024)

  # the output array must have a size equal to or larger than the number of thread blocks
  # in the grid because each block writes to a unique location within the array.
  if length(output) < blocks
    throw(ArgumentError("output array too small, should be at least $blocks elements"))
  end

  @cuda (blocks,threads) reduce_grid(op, input, output, Int32(len))
  @cuda (1,1024) reduce_grid(op, output, output, Int32(blocks))

  return
end
ctx = CuCurrentContext()
dev = device(ctx)
if capability(dev) < v"3.0"
  warn("this example requires a newer GPU")
  exit(0)
end

len = 10^7
input = ones(Int32, len)

# CPU
cpu_val = reduce(+, input)

# CUDAnative
let
  gpu_input = CuArray(input)
  gpu_output = similar(gpu_input)
  gpu_reduce(+, gpu_input, gpu_output)
  gpu_val = Array(gpu_output)[1]
  @assert cpu_val == gpu_val
end


text:
using NPZ
using CSV
using JSON
using Plots
using DataFrames
using Statistics
using CartesianGeneticProgramming

include("utils.jl")

struct Month
    topo0
    topo1
    tp
    hs
    direction
end

gene_json = JSON.parsefile("/home/mn/JuliaProjects/cgp-morphcast/gens/2021-05-19T15:31:14.968/23000/0010.dna")
fitness = gene_json["fitness"]
functions = gene_json["fs"]
chromosome = convert(Array{Float64,1}, gene_json["chromosome"])
xs = gene_json["xs"]
ys = gene_json["ys"]

ind = CGPInd(get_config("topo.yaml"), chromosome)

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

dataset_path = get_dataset_path()
n_months = 236
months = []

for i in 1:n_months
    topo_path = string(dataset_path, '/', i, "_topo.npy")
    topo = npzread(topo_path)
    topo0 = norm_to_range(topo[1, :], -10.0, 10.0, -1, 1)
    topo1 = norm_to_range(topo[2, :], -10.0, 10.0, -1, 1)
    d_topo = topo1 - topo0

    @assert topo0 + d_topo == topo1

    waves_path = string(dataset_path, '/', i, "_waves.csv")
    waves_csv = CSV.read(waves_path, DataFrame)

    tp_month = norm_to_range(waves_csv.tp, 0, 25, -1, 1)
    hs_month = norm_to_range(waves_csv.hs, 0, 10, -1, 1)
    direction_month = norm_to_range(waves_csv.direction, 0, 360, -1, 1)

    input_row = vcat(topo0,
        mean(tp_month), std(tp_month), maximum(tp_month), median(tp_month),
        mean(hs_month), std(hs_month), maximum(hs_month), median(hs_month),
        mean(direction_month), std(direction_month), maximum(direction_month), median(direction_month))
    pred_delta_topo = process(ind, input_row)

    println("final plot:")

    plot_d_topo = plot(d_topo, label="Target Δtopo")
    plot!(plot_d_topo, pred_delta_topo, label="Predicted Δtopo")

    plot_target = plot(norm_to_range(topo0 + pred_delta_topo, -1, 1, -10.0, 10.0))

    display(plot(plot_d_topo, plot_target, layout=(2, 1)))
    readline()
end

ext: jl

max_stars_repo_name: mahmoud-al-najar/cgp-morphcast

max_issues_repo_name: mahmoud-al-najar/cgp-morphcast

max_forks_repo_name: mahmoud-al-najar/cgp-morphcast


text:
using Test, Sym

@testset "tracing" begin
end

@testset "match" begin
end

@testset "simplification" begin
end


text:
using Base.Test
import Deprecations
import Tokenize.Tokens
import CSTParser: KEYWORD, OPERATOR, LITERAL
##########
# isexpr #
##########
let
    expr = CSTParser.parse("if true elseif false end")
    elseif_arg = expr.args[4]
    @test Deprecations.isexpr(elseif_arg, KEYWORD, Tokens.ELSEIF)

    expr = CSTParser.parse("true")
    @test Deprecations.isexpr(expr, LITERAL, Tokens.TRUE)

    expr = CSTParser.parse("+")
    @test Deprecations.isexpr(expr, OPERATOR, Tokens.PLUS)
end

let
    x = CSTParser.parse("foo")
    y = CSTParser.parse(":foo")
    @test !Deprecations.matches_template(x, y)
end

let
    x = CSTParser.parse("foo")
    y = CSTParser.parse("foo")
    @test Deprecations.matches_template(x, y)
end

let
    x = CSTParser.parse("\$name")
    @test Deprecations.is_template_expr(x) == (true, :name, false)

    x = CSTParser.parse("using name...")
    @test Deprecations.is_template_expr(x) == (false, nothing, false)

    x = CSTParser.parse("\$BODY...")
    @test Deprecations.is_template_expr(x) == (true, :BODY, true)

    x = CSTParser.parse("::Type{\$NAME{\$T...}}")
    @test Deprecations.is_template_expr(x) == (false, nothing, false)
end

# VERSIONCOMPARE
let
    x = CSTParser.parse("VERSION")
    y = CSTParser.parse("v\"0.2.0\"")
    @test Deprecations.detect_ver_arguments(x, y) == v"0.2.0"

    x = CSTParser.parse("foo")
    @test Deprecations.is_identifier(x, "foo") == true
end

# OLDCOMPAT
let
    x = CSTParser.parse("f(x::T) where{T}")
    @test Deprecations.is_where_expr(x)
end


text:
# CG preconditioners must be symmetric and positive definite.
# Multigrid cycles will satisfy these requirements for V- and W-cycles if
#   * Restriction is the adjoint of prolongation (up to scaling)
#   * Postsmoothing is performed in reverse order of presmoothing
#   * The coarse solver is either exact or symmetric and positive definite

__DEFAULT_MGPCG_DEPTH(A) = max(2, 1 + floor(Int, log2(minimum(size(A)))-2))

__DEFAULT_MGPCG_PRE(nlevels::Integer) = ntuple(nlevels-1) do l
    Smoothers(
        RedBlackGaussSeidel(ForwardSweep(), iter=1),
        RedBlackGaussSeidel(BoundarySweep(ForwardSweep()), iter=2^l),
    )
end

__DEFAULT_MGPCG_COARSE() = Smoothers(
    RedBlackGaussSeidel(ForwardSweep(), iter=512),
    RedBlackGaussSeidel(BackwardSweep(), iter=512),
)

__DEFAULT_MGPCG_POST(nlevels::Integer) = ntuple(nlevels-1) do l
    Smoothers(
        RedBlackGaussSeidel(BoundarySweep(BackwardSweep()), iter=2^l),
        RedBlackGaussSeidel(BackwardSweep(), iter=1),
    )
end

function mgpcg(
    A::Poisson2{<:AbstractFloat, 3},
    b::AbstractArray{<:AbstractFloat, N};
    kwargs...
) where {N}
    N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
    return mgpcg!(tzero(b), A, b; kwargs...)
end

function mgpcg!(
    x::AbstractArray{<:AbstractFloat, N},
    A::Poisson2{<:AbstractFloat, 3},
    b::AbstractArray{T};
    maxlevel::Integer = __DEFAULT_MGPCG_DEPTH(A),
    presmoother = __DEFAULT_MGPCG_PRE(maxlevel),
    coarsesolver = __DEFAULT_MGPCG_COARSE(),
    postsmoother = __DEFAULT_MGPCG_POST(maxlevel),
    cycle::Cycle = V(),
    ncycles::Integer = 1,
    atol::Real = zero(T),
    rtol::Real = sqrt(eps(T)),
    maxit::Integer = floor(Int, sqrt(length(b))),
    verbose::Bool = false
) where {T<:AbstractFloat, N}
    N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))

    size(x) == size(b) || throw(DimensionMismatch())
    size(A) == size(b)[1:3] || throw(DimensionMismatch())

    M = Multigrid(T, A,
        coarsesolver = coarsesolver,
        presmoother = presmoother,
        postsmoother = postsmoother,
        maxlevel = maxlevel
    )

    q = tzero(b, size(A))
    p = similar(b, size(A))

    M.workspace.x[1] = q

    if N == 3
        M.workspace.b[1] = tcopy(b)
        x = _mgpcg!(x, M, p, cycle, ncycles, atol, rtol, maxit, verbose)

    elseif N == 4
        xt = similar(x, size(A))
        M.workspace.b[1] = similar(b, size(A))

        for t in axes(b, 4)
            _tcopyto!(xt, @view(x[:,:,:,t]))
            _tcopyto!(M.workspace.b[1], @view(b[:,:,:,t]))

            xt = _mgpcg!(xt, M, p, cycle, ncycles, atol, rtol, maxit, verbose)

            _tcopyto!(@view(x[:,:,:,t]), xt)
        end
    end

    return x
end

function _mgpcg!(
    x::AbstractArray{<:AbstractFloat, 3},
    M::Multigrid,
    p::AbstractArray{T, 3},
    cycle::Cycle,
    ncycles::Integer,
    atol::Real,
    rtol::Real,
    maxit::Integer,
    verbose::Bool
) where {T<:AbstractFloat}
    _zero = zero(T)
    ρ = one(T)

    A = M.grids[1].A
    q = M.workspace.x[1]
    r = M.workspace.b[1]
    p = tfill!(p, _zero)

    # q = A * x
    # r -= q
    # res = norm(r)
    r, q, res = _init_r!(r, q, A, x)

    tol = max(rtol * res, atol)

    if verbose
        @printf("==== MGPCG ====\n")
        @printf("Depth: %d\n", length(M.grids))
        @printf("Tolerance: %1.4e\n", tol)
        @printf("%8s %3s %9s\n", "iter", " ", "resnorm")
    end

    if res ≤ tol
        return x
    end

    for i in 1:maxit
        # q = M\r
        tfill!(q, _zero, A)
        for _ in 1:ncycles
            cycle!(M, cycle)
        end

        # ρ = q⋅r
        ρ_prev = ρ
        ρ = _compute_ρ(q, r, A)

        # p = q + β * p
        β = ρ / ρ_prev
        p = _compute_p!(p, q, β, A)

        # q = Ap
        # α = q⋅r / p⋅Ap
        q, pAp = _compute_q!(q, A, p)
        α = ρ / pAp

        # x = x + α * p
        # r = r - α * q
        # res = norm(r)
        x, r, res = _compute_x_r!(x, r, p, q, α, A)

        if verbose
            @printf("[%5d/%d]  %1.4e\n", i, maxit, res)
        end

        if res ≤ tol
            break
        end
    end

    return x
end

function _compute_ρ(q::AbstractArray{T}, r, A::Poisson2) where {T}
    # ρ = q⋅r
    @inbounds @batch threadlocal=zero(T)::T for I in A.RI
        threadlocal = muladd(q[I], r[I], threadlocal)
    end
    return sum(threadlocal::Vector{T})
end

function _compute_p!(p, q, β, A::Poisson2)
    # p = q + β * p
    @inbounds @batch for I in A.RI
        p[I] *= β
        p[I] += q[I]
    end
    return p
end

function _compute_q!(q::AbstractArray{T}, A::Poisson2, p) where {T}
    # q = A * p
    # α = p*A*p
    D = A.D
    idx2 = convert(T, A.idx2[1])
    idy2 = convert(T, A.idx2[2])
    idz2 = convert(T, A.idx2[3])

    @inbounds @batch per=thread threadlocal=zero(T)::T for (I, J, K) in A.R8
        for k in K, j in J, i in I
            z1 = p[i,j,k-1]
            y1 = p[i,j-1,k]
            x1 = p[i-1,j,k]
            x0 = p[i,j,k]
            x2 = p[i+1,j,k]
            y2 = p[i,j+1,k]
            z2 = p[i,j,k+1]

            if A.interior[i,j,k]
                Δz = z1 + z2
                Δy = y1 + y2
                Δx = x1 + x2

                Δ = D * x0
                Δ = muladd(idz2, Δz, Δ)
                Δ = muladd(idy2, Δy, Δ)
                Δ = muladd(idx2, Δx, Δ)

                q[i,j,k] = Δ
                threadlocal = muladd(x0, Δ, threadlocal)
            end
        end
    end

    return q, sum(threadlocal::Vector{T})
end

function _compute_x_r!(x::AbstractArray{T}, r, p, q, α, A::Poisson2) where {T}
    # x = x + α * p
    # r = r - α * q
    # res = norm(r)
    @inbounds @batch threadlocal=zero(T)::T for I in A.RI
        x[I] = muladd( α, p[I], x[I])
        r[I] = muladd(-α, q[I], r[I])
        threadlocal = muladd(r[I], r[I], threadlocal)
    end
    return x, r, sqrt(sum(threadlocal::Vector{T}))
end

function _init_r!(r::AbstractArray{T}, q, A, x) where {T}
    # q = A * x
    # r -= q
    # res = norm(r)
    q = mul!(q, A, x)
    @inbounds @batch threadlocal=zero(T)::T for I in eachindex(r)
        r[I] -= q[I]
        threadlocal = muladd(r[I], r[I], threadlocal)
    end
    return r, q, sqrt(sum(threadlocal::Vector{T}))
end


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

JLLWrappers.@generate_main_file_header("SuperLU")
JLLWrappers.@generate_main_file("SuperLU", UUID("f391fee2-5786-5cc9-9662-0e362ec949e2"))
end  # module SuperLU_jll


text:
using Documenter, myexample

makedocs(;
    modules=[myexample],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/sumidu/myexample.jl/blob/{commit}{path}#L{line}",
    sitename="myexample.jl",
    authors="andré Calero Valdez",
    assets=String[],
)

deploydocs(;
    repo="github.com/sumidu/myexample.jl",
)


text:
using NodeJS

# Wrap in useless try/catch to avoid namespace pollution
try
    isdev = basename(dirname(dirname(dirname(@__FILE__)))) == "dev"
    if haskey(ENV, "CI") || isdev
        cd(joinpath(dirname(@__FILE__), "..", "packages")) do
            npm = NodeJS.npm_cmd()

            run(`$npm install --scripts-prepend-node-path=auto --unsafe-perm`)
            if haskey(ENV, "WEBIO_WEBPACK_ARGS")
                args = [ENV["WEBIO_WEBPACK_ARGS"]]
            else
                args =[]
            end
            run(`$npm run build-prod --scripts-prepend-node-path=auto --unsafe-perm -- $args`)
        end
    else
        @warn(
            "Can't build WebIO JS when not checkout out for development. "
            * "Try running Pkg.dev(\"WebIO\") if you want to rebuild JS."
        )
    end
catch
    rethrow()
end


text:
using FFTW

# Algorithms to choose
abstract type Alg end
"Build SH diversity using the quadratic functions" struct Quadratic <: Alg end
"Build SH diversity using the linear functions" struct Linear <: Alg end

using ..PhaseRetrieval: lineararray

"""
    SHdiversity(xrange, yrange, e₁, e₂, o)

Generate SH diversity phase on array `xrange × yrange` with cell based on orthogonal grid
formed by vectors  `e₁, e₂,` and origin at `o`.
"""
function SHdiversity(size, cellsize, celloffset=0, α=0)
    shd = zeros(size)

    # # first generate indexes
    # ind = orthIndexes(size, cellsize, celloffset, α)

    # # Now find range of indexes
    # irange = extrema(ind[:,:,1])
    # jrange = extrema(ind[:,:,2])
    # indoffset = -first.([irange,jrange]) .+ 1

    # # generate array of cell centres
    # centres = [centre(i, j, cellsize, celloffset, α)
    #             for i in UnitRange(irange...), j in UnitRange(jrange...)]

    (ind, centres, indoffset) = indexes_centres(size, cellsize, celloffset, α)

    # Method through MLA phase and global defocus
    # Now for every point in the array take the distance squared to centre of its cell
    # This will construct MLA phase
    for i in 1:size[1], j in 1:size[2]
        cellind = ind[i,j,:]
        cellcentre = centres[(cellind .+ indoffset)...]
        shd[i,j] = +((([i,j] - cellcentre).^2)...)
    end

    # Now we can subtract the global defocus
    defocus =  quadratic(size)
    shd -= defocus

    return shd
end

SHdiversity(size::Real, cellsize, celloffset, α) = SHdiversity((size, size), cellsize, celloffset, α)

function SHdiversity(size, cellsize, celloffset, α, ::Linear)
    shd = zeros(size)

    (ind, centres, indoffset) = indexes_centres(size, cellsize, celloffset, α)

    # Method through explicit piece-wise linear function
    # - 2 c_i x + c_i^2, with x = i -size/2
    for i in 1:size, j in 1:size
        cellind = ind[i,j,:]
        cellcentre = centres[(cellind .+ indoffset)...]
        shd[i,j] = -2*cellcentre .* (cellind .- size/2) .+ dot(cellcentre, cellcentre)

    end

    return shd
end

SHdiversity(size::Real, cellsize, celloffset, α, method) = SHdiversity((size,size), cellsize, celloffset, α, method)

function SHdiversity(xrange::AbstractRange, yrange::AbstractRange, e₁, e₂, o)
    cellind = indicesmap(xrange, yrange, e₁, e₂,o)
    irange  = UnitRange(extrema(cellind[1])...)
    jrange  = UnitRange(extrema(cellind[2])...)

    cellcentres = centresarray(cellind, e₁, e₂,o)

    s0 = 0.5*(e₁ + e₂) .+ o # centre of the zeroth cell
    p1 = dot(e₁, e₁)
    p2 = dot(e₂, e₂)
    s01 = dot(e₁, s0)
    s02 = dot(e₂, s0)
    s00 = dot(s0, s0)

    ret = lineararray(xrange, yrange, 2*s0, s00) # l0 = (s0,2 x) +s0^2
    l1 = lineararray(xrange, yrange, e₁)
    l2 = lineararray(xrange, yrange, e₂)

    for i in irange
        zi = cellind[1] .== i
        ci = - i^2*p1  - 2i*s01
        ret[zi] += 2i*l1[zi] .+ ci
    end

    for j in jrange
        zj = cellind[2] .== j
        cj = - j^2*p2  - 2j*s02
        ret[zj] += 2j*l2[zj] .+ cj
    end

    return ret
end

function indicesmap(xrange, yrange, e₁, e₂,o)
    iind =fld.(lineararray(xrange, yrange, e₁, - dot(e₁, o)), dot(e₁, e₁))
    jind =fld.(lineararray(xrange, yrange, e₂, - dot(e₂, o)), dot(e₂, e₂))
    return(iind, jind)

end

function centresarray(cellind, e₁, e₂,o)
    ijrange=extrema.(cellind)
    cij = [ (i + 0.5)*e₁ + (j + 0.5)*e₂ .+  o  for i in UnitRange(ijrange[1]...), j in UnitRange(ijrange[2]...)  ]
    return cij
end

# TODO finish SH diversity for a sensor and total size of the phase field
# function SHdiversity(wfs::SHSensor) =
#     SHdiversity(wfs.CameraChip.imagesize, cellsize, celloffset=0, α=0)

"""
    Generate array of indexes of an orthogonal geometry
"""
function orthIndexes(arrsize::Tuple, cellsize, gridorigin=0, α=0)
    ind = zeros(arrsize...,2)
    R = rotationmatrix(-α)
    x = 1- gridorigin[1] : arrsize[1]- gridorigin[1]
    y = reshape(1 - gridorigin[2] : arrsize[2] - gridorigin[2], (1,:))
    # for i in 1:arrsize, j in 1:arrsize
    #     ind[:,i, j] = R * [i, j]
    # end
    ind[:,:,1] = x .* R[1,1] .+ y .* R[1,2]
    ind[:,:,2] = x .* R[2,1] .+ y .* R[2,2]
    ind = fld.(ind , cellsize)
    return Int64.(ind)
end

orthIndexes(arrsize::Tuple, cellsize, gridorigin::Real, α=0) =
    orthIndexes(arrsize, cellsize, (gridorigin, gridorigin), α)

orthIndexes(arrsize::Real, cellsize, gridorigin, α=0) =
    orthIndexes((arrsize,arrsize), cellsize, gridorigin, α)

# this implementation appears to be quite slow
# function orthIndexes(arrsize, cellsize, gridorigin=0, α=0)
#     R = rotationmatrix(-α)
#     ind = [Int.(fld.(R * ( [i, j] .- gridorigin), cellsize))
#             for i in 1:arrsize, j in 1:arrsize]
# end

"""
    Give centre of an orthogonal grid cell of the grid with cellsize,
    centre at grid origin and rotated at angle α counterclockwise.
"""
function centre(i, j, cellsize, gridorigin, α)
    R = rotationmatrix(α)
    R * ([i,j] * cellsize  .+ cellsize / 2 ) .+ gridorigin
end

function indexes_centres(size, cellsize, celloffset, α)

    # first generate indexes
    ind = orthIndexes(size, cellsize, celloffset, α)

    # Now find range of indexes
    irange = extrema(ind[:,:,1])
    jrange = extrema(ind[:,:,2])
    indoffset = -first.([irange,jrange]) .+ 1

    # generate array of cell centres
    centres = [centre(i, j, cellsize, celloffset, α)
                for i in UnitRange(irange...), j in UnitRange(jrange...)]

    return (ind, centres, indoffset)
end

"""
    Micro-lens array structure.

    The main reason to introduce it separable is to handle the case when camera chip is located not exactly
    in the MLA focal plane.
"""
Base.@kwdef struct MLA
    pitch
    focallength
end

"""
    Shack-Hartmann sensor composed from MLA and a camera.

    MLA plane position is defined by distance,  β,  and γ.
    MLA subapertures position is defined by mlaorigin and α.
"""
Base.@kwdef struct SHSensor
    mla::MLA
    cam::CameraChip
    distance = mla.focallength
    mlaorigin = [0.,0.]
    α = 0.
    β = 0.
    γ = 0.
end

struct HPPlan
    wfs::SHSensor
    width
    height

end

struct PFPlan
    xdomain::Domain{2}
    fdomain
    ap
    diversity
    # f
    # l
end

function upscaleFactor(wfs::SHSensor, apertureD, λ )
    upscaleFactor(wfs.cam.pixelsize, wfs.mla.focallength, apertureD, λ )
end

"""
    constructSHPhaseDiversity(wfs::SHSensor, d, λ )

Construct SH phase diversity of minimal size for `wfs` within aperture of diameter `d` and wavelength `λ`.
Return tuple of arrays for the aperture and phase.
"""
function constructSHPhaseDiversity(wfs::SHSensor, apertureD, λ )
    s = wfs.cam.pixelsize
    f = wfs.mla.focallength
    p = wfs.mla.pitch
    o = wfs.mlaorigin
    α = wfs.α
    imsize = wfs.cam.imagesize

    # first check the dimensions
    q = apertureD* s /(f *λ)
    # if upscale=1, we are fine. otherwise we need to upscale
    upscale = upscaleFactor(wfs, apertureD, λ )
    cols,rows = imsize .* upscale
    # size of aperture
    arrsize = ceil(Int,apertureD/s * upscale)
    phasescale = pi/(1*arrsize) * q / upscale
    # or size of the camera chip
    # TODO extend to different scaling
    # arrsize =  imsize .* upscale
    # phasescale = pi/(1*arrsize[1]) * q / upscale
    cellsize = wfs.mla.pitch/s * upscale

    # in pixel coordianates
    # shd = SHdiversity(arrsize, cellsize, o, α) * phasescale
    # ap, apmask = diskmatrix(arrsize[1], q/upscale)

    # in physical coordinates
    # rng = -apertureD/2 : s/upscale : apertureD # wrong, HO size can be bigger
    # xrng= ((1:cols) .- ceil(cols/2)) * 1/(s*cols / upscale) * (λ * f)
    # yrng= ((1:rows) .- ceil(rows/2)) * 1/(s*rows / upscale) * (λ * f)

    # even better, through the domain
    wfsdom = make_centered_domain2D(imsize...,s)
    xrng = dualRange(wfsdom.xrange, upscale) * (λ * f)
    yrng = dualRange(wfsdom.yrange, upscale) * (λ * f)
    basis = rotationmatrix(α) * p
    shd = - π / (f*λ) * SHdiversity(xrng, yrng, basis[:,1], basis[:,2],o  )
    ap, apmask = aperture(xrng,yrng, apertureD)
    return (ap,shd)

end

make_centered_domain2D(wfs::SHSensor) = make_centered_domain2D(wfs.cam.imagesize..., wfs.cam.pixelsize)


text:
using Test
using LinearAlgebra
using SparseArrays
using Statistics

@testset "ConvDiffMIPDECO" begin
@testset "getBoundaryCondition" begin
	include("testGetBoundaryCondition.jl")
end
@testset "getBoundaryFuncVal" begin
	include("testGetBoundaryFuncVal.jl")
end
@testset "getGhostIndices" begin
	include("testGetGhostIndices.jl")
end
@testset "ConvDiffParam" begin
	include("testConvDiffParam.jl")
end
@testset "FEM" begin
	include("testFEM.jl")
end
@testset "FEM3D" begin
	include("testFEM3D.jl")
end
@testset "run2D" begin
	include("testRun2D.jl")
end
@testset "convDiff2D" begin
	include("testConvDiff2D.jl")
end
@testset "convDiff3D" begin
	include("testConvDiff3D.jl")
end
@testset "FEM2Dmms" begin
	include("testFEM2Dmms.jl")
end
@testset "FEM3Dmms" begin
	include("testFEM3Dmms.jl")
end
@testset "getConvDiffFEMConstraintsAMPL" begin
	include("testGetConvDiffFEMConstraintsAMPL.jl")
end

@testset "regularizers" begin
	include("testRegularizers.jl")
end
end

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
################################################################################
#
# TWRI module
#
################################################################################

module TWRIdual

using LinearAlgebra
using JUDI.TimeModeling
using PyCall
using Distributed
using FFTW, DSP
using SetIntersectionProjection
using Caching

const R = Float32; export R
const C = ComplexF32; export C
const SPACE_ORDER = 8
const NB = 40

# TWRI julia sources
include("TWRI_jl/utils.jl")
include("TWRI_jl/gendata.jl")
include("TWRI_jl/data_filter.jl")
include("TWRI_jl/FWIFun.jl")
include("TWRI_jl/TWRIdualFun.jl")

end


text:
function knn(k::Integer, data::DataSet, un_class_fv::FeatureVector, dist_func::Function)
  # get neighbors
  distances = Array((Number,Any),length(data))
  classifications = collect(keys(data))
  i = 1
  for class in classifications
    for fv in data[class].vectors
      distances[i] = dist_func(un_class_fv,fv[2]), class
      i += 1
    end
  end

  neighbors = FeatureVector()
  Base.sort!(distances)
  for x in 1:k
    neighbors[distances[x][2]] += 1
  end

  class_name = ""
  majority = 0
  for class in keys(neighbors)
    if neighbors[class] > majority
      majority = neighbors[class]
      class_name = class
    end
  end

  return class_name
end


text:
using MacroTools: splitstructdef, combinestructdef

macro nothing_macro()
end
@test @expand(@nothing_macro) === nothing

macro splitcombine(fundef) # should be a no-op
    dict = splitdef(fundef)
    esc(MacroTools.combinedef(dict))
end

# Macros for testing that splitcombine doesn't break
# macrocalls in bodies
macro zeroarg()
   :(1)
end
macro onearg(x)
   :(1+$(esc(x)))
end

let
    # Ideally we'd compare the result against :(function f(x)::Int 10 end),
    # but it fails because of :line and :block differences
    @test longdef(:(f(x)::Int = 10)).head == :function
    @test longdef(:(f(x::T) where U where T = 2)).head == :function
    @test shortdef(:(function f(x)::Int 10 end)).head != :function
    @test map(splitarg, (:(f(a=2, x::Int=nothing, y, args...))).args[2:end]) ==
        [(:a, :Any, false, 2), (:x, :Int, false, :nothing),
         (:y, :Any, false, nothing), (:args, :Any, true, nothing)]
    @test splitarg(:(::Int)) == (nothing, :Int, false, nothing)

    @splitcombine foo(x) = x+2
    @test foo(10) == 12
    @splitcombine add(a, b=2; c=3, d=4)::Float64 = a+b+c+d
    @test add(1; d=10) === 16.0
    @splitcombine fparam(a::T) where {T} = T
    @test fparam([]) == Vector{Any}
    struct Orange end
    @splitcombine (::Orange)(x) = x+2
    @test Orange()(10) == 12
    @splitcombine fwhere(a::T) where T = T
    @test fwhere(10) == Int
    @splitcombine manywhere(x::T, y::Vector{U}) where T <: U where U = (T, U)
    @test manywhere(1, Number[2.0]) == (Int, Number)
    @splitcombine fmacro0() = @zeroarg
    @test fmacro0() == 1
    @splitcombine fmacro1() = @onearg 1
    @test fmacro1() == 2

    struct Foo{A, B}
        a::A
        b::B
    end
    # Parametric outer constructor
    @splitcombine Foo{A}(a::A) where A = Foo{A, A}(a,a)
    @test Foo{Int}(2) == Foo{Int, Int}(2, 2)

    @test (@splitcombine x -> x + 2)(10) === 12
    @test (@splitcombine (a, b=2; c=3, d=4) -> a+b+c+d)(1; d=10) === 16
    @test (@splitcombine ((a, b)::Tuple{Int,Int} -> a + b))((1, 2)) == 3
    @test (@splitcombine ((a::T) where {T}) -> T)([]) === Vector{Any}
    @test (@splitcombine ((x::T, y::Vector{U}) where T <: U where U) -> (T, U))(1, Number[2.0]) ==
          (Int, Number)
    @test (@splitcombine () -> @zeroarg)() == 1
    @test (@splitcombine () -> @onearg 1)() == 2
    @test (@splitcombine function (x) x + 2 end)(10) === 12
    @test (@splitcombine function (a::T) where {T} T end)([]) === Vector{Any}
    @test (@splitcombine function (x::T, y::Vector{U}) where T <: U where U
               (T, U)
           end)(1, Number[2.0]) == (Int, Number)
end

@testset "combinestructdef, splitstructdef" begin
    ex = :(struct S end)
    @test ex |> splitstructdef |> combinestructdef |> Base.remove_linenums! ==
        :(struct S <: Any end)

    @test splitstructdef(ex) == Dict(
        :constructors => Any[],
        :mutable => false,
        :params => Any[],
        :name => :S,
        :fields => Any[],
        :supertype => :Any)

    ex = :(mutable struct T end)
    @test splitstructdef(ex)[:mutable] === true
    @test ex |> splitstructdef |> combinestructdef |> Base.remove_linenums! ==
        :(mutable struct T <: Any end)

    ex = :(struct S{A,B} <: AbstractS{B}
                               a::A
                           end)
    @test splitstructdef(ex) == Dict(
        :constructors => Any[],
        :mutable => false,
        :params => Any[:A, :B],
        :name => :S,
        :fields => Any[(:a, :A)],
        :supertype => :(AbstractS{B}),)

    @test ex |> splitstructdef |> combinestructdef |> Base.remove_linenums! ==
        ex |> Base.remove_linenums!

    ex = :(struct S{A} <: Foo; S(a::A) where {A} = new{A}() end)
    @test ex |> splitstructdef |> combinestructdef |>
        Base.remove_linenums! |> MacroTools.flatten ==
        ex |> Base.remove_linenums! |> MacroTools.flatten

    constructors = splitstructdef(ex)[:constructors]
    @test length(constructors) == 1
    @test first(constructors) ==
        :((S(a::A) where A) = new{A}()) |> MacroTools.flatten

    @test_throws ArgumentError splitstructdef(:(call_ex(arg)))
end


text:
################################################################################
# Convenience functions for kinematics
################################################################################

#..oooOOO0Q0Q0OOOooo...oooOOO0Q0Q0OOOooo...oooOOO0Q0Q0OOOooo...oooOOO0Q0Q0OOOooo

"""
	β(γ::Number)

Returns β given γ.
"""
β(γ::Number) = √(γ^2 - 1) / γ

export β
#..oooOOO0Q0Q0OOOooo...oooOOO0Q0Q0OOOooo...oooOOO0Q0Q0OOOooo...oooOOO0Q0Q0OOOooo
"""
    function momentum(E::Quantity,id::ParticleID)

Returns momentum given total energy and rest mass.
"""
function momentum(E::Quantity, id::ParticleID)
    @assert typeof(E) <: Unitful.Energy "E requires units of energy"
    m = mass(id)
    return uconvert(u"GeV/c",√(E^2 - m^2 * Unitful.c^4) / Unitful.c)
end #momentum

"""
	function momentum(E::Quantity,m::Quantity)

Returns momentum given total energy and rest mass.
"""
function momentum(E::Quantity, m::Quantity)
    @assert typeof(E) <: Unitful.Energy "E requires units of energy"
    @assert typeof(m) <: Unitful.Mass "m requires units of mass"
    return uconvert(u"GeV/c",√(E^2 - m^2 * Unitful.c^4) / Unitful.c)
end #momentum

"""
    function momentumsi(E::Float64,id::ParticleID)

Returns momentum given total energy and a ParticleID. SI units used.
"""
function momentum_si(E::Float64, id::ParticleID)
    m = usistrip(mass(id))
    momentum_si(E,m)
end #momentum

"""
	function momentum(E::Number,m::Number)

Returns momentum given total energy and rest mass in SI units.
"""
function momentum_si(E::Number, m::Number)
    return √(E^2 - m^2 * Unitful.c0.val^4) / Unitful.c0.val
end #momentum

export momentum, momentum_si
#..oooOOO0Q0Q0OOOooo...oooOOO0Q0Q0OOOooo...oooOOO0Q0Q0OOOooo...oooOOO0Q0Q0OOOooo

"""
	γFromβ(β::Number)

Returns γ given β
"""
function γFromβ(β::Number)
    @assert β <= 1.0 "β ≦ 1 required"
    return 1 / sqrt(1 - β^2)
end #γFromβ

#..oooOOO0Q0Q0OOOooo...oooOOO0Q0Q0OOOooo...oooOOO0Q0Q0OOOooo...oooOOO0Q0Q0OOOooo

"""
    γFromP(p::Array{T}, m::Quantity) where {T<:Quantity}

Returns gamma given momentum and rest mass.
"""
function γFromP(p::Array{T}, m::Quantity) where {T<:Quantity}
    pnorm = norm(p)
    γFromP(pnorm,m)
end #γFromP

"""
    γFromP(p::Array{T}, m::Number) where {T<:Number}

Returns gamma given momentum and rest mass in SI units.
"""
function γFromP(p::Array{T}, m::Number) where {T<:Number}
    pnorm = norm(p)
    γFromP(pnorm,m)
end #γFromP

"""
γFromP(p::Quantity,m::Quantity)

Returns gamma given momentum and rest mass.
"""
function γFromP(p::Quantity, m::Quantity)
    @assert typeof(p) <: Unitful.Momentum "p requires units of momentum"
    @assert typeof(m) <: Unitful.Mass "m requires units of mass"
    c = Unitful.c
    βγ = upreferred(p / (m * c))
    return upreferred(sqrt(1 + βγ^2))
end #γFromP

"""
	γFromP(p::Number,m::Number)

Returns gamma given momentum and rest mass in SI units.
"""
function γFromP(p::Number, m::Number)
    c = Unitful.c0.val
    return sqrt(1 + (p / (m * c))^2)
end #γFromP

export γFromP
#..oooOOO0Q0Q0OOOooo...oooOOO0Q0Q0OOOooo...oooOOO0Q0Q0OOOooo...oooOOO0Q0Q0OOOooo
"""
	βFromP(p::Quantity,m::Quantity)

Returns β given particle momentum and rest mass.
"""
function βFromP(p::Quantity, m::Quantity)
    @assert typeof(p) <: Unitful.Momentum "p requires units of momentum"
    @assert typeof(m) <: Unitful.Mass "m requires units of mass"
    c = Unitful.c
    η = upreferred(p / (m * c))
    return sqrt(η^2 / (1 + η^2))
end

"""
	βFromP(p::Number,m::Number)

Returns β given particle momentum and rest mass in SI units
"""
function βFromP(p::Number, m::Number)
    c = Unitful.c0.val
    η = p / (m * c)
    return sqrt(η^2 / (1 + η^2))
end

export βFromP
#..oooOOO0Q0Q0OOOooo...oooOOO0Q0Q0OOOooo...oooOOO0Q0Q0OOOooo...oooOOO0Q0Q0OOOooo
"""
    vFromP(p::Array{T}, m::U) where {T <: Quantity, U <: Quantity}

Returns velocity given particle momentum and rest mass in SI units
"""
function vFromP(p::Array{T}, m::U) where {T <: Quantity, U <: Quantity}
    β = βFromP(norm(p),m)
    c = Unitful.c
    phat = p / norm(p)
    return upreferred.(c * β * phat)
end

"""
    vFromP(p::Array{T}, m::T) where {T <: Number}

Returns velocity given particle momentum and rest mass in SI units
"""
function vFromP(p::Array{T}, m::T) where {T <: Number}
    β = βFromP(norm(p),m)
    c = Unitful.c0.val
    phat = p / norm(p)
    return c * β * phat
end

"""
    vFromP(p::Quantity, m::Quantity)

Returns velocity given particle momentum and rest mass in SI units
"""
function vFromP(p::Quantity, m::Quantity)
    β = βFromP(p,m)
    return upreferred.(Unitful.c * β )
end

"""
    vFromP(p::Array{T}, m::T) where {T <: Number}

Returns velocity given particle momentum and rest mass in SI units
"""
function vFromP(p::Number, m::Number)
    β = βFromP(p,m)
    c = Unitful.c0.val
    return c * β
end

export vFromP
#..oooOOO0Q0Q0OOOooo...oooOOO0Q0Q0OOOooo...oooOOO0Q0Q0OOOooo...oooOOO0Q0Q0OOOooo
"""
	gyroradius(p::Quantity,Z::Quantity,B::Quantity)

Returns gyroradius given dimensionally correct momentum, charge, and B field strength
"""
function gyroradius(p::Quantity, Q::Quantity, B::Quantity)
    @assert typeof(p) <: Unitful.Momentum "p must have units of momentum"
    @assert typeof(Q) <: Unitful.Charge "Z must have units of charge"
    @assert typeof(B*1u"m^2") <: Unitful.MagneticFlux "B must have units of magnetic flux density"

    return upreferred(p/(Q*B))
end
# TODO: add functionality to allow for vector p and B

"""
	gyroradius(p::Number, Z::Number, B::Number)

Returns gyroradius (m) given momentum, atomic (charge) number, and B field strength

# Arguments
- `p::Number`: momentum of particle in GeV/c
- `Z::Number`: charge of particle in units of q_e
- `B::Number`: magnetic flux density in gauss
"""
function gyroradius(p::Number, Z::Number, B::Number)
    33356.40951981521 * p / (Z * B)
end

export gyroradius


text:
const sense_alias = Dict(
"max"      => :Max,
"maximize" => :Max,
"maximise" => :Max,
"maximum"  => :Max,
"min"      => :Min,
"minimize" => :Min,
"minimise" => :Min,
"minimum"  => :Min
)

const subject_to_alias = ["subject to", "such that", "st", "s.t."]

# a list of section keywords in lower-case
const KEYWORDS = Dict(
    "max"      => Val{:obj},
    "maximize" => Val{:obj},
    "maximise" => Val{:obj},
    "maximum"  => Val{:obj},
    "min"      => Val{:obj},
    "minimize" => Val{:obj},
    "minimise" => Val{:obj},
    "minimum"  => Val{:obj},

    "subject to" => Val{:constraints},
    "such that"  => Val{:constraints},
    "st"         => Val{:constraints},
    "s.t."       => Val{:constraints},

    "bounds" => Val{:bounds},
    "bound"  => Val{:bounds},

    "gen"      => Val{:integer},
    "general"  => Val{:integer},
    "generals" => Val{:integer},

    "bin"      => Val{:binary},
    "binary"   => Val{:binary},
    "binaries" => Val{:binary},

    "end"      => Val{:quit}
)

const COMMENT_REG = r"(.*?)\\(.*)"
function stripcomment(line::String)
    if contains(line, "\\")
        m = match(COMMENT_REG, line)
        return strip(String(m[1]))
    else
        return strip(line)
    end
end

immutable TempSparseMatrix
    i::Vector{Int}
    j::Vector{Int}
    v::Vector{Float64}
end
TempSparseMatrix() = TempSparseMatrix(Int[], Int[], Float64[])
Base.sparse(m::TempSparseMatrix, nr::Int, nc::Int) = sparse(m.i, m.j, m.v, nr, nc)

newdatastore() = Dict(
    :A => TempSparseMatrix(),
    :collb => Float64[],
    :colub => Float64[],
    :c => Float64[],
    :rowlb => Float64[],
    :rowub => Float64[],
    :sense => :Min,
    :colcat => Symbol[],
    :sos => SOS[],
    :Q => TempSparseMatrix(),
    :modelname => "",
    :colnames => String[],
    :rownames => String[],
    :open_constraint => false
)
setsense!(T, data, line) = nothing
function setsense!(::Type{Val{:obj}}, data, line)
    data[:sense] = sense_alias[lowercase(line)]
end
function read(filename::String)
    data = newdatastore()
    open(filename, "r") do io
        section = :none
        while !eof(io)
            line = stripcomment(readline(io))
            if line == "" # skip blank lines
                continue
            end
            if haskey(KEYWORDS, lowercase(line)) # section has changed
                section = KEYWORDS[lowercase(line)]
                setsense!(section, data, line)
                continue
            end
            parsesection!(section, data, line)
        end
    end
    return sparse(data[:A], length(data[:rownames]), length(data[:colnames])),
        data[:collb],
        data[:colub],
        data[:c],
        data[:rowlb],
        data[:rowub],
        data[:sense],
        data[:colcat],
        data[:sos],
        sparse(data[:Q], length(data[:rownames]), length(data[:colnames])),
        data[:modelname],
        data[:colnames],
        data[:rownames]
end

parsesection!(::Type{Val{:quit}}, data, line) = error("Corrupted LP File. You have the lne $(line) after an end.")

function parsesection!(::Type{Val{:obj}}, data, line)
    # okay so line should be the start of the objective
    if contains(line, ":")
        # throw away name
        m = match(r"(.*?)\:(.*)", line)
        line = String(m[2])
    end
    tokens = tokenize(line)
    if length(tokens) == 0 # no objective
        return
    end
    # tokens should be in order (+/-) (numeric) (variable) ...
    while length(tokens) > 0
        variable = String(pop!(tokens))
        idx = getvariableindex!(data, variable) # catch in here for malformed variables
        coef_token = pop!(tokens)
        try
            coefficient = parse(Float64, coef_token)
        catch
            error("Unable to parse objective. Expected numeric coefficient. Got $(coef_token)")
        end
        if length(tokens) > 0
            _sign = pop!(tokens)
            if _sign == "-"
                coefficient *= -1
            elseif _sign == "+"
            else
                error("Unable to parse objective due to bad operator: $(_sign) $(line)")
            end
        end
        data[:c][idx] = coefficient
    end
end

const constraintsense = Dict(
    "<"  => :le,
    "<=" => :le,
    "="  => :eq,
    "==" => :eq,
    ">"  => :ge,
    ">=" => :ge,
)

function parseconstraintcoefficients!(data, line, tokens, rowidx)
    # tokens should be in order (+/-) (numeric) (variable) ...
    while length(tokens) > 0
        variable = String(pop!(tokens))
        idx = getvariableindex!(data, variable) # catch in here for malformed variables
        coef_token = pop!(tokens)
        try
            coefficient = parse(Float64, coef_token)
        catch
            error("Unable to parse constraint $(line). Expected numeric coefficient. Got $(coef_token)")
        end
        if length(tokens) > 0
            _sign = pop!(tokens)
            if _sign == "-"
                coefficient *= -1
            elseif _sign == "+"
            else
                error("Unable to parse objective due to bad operator: $(_sign) $(line)")
            end
        end
        push!(data[:A].i, rowidx)
        push!(data[:A].j, idx)
        push!(data[:A].v, coefficient)
    end
end

function parsesos!(data, line)
    tokens = tokenize(line)
    if length(tokens) < 3
        error(string("Malformed SOS constraint: ", line))
    end
    if tokens[2] == "S1::"
        order = 1
    elseif tokens[2] == "S2::"
        order = 2
    else
        error("SOS of type $(tokens[2]) not recognised")
    end
    names   = String[]
    weights = Float64[]
    for token in tokens[3:end]
        items = split(token, ":")
        if length(items) != 2
            error(string("Invalid sequence: ", token))
        end
        push!(names, String(items[1]))
        push!(weights, parsefloat(String(items[2])))
    end
    indices = Int[]
    for name in names
        idx = findfirst(data[:colnames], name)
        if idx == 0
            push!(data[:colnames], name)
            push!(indices, length(data[:colnames]))
        else
            push!(indices, idx)
        end
    end
    push!(data[:sos], (order, indices, weights))
end

function parsesection!(::Type{Val{:constraints}}, data, line)
    if match(r" S([0-9]):: ", line) != nothing
        # it's an SOS constraint
        parsesos!(data, line)
        return
    end
    if data[:open_constraint] == false
        push!(data[:rownames], "R$(length(data[:rownames]) + 1)")
        push!(data[:rowlb], -Inf)
        push!(data[:rowub], Inf)
    end
    if contains(line, ":")
        if data[:open_constraint] == true
            error("Malformed constraint $(line). Is the previous one valid?")
        end
        # throw away name
        m = match(r"(.*?)\:(.*)", line)
        data[:rownames][end] = String(m[1])
        line = String(m[2])
    end
    data[:open_constraint] = true

    tokens = tokenize(line)
    if length(tokens) == 0 # no entries
        return
    elseif length(tokens) >= 2 && haskey(constraintsense, tokens[end-1])# test if constraint ends this line
        rhs = parsefloat(pop!(tokens))
        sym = pop!(tokens)
        if constraintsense[sym] == :le
            data[:rowub][end] = rhs
        elseif constraintsense[sym] == :ge
            data[:rowlb][end] = rhs
        elseif constraintsense[sym] == :eq
            data[:rowlb][end] = rhs
            data[:rowub][end] = rhs
        end
        data[:open_constraint] = false # finished
    end
    parseconstraintcoefficients!(data, line, tokens, length(data[:rownames]))
end

function getvariableindex!(data, v)
    i = findfirst(data[:colnames], v)
    if i == 0
        if !verifyname(v)
            error("Invalid variable name $v")
        end
        addnewvariable!(data, v)
        return length(data[:colnames])
    end
    return i
end

function addnewvariable!(data, name)
    push!(data[:collb], -Inf)
    push!(data[:colub],  Inf)
    push!(data[:c], 0)
    push!(data[:colcat], :Cont)
    push!(data[:colnames], name)
end

function parsevariabletype!(data, line, cat)
    items = tokenize(line)
    for v in items
        i = getvariableindex!(data, v)
        data[:colcat][i] = cat
    end
end
parsesection!(::Type{Val{:integer}}, data, line) = parsevariabletype!(data, line, :Int)
parsesection!(::Type{Val{:binary}}, data, line)  = parsevariabletype!(data, line, :Bin)

function parsefloat(val::String)
    if lowercase(val) == "-inf" || lowercase(val) == "-infinity"
        return -Inf
    elseif lowercase(val) == "+inf" || lowercase(val) == "+infinity"
        return Inf
    else
        return parse(Float64, val)
    end
end
function tokenize(line)
    items = String.(split(line, " "))
    items[items .!= ""]
end
bounderror(line) = error("Unable to parse bound: $(line)")
function parsesection!(::Type{Val{:bounds}}, data, line)
    items = tokenize(line)
    v = ""
    lb = -Inf
    ub = Inf
    if length(items) == 5 # ranged bound
        v = items[3]
        if (items[2] == "<=" || items[2] == "<") &&  (items[4] == "<=" || items[4] == "<") # le
            lb = parsefloat(items[1])
            ub = parsefloat(items[5])
        elseif (items[2] == ">=" || items[2] == ">") &&  (items[4] == ">=" || items[4] == ">") # ge
            lb = parsefloat(items[5])
            ub = parsefloat(items[1])
        else
            bounderror(line)
        end
    elseif length(items) == 3 # one sided
        v = items[1]
        if items[2] == "<=" || items[2] == "<" # le
            ub = parsefloat(items[3])
        elseif items[2] == ">=" || items[2] == ">" # ge
            lb = parsefloat(items[3])
        elseif items[2] == "==" || items[2] == "=" # eq
            lb = ub = parsefloat(items[3])
        else
            bounderror(line)
        end
    elseif length(items) == 2 # free
        if items[2] != "free"
            bounderror(line)
        end
        v = items[1]
    else
        bounderror(line)
    end
    i = getvariableindex!(data, v)
    data[:collb][i] = lb
    data[:colub][i] = ub

end


text:
"""
    compose(g::Function,f...)

Returns a new field obtained by composition of function `g` and the fields
`f`. The value of the resulting field at a vector of points `x` is numerically equivalent to

    fx = evaluate_fields(f,x)
    apply_kernel(bcast(g), fx...)

The gradient of the resulting field evaluated at a vector of points `x` is equivalent to

    fx = evaluate_fields(f,x)
    apply_kernel(bcast(gradient(g)), fx...)

Note that it is needed to overload `gradient(::typeof(g))` for the given function `g`
in order to be able to compute the gradient.

As in function [`apply_kernel_to_field`](@ref) if any of the inputs in `f` is a number or an array
instead of a field it will be treated as a "constant field".
"""
function compose(g::Function,f...)
  k = Comp(g)
  apply_kernel_to_field(k,f...)
end

struct Comp{F} <: Kernel
  e::BCasted{F}
  @inline Comp(f::Function) = new{typeof(f)}(BCasted(f))
end

@inline apply_kernel!(cache,k::Comp,x...) = apply_kernel!(cache,k.e,x...)

kernel_cache(k::Comp,x...) = kernel_cache(k.e,x...)

kernel_return_type(k::Comp,x...) = kernel_return_type(k.e,x...)

function apply_kernel_gradient(k::Comp,f...)
  g = gradient(k.e.f)
  compose(g,f...)
end

"""
    compose(g::Function,f::AbstractArray...)

Returns an array of fields numerically equivalent to

    map( (x...)->compose(g,x...), f...)

"""
function compose(g::Function,f::AbstractArray...)
  k = Comp(g)
  apply_to_field_array(k,f...)
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -47.221875796417265 <= q <= 49.11122836315772)

                     add_NL_constraint(m, :(swish(0.42790162605340454 + 0.9948146339173967*swish(0.30429674417798 + 0.01326158862793747*swish(0.7020479917236249 + 0.5157446936496526*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.8506753528903843*swish(0.7594455380381384 + 0.9088497207479063*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.7699605184902198*swish(0.5333438081093944 + -0.9884706408297048*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.9481609482474158*swish(-0.8688273038537231 + -0.5762764679031651*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + 0.5690258167368487*swish(0.4003333530904172 + 0.4565361301167439*swish(0.7020479917236249 + 0.5157446936496526*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.41218949132242644*swish(0.7594455380381384 + 0.9088497207479063*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.6719961553691296*swish(0.5333438081093944 + -0.9884706408297048*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.9454422589116072*swish(-0.8688273038537231 + -0.5762764679031651*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + -0.632632285898699*swish(-0.4740061987673023 + -0.21503151130953757*swish(0.7020479917236249 + 0.5157446936496526*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.7641417558802628*swish(0.7594455380381384 + 0.9088497207479063*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.8796988576372375*swish(0.5333438081093944 + -0.9884706408297048*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.6466469730144713*swish(-0.8688273038537231 + -0.5762764679031651*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + -0.0799939753544896*swish(-0.9533597279323933 + 0.019939841305541606*swish(0.7020479917236249 + 0.5157446936496526*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.8134733780415981*swish(0.7594455380381384 + 0.9088497207479063*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.10348468441851777*swish(0.5333438081093944 + -0.9884706408297048*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.008542712959329624*swish(-0.8688273038537231 + -0.5762764679031651*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + swish(-0.31282642357163937 + 0.400776902433146*swish(0.30429674417798 + 0.01326158862793747*swish(0.7020479917236249 + 0.5157446936496526*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.8506753528903843*swish(0.7594455380381384 + 0.9088497207479063*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.7699605184902198*swish(0.5333438081093944 + -0.9884706408297048*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.9481609482474158*swish(-0.8688273038537231 + -0.5762764679031651*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + -0.8700113611274283*swish(0.4003333530904172 + 0.4565361301167439*swish(0.7020479917236249 + 0.5157446936496526*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.41218949132242644*swish(0.7594455380381384 + 0.9088497207479063*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.6719961553691296*swish(0.5333438081093944 + -0.9884706408297048*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.9454422589116072*swish(-0.8688273038537231 + -0.5762764679031651*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + -0.8173535730312884*swish(-0.4740061987673023 + -0.21503151130953757*swish(0.7020479917236249 + 0.5157446936496526*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.7641417558802628*swish(0.7594455380381384 + 0.9088497207479063*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.8796988576372375*swish(0.5333438081093944 + -0.9884706408297048*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.6466469730144713*swish(-0.8688273038537231 + -0.5762764679031651*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + -0.10383649475667989*swish(-0.9533597279323933 + 0.019939841305541606*swish(0.7020479917236249 + 0.5157446936496526*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.8134733780415981*swish(0.7594455380381384 + 0.9088497207479063*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.10348468441851777*swish(0.5333438081093944 + -0.9884706408297048*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.008542712959329624*swish(-0.8688273038537231 + -0.5762764679031651*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + swish(0.07748802701261592 + -0.48287575350847955*swish(0.30429674417798 + 0.01326158862793747*swish(0.7020479917236249 + 0.5157446936496526*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.8506753528903843*swish(0.7594455380381384 + 0.9088497207479063*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.7699605184902198*swish(0.5333438081093944 + -0.9884706408297048*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.9481609482474158*swish(-0.8688273038537231 + -0.5762764679031651*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + 0.9032398979921337*swish(0.4003333530904172 + 0.4565361301167439*swish(0.7020479917236249 + 0.5157446936496526*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.41218949132242644*swish(0.7594455380381384 + 0.9088497207479063*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.6719961553691296*swish(0.5333438081093944 + -0.9884706408297048*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.9454422589116072*swish(-0.8688273038537231 + -0.5762764679031651*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + 0.8154202752057489*swish(-0.4740061987673023 + -0.21503151130953757*swish(0.7020479917236249 + 0.5157446936496526*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.7641417558802628*swish(0.7594455380381384 + 0.9088497207479063*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.8796988576372375*swish(0.5333438081093944 + -0.9884706408297048*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.6466469730144713*swish(-0.8688273038537231 + -0.5762764679031651*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + 0.8885160528665201*swish(-0.9533597279323933 + 0.019939841305541606*swish(0.7020479917236249 + 0.5157446936496526*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.8134733780415981*swish(0.7594455380381384 + 0.9088497207479063*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.10348468441851777*swish(0.5333438081093944 + -0.9884706408297048*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.008542712959329624*swish(-0.8688273038537231 + -0.5762764679031651*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + swish(0.6389491685818425 + -0.9024035564379238*swish(0.30429674417798 + 0.01326158862793747*swish(0.7020479917236249 + 0.5157446936496526*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.8506753528903843*swish(0.7594455380381384 + 0.9088497207479063*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.7699605184902198*swish(0.5333438081093944 + -0.9884706408297048*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.9481609482474158*swish(-0.8688273038537231 + -0.5762764679031651*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + -0.8688708754259227*swish(0.4003333530904172 + 0.4565361301167439*swish(0.7020479917236249 + 0.5157446936496526*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.41218949132242644*swish(0.7594455380381384 + 0.9088497207479063*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.6719961553691296*swish(0.5333438081093944 + -0.9884706408297048*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.9454422589116072*swish(-0.8688273038537231 + -0.5762764679031651*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + -0.9399903753655909*swish(-0.4740061987673023 + -0.21503151130953757*swish(0.7020479917236249 + 0.5157446936496526*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.7641417558802628*swish(0.7594455380381384 + 0.9088497207479063*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.8796988576372375*swish(0.5333438081093944 + -0.9884706408297048*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.6466469730144713*swish(-0.8688273038537231 + -0.5762764679031651*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + 0.5656846187923494*swish(-0.9533597279323933 + 0.019939841305541606*swish(0.7020479917236249 + 0.5157446936496526*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.8134733780415981*swish(0.7594455380381384 + 0.9088497207479063*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.10348468441851777*swish(0.5333438081093944 + -0.9884706408297048*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.008542712959329624*swish(-0.8688273038537231 + -0.5762764679031651*swish(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*swish(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*swish(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*swish(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using PyPlot
using JLD
nSteps = 200001
X = load(string("../data/sens/solenoid/metrics_",nSteps,"_K.jld"))
dJds_ust = X["dJds_ust"]
dJds_st = X["dJds_st"]
ust_cn_k = X["unstable_contribution_k"]
K = size(ust_cn_k)[1]
fig, ax = subplots()
#ust_pge = abs.((cumsum(ust_cn_k) .- dJds_ust)./dJds_ust)*100
ust_pge = abs.(ust_cn_k./dJds_ust)*100
ax.semilogy(1:K, ust_pge, "ko-", ms=20, lw=2.5, label="% error in sum of k terms")
ax.set_xlabel("k",fontsize=28)
ax.xaxis.set_tick_params(labelsize=28)
ax.yaxis.set_tick_params(labelsize=28)
ax.grid(true)
leg = ax.legend(fontsize=28)
leg.legendHandles[1]._legmarker.set_markersize(10)

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
@testset "MatrixDepot" begin
debit = [ ("hilb", 5),  ("cauchy",4), ("invhilb",4), ("moler",7),("pascal",6),("pei",100),("tridiag",100),("lehmer",100),("poisson",12),("wathen",10)]

for k=1:length(debit)
	A   = matrixdepot(debit[k][1],debit[k][2])
	rhs = randn(size(A,1))
	@printf "%s\tn=[%d,%d]\t" debit[k][1] size(A,1) size(A,2)

	xt  = A\rhs

	x3, = gmres(A,rhs,10,maxIter=50*length(rhs),tol=1e-16,out=0)
	@test norm(xt-x3)/norm(xt) < 1e-3
	@printf " gmres is OK "
	x1, = cg(A,rhs,maxIter=length(rhs),tol=1e-16,out=-1)
	@test norm(xt-x1)/norm(xt) < 1e-3
	@printf "cg is OK\n"
end

end


text:
export getSSORCGFourthOrderRegularizationPreconditioner,SSORCGParamFourthOrder, applySSORCGFourth, setupDoNothing, wFourthOrderSmoothingNodal

function wFourthOrderSmoothingNodal(m::Vector, mref::Vector, M::AbstractMesh; Iact=1.0, C=[])
	dm = m.-mref;
	d2R = wdiffusionRegNodal(m,mref,M,Iact = Iact,C = C)[3];
	clear!(M);
	d2R = d2R'*d2R;
	dR  = d2R*dm;
	Rc  = 0.5*dot(dm,dR);
   return Rc,dR,d2R
end

import jInv.InverseSolve.HessianPreconditioner
type SSORCGParamFourthOrder
	C 			::Vector;
	M			::AbstractMesh;
	Iact
	omega		::Float64
	tol			::Float64
	maxCGIter	::Int64
end

function getSSORCGFourthOrderRegularizationPreconditioner(C,M::AbstractMesh,Iact, omega::Float64=1.0,tol::Float64=1e-2,maxCGIter::Int64=100)
	return HessianPreconditioner(SSORCGParamFourthOrder(C,M,Iact,omega,tol,maxCGIter),applySSORCGFourth,setupDoNothing);
end

function getSAAMGFourthOrderRegularizationPreconditioner(C,M::AbstractMesh,Iact, omega::Float64=1.0,tol::Float64=1e-2,maxCGIter::Int64=100)
	return HessianPreconditioner(SSORCGParamFourthOrder(C,M,Iact,omega,tol,maxCGIter),applySAAMGFourth,setupDoNothing);
end

function applySAAMGFourth(Hs::Function, d2R::SparseMatrixCSC,v::Vector,param)
	d2R = wdiffusionRegNodal(v,v,param.M,Iact = param.Iact, C = param.C)[3]
	MG = getMGparam(5,8,500,1e-8,"SPAI",1.0,1,1,'V',"MUMPS");
	SA_AMGsetup(d2R,MG,eltype(v),true,1,false);
	x = zeros(eltype(v),size(v));
	x = solveCG_MG(d2R,MG,v,x,false)[1];
	v = copy(x); x[:] = 0.0;
	MG.relativeTol = 1e-4;
	x = solveCG_MG(d2R,MG,v,x,false)[1];
	return x;
end

function applySSORCGFourth(Hs::Function, d2R::SparseMatrixCSC,v::Vector,param)
	d2R = wdiffusionRegNodal(v,v,param.M,Iact = param.Iact, C = param.C)[3]
	Dinv = param.omega./diag(d2R);
	aux   = zeros(size(d2R,2));
	SSOR(r) = (aux[:]=0.0; ssorPrecTrans!(d2R,aux,r,Dinv); return aux);
	# println("In prec.")
	# tic();
	x = KrylovMethods.cg(d2R,v,tol=param.tol,maxIter=param.maxCGIter,M=SSOR,out=-1)[1]
	x = KrylovMethods.cg(d2R',copy(x),tol=param.tol,maxIter=param.maxCGIter,M=SSOR,out=-1)[1]
	# println("Done prec in ",toq()," seconds.");
	return x;
end

function setupDoNothing(Hs::Function, d2R::SparseMatrixCSC,param)
	return;
end

export getSSORCGDoubleTVRegularizationPreconditioner,wDoubleTVNodal,applySSORCGDoubleTV

function getSSORCGDoubleTVRegularizationPreconditioner(C,M::AbstractMesh,Iact, omega::Float64=1.0,tol::Float64=1e-2,maxCGIter::Int64=100)
	return HessianPreconditioner(SSORCGParamFourthOrder(C,M,Iact,omega,tol,maxCGIter),applySSORCGDoubleTV,setupDoNothing);
end

function wDoubleTVNodal(m::Vector, mref::Vector, M::AbstractMesh; Iact=1.0, C=[])
	dm = m.-mref;
	d2R = wTVRegNodal(m,mref,M,Iact = Iact,C = C)[3];
	clear(M);
	d2R = d2R'*d2R;
	dR  = d2R*dm;
	Rc  = 0.5*dot(dm,dR);
   return Rc,dR,d2R
end

function applySSORCGDoubleTV(Hs::Function, d2R::SparseMatrixCSC,v::Vector,param)
	d2R = wTVRegNodal(v,v,param.M,Iact = param.Iact, C = param.C)[3]
	Dinv = param.omega./diag(d2R);
	aux   = zeros(size(d2R,2));
	SSOR(r) = (aux[:]=0.0; ssorPrecTrans!(d2R,aux,r,Dinv); return aux);
	x = KrylovMethods.cg(d2R,v,tol=param.tol,maxIter=param.maxCGIter,M=SSOR,out=-1)[1]
	x = KrylovMethods.cg(d2R',copy(x),tol=param.tol,maxIter=param.maxCGIter,M=SSOR,out=-1)[1]
	return x;
end


text:
# Functions for working with irfs

function shift(data::AbstractVector{<:Real},Q::Integer)
    # Shifts the irf by n channels
    nbins = length(data)
    temp = data
    if Q > 0
        nsc = Q
        nec = nbins
        osc = 1
        oec = nbins - Q
        temp[1:(nsc - 1)] .= 0
    elseif Q < 0
        nsc = 1
        nec = nbins - Q
        osc = abs(Q)
        oec = nbins
        temp[nec:nbins] .= 0
    end
    temp[nsc:nec] = data[osc:oec]
end

function shift(irf::TcspcTransient,Q::AbstractFloat)
    if Q != 0
        # Same thing, except rounds Q to nearest tenth
        Q = round(Int,round(Q,digits=1) / 0.1)
        # Linear IRF interpolation
        t = t_vector(irf.t)
        itp = interpolate((t,),irf.data,Gridded(Linear()))
        # Now evaluate at the new grid
        tnew = TcspcTimeAxis(irf.t.tcal / 10,irf.t.nbins * 10,irf.t.ψ)
        tgrid = t_vector(tnew)
        v = itp(tgrid)
        v = shift(v,Q)
        # Now, reverse interpolate back
        itp = interpolate((tgrid,),v,Gridded(Linear()))
        # Evaluate at the old grid
        round(Int,itp(t))
    else
        irf.data
    end

end


text:
# Run this file to download the trained data to `/data`
using HTTP
repos = ["tessdata", "tessdata_fast", "tessdata_best"]
names = ["eng", "equ"]

does_not_exist = [("tessdata_best", "equ")]

data_dir = joinpath(@__DIR__, "..", "data")
isdir(data_dir) || mkdir(data_dir)

for repo in repos
    dir = joinpath(data_dir, repo)
    isdir(dir) || mkdir(dir)

    # Download license
    path = joinpath(dir, "LICENSE")
    url = "https://raw.githubusercontent.com/tesseract-ocr/tessdata/master/LICENSE"
    HTTP.download(url, path)

    # Download `pdf.ttf`
    path = joinpath(dir, "pdf.ttf")
    url = "https://github.com/tesseract-ocr/tessconfigs/raw/master/pdf.ttf"
    HTTP.download(url, path)

    # Download data files
    for name in names
        (repo, name) ∈ does_not_exist && continue
        path = joinpath(dir, name) * ".traineddata"
        url = "https://github.com/tesseract-ocr/$repo/raw/master/$name.traineddata"
        @info "Downloading" repo name
        HTTP.download(url, path)
    end
end


text:
include(joinpath("joint_degree_dist_maxent.jl")) # the degree distribution is obtained from the joint degree distribution of maximum entropy

"""
degree_dist_maxent(S::Int64, L::Int64)
    S: number of species
    L: number of links
Returns the degree distribution of maximum entropy given S and L
"""
function degree_dist_maxent(S::Int64, L::Int64)
    joint_degree_dist = joint_degree_dist_maxent(S, L) # get the joint degree distribution of maximum entropy given S and L

    degree_dist = zeros(S+1) # create vector for the degree distribution (S+1 entries because degrees go from 0 to S)

    # compute probability for all possible degrees (degree k = in-degree + out-degree)
    # there is many combinations of in and out-degrees for a given degree
    # these combinations are on the antidiagonals of k by k submatrices
    for k in 0:S
        degree_dist[k+1] = sum(joint_degree_dist[k-i+1, i+1] for i in 0:k)
    end
    return degree_dist
end


text:
module EquilibriumMeasures
using Base, ClassicalOrthogonalPolynomials, ContinuumArrays, ForwardDiff, IntervalSets, DomainSets, StaticArrays,
    FillArrays, BlockArrays, LazyBandedMatrices

import ForwardDiff: derivative, gradient, jacobian

import LinearAlgebra: dot
import IntervalSets: mean

export equilibriummeasure, _equilibriummeasure, deflation_inner_products, deflation_deriv, deflation_op, deflation_scale

Base.floatmin(::Type{<:ForwardDiff.Dual}) = floatmin(Float64)

function _equilibriummeasure(V, ab::SVector{2})
    a,b = ab
    Typ = float(eltype(ab))
    T = ChebyshevT{Typ}()
    x = Inclusion(a..b)
    y = affine(x, axes(T,1)) # affine map from a..b to -1..1
    T̃ = T[y,:]
    wŨ = Weighted(ChebyshevU{Typ}())[y,:]
    Ũ = ChebyshevU{Typ}()[y,:]
    # Operators
    H = T̃ \ (inv.(x .- x') * wŨ) # Hilbert wU\-> T̃
    D = Ũ \ (Derivative(axes(T̃,1)) * T̃) # Derivative T̃ -> Ũ
    ŨT̃ = Ũ \ T̃  # Converion T̃ -> Ũ
    V_cfs = T̃ \ V.(x)
    Vp_cfs = ŨT̃ \ (D * V_cfs)
    Vp_cfs[1], wŨ * (H[2:end,:] \ Vp_cfs[2:end])/2
end

function _equilibriummeasure(V, ad::SVector{4})
    a,b,c,d = ad
    Typ = float(eltype(ad))
    T1,T2 = chebyshevt(a..b),chebyshevt(c..d)
    U1,U2 = chebyshevu(a..b),chebyshevu(c..d)
    T = PiecewiseInterlace(T1, T2)
    U = PiecewiseInterlace(U1, U2)
    W = PiecewiseInterlace(Weighted(U1), Weighted(U2))
    x = axes(W,1)
    # Operators
    D = U \ (Derivative(x) * T) # Derivative T -> U
    UT = U \ T  # Converion T -> U
    V_cfs = T \ V.(x)
    Vp_cfs = UT \ (D * V_cfs)

    H = T \ (inv.(x .- x') * W) # Hilbert W -> T

    # add extra columns for constants
    H̃ = BlockHcat(Eye((axes(H,1),))[:,Block(1)], 2H)
    c = H̃ \ Vp_cfs
    c[Block(1)], W * c[Block.(2:∞)]
end

# Intentionally hide type for compile time
struct DFunction <: Function
    f
end
(f::DFunction)(x::T) where T = f.f(x)::T
Base.convert(::Type{DFunction}, f::DFunction) = f
Base.convert(::Type{DFunction}, f::Function) = DFunction(f)

struct EquilibriumMeasureMoment
    V::DFunction
end

_logterms(V, μ) = ()
_logterms(V, μ, d) = ()
function _logterms(V, μ, d1, d2)
    x = axes(μ,1)
    z1,z2 = mean(d1),mean(d2)
    (2*(log.(abs.(z1 .- x'))*μ) - V(z1) - 2*(log.(abs.(z2 .- x'))*μ) + V(z2),)
end

function (E::EquilibriumMeasureMoment)(a)
    c0,μ = _equilibriummeasure(E.V, a)
    SVector(c0..., sum(μ) - 1, _logterms(E.V, μ, components(axes(μ,1).domain)...)...)
end

equilibriummeasure(V; a = SVector(-1.0,1.0), maxiterations=1000, knownsolutions=[], power=2, shift=1.0, dampening=1.0, returnendpoint=false) =
    deflate_equilibriummeasure(EquilibriumMeasureMoment(DFunction(V)), a, maxiterations, knownsolutions, power, shift, dampening, returnendpoint)
function deflate_equilibriummeasure(μ, a, maxiterations, knownsolutions, power, shift, dampening, returnendpoint)
    num_found_sols = length(knownsolutions)
    for k=1:maxiterations
        update = - jacobian(μ, a)\μ(a)
        if num_found_sols > 0
            update = deflation_scale(update, a, knownsolutions, power, shift) * update
        end
        an = a + dampening*update
        if an ≈ a
            # improve accuracy a bit more
            for _=1:3
                update = - jacobian(μ, a)\μ(a)
                if num_found_sols > 0
                    update = deflation_scale(update, a, knownsolutions, power, shift) * update
                end
                a = a + dampening*update
            end
            if returnendpoint
                return  _equilibriummeasure(μ.V, a)[2], a
            else
                return  _equilibriummeasure(μ.V, a)[2]
            end
        end
        a = an
    end
    error("Max its")
end

# Deflation code

function  deflation_inner_products(state::AbstractVector, sol)
    num_sols = length(sol)
    T = promote_type(eltype(state), eltype(eltype(sol)))

    inner_products = zeros(T, num_sols)
    for iter = 1:num_sols
        inner_products[iter] = dot(state.-sol[iter], state.-sol[iter])
    end
    inner_products
end

function deflation_op(state::AbstractVector, sol, power::Real, shift::Real)
    num_sols = length(sol)
    T = promote_type(eltype(state), eltype(eltype(sol)))
    power = T(power); shift = T(shift)

    m = one(T)
    inner_products = deflation_inner_products(state, sol)

    for iter = 1:num_sols
        normsq = inner_products[iter]
        factor = one(T)/normsq^(power/2) + shift
        m = m*factor
    end
    m, inner_products
end

function deflation_deriv(update::AbstractVector, state::AbstractVector, sol::Vector, power::Real, shift::Real)
    m, inner_products = deflation_op(state, sol, power, shift)
    num_sols = length(sol)
    T = promote_type(eltype(state), eltype(eltype(sol)))
    power = T(power); shift = T(shift)

    dm = zero(T)
    state_dot_update = dot(state, update)
    for iter = 1:num_sols
        scale = m/deflation_op(state, sol[iter], power, shift)[1]
        deriv = -power*state_dot_update/inner_products[iter]^((power + 1)/2)
        dm += scale*deriv
    end
    m, dm
end

function deflation_scale(update::AbstractVector, state::AbstractVector, sol::Vector, power::Real, shift::Real)
    T = promote_type(eltype(state), eltype(eltype(sol)))
    m, dm = deflation_deriv(update, state, sol, power, shift)
    minv = one(T)/m
    tau = one(T) + minv*dm/(one(T) - minv*dm)
    tau
end

# function equilibriummeasuresupport(V,ab=(-1.0..1.0);maxiterations=100,bounded=:none)
#     if bounded == :left
#         return equilibriummeasuresupportbounded(false,V,ab;maxiterations=maxiterations)
#     elseif bounded == :right
#         return equilibriummeasuresupportbounded(true,V,ab;maxiterations=maxiterations)
#     end

#     a,b = endpoints(ab)

#     # Newton iteration, using dual numbers
#     for k=1:maxiterations
#         F1 = F(V,dual(a,1.),b)
#         p = realpart.(F1)
#         J=[epsilon.(F1) epsilon.(F(V,a,dual(b,1)))]

#         an,bn=sort([a,b]-J\p)
#         if isapprox(an,a) && isapprox(bn,b)
#             return (an..bn)
#         else
#             a,b=an,bn
#         end
#     end

#     warn("maxiterations reached")
#     a..b
# end

# Fbounded(s::Bool,V,a,b) = emmoment(Fun(V,(a..b))'; bounded= s ? (:right) : (:left))-1

# function equilibriummeasuresupportbounded(s::Bool,V,ab=(-1.0..1.0);maxiterations=100)
# a,b=endpoints(ab)
# for k=1:maxiterations
#     F1=s ? Fbounded(s,V,a,dual(b,1.)) : Fbounded(s,V,dual(a,1.),b)

#     p=realpart(F1)
#     J=epsilon(F1)

#     if s
#         bn=b-J\p
#         if isapprox(bn,b)
#             return (a..bn)
#         else
#             b=bn
#         end
#     else
#         an=a-J\p
#         if isapprox(an,a)
#             return (an..b)
#         else
#             a=an
#         end
#     end
# end
# warn("maxiterations reached")
# a..b
# end

# # we assume that the zeroth coefficient is zero
# equilibriummeasure(V,gs...;bounded=:none,opts...) =
# -hilbertinverse(Fun(V,equilibriummeasuresupport(V,gs...;bounded=bounded,opts...))';bounded=bounded)/(2π)

end #module

max_stars_repo_name: ApproxFun/EquilbriumMeasures

max_issues_repo_name: ApproxFun/EquilbriumMeasures

max_forks_repo_name: ApproxFun/EquilbriumMeasures


text:
using MathGL.Capi
mgl = MathGL.Capi

gr = mgl.create_graph(800, 500)
dat = mgl.create_data()
mgl.data_link(dat, 0.8*sin(linspace(-4pi, 4pi, 200)), 200, 1, 1)
mgl.label(gr, 'x', "x", 0, "")
mgl.label(gr, 'y', "y", 0, "")
mgl.box(gr)
mgl.axis(gr, "xyz", "", "")
mgl.axis_grid(gr, "xyz", "H|", "")
#=mgl.set_origin(gr, 0., 0., 0.)=#
mgl.plot(gr, dat, "", "")
mgl.write_frame(gr, "../graphs/sin_C.png", "")


text:
module BenchmarkSets

using BenchmarkTools
export @bench, @benchmarkset

const benchmark_stack = []

macro benchmarkset(title, ex)
    esc(benchmarkset_m(title, ex))
end

macro bench(title, xs...)
    esc(:(var"#suite#"[$title] = @benchmarkable $(xs...)))
end

function benchmarkset_m(title, ex::Expr)
    stack = GlobalRef(BenchmarkSets, :benchmark_stack)
    init = quote
        if isempty($stack)
            push!($stack, $BenchmarkGroup())
        end
    end
    exec = quote
        if length($stack) == 1
            pop!($stack)
        end
    end
    return if ex.head === :block
        quote
            $init
            $(benchmarkset_block(title, ex))
            $exec
        end
    elseif ex.head === :for
        quote
            $init
            $(Expr(ex.head, ex.args[1], benchmarkset_block(title, ex.args[2])))
            $exec
        end
    end
end

function benchmarkset_block(title, ex::Expr)
    stack = GlobalRef(BenchmarkSets, :benchmark_stack)
    quote
        let var"#root#" = last($stack)
            var"#root#"[$title] = var"#suite#" = BenchmarkGroup()
            push!($stack, var"#suite#")
            $ex
            pop!($stack)
        end
    end
end

end


text:
using Distributions

const GA_DefaultOptions = ParamsDict(
    :CrossoverOperator => MultiPointCrossover(2),
    :MutationOperator => BinaryFlipMutation(.01),
    :EmbeddingOperator => NoEmbedding(),
    :Selector => SimpleSelector()
)

mutable struct GeneticAlgorithm{P<:Population,
                                S<:IndividualsSelector,
                                X<:CrossoverOperator,
                                M<:MutationOperator,
                                E<:EmbeddingOperator} <: EvolutionaryOptimizer
    name::String
    population::P
    so::S
    xo::X        # xover operator
    mo::M        # mutation operator
    eo::E         # embedding operator

    function GeneticAlgorithm(name, pop::P, so::S, xo::X, mo::M, eo::E) where {P,S,X,M,E}
        return new{P,S,X,M,E}(name, pop, so, xo, mo, eo)
    end
end

function ask(ga::GeneticAlgorithm)
    selected_indices = select(ga.so, population(ga), numparents(ga.xo) + numchildren(ga.xo))

    @inbounds parent_indices = selected_indices[1:numparents(ga.xo)]
    @inbounds target_indices = selected_indices[(numparents(ga.xo) + 1):end]

    targets = map( ti -> acquire_candi(population(ga), ti), target_indices)
    trials = map( t -> acquire_candi(population(ga), t), targets)
    trial_params = map( t -> view(t.params, :), trials)

    # xover
    apply!(ga.xo, trial_params, target_indices, population(ga), parent_indices)
    # mutations
    foreach(idx -> apply!(ga.mo, trial_params[idx], target_indices[idx]), eachindex(target_indices))

    return post_process(ga, targets, trials, ga.xo)
end

# implementation and integration into bbo

function ga(problem::OptimizationProblem,
                            options::Parameters = EMPTY_PARAMS,
                            name = "GA/ga")
    opts = chain(GA_DefaultOptions, options)
    pop = get(opts, :Population, population(problem, options))
    so = opts[:Selector]
    xo = opts[:CrossoverOperator]
    mo = opts[:MutationOperator]
    eo = opts[:EmbeddingOperator]

    return GeneticAlgorithm(name, pop, so, xo, mo, eo)
end