text:
module StringMethod

using LinearAlgebra
using Dierckx
using Printf

# types
include("types.jl")

# utility functions
include("utils.jl")
export StringOptions, SaddleOptions, periodic_dist, linear_string, upwind_tangent

# integration routines
include("integrate.jl")
export stepEuler!, stepRK4!, stepPCEuler!

# reparametrization routines
include("reparameterize.jl")
export  spline_reparametrize!

# string method routines
include("string.jl")
include("pcstring.jl")
export SimplifiedString, PCSimplifiedString, simplified_string, simplified_string!

# saddle point method routines
include("saddle.jl")
export ClimbingImage, climbing_image, climbing_image!
end #end module


text:
using Test
using BenchmarkTools

using iLQGames:
    iLQGames,
    Unicycle4D,
    NPlayerUnicycleCost,
    quadraticize,
    _quadraticize_ad,
    linearize_discrete,
    n_states,
    n_controls,
    dynamics,
    player_costs,
    SystemTrajectory,
    lq_approximation,
    iLQSolver,
    solve!,
    AffineStrategy,
    generate_nplayer_navigation_game,
    samplingtime,
    horizon

using iLQGames.TestUtils

using StaticArrays
using LinearAlgebra

# generate a game
T_horizon = 10.
ΔT = 0.1
# initial conditions:
# x = (x, y, phi, β, v)
x01 = @SVector [-3., 0., 0., 0.]
x02 = @SVector [0.,  3., -pi/2, 0.]
x0 = vcat(x01, x02)
# goal states (goal position of other player with opposite orientation)
xg1 = @SVector [3., 0., 0., 0.]
xg2 = @SVector [0., -3., -pi/2, 0.]
# generate game
g = generate_nplayer_navigation_game(Unicycle4D, NPlayerUnicycleCost, T_horizon, ΔT, xg1, xg2)

# unpack for testing
dyn = dynamics(g)
nx = n_states(dyn)
nu = n_controls(dyn)
c1, c2 = player_costs(g)

# test quadraticization of the cost and quadratization of cost
x = @SVector zeros(nx)
u = @SVector zeros(nu)
t = 0.
# solve the lq game
solver = iLQSolver(g)

quad_sanity_check(g)

# test the lq approximation:
# generate an operating point
h = Int(T_horizon/ΔT)
zero_op = zero(SystemTrajectory{h, ΔT, nx, nu})
lqg = lq_approximation(solver, g, zero_op)

# the lqg approximation evaluated at zero should be approximate the true cost:

# - setup initial_strategy
steer_init(k::Int) = cos(k/h*pi) * deg2rad(0)
acc_init(k::Int) = -cos(k/h*pi)*0.3
γ_init = SizedVector{h}([AffineStrategy((@SMatrix zeros(nu, nx)),

                                 (@SVector [steer_init(k), 0.7*acc_init(k),
                                            steer_init(k), acc_init(k)])) for k in 1:h])
# generate initial operating point from simulating initial strategy
# solve the game
@benchmark(solve!(copy(zero_op), copy(γ_init), $g, $solver, $x0)) |> display

_, traj, _ = solve!(copy(zero_op), copy(γ_init), g, solver, x0)
@testset "solution trajectory sanity checks" begin
    @test horizon(traj) == h
    @test samplingtime(traj) == ΔT
end

ext: jl


text:
using PQEcolaPoint
using Dates

# setup DTE -
T = 45.0
DTE = range(T, stop=0.0,step=-1.0) |> collect

# setup underlying range -
Sₒ = 79.74
IV = 0.4456
σ = Sₒ*(IV)*(√(T/365))
L = round(Sₒ - 2.96*σ)
U = round(Sₒ + 2.96*σ)
S = range(L,stop=U, step=1.0) |> collect

# setup IV array -
IV_array = range(0.3, stop = 0.5, step = 0.1) |> collect

# setup strike coordinates -
K = range(80.0, stop = 125.0, step = 5.0) |> collect

# compute the grid -
lattice_model = build(CRRContractPremiumLatticeModel, PutContractModel, S, K, IV_array, DTE; number_of_levels = 80)

# # try to save -
path_to_save = "./test/tmp/Test.bson"
save(path_to_save, lattice_model)


text:
"""
    Options

Options configurable by the user with recommended default values.
"""
@with_kw struct Options
    num_start_samples::Int = 5
    trace::Bool = true
    sampling_plan_opt_gens::Int = 5000
    rippa::Bool = true
    kerns = [ScatteredInterpolation.Gaussian, ScatteredInterpolation.InverseQuadratic,
             ScatteredInterpolation.InverseMultiquadratic]
    rbf_opt_gens::Int = 50_000
    rbf_opt_pop::Int = 50
    rbf_opt_method::Symbol = :de_rand_1_bin_radiuslimited
    rbf_dist_metric = Distances.Euclidean()
    variable_kernel_width::Bool = true
    variable_dim_scaling::Bool = true
    cond_max::Float64 = 1e4
    max_rbf_width::Float64 = 1000.0
    min_rbf_width::Float64 = 1e-4
    max_scale::Float64 = 10.0
    min_scale::Float64 = 1e-4
    num_interpolants::Int = 20
    smooth = :single
    max_smooth::Float64 = 1.0
    smooth_user::Float64 = 0.0
    iterations::Int64 = 10
    num_infill_points::Int64 = 1
    parallel_surrogate::Bool = true
    infill_funcs::Array{Symbol,1} = [:median,:std]
    infill_iterations::Int64 = 25_000

    @assert ((smooth == false) || (smooth == :variable) ||
    (smooth == :single) || (smooth == :single_user)) "Not supported option for smooth"
end

"""
    RBFHypers(width::S,kernelFunc,scaling::U,fitness::Float64)

Data structure to store results from the optimization of an RBF interpolation kernel
"""
struct RBFHypers{T,U}
    kernelFunc
    scaling::T
    smooth::U
end

"""
    SurrogateResult

Data structure to store the results of an optimization task.
"""
struct SurrogateResult
    lhc_samples::Array{Float64,2}
    lhc_plan::Array{Float64,2}
    sm_interpolant
    sm_interpolant_settings::Array{SurrogateModelOptim.RBFHypers{T,U},1} where U where T
    infill_samples::Array{Float64,2}
    infill_type::Array{Symbol,1}
    infill_plan::Array{Float64,2}
    infill_prediction::Array{Float64,1}
    options::Options
end


text:
type LSTMLayer # a single layer of a multi-layer RNN
        # cell's input gate params
        wix::NNMatrix
        wih::NNMatrix
        bi::NNMatrix
        # cell's forget gate parameters
        wfx::NNMatrix
        wfh::NNMatrix
        bf::NNMatrix
        # cell's out gate parameters
        wox::NNMatrix
        woh::NNMatrix
        bo::NNMatrix
        # cell's write parameters
        wcx::NNMatrix
        wch::NNMatrix
        bc::NNMatrix
    function LSTMLayer(prevsize::Int, hiddensize::Int, std::Float64)
        ###  gate parameters ###
        # cell's input gate params
        wix = randNNMat(hiddensize, prevsize, std)
        wih = randNNMat(hiddensize, hiddensize, std)
        bi = NNMatrix(hiddensize, 1, zeros(hiddensize,1), zeros(hiddensize,1))

        # cell's forget gate parameters
        wfx = randNNMat(hiddensize, prevsize, std)
        wfh = randNNMat(hiddensize, hiddensize, std)
        bf = NNMatrix(hiddensize, 1, zeros(hiddensize,1), zeros(hiddensize,1))

        # cell's out gate parameters
        wox = randNNMat(hiddensize, prevsize, std)
        woh = randNNMat(hiddensize, hiddensize, std)
        bo = NNMatrix(hiddensize, 1, zeros(hiddensize,1), zeros(hiddensize,1))

        ### cell's write parameters ###
        wcx = randNNMat(hiddensize, prevsize, std)
        wch = randNNMat(hiddensize, hiddensize, std)
        bc = NNMatrix(hiddensize, 1, zeros(hiddensize,1), zeros(hiddensize,1))

        new(wix,wih,bi,　wfx,wfh,bf,　wox,woh,bo,　wcx,wch,bc)
    end
end

type LSTM <: Model
    hdlayers::Array{LSTMLayer,1} # holds variable number of hidden layers
    whd::NNMatrix # hidden to decoder weights & gradients
    bd::NNMatrix  # bias of hidden to decoder layer
    matrices::Array{NNMatrix,1} # used by solver - holds references to each of matrices in model
    hiddensizes::Array{Int,1}
    function LSTM(inputsize::Int, hiddensizes::Array{Int,1}, outputsize::Int, std::Float64=0.08)
        hdlayers = Array(LSTMLayer, length(hiddensizes))
        matrices = Array(NNMatrix, 0)
        for d in 1:length(hiddensizes)
            prevsize = d == 1? inputsize : hiddensizes[d-1]
            layer = LSTMLayer(prevsize, hiddensizes[d], std)
            hdlayers[d] = layer
            # input gate
            push!(matrices, layer.wix)
            push!(matrices, layer.wih)
            push!(matrices, layer.bi)
            # forget gate
            push!(matrices, layer.wfx)
            push!(matrices, layer.wfh)
            push!(matrices, layer.bf)
            # output gate
            push!(matrices, layer.wox)
            push!(matrices, layer.woh)
            push!(matrices, layer.bo)
            # cell params
            push!(matrices, layer.wcx)
            push!(matrices, layer.wch)
            push!(matrices, layer.bc)
        end
        whd = randNNMat(outputsize, hiddensizes[end], std)
        bd = NNMatrix(outputsize, 1, zeros(outputsize,1), zeros(outputsize,1))
        push!(matrices, whd) # bias of hidden to hidden
        push!(matrices, bd) # bias of hidden to hidden
#         println("bh.w size($(size(bd.w)))")
        new(hdlayers, whd, bd, matrices, hiddensizes)
    end
end

function forwardprop!(g::Graph, model::LSTM, x, prev)

    # forward prop for a single tick of LSTM
    # g is graph to append ops to
    # model contains LSTM parameters
    # x is 1D column vector with observation
    # prev is a struct containing hidden and cell from previous iteration

    prevhd, prevcell, _ = prev # extract previous hidden and cell states from the tuple
    hiddenprevs = Array(NNMatrix,0)
    cellprevs = Array(NNMatrix,0)
    if length(prevhd) == 0
        for hdsize in model.hiddensizes
            push!(hiddenprevs, NNMatrix(hdsize,1))
            push!(cellprevs, NNMatrix(hdsize,1))
        end
    else
        hiddenprevs = prevhd
        cellprevs = prevcell
    end

    hidden = Array(NNMatrix,0)
    cell = Array(NNMatrix,0)
    for d in 1:length(model.hiddensizes) # for each hidden layer

        input = d == 1 ? x : hidden[d-1]
        hdprev = hiddenprevs[d]
        cellprev = cellprevs[d]
        # cell's input gate params
        wix = model.hdlayers[d].wix
        wih = model.hdlayers[d].wih
        bi = model.hdlayers[d].bi
        # cell's forget gate parameters
        wfx = model.hdlayers[d].wfx
        wfh = model.hdlayers[d].wfh
        bf = model.hdlayers[d].bf
        # cell's out gate parameters
        wox = model.hdlayers[d].wox
        woh = model.hdlayers[d].woh
        bo = model.hdlayers[d].bo
        # cell's write parameters
        wcx = model.hdlayers[d].wcx
        wch = model.hdlayers[d].wch
        bc = model.hdlayers[d].bc

        # input gate
        h0 = mul!(g, wix, input)
        h1 = mul!(g, wih, hdprev)
        inputgate = sigmoid!(g, add!(g, h0, h1, bi))

        # forget gate
        h2 = mul!(g, wfx, input)
        h3 = mul!(g, wfh, hdprev)
        forgetgate = sigmoid!(g, add!(g, h2, h3, bf))

        # output gate
        h4 =mul!(g, wox, input)
        h5 = mul!(g, woh, hdprev)
        outputgate = sigmoid!(g, add!(g, h4, h5, bo))

        # write operation on cells
        h6 = mul!(g, wcx, input)
        h7 = mul!(g, wch, hdprev)
        cellwrite = tanh!(g, add!(g, h6, h7, bc))

        # compute new cell activation
        retaincell = eltmul!(g, forgetgate, cellprev) # what do we keep from cell
        writecell = eltmul!(g, inputgate, cellwrite)  # what do we write to cell
        cell_d = add!(g, retaincell, writecell)        # new cell contents

        # compute hidden state as gated, saturated cell activations
        hidden_d = eltmul!(g, outputgate, tanh!(g, cell_d))

        push!(hidden,hidden_d)
        push!(cell, cell_d)
    end

    # one decoder to outputs at end
    output = add!(g, mul!(g, model.whd, hidden[end]),model.bd)

    # return cell memory, hidden representation and output
    return hidden, cell, output
end

max_forks_repo_name: haroldsoh/RecurrentNN


text:
include("src/model.jl")
using Printf
using PyPlot
using PyCall
using JLD2,FileIO

K = round.(Int, 10 .^range(log10(500), log10(20000), length=10))[1:7]
scalegs = [-0.25, -0.5]

ϕ̄ = zeros(length(K),2)
f = zeros(length(K),2)
μ = zeros(length(K),2)

for k = 1:2
   v.dp.scaleg = scalegs[k]
   if k==2
      v.dp.J_ii = 1.05
   end
   for i = 1:length(K)
      v.dp.K = K[i]
      v.dp.N = K[i]
      v.dp.rcNsamp = min(1000, v.dp.N)

      filename = v.genfilename(v.dp)

      if !isfile(filename)
         D,F = v.runsim()
         v.writefile(D, v.dp)
      else
         D = v.readfile(v.dp)
      end

      ϕ̄[i,k] = D[:m]
      f[i,k] = D[:f̄]
      μ[i,k] = D[:μ̄]
   end
end

lowvar  = "#b2581e"
highvar = "#1e78b2"

fig, ax = plt.subplots(1,3,figsize=(19/2.54,7/2.54))
mticker = pyimport("matplotlib.ticker")

ax[1].plot(K, ϕ̄[:,1], color=highvar, label="high var")
ax[1].plot(K, ϕ̄[:,2], color=lowvar, linestyle="dashed", label="low var")
ax[1].set(ylabel=L"mean response ($\overline{\phi}$)")
ax[1].legend(frameon=false)

ax[2].plot(K, f[:,1], color=highvar)
ax[2].plot(K, f[:,2], color=lowvar)
ax[2].set(ylabel=L"fraction active ($f$)")

ax[3].plot(K, μ[:,1], color=highvar)
ax[3].plot(K, μ[:,2], color=lowvar)
ax[3].set(ylabel=L"mean active response ($\mu$)")

ax[1].set_xscale("log");ax[1].set_yscale("log")
ax[2].set_xscale("log");ax[2].set_yscale("log")
ax[3].set_xscale("log");ax[3].set_yscale("log")

ax[1].yaxis.set_minor_formatter(mticker.NullFormatter())
ax[2].yaxis.set_minor_formatter(mticker.NullFormatter())
ax[3].yaxis.set_minor_formatter(mticker.NullFormatter())

ax[1].set_yticks([0.01,0.03])
ax[2].set_yticks([0.1,0.3])
ax[3].set_yticks([0.01,0.1])

ax[1].get_yaxis().set_major_formatter(mticker.ScalarFormatter())
ax[2].get_yaxis().set_major_formatter(mticker.ScalarFormatter())
ax[3].get_yaxis().set_major_formatter(mticker.ScalarFormatter())

for i in 1:length(ax)
   ax[i].set_xticks([1000,10000])
   ax[i].set(xlabel=L"in-degree ($K$)")
   ax[i].spines["right"].set_visible(false)
   ax[i].spines["top"].set_visible(false)
   ax[i].xaxis.set_ticks_position("bottom")
   ax[i].yaxis.set_ticks_position("left")
end

tight_layout()
savefig("plots/Kscalings.pdf")

###############################################################
fig, ax = plt.subplots(1,3,figsize=(19/2.54,7/2.54))

v.dp.scaleg = -0.5
v.dp.J_ii   = 1.05
v.dp.N = 1000
v.dp.K = 1000
v.dp.T = 500
D, F = v.runsim()

times = 0 : v.dp.rcDt : (v.dp.T-v.dp.rcDt)
ax[1].plot(times, F[:ϕs][:,1:7], alpha=0.4)
ax[3].hist(v.unroll(F[:xs]), 100, alpha=0.7, density=true, color=lowvar, label="low var")
ax[3].axvline(x=D[:x̄], color=lowvar, linestyle="dashed", alpha=0.7)

v.dp.scaleg = -0.25
v.dp.J_ii   = 2.0
v.dp.N = 1000
v.dp.K = 1000
v.dp.T = 500
D, F = v.runsim()

ax[2].plot(times, F[:ϕs][:,1:7], alpha=0.4)
ax[3].hist(v.unroll(F[:xs]), 100, alpha=0.7, density=true, color=highvar, label="high var")
ax[3].axvline(x=0,color="gray")
ax[3].axvline(x=D[:x̄], color=highvar, linestyle="dashed", alpha=0.7)

ax[1].set(title="low var", xlabel="time",ylabel="response")
ax[2].set(title="high var",xlabel="time",ylabel="response")

ax[1].set_ylim([0,0.6])
ax[2].set_ylim([0,0.6])
ax[3].set_xlim([-4,1])
ax[3].legend(frameon=false)

for i in 1:length(ax)
   ax[i].spines["right"].set_visible(false)
   ax[i].spines["top"].set_visible(false)
   ax[i].xaxis.set_ticks_position("bottom")
   ax[i].yaxis.set_ticks_position("left")
end

tight_layout()
savefig("plots/activities.pdf")

max_stars_repo_name: raminkhajeh/sparse-balance

max_issues_repo_name: raminkhajeh/sparse-balance

max_forks_repo_name: raminkhajeh/sparse-balance


text:
using Parameters
@with_kw struct InterpRange{T}
    xmin::T
    xmax::T
    dx_inv::T
    ymin::T
    ymax::T
    dy_inv::T
end

function interp(x::T,y::T, tbl::Matrix{T}, par::InterpRange{T}) where {T}
    @unpack xmin,xmax,dx_inv,ymin,ymax,dy_inv = par
    x = x >= xmax ? (1-1e-7)*xmax : x
    x = x <= xmin ? xmin : x
    y = y >= ymax ? (1-1e-7)*ymax : y
    y = y <= ymin ? ymin : y
    xi = (x - xmin) * dx_inv
    x0 = Int64(floor(xi))
    yi = (y - ymin) * dy_inv
    y0 = Int64(floor(yi))
    wx1 = xi - x0
    wy1 = yi - y0
    wx0 = 1.0 - wx1
    wy0 = 1.0 - wy1
    f00 = tbl[x0+1,y0+1] #julia array is 1-based
    f10 = tbl[x0+2,y0+1]
    f01 = tbl[x0+1,y0+2]
    f11 = tbl[x0+2,y0+2]

    fi = wx0*wy0 * f00 + wx0*wy1 * f01 + wx1*wy0 * f10 + wx1*wy1 * f11
end


text:
corr_from_chol(w) = w*w'

function chol_from_vec(x, K)
    w = vec_to_lt(x, K)
    for i in 1:K
        if sum(w[i,:].^2) > 1
            w[i,i] = 0.0
        else
            w[i,i] = sqrt(1 - sum(w[i,:].^2))
        end
    end
    LowerTriangular(w)
end

## transform cholesky factor (vector repr) to KxK corr matrix (vector repr)
function chol_vec_to_corr_vec(y, K)
    lt_to_vec(LowerTriangular(corr_from_chol(chol_from_vec(y, K))))
end

function log_jacobian_det_chol_to_corr(w::LowerTriangular)
    K = size(w,1)
    sum((K .- (2:K)).*log.(diag(w)[2:K]))
end

function log_jacobian_det_free_to_corr(y, K)
    w = free_to_chol(y, K)
    log_jacobian_det_free_to_chol(y, K) + log_jacobian_det_chol_to_corr(w)
end

function lt_to_vec(w::LowerTriangular)
    K = size(w)[1]
    x = Vector{eltype(w)}(undef, Int(K*(K-1)/2))
    k = 1
    for i in 2:K
        for j in 1:(i-1)
            x[k] = w[i,j]
            k += 1
        end
    end
    x
end

function vec_to_lt(y, K)
    w = zeros(eltype(y), K, K)
    k = 1
    for i in 2:K
        for j in 1:(i-1)
            w[i,j] = y[k]
            k += 1
        end
    end
    LowerTriangular(w)
end

ut_to_vec(w::UpperTriangular) = lt_to_vec(LowerTriangular(w'))
vec_to_ut(y, K) = UpperTriangular(vec_to_lt(y, K)')

function corr_to_free(w::UpperTriangular)
    K = size(w,1)

    z = UpperTriangular(zeros(eltype(w), K, K))
    z[1,2:end] = w[1,2:end]

    for i in 2:K
        for j in (i+1):K
            ##z[i,j] = w[i,j] / sqrt(prod(1 .- z[1:(i-1),j].^2))
            z[i,j] = w[i,j] / sqrt(1 - sum(w[1:(i-1),j].^2))
        end
    end
    atanh.(z)
end

function free_to_corr(y::UpperTriangular)
    z = tanh.(y)

    K = size(z)[1]
    w = zeros(eltype(y), K, K)

    w[1,1] = 1
    w[1,2:end] = z[1,2:end]

    for i in 2:K
        w[i,i] = sqrt(prod(1 .- z[1:(i-1),i].^2))
        for j in (i+1):K
            w[i,j] = z[i,j]*sqrt(prod(1 .- z[1:(i-1),j].^2))
        end
    end

    UpperTriangular(w) ##w'*w
end

## p_y(y) \propto |J|p_w(w) \propto |J'|p_x(x)
function logpdf_free(y::Vector, K::Int, eta::Float64)
    w = free_to_corr(y, K)
    d = LKJcorrChol(K, eta)
    logpdf(d, w) + log_jacobian_det_free_to_chol(y, K)
end


text:
mutable struct pbtfProblem{T<:Real}
    y       :: Vector{T}               # response
    x       :: Vector{T}               # predictor
    k       :: Int                     # order of desired polynomial, only 1 or 2 is allowed
    λ       :: Union{T,Nothing}        # Moreau-Yosida envelope parameter
    s       :: T                       # IG(s, r) prior for σ²
    r       :: T
    s₂      :: Union{T,Nothing}        # β'(n-k,s₂) prior for α
    m       :: Int                     # total number of observations
    n       :: Int                     # number of grid locations
    xgrid   :: Vector{T}               # strictly increasing grid locations (length n)
    w       :: Vector{T}               # number of observations at each location (length n)
    ybar    :: Vector{T}               # the average of y at each location (length n)
    sse     :: T                       # sum of squared errors
    TM      :: SparseMatrixCSC{T,Int}  # n×n transformation matrix
    TMᵀ     :: SparseMatrixCSC{T,Int}  # n×n transformation matrix transposed
    D       :: SparseMatrixCSC{T,Int}  # (n-k-1)×n difference matrix
    inv_DDᵀ_D :: Matrix{T}             # used to determine the bisection interval
    θ_proj  :: Vector{T}               # the following arrays are pre-created containers to avoid creating
    β       :: Vector{T}               # new arrays during the course of MCMC sampling
    res     :: Vector{T}
end

function pbtfProblem(y::AbstractArray{T}, x::AbstractArray{T}, k::Int;
                     thinning::Bool=false,nbins::Int=100,λ::Union{T,Nothing}=nothing,
                     s::T=1e-3,r::T=1e-3,s₂::Union{T,Nothing}=nothing) where T<: Real
    if (k!=1) & (k!=2)
        throw(ArgumentError("Currenly only k=1 and k=2 are accepted for pbtf."))
    end
    m = length(x)
    xgrid = sort(unique(x))
    n = length(xgrid)
    if ((k==1) & (n>1000)) || ((k==2) & (n>200))
        @warn "Transformation matrix is likely too ill-conditioned, setting thinning=true is recommended."
    end
    w = zeros(T,n)
    ybar = zeros(T,n)
    sse = 0.0
    for i in 1:n
        w[i] = length(x[x.==xgrid[i]])
        ybar[i] = mean(y[x.==xgrid[i]])
        sse += sum(abs2, y[x.==xgrid[i]].-ybar[i])
    end
    if thinning
        xgrid,w,ybar,sse = thin(xgrid,w,ybar,sse,nbins)
        n = nbins
    end
    if λ==nothing
        λ = min(1/n^2,1e-4*var(y))
    end
    if s₂==nothing
        s₂ = sqrt(n)
    end
    TM = getTM(k,n,xgrid)
    TMᵀ = SparseMatrixCSC(TM')
    D1 = getD(0,n-k)
    inv_DDᵀ_D = inv(Matrix(D1*D1'))*D1
    D = getD(k,n,xgrid)
    θ_proj = zeros(T,n)
    β = zeros(T,n)
    res = zeros(T,n)
    pbtfProblem{T}(y,x,k,λ,s,r,s₂,m,n,xgrid,w,ybar,sse,TM,TMᵀ,D,inv_DDᵀ_D,θ_proj,β,res)
end

dimension(problem::pbtfProblem) = problem.n + 2
capabilities(::Type{<:pbtfProblem}) = LogDensityOrder{1}()

function logdensity_and_gradient(problem::pbtfProblem, z)
    @unpack y,x,k,λ,s,r,s₂,m,n,xgrid,w,ybar,sse,TM,TMᵀ,D,inv_DDᵀ_D,θ_proj,β,res=problem
    θ         = view(z,1:n)
    logσ²     = z[n+1]
    σ²        = exp(logσ²)
    logα      = z[n+2]
    α         = exp(logα)
    α_proj    = 0.0
    if (k==1) & (n<=200)
        α_proj = absepi_projection!(θ_proj,k+1,θ,α)
        forward_solve!(β,k+2,TM,θ)
    else
        α_proj = diffepi_projection!(θ_proj,k,θ,α,inv_DDᵀ_D)
        forward_solve!(β,k+1,TM,θ)
    end
    res      .= ybar .- β
    ssm       = sum(res .* w .* res)
    qf        = (ssm+sse+2*r)/(2*σ²)
    dist²     = norm(θ.-θ_proj,2)^2+(α-α_proj)^2
    logl      = -dist²/(2*λ)-qf-(m/2+s)*logσ²
    ∇         = similar(z)
    if (k==1) & (n<=200)
        back_solve!(view(∇,1:n),k+2,TMᵀ,w.*res)
    else
        back_solve!(view(∇,1:n),k+1,TMᵀ,w.*res)
    end
    logl += logα-(n-k+s₂)*log(1+α)
    ∇[n+2]  = 1 -α/λ*(α-α_proj)-(n-k+s₂)*α/(1+α)
    ∇[1:n]  ./= σ²
    ∇[1:n]  .-= (θ.-θ_proj)./λ
    ∇[n+1]    = qf-(m/2+s)
    logl, ∇
end


text:
abstract type ExtensibleUnion end

members(::Type{U}) where {U<:ExtensibleUnion} =
    collect(ExtensibleUnions._registry_extensibleunion_to_members[U])[2:end]

function Base.show(io::IO, ::Type{U}) where {U <: ExtensibleUnion}
    membs = members(U)
    print(io, "Extensible Union $((string ∘ nameof)(U)) with ")
    if isempty(membs)
        print(io, "no members")
    else
        print(io, "members: $(membs[1])")
        for memb in membs[2:end]
            print(io, ", $(memb)")
        end
    end
    print(io, ".")
end

macro ExtensibleUnion(T)
    quote
        struct $T <: ExtensibleUnion end
        extensibleunion!($T)
    end |> esc
end

function extensibleunion!(@nospecialize(u))
    global _registry_extensibleunion_to_genericfunctions
    global _registry_extensibleunion_to_members
    if !isconcretetype(u)
        throw(ArgumentError("The provided type must be a concrete type"))
    end
    if !isstructtype(u)
        throw(ArgumentError("The provided type must be a struct type"))
    end
    if length(fieldnames(u)) > 0
        throw(ArgumentError("The provided type must have no fields"))
    end
    if u.mutable
        throw(ArgumentError("The provided type must be an immutable type"))
    end
    if !(supertype(u) === ExtensibleUnion)
        throw(ArgumentError(
            "The immediate supertype of the provided type must be Any"))
    end
    if !haskey(_registry_extensibleunion_to_members, u)
        _registry_extensibleunion_to_genericfunctions[u] = Set{Any}()
        _registry_extensibleunion_to_members[u] = Set{Any}([u])
    end
    _update_all_methods_for_extensibleunion!(u)
    return u
end

function isextensibleunion(@nospecialize(u))
    global _registry_extensibleunion_to_members
    return haskey(_registry_extensibleunion_to_members, u)
end

function addtounion!(@nospecialize(u), varargs...)
    return addtounion!(u, varargs)
end

function addtounion!(@nospecialize(u), @nospecialize(varargs::Tuple))
    global _registry_extensibleunion_to_members
    if isextensibleunion(u)
        old_members_set = deepcopy(_registry_extensibleunion_to_members[u])
        for i = 1:length(varargs)
            push!(_registry_extensibleunion_to_members[u], varargs[i])
        end
        new_members_set = deepcopy(_registry_extensibleunion_to_members[u])
        old_members_union = _set_to_union(old_members_set)
        new_members_union = _set_to_union(new_members_set)
        _update_all_methods_for_extensibleunion!(u, old_members_union =>
                                                    new_members_union)
    else
        throw(ArgumentError(
            "First argument must be a registered extensible union."))
    end
end

function unioncurrentlycontains(@nospecialize(u), @nospecialize(t))
    global _registry_extensibleunion_to_members
    if isextensibleunion(u)
        return t in _registry_extensibleunion_to_members[u]
    else
        throw(ArgumentError(
            "First argument must be a registered extensible union."))
    end
end


text:
module HomotopyOpt

import HomotopyContinuation
import LinearAlgebra
import ImplicitPlots: implicit_plot
import Plots
import Statistics
import Implicit3DPlotting: plot_implicit_surface, plot_implicit_surface!, plot_implicit_curve, plot_implicit_curve!#, GLMakiePlottingLibrary
import Implicit3DPlotting.GLMakiePlottingLibrary as GLMakiePlottingLibrary
import ForwardDiff

export ConstraintVariety,
       findminima,
       watch,
       draw,
	   addSamples!,
	   setEquationsAtp!

# this code modifies `ConstrainedOptimizationByParameterHomotopy.jl`
# so that instead of an ObjectiveFunction, you specify a function called `evaluateobjectivefunctiongradient`
# this makes more sense. You have to define this function yourself,
# but now it does not depend on symbolic algebra from HomotopyContinuation.jl

# TODO add basepoint p+εv to the struct TrackerWithStartSolution, so we only have to solve the ED linear system once
# TODO reduce redundancy

#=
 Equips a HomotopyContinuation.Tracker with a start Solution that can be changed on the fly
=#
mutable struct TrackerWithStartSolution
	tracker
	startSolution
	#basepoint

	function TrackerWithStartSolution(T::HomotopyContinuation.Tracker, startSol::Vector)
		new(T,startSol)
	end
end

function setStartSolution(T::TrackerWithStartSolution, startSol::Vector)
	setfield!(T, :startSolution, startSol)
end

#=
 An object that describes a constraint variety by giving its generating equations, coordinate variables, its dimension and its jacobian.
 Additionally, it contains the system describing the Euclidian Distance Problem and samples from the variety.
=#
mutable struct ConstraintVariety
    variables
    equations
	fullequations
    jacobian
    ambientdimension
	dimensionofvariety
    samples
    implicitequations
	EDTracker

	# Given variables and HomotopyContinuation-based equations, sample points from the variety and return the corresponding struct
	function ConstraintVariety(varz, eqnz, N::Int, d::Int, numsamples::Int)
        dg = HomotopyContinuation.differentiate(eqnz, varz)
		impliciteq = [p->eqn(varz=>p) for eqn in eqnz]
        randL = nothing
		randresult=nothing
        Ωs = []
		if numsamples > 0
			randL = HomotopyContinuation.rand_subspace(N; codim=d)
			randResult = HomotopyContinuation.solve(eqnz; target_subspace = randL, variables=varz, show_progress = true)
		end
        for _ in 1:numsamples
            newΩs = HomotopyContinuation.solve(
                    eqnz,
                    HomotopyContinuation.solutions(randResult);
                    variables = varz,
                    start_subspace = randL,
                    target_subspace = HomotopyContinuation.rand_subspace(N; codim = d, real = true),
                    transform_result = (R,p) -> HomotopyContinuation.real_solutions(R),
                    flatten = true,
					show_progress = true
            )
            realsols = HomotopyContinuation.real_solutions(newΩs)
            push!(Ωs, realsols...)
        end
		fulleqnz = eqnz
		if length(eqnz) + d > N
			eqnz = randn(Float64, N-d, length(eqnz))*eqnz
		end

		HomotopyContinuation.@var u[1:N]
		HomotopyContinuation.@var λ[1:length(eqnz)]
		Lagrange = Base.sum((varz-u).^2) + sum(λ.*eqnz)
		∇Lagrange = HomotopyContinuation.differentiate(Lagrange, vcat(varz,λ))
		EDSystem = HomotopyContinuation.System(∇Lagrange, variables=vcat(varz,λ), parameters=u)
		p0 = HomotopyContinuation.randn(Float64, N)
		H = HomotopyContinuation.ParameterHomotopy(EDSystem, start_parameters = p0, target_parameters = p0)
		EDTracker = TrackerWithStartSolution(HomotopyContinuation.Tracker(H),[])
        new(varz,eqnz,fulleqnz,dg,N,d,Ωs,impliciteq,EDTracker)
    end

	# Given implicit equations, sample points from the corresponding variety and return the struct
    function ConstraintVariety(eqnz::Function, N::Int, d::Int, numsamples::Int)
        HomotopyContinuation.@var varz[1:N]
        algeqnz = [eqn(varz) for eqn in eqnz]
		ConstraintVariety(varz, algeqnz, N::Int, d::Int, numsamples::Int)
    end

	# Implicit Equations, no sampling
    function ConstraintVariety(eqnz,N::Int,d::Int)
		ConstraintVariety(eqnz::Function, N::Int, d::Int, 0)
	end

    # HomotopyContinuation-based expressions and variables, no sanples
    function ConstraintVariety(varz,eqnz,N::Int,d::Int)
		ConstraintVariety(varz, eqnz, N::Int, d::Int, 0)
    end
end

#=
Add Samples to an already existing ConstraintVariety
=#
function addSamples!(G::ConstraintVariety, newSamples)
	setfield!(G, :samples, vcat(newSamples, G.samples))
end

#=
Add Samples to an already existing ConstraintVariety
=#
function setEquationsAtp!(G::ConstraintVariety, p; tol=1e-6)
	jacobianRank = LinearAlgebra.rank(HomotopyContinuation.evaluate(G.jacobian, G.variables=>p); atol=tol)
	eqnz = G.fullequations
	if length(eqnz) + (G.ambientdimension-jacobianRank) > G.ambientdimension
		eqnz = randn(Float64, jacobianRank, length(eqnz))*eqnz
	end
	setfield!(G, :equations, eqnz)
	setfield!(G, :dimensionofvariety, (G.ambientdimension-jacobianRank))

	HomotopyContinuation.@var u[1:G.ambientdimension]
	HomotopyContinuation.@var λ[1:length(eqnz)]
	Lagrange = Base.sum((G.variables-u).^2) + sum(λ.*eqnz)
	∇Lagrange = HomotopyContinuation.differentiate(Lagrange, vcat(G.variables,λ))
	EDSystem = HomotopyContinuation.System(∇Lagrange, variables=vcat(G.variables,λ), parameters=u)
	H = HomotopyContinuation.ParameterHomotopy(EDSystem, start_parameters = p, target_parameters = p)
	EDTracker = TrackerWithStartSolution(HomotopyContinuation.Tracker(H),[])
	setfield!(G, :EDTracker, EDTracker)
end

#=
Compute the system that we need for the onestep and twostep method
=#
function computesystem(p, G::ConstraintVariety,
                evaluateobjectivefunctiongradient::Function)

    dgp = HomotopyContinuation.ModelKit.evaluate(G.jacobian, G.variables => p)
    Up,_ = LinearAlgebra.qr( LinearAlgebra.transpose(dgp) )
    Np = Up[:, 1:(G.ambientdimension - G.dimensionofvariety)] # gives ONB for N_p(G) normal space

    # we evaluate the gradient of the obj fcn at the point `p`
    ∇Qp = evaluateobjectivefunctiongradient(p)
    #display("evaluated the gradient of the objective fcn and got: $∇Qp")

    w = -∇Qp # direction of decreasing energy function
    v = w - Np * (Np' * w) # projected gradient -∇Q(p) onto the tangent space, subtract the normal components
	g = G.equations
	#display("variables $(G.variables) and equations $g")

    if G.dimensionofvariety > 1 # Need more linear equations when tangent space has dim > 1
        A,_ = LinearAlgebra.qr( hcat(v, Np))
        A = A[:, (G.ambientdimension - G.dimensionofvariety + 1):end] # basis of the orthogonal complement of v inside T_p(G)
        L = A' * G.variables - A' * p # affine linear equations through p, containing v, give curve in variety along v
		u = v/LinearAlgebra.norm(v)
        S = u' * G.variables - u' * (p + HomotopyContinuation.Variable(:ε)*u) # create and use the variable ε here.
        F = HomotopyContinuation.System( vcat(g,L,S); variables=G.variables,
                                        parameters=[HomotopyContinuation.Variable(:ε)])
        return F
    else
        u = LinearAlgebra.normalize(v)
        S = u' * G.variables - u' * (p + HomotopyContinuation.Variable(:ε)*u) # create and use the variable ε here.
        F = HomotopyContinuation.System( vcat(g,S); variables=G.variables,
                                        parameters=[HomotopyContinuation.Variable(:ε)])
        return F
    end
end

#=
If we are at a point of slow progression / singularity we blow the point up to a sphere and check the intersections (witness sets) with nearby components
for the sample with lowest energy
=#
function resolveSingularity(p, G::ConstraintVariety, Q::Function, evaluateobjectivefunctiongradient, whichstep; initialtime = Base.time(), maxseconds = 50)
	if length(p)>10
		q = gaussnewtonstep(G, p, 1e-2, -evaluateobjectivefunctiongradient(p); initialtime=initialtime, maxseconds=maxseconds)[1]
		if Q(q) < Q(p)
			return q, false
		else
			return p, true
		end
	end

	eqn = G.fullequations
	var = G.variables
	d = G.dimensionofvariety
	sphereAtPoint = sum((var.-p).^2)-0.001
	samples = []
	try
		F = HomotopyContinuation.System(vcat(eqn,[sphereAtPoint]))
		rel = HomotopyContinuation.solve(F; show_progress=false)
		samples = HomotopyContinuation.real_solutions(rel)
	catch e
		println("dimension -1: ", e)
	end
	for j in 1:d-1
		try
			for _ in 1:3
				a = rand(Float64, length(var), j)
				L = a'*var-a'*p
				F = HomotopyContinuation.System(vcat(eqn,[sphereAtPoint],L))
				append!(samples, HomotopyContinuation.real_solutions(HomotopyContinuation.solve(F; show_progress=false)))
			end
		catch e
			println("dimension -$(j+1): ", e)
		end
	end

	#TODO Alternative for too large to sample varieties.
	#Random directions? Sampling via Gaussnewton? Gaussnewtonstep altogether?
	minimumvalue = Q(p)
	q = Base.copy(p)
	for sol in samples
		if Q(sol)<minimumvalue
			minimumvalue = Q(sol)
			q = sol
		end
	end
	Basenormal, _, basegradient = getNandTandv(q, G, evaluateobjectivefunctiongradient)

	if q==p && !isempty(samples)
		#In this case, the singularity is optimal in a sense
		return(p,true)
	else
		return(q,false)
	end
end

#=
 We predict in the projected gradient direction and correct by using the Gauss-Newton method
=#
function gaussnewtonstep(G::ConstraintVariety, p, stepsize, v; tol=1e-8, initialtime, maxseconds)
	q = p+stepsize*v
	jac = Base.hcat([HomotopyContinuation.differentiate(eq, G.variables) for eq in G.fullequations]...)
	while(LinearAlgebra.norm([eq(G.variables=>q) for eq in G.fullequations]) > tol)
		J = jac(G.variables=>q)
		q .= q .- LinearAlgebra.pinv(J*LinearAlgebra.transpose(J))*J*[eq(G.variables=>q) for eq in G.fullequations]
		if time()-initialtime > maxseconds
			return p, false
		end
	end
	return q, true
end

#=
We predict in the projected gradient direction and correct by solving a Euclidian Distance Problem
=#
function EDStep(ConVar, p, stepsize, v)
	q = p+stepsize*v
	HomotopyContinuation.target_parameters!(ConVar.EDTracker.tracker, q)
	tracker = HomotopyContinuation.track(ConVar.EDTracker.tracker, ConVar.EDTracker.startSolution)
	res = HomotopyContinuation.solution(tracker)
	if all(point->Base.abs(point.im)<1e-4, res)
		# TODO set start parameter to be q? Maybe also at a different point?
		return [point.re for point in res[1:length(p)]], true
	else
		return p, false
	end
end

#=
 Move a line along the projected gradient direction for the length stepsize and calculate the resulting point of intersection with the variety
=#
function onestep(F, p, stepsize)
    # we want parameter homotopy from 0.0 to stepsize, so we take two steps
    # first from 0.0 to a complex number parameter, then from that parameter to stepsize.
    solveresult = HomotopyContinuation.solve(F, [p]; start_parameters=[0.0], target_parameters=[stepsize],
                                                     show_progress=false)
    sol = HomotopyContinuation.real_solutions(solveresult)
    success = false
    if length(sol) > 0
        q = sol[1] # only tracked one solution path, thus there should only be one solution
        success = true
    else
        q = p
    end
    return q, success
end

#=
 Similar to onestep. However, we take an intermediate, complex step to avoid singularities
=#
function twostep(F, p, stepsize)
    # we want parameter homotopy from 0.0 to stepsize, so we take two steps
    # first from 0.0 to a complex number parameter, then from that parameter to stepsize.
    midparam = stepsize/2 + stepsize/2*1.0im # complex number *midway* between 0 and stepsize, but off real line
    solveresult = HomotopyContinuation.solve(F, [p]; start_parameters=[0.0 + 0.0im], target_parameters=[midparam],
                                                    show_progress=false)
    midsols = HomotopyContinuation.solutions(solveresult)
    success = false
    if length(midsols) > 0
        midsolution = midsols[1] # only tracked one solution path, thus there should only be one solution
        solveresult = HomotopyContinuation.solve(F, [midsolution]; start_parameters=[midparam],
                                                    target_parameters=[stepsize + 0.0im],
                                                    show_progress=false)
        realsols = HomotopyContinuation.real_solutions(solveresult)
        if length(realsols) > 0
            q = realsols[1] # only tracked one solution path, thus there should only be one solution
            success = true
        else
            q = p
        end
    else
        q = p
    end
    return q, success
end

#=
 Checks, whether p is a local minimum of the objective function Q w.r.t. the tangent space Tp
=#
function isMinimum(G::ConstraintVariety, Q::Function, evaluateobjectivefunctiongradient, Tp, v, p::Vector; tol=1e-4, criticaltol=1e-3)
	H = ForwardDiff.hessian(Q, p)
	HConstraints = [HomotopyContinuation.evaluate(HomotopyContinuation.differentiate(HomotopyContinuation.differentiate(eq, G.variables), G.variables), G.variables=>p) for eq in G.fullequations]
	Qalg = Q(p)+(G.variables-p)'*ForwardDiff.gradient(Q,p)+0.5*(G.variables-p)'*H*(G.variables-p) # Taylor Approximation of x, since only the Hessian is of interest anyway
	HomotopyContinuation.@var λ[1:length(G.fullequations)]
	L = Qalg+λ'*G.fullequations
	∇L = HomotopyContinuation.differentiate(L, vcat(G.variables, λ))
	gL = Matrix{Float64}(HomotopyContinuation.differentiate(HomotopyContinuation.evaluate(∇L,G.variables=>p), λ))
	bL = -HomotopyContinuation.evaluate(HomotopyContinuation.evaluate(∇L,G.variables=>p), λ=>[0 for _ in 1:length(λ)])
	λ0 = gL\bL
	for i in 1:length(λ0)
		if λ0[i]==NaN || λ0[i]==Inf
			λ0[i]=1
		end
	end
	Htotal = H+λ0'*HConstraints
	projH = Tp'*Htotal*Tp
	projEigvals = real(LinearAlgebra.eigvals(projH)) #projH symmetric => all real eigenvalues
	indices = filter(i->LinearAlgebra.abs(projEigvals[i])<=tol, 1:length(projEigvals))
	projEigvecs = real(LinearAlgebra.eigvecs(projH))[:, indices]
	projEigvecs = Tp*projEigvecs
	precision=5e-2
	if all(q-> q>=tol, projEigvals) && LinearAlgebra.norm(v) <= criticaltol
		return true
	elseif any(q-> q<=-tol, projEigvals) || LinearAlgebra.norm(v) > criticaltol
		return false
		#TODO Third derivative at x_0 at proj hessian sing. vectors not 0?!
	else
		q, optimality = resolveSingularity(p, G, Q, evaluateobjectivefunctiongradient, "gaussnewtonstep")
		return optimality

		#=
		A = ∇Qp'*G.variables - ∇Qp'*p + ∇Qp'*∇Qp*t
		algObjective = Q(p)+ForwardDiff.gradient(Q,p)'*(G.variables-p)+0.5*(G.variables-p)'*ForwardDiff.hessian(Q,p)*(G.variables-p)
		L = algObjective + λ'*vcat(G.equations, A)
		∇L = HomotopyContinuation.differentiate(L, vcat(G.variables, λ))
		F = HomotopyContinuation.System(∇L, variables=vcat(G.variables, λ), parameters=vcat(t,r))
		F = HomotopyContinuation.System(vcat(A,G.equations), variables=vcat(G.variables, λ), parameters=vcat(t,r))
		=#
		#randsummand = randn(Float64,length(r))
		#=g = HomotopyContinuation.evaluate(HomotopyContinuation.differentiate(∇L, λ), vcat(G.variables, λ, t, r)=>vcat(p,[0 for _ in 1:length(λ)],0, randsummand))
		b = - HomotopyContinuation.evaluate(∇L, vcat(G.variables, λ, t, r)=>vcat(p,[0 for _ in 1:length(λ)],0, randsummand))
		λ0 = g\b
		HomotopyContinuation.@var a[1:G.ambientdimension-(G.dimensionofvariety+1),1:length(G.variables)] b[1:G.ambientdimension-(G.dimensionofvariety+1)]
		LinSp = a*G.variables+b
		realsols=[]
		if length(LinSp)==0
			Sys = HomotopyContinuation.System(vcat(A,G.equations), variables = G.variables, parameters=[t])
			res = HomotopyContinuation.solve(Sys,[p]; start_parameters=[0], target_parameters=[precision])
			reals = HomotopyContinuation.real_solutions(res)
			if length(reals)==0
				return true
			else
				push!(realsols, reals[Base.argmin([LinearAlgebra.norm(sol-p) for sol in reals])])
			end
		else
			Sys = HomotopyContinuation.System(vcat(A,G.equations,LinSp), variables = G.variables, parameters=vcat(Base.collect(Iterators.flatten(a)),b,t))
			for _ in 1:15
				res =  HomotopyContinuation.solve(Sys, target_parameters=vcat(randn(Float64, length(vcat(Base.collect(Iterators.flatten(a)),b))), precision))
				reals = HomotopyContinuation.real_solutions(res)
				if length(reals)!=0
					reals = realsols[Base.argmin([LinearAlgebra.norm(sol-p) for sol in reals])]
					reals = length(reals)==1 ? [reals] : reals
				end
				append!(realsols,reals)
			end
		end
		try
			if any(sol->Q(sol) < Q(p), realsols)
				# No strict minimum: We found a point that yields a smaller energy function
				return false
			else
				return true
			end
		catch e
			# If an error occurs in the tracking it means that we cannot continue aby path in the negative gradient direction from p. Thus, p is a minimum
			return true
		end
		=#
	end
end

#=
Determines, which optimization algorithm to use
=#
function stepchoice(F, ConstraintVariety, whichstep, stepsize, p, v; initialtime, maxseconds)
	if(whichstep=="twostep")
		return(twostep(F, p, stepsize))
	elseif whichstep=="onestep"
		return(onestep(F, p, stepsize))
	elseif whichstep=="gaussnewtonstep"
		return(gaussnewtonstep(ConstraintVariety, p, stepsize, v; initialtime, maxseconds))
	elseif whichstep=="EDStep"
		return(EDStep(ConstraintVariety, p, stepsize, v))
	else
		throw(error("A step method needs to be provided!"))
	end
end

# WARNING This one is worse than backtracking_linesearch
function alternative_backtracking_linesearch(Q::Function, F::HomotopyContinuation.ModelKit.System, G::ConstraintVariety, evaluateobjectivefunctiongradient::Function, p0::Vector, stepsize::Float64; maxstepsize=100.0, τ=0.4, r=1e-4, s=0.9, whichstep="twostep")
    α=Base.copy(stepsize)
    p=Base.copy(p0)

	Basenormal, _, basegradient = getNandTandv(p0, G, evaluateobjectivefunctiongradient)
	if whichstep=="EDStep"
		q0 = p+1e-3*Basenormal[:,1]
		HomotopyContinuation.start_parameters!(G.EDTracker.tracker, q0)
		A = HomotopyContinuation.evaluate(HomotopyContinuation.differentiate(G.EDTracker.tracker.homotopy.F.interpreted.system.expressions, G.EDTracker.tracker.homotopy.F.interpreted.system.variables[length(p)+1:end]), G.variables => p)
		λ0 = A\(-HomotopyContinuation.evaluate(HomotopyContinuation.evaluate(HomotopyContinuation.evaluate(G.EDTracker.tracker.homotopy.F.interpreted.system.expressions, G.EDTracker.tracker.homotopy.F.interpreted.system.variables[length(p)+1:end] => [0 for _ in length(p)+1:length(G.EDTracker.tracker.homotopy.F.interpreted.system.variables)]), G.variables => p),  G.EDTracker.tracker.homotopy.F.interpreted.system.parameters=>q0))
		setStartSolution(G.EDTracker, vcat(p,λ0))
	end
    while(true)
		q, success = stepchoice(F, G, whichstep, α, p0, basegradient)
        success ? p=q : nothing
        Nq, Tq, vq = getNandTandv(p, G, evaluateobjectivefunctiongradient)
        # Proceed until the Wolfe condition is satisfied or the stepsize becomes too small. First we quickly find a lower bound, then we gradually increase this lower-bound
		if (Q(p0)-Q(p) >= r*α*Base.abs(basegradient'*evaluateobjectivefunctiongradient(p0)) && vq'*basegradient >= 0 && success)

			for indicator in 1:12
                αsub = α*1.1
				q, success = stepchoice(F, G, whichstep, αsub, p0, basegradient)
                if(!success)
                    return(p, Nq, Tq, vq, false, α)
                end
                Nqsub, Tqsub, vqsub = getNandTandv(q, G, evaluateobjectivefunctiongradient)
                if Q(p0)-Q(q) < r*αsub*Base.abs(basegradient'*evaluateobjectivefunctiongradient(p0)) || vqsub'*basegradient < 0
                    return(p, Nq, Tq, vq, true, α)
                elseif Base.abs(basegradient'*evaluateobjectivefunctiongradient(q)) <= Base.abs(basegradient'*evaluateobjectivefunctiongradient(p0))*s
                    return(q, Nqsub, Tqsub, vqsub, true, αsub)
				elseif αsub>maxstepsize
					return(p, Nq, Tq, vq, false, α*τ)
                else
                    p=q; α=αsub; vq=vqsub; Tq=Tqsub; Nq=Nqsub;
                end
            end
			return(p, Nq, Tq, vq, true, α)

		elseif α<1e-6
	    	return(q, Nq, Tq, vq, false, stepsize)
        else
            α=τ*α
        end
    end
end

#=
Use line search with the strong Wolfe condition to find the optimal step length.
This particular method can be found in Nocedal & Wright: Numerical Optimization
=#
function backtracking_linesearch(Q::Function, F::HomotopyContinuation.ModelKit.System, G::ConstraintVariety, evaluateobjectivefunctiongradient::Function, p0::Vector, stepsize::Float64; maxstepsize=100.0, r=1e-3, s=0.9, whichstep="EDStep", initialtime, maxseconds)
	Basenormal, Basetangent, basegradient = getNandTandv(p0, G, evaluateobjectivefunctiongradient)
	α0 = 0
	α = [0, stepsize]
	p = Base.copy(p0)
	if whichstep=="EDStep"
		q0 = p+1e-4*Basenormal[:,1]
		HomotopyContinuation.start_parameters!(G.EDTracker.tracker, q0)
		A = HomotopyContinuation.evaluate(HomotopyContinuation.differentiate(G.EDTracker.tracker.homotopy.F.interpreted.system.expressions, G.EDTracker.tracker.homotopy.F.interpreted.system.variables[length(p)+1:end]), G.variables => p)
		λ0 = A\(-HomotopyContinuation.evaluate(HomotopyContinuation.evaluate(HomotopyContinuation.evaluate(G.EDTracker.tracker.homotopy.F.interpreted.system.expressions, G.EDTracker.tracker.homotopy.F.interpreted.system.variables[length(p)+1:end] => [0 for _ in length(p)+1:length(G.EDTracker.tracker.homotopy.F.interpreted.system.variables)]), G.variables => p),  G.EDTracker.tracker.homotopy.F.interpreted.system.parameters => q0))
		setStartSolution(G.EDTracker, vcat(p, λ0))
	end
    while true
		q, success = stepchoice(F, G, whichstep, α[end], p0, basegradient; initialtime, maxseconds)
		if time()-initialtime > maxseconds
			Nq, Tq, vq = getNandTandv(q, G, evaluateobjectivefunctiongradient)
			return q, Nq, Tq, vq, success, α[end]
		end
        Nq, Tq, vq = getNandTandv(q, G, evaluateobjectivefunctiongradient)
		if ( ( Q(q) > Q(p0) + r*α[end]*basegradient'*basegradient || (Q(q) > Q(p0) && q!=p0) ) && success)
			helper = zoom(α[end-1], α[end], Q, evaluateobjectivefunctiongradient, F, G, whichstep, p0, basegradient, r, s; initialtime, maxseconds)
			Nq, Tq, vq = getNandTandv(helper[1], G, evaluateobjectivefunctiongradient)
			return helper[1], Nq, Tq, vq, helper[2], helper[end]
		end
		if ( Base.abs(basegradient'*vq) <= s*Base.abs(basegradient'*basegradient) ) && success
			return q, Nq, Tq, vq, success, α[end]
		end
		if basegradient'*vq <= 0 && success
			helper = zoom(α[end], α[end-1], Q, evaluateobjectivefunctiongradient, F, G, whichstep, p0, basegradient, r, s; initialtime, maxseconds)
			Nq, Tq, vq = getNandTandv(helper[1], G, evaluateobjectivefunctiongradient)
			return helper[1], Nq, Tq, vq, helper[2], helper[end]
		end
		if (success)
			push!(α, (α[end]))
			p = q
		else
			Np, Tp, vp = getNandTandv(p, G, evaluateobjectivefunctiongradient)
			return p, Np, Tp, vp, success, α[end]
		end
		deleteat!(α, 1)
		if α[end] > maxstepsize
			return q, Nq, Tq, vq, success, maxstepsize
		end
    end
end

#=
Zoom in on the step lengths between αlo and αhi to find the optimal step size here. This is part of the backtracking line search
=#
function zoom(αlo, αhi, Q, evaluateobjectivefunctiongradient, F, G, whichstep, p0, basegradient, r, s; initialtime, maxseconds)
	qlo, suclo = stepchoice(F, G, whichstep, αlo, p0, basegradient; initialtime, maxseconds)
	# To not get stuck in the iteration, we use a for loop instead of a while loop
	# TODO Add a more meaningful stopping criterion
	for i in 1:15
		global α = 0.5*(αlo+αhi)
		global q, success = stepchoice(F, G, whichstep, α, p0, basegradient; initialtime, maxseconds)
		Nq, Tq, vq = getNandTandv(q, G, evaluateobjectivefunctiongradient)
		if !success || time()-initialtime > maxseconds
			return q, success, α
		end

		if  Q(q) > Q(p0) + r*α*basegradient'*basegradient || Q(q) >= Q(qlo)
			αhi = α
		else
			if Base.abs(basegradient'*vq) <= Base.abs(basegradient'*basegradient)*s
				return q, success, α
			end
			if basegradient'*vq*(αhi-αlo) >= 0
				αhi = αlo
			end
			αlo = α
			qlo, suclo = q, success
		end
	end
	return q, success, α
end

#=
 Get the tangent and normal space of a ConstraintVariety at a point q
=#
function getNandTandv(q, G::ConstraintVariety,
                    evaluateobjectivefunctiongradient::Function)
    dgq = HomotopyContinuation.ModelKit.evaluate(G.jacobian, G.variables => q)
    Qq,_ = LinearAlgebra.qr(transpose(dgq))
	#TODO This approach does not work for hypersurfaces I guess?
	#index = count(p->p>1e-8, S)
    Nq = Qq[:, 1:(G.ambientdimension - G.dimensionofvariety)] # O.N.B. for the normal space at q
    Tq = Qq[:, (G.ambientdimension - G.dimensionofvariety + 1):end] # O.N.B. for tangent space at q
    # we evaluate the gradient of the obj fcn at the point `q`
    ∇Qq = evaluateobjectivefunctiongradient(q)

    w = -∇Qq # direction of decreasing energy function
    vq = w - Nq * (Nq' * w) # projected gradient -∇Q(p) onto the tangent space, subtract the normal components
    return Nq, Tq, vq
end

#=
 Parallel transport the vector vj from the tangent space Tj to the tangent space Ti
=#
function paralleltransport(vj, Tj, Ti)
    # transport vj ∈ Tj to become a vector ϕvj ∈ Ti
    # cols(Tj) give ONB for home tangent space, cols(Ti) give ONB for target tangent space
    U,_,Vt = LinearAlgebra.svd( Ti' * Tj )
    Oij = U * Vt # closest orthogonal matrix to the matrix (Ti' * Tj) comes from svd, remove \Sigma
    ϕvj = Ti * Oij * (Tj' * vj)
    return ϕvj
end

#=
An object that contains the iteration's information like norms of the projected gradient, step sizes and search directions
=#
struct LocalStepsResult
    initialpoint
    initialstepsize
    allcomputedpoints
    allcomputedprojectedgradientvectors
    allcomputedprojectedgradientvectornorms
    newsuggestedstartpoint
    newsuggestedstepsize
    converged
    timesturned
    valleysfound

    function LocalStepsResult(p,ε0,qs,vs,ns,newp,newε0,converged,timesturned,valleysfound)
        new(p,ε0,qs,vs,ns,newp,newε0,converged,timesturned,valleysfound)
    end
end

#= Take `maxsteps` steps to try and converge to an optimum. In each step, we use backtracking linesearch
to determine the optimal step size to go along the search direction
WARNING This is redundant and can be merged with findminima
=#
function takelocalsteps(p, ε0, tolerance, G::ConstraintVariety,
                objectiveFunction::Function,
                evaluateobjectivefunctiongradient::Function;
                maxsteps, decreasefactor=2, initialtime, maxseconds, whichstep="twostep", maxstepsize=200.0)
    timesturned, valleysfound, F = 0, 0, HomotopyContinuation.System([G.variables[1]])
    _, Tp, vp = getNandTandv(p, G, evaluateobjectivefunctiongradient)
    Ts = [Tp] # normal spaces and tangent spaces, columns of Np and Tp are orthonormal bases
    qs, vs, ns = [p], [vp], [LinearAlgebra.norm(vp)] # qs=new points on G, vs=projected gradients, ns=norms of projected gradients
    stepsize = Base.copy(ε0)
    for count in 1:maxsteps
        if Base.time() - initialtime > maxseconds
			break;
        end
		if whichstep=="onestep" || whichstep=="twostep"
        	F = computesystem(qs[end], G, evaluateobjectivefunctiongradient)
		end
        q, Nq, Tq, vq, success, stepsize = backtracking_linesearch(objectiveFunction, F, G, evaluateobjectivefunctiongradient, qs[end], Float64(stepsize); whichstep, maxstepsize, initialtime, maxseconds)
		push!(qs, q)
        push!(Ts, Tq)
		length(Ts)>3 ? deleteat!(Ts, 1) : nothing
        push!(ns, LinearAlgebra.norm(vq))
		push!(vs, vq)
		length(vs)>3 ? deleteat!(vs, 1) : nothing
        if ns[end] < tolerance
            return LocalStepsResult(p,ε0,qs,vs,ns,q,stepsize,true,timesturned,valleysfound)
		# TODO I believe parallel transport is redundant. I think it is already covered by the backtracking linesearch. Not sure though.
        elseif ((ns[end] - ns[end-1]) > 0.0)
            if length(ns) > 2 && ((ns[end-1] - ns[end-2]) < 0.0)
                # projected norms were decreasing, but started increasing!
                # check parallel transport dot product to see if we should slow down
                valleysfound += 1
                ϕvj = paralleltransport(vs[end], Ts[end], Ts[end-2])
                if ((vs[end-2]' * ϕvj) < 0.0)
                    # we think there is a critical point we skipped past! slow down!
                    return LocalStepsResult(p,ε0,qs,vs,ns,qs[end-2],stepsize/decreasefactor,false,timesturned+1,valleysfound)
                end
            end
        end
        # The next (initial) stepsize is determined by the previous step and how much the energy function changed - in accordance with RieOpt.
		stepsize = Base.minimum([ Base.maximum([ success ? 1.5*stepsize*vs[end-1]'*evaluateobjectivefunctiongradient(qs[end-1])/(vs[end]'*evaluateobjectivefunctiongradient(qs[end]))  : 0.1*stepsize, 1e-4]), maxstepsize])
    end
    return LocalStepsResult(p,ε0,qs,vs,ns,qs[end],stepsize,false,timesturned,valleysfound)
end

#=
 Output object of the method `findminima`
=#
struct OptimizationResult
    computedpoints
    initialpoint
    initialstepsize
    tolerance
    converged
    lastlocalstepsresult
    constraintvariety
    objectivefunction
	lastpointisminimum

    function OptimizationResult(ps,p0,ε0,tolerance,converged,lastLSResult,G,Q,lastpointisminimum)
        new(ps,p0,ε0,tolerance,converged,lastLSResult,G,Q,lastpointisminimum)
    end
end

#=
 The main function of this package. Given an initial point, a tolerance, an objective function and a constraint variety,
 we try to find the objective function's closest local minimum to the initial guess.
=#
function findminima(p0, tolerance,
                G::ConstraintVariety,
                objectiveFunction::Function;
                maxseconds=100, maxlocalsteps=3, initialstepsize=1.0, whichstep="twostep", initialtime = Base.time())
	#TODO Rework minimality: We are not necessarily at a minimality, if resolveSingularity does not find any better point. => first setequations, then ismin
	setEquationsAtp!(G,p0)
    p = copy(p0) # initialize before updating `p` below
    ps = [p0] # record the *main steps* from p0, newp, newp, ... until converged
	jacobianG = HomotopyContinuation.evaluate(HomotopyContinuation.differentiate(G.equations, G.variables), G.variables=>p0)
	jacRank = LinearAlgebra.rank(jacobianG; atol=tolerance)
    evaluateobjectivefunctiongradient = x -> ForwardDiff.gradient(objectiveFunction, x)
    _, Tq, v = getNandTandv(p0, G, evaluateobjectivefunctiongradient) # Get the projected gradient at the first point
	# initialize stepsize. Different to RieOpt! Logic: large projected gradient=>far away, large stepsize is admissible.
	ε0 = 2*initialstepsize
    lastLSR = LocalStepsResult(p,ε0,[],[],[],p,ε0,false,0,0)
    while (Base.time() - initialtime) <= maxseconds
        # update LSR, only store the *last local run*
        lastLSR = takelocalsteps(p, ε0, tolerance, G, objectiveFunction, evaluateobjectivefunctiongradient; maxsteps=maxlocalsteps, initialtime, maxseconds, whichstep)
		push!(ps, lastLSR.allcomputedpoints[end])
		println(lastLSR.allcomputedprojectedgradientvectornorms[end])
        if lastLSR.converged
			jacobian = HomotopyContinuation.evaluate(HomotopyContinuation.differentiate(G.equations, G.variables), G.variables=>lastLSR.newsuggestedstartpoint)
			jR = LinearAlgebra.rank(jacobian; atol=tolerance)
			# if we are in a singularity do a few steps again - if we revert back to the singularity, it is optiomal
			if jR != jacRank || LinearAlgebra.norm(ps[end-1]-ps[end]) < tolerance^2
				p, optimality = resolveSingularity(lastLSR.allcomputedpoints[end], G, objectiveFunction, evaluateobjectivefunctiongradient, whichstep; initialtime=initialtime, maxseconds=maxseconds)
				setEquationsAtp!(G, p; tol=tolerance)
				if !optimality
					maxseconds = maxseconds+5
					optRes = findminima(p, tolerance, G, objectiveFunction; maxseconds = maxseconds, maxlocalsteps=maxlocalsteps, initialstepsize=initialstepsize, whichstep=whichstep, initialtime=initialtime)
					return OptimizationResult(vcat(ps, optRes.computedpoints),p0,lastLSR.newsuggestedstepsize,tolerance,optRes.lastlocalstepsresult.converged,optRes.lastlocalstepsresult,G,evaluateobjectivefunctiongradient,optRes.lastpointisminimum)
				end
				return OptimizationResult(ps,p0,initialstepsize,tolerance,true,lastLSR,G,evaluateobjectivefunctiongradient,optimality)
			else
				_, Tq, v = getNandTandv(ps[end], G, evaluateobjectivefunctiongradient)
				setEquationsAtp!(G, p; tol=tolerance)
				optimality = isMinimum(G, objectiveFunction, evaluateobjectivefunctiongradient, Tq, v, ps[end]; criticaltol=tolerance)
				if !optimality
					optRes = findminima(ps[end], tolerance, G, objectiveFunction; maxseconds = maxseconds, maxlocalsteps=maxlocalsteps, initialstepsize=initialstepsize, whichstep=whichstep, initialtime=initialtime)
					return OptimizationResult(vcat(ps, optRes.computedpoints),p0,lastLSR.newsuggestedstepsize,tolerance,optRes.lastlocalstepsresult.converged,optRes.lastlocalstepsresult,G,evaluateobjectivefunctiongradient,optRes.lastpointisminimum)
				end
				return OptimizationResult(ps,p0,initialstepsize,tolerance,true,lastLSR,G,evaluateobjectivefunctiongradient,optimality)
			end
        else
			jacobian = HomotopyContinuation.evaluate(HomotopyContinuation.differentiate(G.equations, G.variables), G.variables=>lastLSR.newsuggestedstartpoint)
			jR = LinearAlgebra.rank(jacobian; atol=tolerance)
			# If we are in a point of slow progress or jacobian rank change, we search the neighborhood
			if jR != jacRank || LinearAlgebra.norm(ps[end-1]-ps[end]) < tolerance^3
				p, optimality = resolveSingularity(lastLSR.allcomputedpoints[end], G, objectiveFunction, evaluateobjectivefunctiongradient, whichstep; initialtime=initialtime, maxseconds=maxseconds)
				setEquationsAtp!(G, p; tol=tolerance)
				if optimality
					return OptimizationResult(ps,p0,initialstepsize,tolerance,true,lastLSR,G,evaluateobjectivefunctiongradient,optimality)
				end
				jacobian = HomotopyContinuation.evaluate(HomotopyContinuation.differentiate(G.equations, G.variables), G.variables=>p)
				jR = LinearAlgebra.rank(jacobian; atol=tolerance)
			else
				p = lastLSR.newsuggestedstartpoint
			end
			jacRank = jR
            ε0 = lastLSR.newsuggestedstepsize # update and try again!
        end
    end

	display("We ran out of time... Try setting `maxseconds` to a larger value than $(maxseconds)")
	p, optimality = resolveSingularity(ps[end], G, objectiveFunction, evaluateobjectivefunctiongradient, whichstep; initialtime=initialtime, maxseconds=maxseconds)
	setEquationsAtp!(G, p; tol=tolerance)
	_, Tq, v = getNandTandv(ps[end], G, evaluateobjectivefunctiongradient)
    return OptimizationResult(ps,p0,ε0,tolerance,lastLSR.converged,lastLSR,G,evaluateobjectivefunctiongradient,optimality)
end

# Below are functions `watch` and `draw`
# to visualize low-dimensional examples
function watch(result::OptimizationResult; totalseconds=5.0, kwargs...)
    ps = result.computedpoints
    samples = result.constraintvariety.samples
	# TODO centroid approach rather than mediannorm and then difference from centroid.
    mediannorm = Statistics.median([LinearAlgebra.norm(p) for p in samples])
    samples = filter(x -> LinearAlgebra.norm(x) < 2*mediannorm+0.5, samples)
    initplt = Plots.plot() # initialize
    M = length(ps)
    framespersecond = M / totalseconds
    if framespersecond > 45
        framespersecond = 45
    end
    startingtime = Base.time()
    dim = length(ps[1])
    anim = Plots.Animation()
    if dim == 2
        fullx = [minimum([q[1] for q in vcat(samples, ps)]) - 0.05, maximum([q[1] for q in vcat(samples, ps)]) + 0.05]
        fully = [minimum([q[2] for q in vcat(samples, ps)]) - 0.05, maximum([q[2] for q in vcat(samples, ps)]) + 0.05]
        g1 = result.constraintvariety.equations[1] # should only be a curve in ambient R^2
        initplt = implicit_plot(g1, xlims=fullx, ylims=fully, legend=false)
		initplt = Plots.scatter!(initplt, [ps[end][1]], [ps[end][2]], legend=false, markersize=5, color=:red, xlims=fullx, ylims=fully)
        Plots.frame(anim)
        for p in ps
            # BELOW: only plot next point, delete older points during animation
            # plt = scatter!(initplt, [p[1]], [p[2]], legend=false, color=:black, xlims=fullx, ylims=fully)
            # BELOW: keep old points during animation.
			initplt = Plots.scatter!(initplt, [p[1]], [p[2]], legend=false, markersize=3.5, color=:black, xlims=fullx, ylims=fully)
            Plots.frame(anim)
        end
        return Plots.gif(anim, "watch$startingtime.gif", fps=framespersecond)
    elseif dim == 3
        fullx = [minimum([q[1] for q in vcat(samples, ps)]) - 0.05, maximum([q[1] for q in vcat(samples, ps)]) + 0.05]
        fully = [minimum([q[2] for q in vcat(samples, ps)]) - 0.05, maximum([q[2] for q in vcat(samples, ps)]) + 0.05]
        fullz = [minimum([q[3] for q in vcat(samples, ps)]) - 0.05, maximum([q[3] for q in vcat(samples, ps)]) + 0.05]
        g1 = result.constraintvariety.implicitequations[1]
        if(length(result.constraintvariety.implicitequations)>1)
            # should be space curve
            g2 = result.constraintvariety.implicitequations[2]
            initplt = plot_implicit_curve(g1,g2; xlims = (fullx[1], fullx[2]), ylims = (fully[1], fully[2]), zlims = (fullz[1], fullz[2]), kwargs...)
        else
            #should be surface
            initplt = plot_implicit_surface(g1;  xlims = (fullx[1], fullx[2]), ylims = (fully[1], fully[2]), zlims = (fullz[1], fullz[2]), kwargs...)
        end
        pointsys=[GLMakiePlottingLibrary.Point3f0(p) for p in ps]
		GLMakiePlottingLibrary.scatter!(initplt, pointsys[end];
					color=:red, markersize=40.0)
        GLMakiePlottingLibrary.record(initplt, "watch$startingtime.gif", 1:length(pointsys); framerate = Int64(round(framespersecond))) do i
			GLMakiePlottingLibrary.scatter!(initplt, pointsys[i];
	                    color=:black, markersize=30.0)
        end
        return(initplt)
    end
end

function draw(result::OptimizationResult; kwargs...)
    dim = length(result.computedpoints[1]) # dimension of the ambient space
	ps = result.computedpoints
	samples = result.constraintvariety.samples
	mediannorm = Statistics.median([LinearAlgebra.norm(p) for p in samples])
	# TODO centroid approach rather than mediannorm and then difference from centroid.
	samples = filter(x -> LinearAlgebra.norm(x) < 2*mediannorm+0.5, samples)
    if dim == 2
		fullx = [minimum([q[1] for q in vcat(samples, ps)]) - 0.05, maximum([q[1] for q in vcat(samples, ps)]) + 0.05]
        fully = [minimum([q[2] for q in vcat(samples, ps)]) - 0.05, maximum([q[2] for q in vcat(samples, ps)]) + 0.05]
        g1 = result.constraintvariety.equations[1] # should only be a curve in ambient R^2
        plt1 = implicit_plot(g1, xlims=fullx, ylims=fully, legend=false)
        #f(x,y) = (x^4 + y^4 - 1) * (x^2 + y^2 - 2) + x^5 * y # replace this with `curve`
        #plt1 = implicit_plot(curve; xlims=(-2,2), ylims=(-2,2), legend=false)
        #plt2 = implicit_plot(curve; xlims=(-2,2), ylims=(-2,2), legend=false)
        localqs = result.lastlocalstepsresult.allcomputedpoints
        zoomx = [minimum([q[1] for q in localqs]) - 0.05, maximum([q[1] for q in localqs]) + 0.05]
        zoomy = [minimum([q[2] for q in localqs]) - 0.05, maximum([q[2] for q in localqs]) + 0.05]
		plt2 = implicit_plot(g1, xlims=zoomx, ylims=zoomy, legend=false)
        for q in ps
            plt1 = Plots.scatter!(plt1, [q[1]], [q[2]], legend=false, color=:black, xlims=fullx, ylims=fully)
        end
        for q in localqs
            plt2 = Plots.scatter!(plt2, [q[1]], [q[2]], legend=false, color=:blue, xlims=zoomx, ylims=zoomy)
        end
        vnorms = result.lastlocalstepsresult.allcomputedprojectedgradientvectornorms
        pltvnorms = Plots.plot(vnorms, legend=false, title="norm(v) for last local steps")
        plt = Plots.plot(plt1,plt2,pltvnorms, layout=(1,3), size=(900,300) )
        return plt
    elseif dim == 3
		fullx = [minimum([q[1] for q in vcat(samples, ps)]) - 0.05, maximum([q[1] for q in vcat(samples, ps)]) + 0.05]
        fully = [minimum([q[2] for q in vcat(samples, ps)]) - 0.05, maximum([q[2] for q in vcat(samples, ps)]) + 0.05]
        fullz = [minimum([q[3] for q in vcat(samples, ps)]) - 0.05, maximum([q[3] for q in vcat(samples, ps)]) + 0.05]
		localqs = result.lastlocalstepsresult.allcomputedpoints
		zoomx = [minimum([q[1] for q in ps]) - 0.05, maximum([q[1] for q in ps]) + 0.05]
        zoomy = [minimum([q[2] for q in ps]) - 0.05, maximum([q[2] for q in ps]) + 0.05]
        zoomz = [minimum([q[3] for q in ps]) - 0.05, maximum([q[3] for q in ps]) + 0.05]

        fig = GLMakiePlottingLibrary.Figure(resolution = (1450, 550))
        ax1 = fig[1, 1] = GLMakiePlottingLibrary.AbstractPlotting.MakieLayout.LScene(fig, width=500, height=500, camera = GLMakiePlottingLibrary.cam3d!, raw = false, limits=GLMakiePlottingLibrary.FRect((fullx[1], fully[1], fullz[1]), (fullx[2]-fullx[1], fully[2]-fully[1], fullz[2]-fullz[1])))
        ax2 = fig[1, 2] = GLMakiePlottingLibrary.AbstractPlotting.MakieLayout.LScene(fig, width=500, height=500, camera = GLMakiePlottingLibrary.cam3d!, raw = false, limits=GLMakiePlottingLibrary.FRect((zoomx[1], zoomy[1], zoomz[1]), (zoomx[2]-zoomx[1], zoomy[2]-zoomy[1], zoomz[2]-zoomz[1])))
        ax3 = fig[1, 3] = GLMakiePlottingLibrary.AbstractPlotting.MakieLayout.Axis(fig, width=300, height=450, title="norm(v) for last local steps")
        g1 = result.constraintvariety.implicitequations[1]
        if(length(result.constraintvariety.implicitequations)>1)
            # should be space curve
            g2 = result.constraintvariety.implicitequations[2]
            plot_implicit_curve!(ax1,g1,g2; xlims=(fullx[1],fullx[2]), ylims=(fully[1],fully[2]), zlims=(fullz[1],fullz[2]), kwargs...)
            plot_implicit_curve!(ax2,g1,g2; xlims=(zoomx[1],zoomx[2]), ylims=(zoomy[1],zoomy[2]), zlims=(zoomz[1],zoomz[2]), kwargs...)
        else
            plot_implicit_surface!(ax1,g1; xlims=(fullx[1],fullx[2]), ylims=(fully[1],fully[2]), zlims=(fullz[1],fullz[2]), kwargs...)
            plot_implicit_surface!(ax2,g1; xlims=(zoomx[1],zoomx[2]), ylims=(zoomy[1],zoomy[2]), zlims=(zoomz[1],zoomz[2]), kwargs...)
        end

        for q in ps
            GLMakiePlottingLibrary.scatter!(ax1, GLMakiePlottingLibrary.Point3f0(q);
                legend=false, color=:black, markersize=15)
            GLMakiePlottingLibrary.scatter!(ax2, GLMakiePlottingLibrary.Point3f0(q);
                legend=false, color=:black)
        end

        vnorms = result.lastlocalstepsresult.allcomputedprojectedgradientvectornorms
        GLMakiePlottingLibrary.plot!(ax3,vnorms; legend=false)

        return fig
    end
end

end


text:
function hello()

end


text:
import Base.circshift!

struct QNCache{T1, T2}
    y::T1 # change in successive gradients
    d::T2 # search direction
    s::T2 # change in x
end
function QNCache(x, g)
    QNCache(copy(g), copy(x), copy(x))
end
function preallocate_qn_caches_inplace(x0)
    # Maintain gradient and state pairs in QNCache
    cache = QNCache(x0, x0)
    return cache
end

function update_obj!(objective, d, s, y, ∇fx, z, ∇fz, B, scheme, scale=nothing)
    fz, ∇fz = objective(∇fz, z)
    # add Project gradient

    # Update y
    @. y = ∇fz - ∇fx

    # Update B
    if scale == nothing
        if isa(scheme.approx, DirectApprox)
            Badj = dot(y, d)/dot(y, y).*B
        else
            Badj = dot(y, y)/dot(y, d).*B
        end
    else
        Badj = B
    end

    # Quasi-Newton update
    B = update!(Badj, s, y, scheme)

    return fz, ∇fz, B
end

function update_obj!(objective, d, s, y, ∇fx, z, ∇fz, B, scheme::Newton, scale=nothing)
    fz, ∇fz, B = objective(B, ∇fz, z)

    return fz, ∇fz, B
end

function update_obj(objective, d, s, ∇fx, z, ∇fz, B, scheme, scale=nothing)
    fz, ∇fz = objective(∇fz, z)
    # add Project gradient

    # Update y
    y = ∇fz - ∇fx

    # Update B
    if scale == nothing
        if isa(scheme.approx, DirectApprox)
            Badj = dot(y, d)/dot(y, B* y)*B
        else
            Badj = dot(y, B* y)/dot(y, d)*B
        end
    else
         Badj = B
    end

    # Quasi-Newton update
    B = update(Badj, s, y, scheme)

    return fz, ∇fz, B
end

function update_obj(objective, d, s, ∇fx, z, ∇fz, B, scheme::Newton, is_first=nothing)
    fz, ∇fz, B = objective(B, ∇fx, z)

    return fz, ∇fz, B
end


text:
using Mass2
using ReferenceMicrocalFiles
using ZMQ
using Base.Test

"CountRateHistogram contains a RunningVector of Histograms, where hists[i] contains
all counts observed with `first(posix_timestamp_usec_bin_edges)<timestamp<=posix_timestamp_usec_bin_edges[i]`.
Count rates over ROIS can easily be calculated by differences of histograms."
type CountRateHistogram <: AbstractHistogram
  posix_timestamp_usec_bin_edges::Range{Int}
  hists::RunningVector{Histogram}
  x::FloatRange{Float64} # bins of the histograms in hists
  countedname::ASCIIString
end

"`posix_timestamp_usec_bin_edges` should always start at 0 (this should be an inner constructor?)"
function CountRateHistogram(posix_timestamp_usec_bin_step,x, countedname)
  posix_timestamp_usec_bin_edges = 0:posix_timestamp_usec_bin_step:typemax(Int64)
  h = CountRateHistogram(posix_timestamp_usec_bin_edges, RunningVector(Vector{Histogram}(), 0), x, countedname)
end
"update_histogram!(h::CountRateHistogram, selection, events, timestamps) adds the values in `events[selection]`
to the `h`, adding to the appropriate entry in `hists` based on `timestamps`. returns nothing.
timestamps must be sorted!!"
function update_histogram!(h::CountRateHistogram, selection, events, timestamps, last_timestamp_checked)
  @assert(length(selection)==length(events)==length(timestamps))
  if length(selection) == 0 # if the input arrays are empty, return early
    return donethru(h)
  end
  #last_timestamp_checked = max(last_timestamp_checked, last(timestamps)
  last_timestamp_checked = last(timestamps)
  j_hi = 0
  while true
    if length(available(h.hists)) >=1 # figure out what timestamp range we're interested in
      i_lo, i_hi = length(h.hists), length(h.hists)+1 # indexes into h.posix_timestamp_usec_bin_edges
    else
      i_range = searchsorted(h.posix_timestamp_usec_bin_edges, first(timestamps))
      i_lo, i_hi = last(i_range), first(i_range)
      i_lo = i_lo == i_hi ? i_lo-1 : i_lo # if they are equal, it means the fist timestamp is exactly equal to h.posix_timestamp_usec_bin_edges[i_hi]
      freeuntil!(h.hists, i_lo)
    end
    i_lo, i_hi
    t_lo, t_hi = h.posix_timestamp_usec_bin_edges[i_lo], h.posix_timestamp_usec_bin_edges[i_hi]
    if last_timestamp_checked<t_hi
      return donethru(h)
    end # only add histogram for which we have seen all events
    # t_lo and t_hi are the edges of the time bin to be used for hists[i_hi]
    j_lo = j_hi+1
    while j_hi+1 <= length(timestamps) && timestamps[j_hi+1]<=t_hi
      j_hi+=1
    end
    if j_lo<=j_hi
      @assert(timestamps[j_lo]>t_lo)
    else
      return donethru(h)
    end
    newhist = Histogram(h.x)
    update_histogram!(newhist, selection[j_lo:j_hi], events[j_lo:j_hi])
    if i_lo in available(h.hists)
      newhist = newhist+h.hists[i_lo]
    end
    push!(h.hists, newhist)
  end
  return donethru(h)
end

function donethru(h::CountRateHistogram)
  if length(available(h.hists)) == 0
    0
  else
    h.hists[end].seen
  end
end

tstep = 1000000 # usec
Nstep = 10
Npoints = 1000000
tfirst = 8000000000
crh = CountRateHistogram(tstep, 0:1:1000, "test")

events = ones(Npoints);
selection = trues(Npoints);
tlast = tfirst+Nstep*tstep
timestamps = round(Int,linspace(tfirst+1, tlast, Npoints))

update_histogram!(crh, selection, events, timestamps, last(timestamps))
donethru(crh) == Npoints
length(available(crh.hists)) == Nstep
for (i,j) in enumerate(available(crh.hists))
  @assert counts(crh.hists[j])[1]==div(Npoints,Nstep)*i
end

using ZMQ
port = 5555
context = ZMQ.Context()
socket = ZMQ.Socket(context,ZMQ.PUB)
ZMQ.bind(socket, "tcp://*:$port")
function sendtest()
  ZMQ.send(socket,"test")
end
"(channel_number::Int, crh, i) send CountRateHistogram `crh` over a zmq socket with format:
channel number; bin edges; counts; total counts seen; timestamps of start and end of new observations"
function send_crh(channel_number::Integer, crh, i)
  h = crh.hists[i]
  # last argument is SNDMORE, if true, it builds up a multipart message
  ZMQ.send(socket, crh.countedname, true)
  ZMQ.send(socket, "$channel_number", true)
  ZMQ.send(socket, repr(bin_edges(h)), true)
  ZMQ.send(socket, repr(counts(h)), true)
  ZMQ.send(socket, repr(seen(h)), true)
  ZMQ.send(socket, repr(crh.posix_timestamp_usec_bin_edges[[i-1,i]]))
end
function send_and_free_crh(channel_number::Integer, crh, i)
  send_crh(channel_number, crh, i)
  freeuntil!(crh.hists,i)
end
sub_socket = ZMQ.Socket(context, ZMQ.SUB)
ZMQ.subscribe(sub_socket)
ZMQ.connect(sub_socket, "tcp://localhost:5555")

function recv_multipart(socket)
  out = []
  push!(out, recv(socket))
  while ZMQ.ismore(socket)
    push!(out, recv(socket))
  end
  out
end

# @schedule for i = 1:100000 println(i);sleep(0.001); println(map(bytestring,recv_multipart(sub_socket))) end

# sendtest()
# send_and_free_crh(1, crh, 8011)

function update_histogram_send!(h::CountRateHistogram, selection, events, timestamps, last_timestamp_checked, channel_number)
  update_histogram!(h::CountRateHistogram, selection, events, timestamps, last_timestamp_checked)
  if length(available(h.hists))>=2
    i = available(h.hists)[end-1]
    send_and_free_crh(channel_number, h, i)
  end
  return nothing
end

"setup_channel(ljh_filename, noise_filename) should return a `MassChannel`"
function setup_channel(ljh_filename, noise_filename)
  # TODO:
  # pass in HDF5 group, need to close it on exit
  # pass in "init dictionary"?

    mc=MassChannel()
    mc[:ljh] = LJHGroup(ljh_filename)
    mc[:pretrig_nsamples] = LJH.pretrig_nsamples(mc[:ljh])
    mc[:samples_per_record] = LJH.record_nsamples(mc[:ljh])
    mc[:frametime] = LJH.frametime(mc[:ljh])
    mc[:channel_number] = LJH.channel(mc[:ljh])
    mc[:last_timestamp_checked] = 0.0 # this is updated via a total HACK, its hardcoded into GetPulsesStep

    # look for an init_hdf5 file, extract some values into init_dict
    # if init_dict is empty, get(mc.task) should show fail with a somewhat descriptive error
    init_dict = get_init_dict(mc[:channel_number])
    if isempty(init_dict)
	return mc
    end
    for (k,v) in init_dict
	mc[k]=v
    end
    mc[:pretrig_mean] = RunningVector(Float32)
    mc[:pretrig_rms] = RunningVector(Float32)
    mc[:pulse_average] = RunningVector(Float32)
    mc[:pulse_rms] = RunningVector(Float32)
    mc[:rise_time] = RunningVector(Float32)
    mc[:postpeak_deriv] = RunningVector(Float32)
    mc[:peak_index] = RunningVector(UInt16)
    mc[:peak_value] = RunningVector(UInt16)
    mc[:min_value] = RunningVector(UInt16)
    mc[:filt_value] = RunningVector(Float32)
    mc[:filt_phase] = RunningVector(Float32)
    mc[:energy] = RunningVector(Float32)
    mc[:energy_from_pulse_rms] = RunningVector(Float32)
    mc[:pulse] = RunningVector(Vector{UInt16})
    mc[:rowcount] = RunningVector(Int)
    mc[:timestamp_posix_usec] = RunningVector(Int)
    mc[:selection_good] = RunningSumBitVector()
    mc[:filt_value_hist] = CountRateHistogram(10_000_000, 0:0.5:1500,"filt_value")
    mc[:energy_hist] = CountRateHistogram(10_000_000, 0:0.5:1500,"energy_from_filt_value")
    mc[:energy_from_pulse_rms_hist] = CountRateHistogram(10_000_000, 0:0.5:1500,"energy_from_pulse_rms")

    mc[:noise_filename]=ascii(noise_filename)
    mc[:ljh_filename]=ascii(ljh_filename) # h5py can't read UFT8String written by julia
    mc[:name] = "summarize and filter test"
    mc[:hdf5_filename] = "$(splitext(ljh_filename)[1])_jl.hdf5"
    mc[:init_hdf5_filename] = "$(splitext(ljh_filename)[1])_init_jl.hdf5"

    mc[:oncleanfinish] = markhdf5oncleanfinish
    isfile(mc[:hdf5_filename]) && rm(mc[:hdf5_filename])

    steps = AbstractStep[
    GetPulsesStep(mc[:ljh], [:pulse, :rowcount, :timestamp_posix_usec], 0,100, 0)
    @perpulse pretrig_mean, pretrig_rms, pulse_average, pulse_rms, rise_time, postpeak_deriv, peak_index, peak_value, min_value = compute_summary(pulse, pretrig_nsamples, frametime)
    @threshold calibration_pulse_rms = calibrate_xy(calibration_pulse_rms_pulseheights, calibration_pulse_rms_energies) when length(peak_index)>0
    @threshold calibration_filt_value = calibrate_xy(calibration_filt_value_dc_pulseheights, calibration_filt_value_dc_energies) when length(peak_index)>0 # calibration from filt value dc should work fine for filt value, probably still change this in future
    @perpulse selection_good = selectfromcriteria(pretrig_rms, pretrig_rms_criteria, peak_index, peak_index_criteria, postpeak_deriv, postpeak_deriv_criteria)
    @perpulse filt_value, filt_phase = filter5lag(filter, pulse)
    @perpulse energy = apply_calibration(calibration_filt_value, filt_value)
    @perpulse energy_from_pulse_rms = apply_calibration(calibration_pulse_rms, pulse_rms)
    @histogram update_histogram_send!(filt_value_hist, selection_good, filt_value, timestamp_posix_usec, last_timestamp_checked, channel_number)
    @histogram update_histogram_send!(energy_hist, selection_good, energy, timestamp_posix_usec, last_timestamp_checked, channel_number)
    ToJLDStep([:filt_value, :energy, :filt_phase, :pretrig_rms, :postpeak_deriv, :rise_time, :peak_index,
    :pretrig_mean, :pulse_average, :pulse_rms, :peak_value, :min_value, :rowcount],
    Pair[:filter=>"filter/filter", :f_3db=>"filter/f_3db", :frametime=>"filter/frametime", :noise_autocorr=>"filter/noise_autocorr", :average_pulse=>"filter/average_pulse",
    :average_pulse=>"average_pulse",
    :samples_per_record=>"samples_per_record", :frametime=>"frametime", :pretrig_nsamples=>"pretrig_nsamples",
    :ljh_filename=>"ljh_filename", :noise_filename=>"noise_filename",
  	:peak_index_criteria=>"selection_criteria/peak_index", :pretrig_rms_criteria=>"selection_criteria/pretrig_rms", :postpeak_deriv_criteria=>"selection_criteria/postpeak_deriv",
  	:pretrigger_mean_drift_correction_params=>"pretrigger_mean_drift_correction_params"],
    mc[:hdf5_filename])
    MemoryLimitStep(Int(4e6)) # throw error if mc uses more than 4 MB
    FreeMemoryStep()
    ] # end steps
    setsteps!(mc, steps)
    preptasks!(mc)
    mc
end
get_init_hdf5_filename_groupname(channel_number) = "ssrl_realtime_init_20160602.hdf5","chan$channel_number"
function get_init_dict(channel_number::Integer)
	init_hdf5_filename, groupname = get_init_hdf5_filename_groupname(channel_number)
	d = Dict{Symbol, Any}()
	h5open(init_hdf5_filename,"r") do h5
	    if !exists(h5,groupname)
		return nothing
	    end
	    grp = h5[groupname]
	    d[:calibration_pulse_rms_energies] = read(grp["calibrations/calibration_pulse_rms_energies"])
	    d[:calibration_pulse_rms_pulseheights] = read(grp["calibrations/calibration_pulse_rms_pulseheights"])
	    d[:calibration_filt_value_dc_energies] = read(grp["calibrations/calibration_filt_value_dc_energies"])
	    d[:calibration_filt_value_dc_pulseheights] = read(grp["calibrations/calibration_filt_value_dc_pulseheights"])
	    d[:filter] = read(grp["filter/filter"])
	    d[:pretrigger_mean_drift_correction_params] = LinearDriftCorrect(read(grp["pretrigger_mean_drift_correction_params"])...)
	    d[:peak_index_criteria] = read(grp["selection_criteria"]["peak_index"])
	    d[:pretrig_rms_criteria] = read(grp["selection_criteria"]["pretrig_rms"])
	    d[:postpeak_deriv_criteria] = read(grp["selection_criteria"]["postpeak_deriv"])
	    return nothing
	end # do
	return d
end # get_init_dict

if true # change this to false if you want to type in the filename
  ljhname, writingbool = LJHUtil.matter_writing_status()
else
  ljhname = "put filename here"
  writingbool = false
end
println("***********************")
@show ljhname
@show maxchannels = 240
masschannels = Dict()
t0 = time()
channums = LJHUtil.allchannels(ljhname)
channums = channums[1:min(maxchannels,length(channums))]
ljh_filenames = [LJHUtil.fnames(ljhname,channum) for channum in channums]
noise_filenames = ljh_filenames
for (i, channum) in enumerate(channums)
  masschannels[channums[i]] = setup_channel(ljh_filenames[i], noise_filenames[i])
end
tf = time()
println("setup $(length(channums)) channels in $(tf-t0) seconds")
tschedule = @elapsed map(schedule, values(masschannels))
println("launched $(length(channums)) channels in $tschedule seconds")
# stop tasks for previous file once they've finished their work
println("allowing analysis tasks for $oldljhname to end")

while writingbool
  watch(watch_file(LJHUtil.sentinel_file_path,10))
  ljhname, writingbool = LJHUtil.matter_writing_status()
end

map(plantoend, values(masschannels)) # tell each channel it can end


text:
using SearchBenchmarks
using Test
import Random

function module_functions(modname)
    list = Symbol[]
    for nm in names(modname)
        typeof(eval(nm)) == Function && push!(list,nm)
    end
    return list
end

@testset "SearchBenchmarks.jl" begin
    # test constants
    test_textsize = 10000 # length of searchable array
    test_highestbyte = 0x7f # highest value in alphabet (must be UInt8)
    test_patterns = 8 # no.of pattern sizes to test
    test_seed = 2345 # random generator seed
    # get search functions
    searches = all_search_functions()
    @info "search functions are: ",searches

    # compare search results of all functions
    f0 = pop!(searches)
    s = Vector{UInt8}(undef,test_textsize)
    Random.seed!(test_seed)
    Random.rand!(s,0x00:test_highestbyte) # array to search in
    psize = 1
    for pi in 1:test_patterns
        #@info "pattern size", psize
        t = s[test_textsize-psize+1:test_textsize]
        pos0 = 0
        while true
            pos = pos0 + 1
            pos0 = f0(s,t,pos)
            for f in searches
                posf = f(s,t,pos)
                if pos0!=posf
                    # error! report params
                    @error "test fails on pattern size $psize:  $f0(s,t,$pos) == $f(s,t,$pos)"
                end
                @test pos0 == posf
            end
            if pos0 == 0
                break
            end
        end
        psize = psize*3-1 # sequence 1 2 5 14 41 122 ...
    end
    # test benchmarking
    benchmark(5, 100, "0x7f,1234", searches, "SearchBenchmarks_test_result",false)
    benchmark(5, 100, "SearchBenchmarks_test_result.png", searches, "SearchBenchmarks_test_result2",true)
end


text:
module TestDoctest

using Documenter
using Try
using TryExperimental

function test_try()
    VERSION < v"1.7" && return
    doctest(Try; manual = false)
end

function test_tryexperimental()
    VERSION < v"1.7" && return
    doctest(TryExperimental; manual = false)
end

end  # module


text:
# Julia wrapper for header: /usr/include/scip/pub_reopt.h
# Automatically generated using Clang.jl wrap_c

function SCIPreoptnodeGetNVars(reoptnode)
    ccall((:SCIPreoptnodeGetNVars, libscip), Cint, (Ptr{SCIP_REOPTNODE},), reoptnode)
end

function SCIPreoptnodeGetNConss(reoptnode)
    ccall((:SCIPreoptnodeGetNConss, libscip), Cint, (Ptr{SCIP_REOPTNODE},), reoptnode)
end

function SCIPreoptnodeGetNDualBoundChgs(reoptnode)
    ccall((:SCIPreoptnodeGetNDualBoundChgs, libscip), Cint, (Ptr{SCIP_REOPTNODE},), reoptnode)
end

function SCIPreoptnodeGetNChildren(reoptnode)
    ccall((:SCIPreoptnodeGetNChildren, libscip), Cint, (Ptr{SCIP_REOPTNODE},), reoptnode)
end

function SCIPreoptnodeGetLowerbound(reoptnode)
    ccall((:SCIPreoptnodeGetLowerbound, libscip), Cdouble, (Ptr{SCIP_REOPTNODE},), reoptnode)
end

function SCIPreoptnodeGetType(reoptnode)
    ccall((:SCIPreoptnodeGetType, libscip), SCIP_REOPTTYPE, (Ptr{SCIP_REOPTNODE},), reoptnode)
end

function SCIPreoptnodeGetSplitCons(reoptnode, vars, vals, constype, conssize, nvars)
    ccall((:SCIPreoptnodeGetSplitCons, libscip), Cvoid, (Ptr{SCIP_REOPTNODE}, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Ptr{REOPT_CONSTYPE}, Cint, Ptr{Cint}), reoptnode, vars, vals, constype, conssize, nvars)
end

function SCIPreoptnodeGetConss(reoptnode, vars, bounds, boundtypes, mem, nconss, nvars)
    ccall((:SCIPreoptnodeGetConss, libscip), Cvoid, (Ptr{SCIP_REOPTNODE}, Ptr{Ptr{Ptr{SCIP_VAR}}}, Ptr{Ptr{Cdouble}}, Ptr{Ptr{SCIP_BOUNDTYPE}}, Cint, Ptr{Cint}, Ptr{Cint}), reoptnode, vars, bounds, boundtypes, mem, nconss, nvars)
end

function SCIPreoptnodeSetParentID(reoptnode, parentid)
    ccall((:SCIPreoptnodeSetParentID, libscip), Cvoid, (Ptr{SCIP_REOPTNODE}, UInt32), reoptnode, parentid)
end

function SCIPreoptGetNRestartsGlobal(reopt)
    ccall((:SCIPreoptGetNRestartsGlobal, libscip), Cint, (Ptr{SCIP_REOPT},), reopt)
end

function SCIPreoptGetNRestartsLocal(reopt)
    ccall((:SCIPreoptGetNRestartsLocal, libscip), Cint, (Ptr{SCIP_REOPT},), reopt)
end

function SCIPreoptGetNTotalRestartsLocal(reopt)
    ccall((:SCIPreoptGetNTotalRestartsLocal, libscip), Cint, (Ptr{SCIP_REOPT},), reopt)
end

function SCIPreoptGetFirstRestarts(reopt)
    ccall((:SCIPreoptGetFirstRestarts, libscip), Cint, (Ptr{SCIP_REOPT},), reopt)
end

function SCIPreoptGetLastRestarts(reopt)
    ccall((:SCIPreoptGetLastRestarts, libscip), Cint, (Ptr{SCIP_REOPT},), reopt)
end

function SCIPreoptGetNFeasNodes(reopt)
    ccall((:SCIPreoptGetNFeasNodes, libscip), Cint, (Ptr{SCIP_REOPT},), reopt)
end

function SCIPreoptGetNTotalFeasNodes(reopt)
    ccall((:SCIPreoptGetNTotalFeasNodes, libscip), Cint, (Ptr{SCIP_REOPT},), reopt)
end

function SCIPreoptGetNPrunedNodes(reopt)
    ccall((:SCIPreoptGetNPrunedNodes, libscip), Cint, (Ptr{SCIP_REOPT},), reopt)
end

function SCIPreoptGetNTotalPrunedNodes(reopt)
    ccall((:SCIPreoptGetNTotalPrunedNodes, libscip), Cint, (Ptr{SCIP_REOPT},), reopt)
end

function SCIPreoptGetNCutoffReoptnodes(reopt)
    ccall((:SCIPreoptGetNCutoffReoptnodes, libscip), Cint, (Ptr{SCIP_REOPT},), reopt)
end

function SCIPreoptGetNTotalCutoffReoptnodes(reopt)
    ccall((:SCIPreoptGetNTotalCutoffReoptnodes, libscip), Cint, (Ptr{SCIP_REOPT},), reopt)
end

function SCIPreoptGetNInfNodes(reopt)
    ccall((:SCIPreoptGetNInfNodes, libscip), Cint, (Ptr{SCIP_REOPT},), reopt)
end

function SCIPreoptGetNTotalInfNodes(reopt)
    ccall((:SCIPreoptGetNTotalInfNodes, libscip), Cint, (Ptr{SCIP_REOPT},), reopt)
end


text:
tol = 1e-10
#
# Sanity check.
#
@test size(Planets.mr_table) == (91,42)
@test haskey(Planets.mr_table,:Mearth)
@test haskey(Planets.mr_table,:rocky)
@test haskey(Planets.mr_table,Symbol("5%fe"))
@test haskey(Planets.mr_table,Symbol("5%h2o"))

@test isa(Planets.mr_table[:Mearth][1], Float64)
@test Planets.mr_table[:Mearth][1] == 0.0625

#
# Test grid points, including edge cases.
#
@test abs(core_radius(0.07179, fe=0.0 )  -  0.4888) < tol
@test abs(core_radius(0.07179, h2o=0.0)  -  0.4888) < tol

@test abs(core_radius(0.07179, fe=1.0 )  -  0.3733) < tol
@test abs(core_radius(0.07179, fe=0.5 )  -  0.4394) < tol
@test abs(core_radius(0.07179, fe=0.05)  -  0.4841) < tol
@test abs(core_radius(0.07179, h2o=.05)  -  0.5041) < tol
@test abs(core_radius(0.07179, h2o=1.0)  -  0.6588) < tol

@test abs(core_radius(1.0, fe=1.0)  -  0.8228) < tol
@test abs(core_radius(1.0, fe=0.0)  -  1.0667) < tol
@test abs(core_radius(1.0, fe=.05)  -  1.0568) < tol
@test abs(core_radius(1.0, fe=0.1)  -  1.0466) < tol
@test abs(core_radius(1.0, fe=0.2)  -  1.0260) < tol
@test abs(core_radius(1.0, fe=0.3)  -  1.0050) < tol
@test abs(core_radius(1.0, fe=0.4)  -  0.9834) < tol
@test abs(core_radius(1.0, fe=.45)  -  0.9723) < tol
@test abs(core_radius(1.0, fe=.35)  -  0.9943) < tol
@test abs(core_radius(1.0, fe=1.0)  -  0.8228) < tol
#
# Test interpolation across columns.
#
@test abs(core_radius(1.0, fe=0.975)  -  0.83350) < tol
@test abs(core_radius(1.0, fe=0.025)  -  1.06175) < tol
#
# Test interpolation across rows.
#
@test abs(core_radius(1.0359, fe=1.0)  -  0.8307) < tol
#
# Test interpolation across columns AND rows.
#
@test abs(core_radius(1.035900, fe =0.975)  -  0.841525) < tol
@test abs(core_radius(0.064745, h2o=0.025)  -  0.480450) < tol
#
# Test extrapolation of mass.
#
M_small = 0.0625
R_small = core_radius(M_small, fe=0.3)
M_large = 32.0
R_large = core_radius(M_large, fe=0.3)
@test abs(core_radius(0.01, fe=0.3)  -  R_small * (0.01/M_small)^(1/3.7)) < tol
@test abs(core_radius(35.0, fe=0.3)  -  R_large * (35.0/M_large)^(1/3.7)) < tol