text:
# # Example 002: Solid Body Rotation
#
#md # !!! jupyter
#md #     This example is also available as a Jupyter notebook:
#md #     [`ex_002_solid_body_rotation.ipynb`](@__NBVIEWER_ROOT_URL__examples/DGmethods/generated/ex_002_solid_body_rotation.html)
#
# Key ideas of this tutorial:
#   - Setting up auxiliary state variables
#   - Defining the boundary condition treatment
#
# ## Introduction
#
# In this example we will solve the variable coefficient advection equation. The
# velocity field used is solid body rotation where the domain is the square or
# domain $\Omega = [-1, 1]^{d}$ where $d=2$ or $3$.
#
# The partial differential equation we wish to solve is
# ```math
# \frac{\partial q}{\partial t} + \nabla \cdot (\vec{u} q) = 0,
# ```
# where $q$ is the advected field and the velocity field is $\vec{u} = 2\pi r
# (-\sin(\theta), \cos(\theta), 0)^{T}$ with $r = \sqrt{x^2 + y^2}$ and $\theta
# = \arctan(y / x)$.
#
# The quantity $\vec{u} q$ is more generally called the flux and denoted in the
# tutorial below as $\boldsymbol{F}(q) = \vec{u} q$.
#
# Below is a program interspersed with comments.
#md # The full program, without comments, can be found in the next
#md # [section](@ref ex_002_solid_body_rotation-plain-program).
#
# ## Commented Program

#------------------------------------------------------------------------------

# ### Preliminaries
# Load in modules needed for solving the problem
using MPI
using CLIMA.Topologies
using CLIMA.Grids
using CLIMA.DGBalanceLawDiscretizations
using CLIMA.MPIStateArrays
using CLIMA.LowStorageRungeKuttaMethod
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using CLIMA.Vtk
using LinearAlgebra
using Logging
using Dates
using Printf
using StaticArrays
# Start up MPI if this has not already been done
MPI.Initialized() || MPI.Init()
#md nothing # hide

#------------------------------------------------------------------------------

# ### Initializing the Velocity Field
# The key difference between this example and [example
# 001](ex_001_periodic_advection.html) is that in this case we have a
# non-constant velocity field. In the balance law solver, in addition to the
# PDE state at every degree of freedom we can also define a constant in time
# *auxiliary state,* and it is in this auxiliary state that we will store the
# velocity field.
#
# Initialization of the auxiliary state can happen in several ways, but here we
# will use the default `DGBalanceLaw` initialization interface which requires
# the user-defined function which given `x`, `y`, and `z` defines the auxiliary
# state
const num_aux_states = 3
function velocity_initilization!(uvec::MVector{num_aux_states, DFloat},
                                 x, y, z) where DFloat
  @inbounds begin
    r = hypot(x, y)
    θ = atan(y, x)
    uvec .= 2DFloat(π) * r .* (-sin(θ), cos(θ), 0)
  end
end
#md nothing # hide
# Note: We have caught the type of the elements in order the properly cast $\pi$
# since `2pi` would be default be a `Float64`.
#
# Remark: Though not needed for this problem, if the user wishes to have access
# to the coordinate points during the simulation these should be stored in the
# auxiliary state which would increase the size of the auxiliary state.

#------------------------------------------------------------------------------

# ### Physical Flux
#md # Now we define a function which given a value for $q$ computes the physical
#md # flux $\boldsymbol{F} = \vec{u} q$. The balance law solver will will pass
#md # user-defined auxiliary state at a degree of freedom through to the flux
#md # function as the fourth argument; the third and fifth arguments which are
#md # not needed for this example is the viscous state and simulation time).
function advectionflux!(F, state, _, uvec, _)
  DFloat = eltype(state) # get the floating point type we are using
  @inbounds begin
    q = state[1]
    F[:, 1] = uvec * q
  end
end
#md nothing # hide

#------------------------------------------------------------------------------

# ### Numerical Flux
# As in [example 001](ex_001_periodic_advection.html) we will use an upwind
# numerical flux; more discussion of this can be seen in the [numerical
# flux](ex_001_periodic_advection.html#numerical_flux-1) section of example 001.
#
# The auxiliary state for the minus and plus sides of the interface will be
# passed in through arguments 4 and 6 of the numerical flux callback. Since the
# two sides of the interface are collocated the auxiliary state on the two sides
# should be the same.
function upwindflux!(fs, nM, stateM, viscM, uvecM, stateP, viscP, uvecP, t)
  DFloat = eltype(fs)
  @inbounds begin
    ## determine the advection speed and direction
    un = dot(nM, uvecM)
    qM = stateM[1]
    qP = stateP[1]
    ## Determine which state is "upwind" of the minus side
    fs[1] = un ≥ 0 ? un * qM : un * qP
  end
end

#------------------------------------------------------------------------------

# ### Boundary Numerical Flux
# Since we will not assume that the domain is periodic, we also need to define a
# boundary numerical flux which will be used to define the boundary conditions.
# Generally speaking, boundary conditions for purely hyperbolic problems should
# be imposed by relating incoming characteristic variables to outgoing
# characteristics. In this case we will use characteristic outflow boundary
# conditions and zero inflow boundary conditions.
#
# The syntax of the boundary flux is almost the same as the numerical flux
# except that the boundary condition type is passed into the function, though in
# this case we can neglect the value of the boundary condition flag. In the case
# of boundary conditions the plus state is set to the minus side state; this is
# done since in the case of model coupling this could be set to some values
# derived from the neighbouring model.
#
# In the case of advection with the outflow boundary condition and zero inflow,
# the boundary numerical flux is the same as the upwind flux except with $q^{+}$
# set to zero; more complicated PDES and boundary conditions would require more
# complex constructions.
function upwindboundaryflux!(fs, nM, stateM, viscM, uvecM, stateP, viscP, uvecP,
                             bctype, t)
  DFloat = eltype(fs)
  @inbounds begin
    ## determine the advection speed and direction
    un = dot(nM, uvecM)
    qM = stateM[1]
    ## Determine which state is "upwind" of the minus side
    fs[1] = un ≥ 0 ? un * qM : 0
  end
end
#md nothing # hide

#------------------------------------------------------------------------------

# ### Initial Condition
# In this example we take the initial condition to be
# ```math
# q(\vec{x}, t=0) =
# \exp\left(-\left(8\left\|\vec{x}-\frac{1}{2}\vec{e}_{1}\right\|_2\right)^2\right)
# ```
# where $\vec{e}_{1} = (1, 0, 0)^{T}$.
#
# Note: The initial condition will always be called as though the dimensionality
# of the problem is 3. For the domain used below `z = 0` when the problem is
# actually two-dimensional and thus the `hypot` call before is not effected by
# `z`
#
# Note: When the balance law solver calls the initial condition function the
# auxiliary state, in this case the velocity field, will also be included since
# the number of auxiliary variables is greater than zero.
function initialcondition!(Q, x, y, z, _)
  @inbounds Q[1] = exp(-(8 * hypot(x - 1//2, y, z))^2)
end
#md nothing # hide

#------------------------------------------------------------------------------

# ### Exact Solution
# For solid body rotation the exact solution is computed by tracing back the
# rotation to the initial state.
#
# Note: `uvec` is included to match calling convention of `initialcondition!`
function exactsolution!(Q, t, x, y, z, uvec)
  @inbounds begin
    DFloat = eltype(Q)

    r = hypot(x, y)
    θ = atan(y, x) - 2DFloat(π) * t

    x, y = r * cos(θ), r * sin(θ)

    initialcondition!(Q, x, y, z, uvec)
  end
end
#md nothing # hide

#------------------------------------------------------------------------------

# ### Initialize the DG Method
# The initialization of the DG method is largely the same as the
# [intialization](ex_001_periodic_advection.html#Initial-Condition-1) discussion
# of [ex 001](ex_001_periodic_advection.html).
function setupDG(mpicomm, dim, Ne, polynomialorder, DFloat=Float64,
                 ArrayType=Array)

  @assert ArrayType === Array

  brickrange = (range(DFloat(-1); length=Ne+1, stop=1),
                range(DFloat(-1); length=Ne+1, stop=1),
                range(DFloat(-1); length=Ne+1, stop=1))

  # By default the `BrickTopology` is not periodic, so unlike ex 001, we do not
  # need to specify the periodicity
  topology = BrickTopology(mpicomm, brickrange[1:dim])

  grid = DiscontinuousSpectralElementGrid(topology; polynomialorder =
                                          polynomialorder, FloatType = DFloat,
                                          DeviceArray = ArrayType,)

  # Note the additional keyword arguments: `numerical_boundary_flux!`
  # which is used to pass the numerical flux function that implements the
  # boundary condition, `auxiliary_state_length` which defines the number of
  # auxiliary state fields at each degree of freedom, and
  # `auxiliary_state_initialization!` which initializes the auxiliary state.
  spatialdiscretization = DGBalanceLaw(grid = grid, length_state_vector = 1,
                                       flux! = advectionflux!,
                                       numerical_flux! = upwindflux!,
                                       numerical_boundary_flux! =
                                       upwindboundaryflux!,
                                       auxiliary_state_length = num_aux_states,
                                       auxiliary_state_initialization! =
                                       velocity_initilization!)

end
#md nothing # hide

#------------------------------------------------------------------------------

# ### Initializing and run the DG method
# This `let` statement is largely the same as the [Using ODE solver callback
# functions](ex_001_periodic_advection.html#Using-ODE-solver-callback-functions-1)
# block from ex 001. Difference are highlighted.
let
  mpicomm = MPI.COMM_WORLD
  mpi_logger = ConsoleLogger(MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull)
  dim = 2
  Ne = 20
  polynomialorder = 4
  spatialdiscretization = setupDG(mpicomm, dim, Ne, polynomialorder)
  Q = MPIStateArray(spatialdiscretization, initialcondition!)
  filename = @sprintf("initialcondition_mpirank%04d", MPI.Comm_rank(mpicomm))
  writevtk(filename, Q, spatialdiscretization,
                                       ("q",))

  h = 1 / Ne
  # Since we are on the $[-1, 1]^d$ domain, the maximum velocity will by $2\pi$,
  # thus this defines the CFL restriction
  CFL = h / (2π)
  dt = CFL / polynomialorder^2
  lsrk = LowStorageRungeKutta(spatialdiscretization, Q; dt = dt, t0 = 0)
  finaltime = 1.0
  if (parse(Bool, lowercase(get(ENV,"TRAVIS","false")))       #src
      && "Test" == get(ENV,"TRAVIS_BUILD_STAGE_NAME","")) ||  #src
    parse(Bool, lowercase(get(ENV,"APPVEYOR","false")))       #src
    finaltime = 2dt                                           #src
  end                                                         #src

  # For simplicity we only include the vtk callback

  vtk_step = 0
  mkpath("vtk")
  cb_vtk = GenericCallbacks.EveryXSimulationSteps(20) do
    vtk_step += 1
    filename = @sprintf("vtk/solid_body_rotation_mpirank%04d_step%04d",
                         MPI.Comm_rank(mpicomm), vtk_step)
    writevtk(filename, Q, spatialdiscretization,
                                         ("q",))
    nothing
  end

  solve!(Q, lsrk; timeend = finaltime, callbacks = (cb_vtk, ))

  filename = @sprintf("finalsolution_mpirank%04d", MPI.Comm_rank(mpicomm))
  writevtk(filename, Q, spatialdiscretization,
                                       ("q",))

  # As with the initial condition, we need to catch the auxiliary state `uvec`
  # in this initialization call.
  Qe = MPIStateArray(spatialdiscretization) do Qin, x, y, z, uvec
    exactsolution!(Qin, finaltime, x, y, z, uvec)
  end

  error = euclidean_distance(Q, Qe)
  with_logger(mpi_logger) do
    @info @sprintf("""Run with
                   dim              = %d
                   Ne               = %d
                   polynomial order = %d
                   error            = %e
                   """, dim, Ne, polynomialorder, error)
  end
end
#md nothing # hide

#------------------------------------------------------------------------------

# ### Computing rates and errors
# As with ex 001, since the analytic solution is known we can compute the rate
# of convergence of the scheme
let
  mpicomm = MPI.COMM_WORLD
  mpi_logger = ConsoleLogger(MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull)

  dim = 2
  polynomialorder = 4
  finaltime = 1.0

  with_logger(mpi_logger) do
    @info @sprintf("""Running with
                   dim              = %d
                   polynomial order = %d
                   """, dim, polynomialorder)
  end

  base_Ne = 5
  lvl_error = zeros(4) # number of levels to compute is length(lvl_error)
  for lvl = 1:length(lvl_error)
    ## `Ne` for this mesh level
    Ne = base_Ne * 2^(lvl-1)
    spatialdiscretization = setupDG(mpicomm, dim, Ne, polynomialorder)

    Q = MPIStateArray(spatialdiscretization, initialcondition!)
    h = 1 / Ne
    CFL = h / (2π)
    dt = CFL / polynomialorder^2
    if (parse(Bool, lowercase(get(ENV,"TRAVIS","false")))       #src
        && "Test" == get(ENV,"TRAVIS_BUILD_STAGE_NAME","")) ||  #src
      parse(Bool, lowercase(get(ENV,"APPVEYOR","false")))       #src
      finaltime = 2dt                                           #src
    end                                                         #src
    lsrk = LowStorageRungeKutta(spatialdiscretization, Q; dt = dt, t0 = 0)

    solve!(Q, lsrk; timeend = finaltime)

    Qe = MPIStateArray(spatialdiscretization) do Qin, x, y, z, uvec
      exactsolution!(Qin, finaltime, x, y, z, uvec)
    end

    lvl_error[lvl] = euclidean_distance(Q, Qe)
    msg =  @sprintf   "Level      = %d" lvl
    msg *= @sprintf "\nNe               = %d" Ne
    msg *= @sprintf "\nerror            = %.4e" lvl_error[lvl]
    if lvl > 1
      rate = log2(lvl_error[lvl-1]) - log2(lvl_error[lvl])
      msg *= @sprintf "\nconvergence rate = %.4e" rate
    end
    with_logger(mpi_logger) do
      @info msg
    end
  end
end
#md nothing # hide

#------------------------------------------------------------------------------

# ### Finalizing MPI (if necessary)
Sys.iswindows() || MPI.finalize_atexit()
Sys.iswindows() && !isinteractive() && MPI.Finalize()
#md nothing # hide

#md # ## [Plain Program](@id ex_002_solid_body_rotation-plain-program)
#md #
#md # Below follows a version of the program without any comments.
#md # The file is also available here: [ex\_002\_solid\_body\_rotation\_periodic\_advection.jl](ex_002_solid_body_rotation.jl)
#md #
#md # ```julia
#md # @__CODE__
#md # ```

max_stars_repo_name: jkozdon/CLIMA

max_issues_repo_name: jkozdon/CLIMA

max_forks_repo_name: jkozdon/CLIMA


text:
# using Revise
using PiecewiseDeterministicMarkovProcesses, LinearAlgebra, Random, DifferentialEquations, Sundials
const PDMP = PiecewiseDeterministicMarkovProcesses

function AnalyticalSample(xc0,xd0,ti,nj::Int64)
	xch = [xc0[1]]
	xdh = [xd0[1]]
	th  = [ti]
	t = ti
	while length(th)<nj
		xc = xch[end]
		xd = xdh[end]
		S = -log(rand())
		a = mod(xd[1],2)==0 ? -1 : 1
		dt = (exp(a*S)-1)*exp(-a*S)/(a*xc)
		t += dt
		push!(th, t)
		push!(xch, xc * exp(a*S) )
		push!(xdh, xd .+ 1 )
		S = -log(rand())
	end
	return th,xch,xdh
end

function F_tcp!(ẋ, xc, xd, parms, t)
	# vector field used for the continuous variable
	if mod(xd[1], 2)==0
		 ẋ[1] = 1.
	else
		 ẋ[1] = -1.
	end
	nothing
end

rate_tcp(x) = 1/x

function R_tcp!(rate, xc, xd, parms, t, issum::Bool)
	if issum==false
		rate[1] = rate_tcp(xc[1])
		rate[2] = 0.0
		return 0., 100.
	else
		return rate_tcp(xc[1]), 100.
	end
end

xc0 = [ 1.0 ]
xd0 = [0, 1]

nu_tcp = [1 0;0 -1]
parms = [0.0]
tf = 100000.
nj = 100

Random.seed!(1234)
	res_a = AnalyticalSample(xc0,xd0,0.,nj)

errors = Float64[]

println("\n\nComparison of solvers")
	for ode in [(:cvode,"cvode"),(:lsoda,"lsoda"),(CVODE_BDF(),"CVODEBDF"),(CVODE_Adams(),"CVODEAdams"),(Rosenbrock23(),"RS23"),(Tsit5(),"tsit5"),(Rodas4P(autodiff=true),"rodas4P-AutoDiff"),(Rodas5(),"rodas5"),(AutoTsit5(Rosenbrock23()),"AutoTsit5RS23")]
		Random.seed!(1234)
		problem = PDMP.PDMPProblem(F_tcp!, R_tcp!, nu_tcp, xc0, xd0, parms, (0.0, tf))
		res =  @time PDMP.solve(problem, CHV(ode[1]); n_jumps = nj)
		printstyled(color=:green,"--> norm difference = ", norm(res.time - res_a[1], Inf64), "  - solver = ",ode[2],"\n\n")
		push!(errors,norm(res.time - res_a[1],Inf64))
end

# plot!(res.time,res.xc[:,1])
# case with no allocations  0.000721 seconds (330 allocations: 26.266 KiB)
# Random.seed!(1234)
# 	problem = PDMP.PDMPProblem(F_tcp!, R_tcp!, nu_tcp, xc0, xd0, parms, (0.0, 1e19))
# 	res =  @time PDMP.solve(problem, CHV(Tsit5()); n_jumps = 4nj, save_positions = (false, false))

# res =  @timed PDMP.solve(problem, CHV(Tsit5()); n_jumps = nj, save_positions = (false, false))
# res[end].poolalloc

# # Random.seed!(1234)
# #	 using PiecewiseDeterministicMarkovProcesses
# #	 println("\n"*"+"^40)
# #	 res =  @time PiecewiseDeterministicMarkovProcesses.pdmp!(xc0, xd0, F_tcp!, R_tcp!, nu_tcp, parms, 0.0, tf, n_jumps = 10,   ode =Tsit5(), algo=:rejection, verbose=true)
# #
# #	 res.time |> println

ext: jl


text:
#############################################################################
## pretty_print_xml()
function pretty_print_xml(xmlstr::String)
    # convert a utf-8 string `xmlstr` to a vector of unicode chars so we can
    # index into it easily
    xml_chars = Char[]
    i = 1
    while i <= ncodeunits(xmlstr)
        push!(xml_chars, xmlstr[i])
        i = nextind(xmlstr, i)
    end

    # print the tree
    ele_stack = [("string", -1, -1)]
    i = 1
    while i <= length(xml_chars)
        #skip whitespaces
        while isspace(xml_chars[i]) i += 1 end
        xml_chars[i] != '<' && error("Unknown token: ", xml_chars[i])
        i = pretty_print_element(xml_chars, i, 0, ele_stack)
    end
end

const tbsz = 3
function pretty_print_element(
    xml::Vector{Char}, start::Int, indent::Int, ele_stack::Vector{Tuple{String, Int, Int}}
)
    # mml[start] must be '<'
    t_start, t_end, i, ele_name = start, -1, start + 1, ""
    hasChild::Bool = false
    # get starting tag
    while xml[i] != '>' && xml[i] != ' '
        ele_name *= xml[i]
        i += 1
    end
    # move the index to the position right behind '>'
    if xml[i] == '>'
        t_end = i
        i += 1
    else
        # xml[i] == ` `, there are attributes following the tag name
        while xml[i] != '>'; i += 1; end
        t_end = i
        i += 1
    end

    push!(ele_stack, (ele_name, t_start, t_end))
    #print(" "^indent*"<", ele_name, '>')
    if xml[i - 2] == '/'
        # simple leaf elements <ele ... /> are printed in one line
        print(" "^indent*String(xml[t_start:t_end]))
        #while xml[i] != '<'
        #    print(mml[i])
        #    i += 1
        #end
        # Now mml[i] == '<' in the ending tag of <mi>, <mn>, <mo>, <mspace>
        #while mml[i] != '>' i += 1 end
        #i += 1
        i <= length(xml) && while isspace(xml[i]) i += 1 end
        #print("</", ele_name, ">")
        #println()
        pop!(ele_stack)
        return i
    else
        # the starting tag name, we do not know yet if this element has a child
        print(" "^indent*"<", ele_name, '>')
        # keep finding the next tag and print
        while true
            while isspace(xml[i]) i += 1 end
            i > length(xml) && break

            if xml[i] == '<' && xml[i+1] != '/'
                # a child element
                println() # printing child starts a new line
                hasChild = true
                i = pretty_print_element(xml, i, indent + tbsz, ele_stack)
            elseif xml[i] == '<' && xml[i+1] == '/'
                # we encounter the closing tag of an element
                i += 2
                en2 = ""
                while xml[i] != '>'
                    en2 *= xml[i]
                    i += 1
                end
                en = pop!(ele_stack)
                # we could have sth like <mspace width=...> </mspace>
                !startswith(en[1], en2) && begin
                    println(ele_stack)
                    error("tag not match: ", en, " -- ", en2)
                end
                # xml[i] == '>'
                if hasChild
                    println()
                    println(" "^indent*"</", en2, '>')
                else
                    println("</", en2, '>')
                end
                i += 1
                #i > length(xml) && return i
                #while isspace(xml[i]) i += 1 end
                return i
            else
                # xml[i] != `<`, we encounter the content of an element
                content = ""
                while xml[i] != '<'
                    content *= xml[i]
                    i += 1
                end
                if hasChild # some child element is already printed
                    print(" "^indent*content)
                else
                    print(content)
                end
            end
        end
    end
end


text:
abstract type ConvexMinimumDistanceMethod <: EBayesMethod end

"""
    NPMLE(convexclass, solver) <: Empirikos.EBayesMethod

Given ``n`` independent samples ``Z_i`` from the empirical Bayes problem with prior ``G`` known to lie
in the `convexclass` ``\\mathcal{G}``, estimate ``G`` by Nonparametric Maximum Likelihood (NPMLE)
```math
\\widehat{G}_n \\in \\operatorname{argmax}_{G \\in \\mathcal{G}}\\left\\{\\sum_{i=1}^n \\log( f_{i,G}(Z_i)) \\right\\},
```
where ``f_{i,G}(z) = \\int p_i(z \\mid \\mu) dG(\\mu)`` is the marginal density of the ``i``-th sample.
The optimization is conducted by a JuMP compatible `solver`.
"""
struct NPMLE{C} <: ConvexMinimumDistanceMethod
    convexclass::C
    solver::Any
    dict::Any
end

NPMLE(convexclass, solver; kwargs...) = NPMLE(convexclass, solver, kwargs)
NPMLE(;convexclass, solver, kwargs...) = NPMLE(convexclass, solver; kwargs)

function Base.show(io::IO, npmle::NPMLE)
    print(io, "NPMLE with ")
    show(io, npmle.solver)
    println(io, " and 𝒢:")
    show(io, npmle.convexclass)
end

struct FittedConvexMinimumDistance{D,N<:ConvexMinimumDistanceMethod}
    prior::D
    method::N
    model::Any # add status?
end
Base.broadcastable(fitted_method::FittedConvexMinimumDistance) = Ref(fitted_method)

function Base.show(io::IO, fittedmethod::FittedConvexMinimumDistance)
    print(io, "Fitted ")
    show(io, fittedmethod.method)
end

# TODO: replace these types of methods by a fit_if_not_fitted
StatsBase.fit(fitted_method::FittedConvexMinimumDistance, args...; kwargs...) =
    fitted_method

marginalize(Z, fitted_method::FittedConvexMinimumDistance) =
    marginalize(Z, fitted_method.prior)

function (target::EBayesTarget)(fitted_method::FittedConvexMinimumDistance, args...)
    target(fitted_method.prior, args...)
end

# TODO: macro to all of pdf, cdf,...
Distributions.pdf(fitted_method::FittedConvexMinimumDistance, Z) =
    Distributions.pdf(fitted_method.prior, Z)
Distributions.logpdf(fitted_method::FittedConvexMinimumDistance, Z) =
    Distributions.logpdf(fitted_method.prior, Z)

# seems like template that could be useful more generally..
function StatsBase.fit(method::ConvexMinimumDistanceMethod, Zs; kwargs...)
    Zs = summarize_by_default(Zs) ? summarize(Zs) : Zs
    method = set_defaults(method, Zs; kwargs...)
    _fit(method, Zs)
end

function _fit(npmle::NPMLE, Zs)
    @unpack convexclass, solver = npmle
    model = Model(solver)

    π = Empirikos.prior_variable!(model, convexclass)
    f = pdf.(π, Zs)

    _mult = multiplicity(Zs)
    n = length(_mult)

    @variable(model, u)

    @constraint(model, vcat(u, f, _mult) in MathOptInterface.RelativeEntropyCone(2n + 1))
    @objective(model, Min, u)
    optimize!(model)
    estimated_prior = π()
    FittedConvexMinimumDistance(estimated_prior, npmle, model)
end

"""
    KolmogorovSmirnovMinimumDistance(convexclass, solver) <: Empirikos.EBayesMethod

Given ``n`` i.i.d. samples from the empirical Bayes problem with prior ``G`` known to lie
in the `convexclass` ``\\mathcal{G}`` , estimate ``G`` as follows:
```math
\\widehat{G}_n \\in \\operatorname{argmin}_{G \\in \\mathcal{G}}\\{\\sup_{t \\in \\mathbb R}\\lvert F_G(t) - \\widehat{F}_n(t)\\rvert\\},
```
where ``\\widehat{F}_n`` is the ECDF of the samples. The optimization is conducted by a JuMP compatible `solver`.
"""
struct KolmogorovSmirnovMinimumDistance{C} <: ConvexMinimumDistanceMethod
    convexclass::C
    solver::Any
    dict::Any
end

function KolmogorovSmirnovMinimumDistance(convexclass, solver)
    KolmogorovSmirnovMinimumDistance(convexclass, solver, nothing)
end

function KolmogorovSmirnovMinimumDistance(; convexclass, solver)
    KolmogorovSmirnovMinimumDistance(convexclass, solver, nothing)
end

function Base.show(io::IO, ks::KolmogorovSmirnovMinimumDistance)
    print(io, "Kolmogorov-Smirnov Minimum Distance with ")
    show(io, ks.solver)
    println(io, " and 𝒢:")
    show(io, ks.convexclass)
end

function _fit(method::KolmogorovSmirnovMinimumDistance, Zs)
    @unpack convexclass, solver = method
    model = Model(solver)

    π = Empirikos.prior_variable!(model, convexclass)

    dkw = fit(DvoretzkyKieferWolfowitz(;max_constraints=Inf), Zs)

    F = cdf.(π, keys(dkw.summary))
    Fhat = collect(values(dkw.summary))

    @variable(model, u)

    @constraint(model, F - Fhat .<= u)
    @constraint(model, F - Fhat .>= -u)

    @objective(model, Min, u)
    optimize!(model)
    estimated_prior = π()
    FittedConvexMinimumDistance(estimated_prior, method, model)
end


text:
"""
    FieldModule

Module for abstract fields.
"""
module FieldModule

__precompile__(true)

using ..FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import Base.copyto!

"""
	AbstractField

Abstract field.

Expected  attributes:
  + `values::FMat{T}`: Array of degree of freedom parameters,  indexed by entity number
  + `dofnums::FIntMat`: Array of degree of freedom numbers, indexed by entity number
  + `is_fixed::Matrix{Bool}`: Array of Boolean flags, indexed by entity number
  + `fixed_values::FMat{T}`: Array of fixed values, indexed by entity number
  + `nfreedofs::FInt`: Total number of free degrees of freedom

See also: [`@add_Field_fields()`](@ref) .
"""
abstract type AbstractField end

"""
    add_Field_fields()

Generate the attributes (i. e. fields) of a `Field`. The methods defined for
the abstract type depend on these attributes to be present.
"""
macro add_Field_fields()
    return esc(:(
    values::Array{T, 2};
    dofnums::Array{FInt, 2};
    is_fixed::Matrix{Bool};
    fixed_values::Array{T, 2};;
    nfreedofs::FInt;
    )
    )
end

"""
    ndofs{F<:AbstractField}(self::T)

Dimension of the degree of freedom parameters (i. e. how many degrees of
freedom per entity).
"""
ndofs(self::AbstractField)  = size(self.values, 2)

"""
    nents{F<:AbstractField}(self::F)

Number of nodes associated with the field.
"""
nents(self::AbstractField)  = size(self.values, 1)

"""
    copyto!(DEST::F,  SRC::F) where {F<:AbstractField}

Copy data from one field to another.
"""
function copyto!(DEST::F,  SRC::F) where {F<:AbstractField}
    copyto!(DEST.values, SRC.values)
    copyto!(DEST.dofnums, SRC.dofnums)
    copyto!(DEST.is_fixed, SRC.is_fixed)
    copyto!(DEST.fixed_values, SRC.fixed_values)
    DEST.nfreedofs = SRC.nfreedofs
    return  DEST
end

"""
    wipe!(self::AbstractField)

Wipe all the data from the field.

This includes values, prescribed values, degree of freedom numbers, and "is
fixed" flags. The number of free degrees of freedom is set to zero.
"""
function wipe!(self::AbstractField)
    Zer = zero(eltype(self.fixed_values[1]))
    self.nfreedofs = 0
    fill!(self.dofnums, 0)
    fill!(self.is_fixed, false)
    fill!(self.values, Zer)
    fill!(self.fixed_values, Zer)
    return self
end

"""
    gathersysvec{F<:AbstractField}(self::F)

Gather values from the field for the whole system vector.
Return a vector.
"""
function gathersysvec(self::F) where {F<:AbstractField}
    nents,dim = size(self.values)
    vec = zeros(typeof(self.values[1,1]),self.nfreedofs)
    for i = 1:nents
        for j = 1:dim
            en = self.dofnums[i,j]
            if (en > 0) && (en <= self.nfreedofs)
                vec[en] = self.values[i,j]
            end
        end
    end
    return vec
end

"""
    gathersysvec!(self::F, vec::FVec{T}) where {F<:AbstractField, T}

Gather values from the field for the whole system vector.
"""
function gathersysvec!(self::F, vec::FVec{T}) where {F<:AbstractField, T}
    nents,dim = size(self.values)
    @assert length(vec) == self.nfreedofs
    for i = 1:nents
        for j = 1:dim
            en = self.dofnums[i,j]
            if (en > 0) && (en <= self.nfreedofs)
                vec[en] = self.values[i,j]
            end
        end
    end
    return vec
end

"""
    gathervalues_asvec!(self::AbstractField, dest::AbstractArray{T, 1},
        conn::CC) where {CC, T}

Gather values from the field into a vector.

The order is: for each node  in the connectivity, copy into the buffer all the
degrees of freedom,  then the next node and so on.

`dest` = destination buffer: overwritten  inside,  must be preallocated
in the correct size
"""
function gathervalues_asvec!(self::AbstractField, dest::AbstractArray{T, 1}, conn::CC) where {CC, T}
    # The order of the loops matters, first i, then j
    en::FInt = 1;
    for i = 1:length(conn)
        for j = 1:size(self.values,2)
            dest[en] = self.values[conn[i],j];
            en = en + 1;
        end
    end
    return dest
end

"""
    gathervalues_asmat!(self::AbstractField, dest::AbstractArray{T, 2},
        conn::CC) where {CC, T}

Gather values from the field into a two-dimensional array.

The order is: for each node  in the connectivity, copy into the corresponding
row of the buffer all the degrees of freedom,  then the next node into the next
row and so on.

`dest` = destination buffer: overwritten  inside,  must be preallocated
in the correct size
"""
function gathervalues_asmat!(self::AbstractField, dest::AbstractArray{T, 2},    conn::CC) where {CC, T}
    @inbounds for j in 1:size(self.values,2)
        @inbounds for i in 1:length(conn)
            dest[i, j] = self.values[conn[i], j];
        end
    end
    return dest
end

"""
    gatherfixedvalues_asvec!(self::AbstractField, dest::AbstractArray{T, 1},
        conn::CC) where {CC, T}

Gather FIXED values from the field into a vector.

The order is: for each node  in the connectivity, copy into the buffer all the
fixed degrees of freedom,  then the next node and so on. If a degree of freedom
is NOT fixed, the corresponding entry is  set to zero.

`dest` = destination buffer: overwritten  inside,  must be preallocated
in the correct size
"""
function gatherfixedvalues_asvec!(self::AbstractField, dest::AbstractArray{T, 1},    conn::CC) where {CC, T}
    # The order of the loops matters here! It must be i, j
    en::FInt = 1;
    for i = 1:length(conn)
        for j = 1:size(self.fixed_values,2)
            if self.is_fixed[conn[i],j] # free degree of freedom
                dest[en] = self.fixed_values[conn[i], j];
            else
                dest[en] = 0.0
            end
            en = en + 1;
        end
    end
    return dest
end

"""
    gatherfixedvalues_asmat!(self::AbstractField, dest::AbstractArray{T, 2},
        conn::CC) where {CC, T}

Gather FIXED values from the field into a two-dimensional array.

The order is: for each node  in the connectivity, copy into the corresponding
row of the buffer all the degrees of freedom,  then the next node into the next
row and so on.  If a degree of freedom
is NOT fixed, the corresponding entry is  set to zero.

`dest` = destination buffer: overwritten  inside,  must be preallocated
in the correct size
"""
function gatherfixedvalues_asmat!(self::AbstractField, dest::AbstractArray{T, 2},    conn::CC) where {CC, T}
    for j = 1:size(self.fixed_values,2)
        for i = 1:length(conn)
            if self.is_fixed[conn[i],j] # fixed degree of freedom
                dest[i, j] = self.fixed_values[conn[i], j];
            else
                dest[i, j] = 0.0
            end
        end
    end
    return dest
end

function anyfixedvaluenz(self::AbstractField, conn::CC) where {CC}
    for i = 1:length(conn)
        for j = 1:size(self.fixed_values,2)
            if self.is_fixed[conn[i],j] # free degree of freedom
                if  abs(self.fixed_values[conn[i], j]) > 0.0
                    return true
                end
            end
        end
    end
    return false
end

"""
    gatherdofnums!(self::AbstractField, dest::A, conn::CC) where {A, CC}

Gather dofnums from the field.

The order is: for each node  in the connectivity, copy into the buffer all the degrees of
freedom for that node,  then the next node  and so on.
"""
function gatherdofnums!(self::AbstractField, dest::A, conn::CC) where {A, CC}
    en::FInt = 1;
    for i = 1:length(conn)
        for j = 1:size(self.dofnums,2)
            dest[en] = self.dofnums[conn[i],j];
            en = en+1;
        end
    end
    return dest
end

"""
    numberdofs!(self::AbstractField)

Number the degrees of freedom.

The free components in the field are numbered consecutively. No effort is
made to optimize the numbering in any way. If you'd like to optimize the
numbering of the degrees of freedom, use a form that sets the
permutation of the degrees of freedom, or the permutation of the nodes.
"""
function numberdofs!(self::AbstractField)
    return numberdofs!(self, 1:nents(self))
end

"""
    numberdofs!(self::AbstractField, entperm)

Number the degrees of freedom.

The free components in the field are numbered consecutively. The sequence of the
entities is given by the `entperm` permutation (array or range).
"""
function numberdofs!(self::AbstractField, entperm)
    fixed_dofnum::FInt = 0
    nents,dim = size(self.values)
    self.nfreedofs::FInt =0
    for i in entperm
        for j=1:dim
            if !self.is_fixed[i,j] # free degree of freedom
                self.nfreedofs = self.nfreedofs + 1
                self.dofnums[i,j] = self.nfreedofs
            else # fixed degree of freedom: no equation
                self.dofnums[i,j] = fixed_dofnum
            end
        end
    end
    return  self
end

function _setebc!(self::AbstractField, fenid::FInt, is_fixed::Bool, comp::FInt, val::T) where {T<:Number}
    self.is_fixed[fenid,comp] = is_fixed;
    if self.is_fixed[fenid,comp]
    	self.fixed_values[fenid,comp] = val;
    else
    	self.fixed_values[fenid,comp] = zero(T)
    end
    return  self
end

"""
    setebc!(self::AbstractField, fenid::FInt, is_fixed::Bool, comp::FInt, val::T) where {T<:Number}

Set the EBCs (essential boundary conditions).

`fenids`         - array of N node identifiers
`is_fixed` = scaler Boolean: are the degrees of freedom being fixed (true)
             or released (false),
`comp` = integer, which  degree of freedom (component),
`val` = array of N values of type T

Note:  Any call to `setebc!()` potentially changes the current assignment
which degrees of freedom are free and which are fixed and therefore is
presumed to invalidate the current degree-of-freedom numbering. In such a case
this method sets `nfreedofs = 0`; and  `dofnums=0`.
"""
function setebc!(self::AbstractField, fenid::FInt, is_fixed::Bool, comp::FInt, val::T) where {T<:Number}
    @assert 1 <= comp <= size(self.values,2) "Requested  nonexistent  degree of freedom"
    @assert 1 <= fenid <= size(self.values,1) "Requested nonexistent node"
    _setebc!(self, fenid, is_fixed, comp, val)
    self.nfreedofs = 0
    fill!(self.dofnums, 0)
    return  self
end

"""
    setebc!(self::AbstractField, fenids::FIntVec, is_fixed::Bool, comp::FInt,
      val::FVec{T}) where {T<:Number}

Set the EBCs (essential boundary conditions).

`fenids`         - array of N node identifiers
`is_fixed` = scaler Boolean: are the degrees of freedom being fixed (true)
             or released (false),
`comp` = integer, which  degree of freedom (component),
`val` = array of N values of type T

Note:  Any call to `setebc!()` potentially changes the current assignment
which degrees of freedom are free and which are fixed and therefore is
presumed to invalidate the current degree-of-freedom numbering. In such a case
this method sets `nfreedofs = 0`; and  `dofnums=0`.
"""
function setebc!(self::AbstractField, fenids::FIntVec, is_fixed::Bool, comp::FInt, val::FVec{T}) where {T<:Number}
    @assert comp <= size(self.values,2) "Requested  nonexistent  degree of freedom"
    @assert maximum(fenids) <= size(self.values,1) "Requested nonexistent node"
    @assert size(fenids) == size(val) "Arrays of mismatched sizes"
    for  j = 1:length(fenids)
        _setebc!(self, fenids[j], is_fixed, comp, val[j])
    end
    self.nfreedofs = 0
    fill!(self.dofnums, 0)
    return  self
end

"""
    setebc!(self::AbstractField, fenids::FIntVec, is_fixed::Bool, comp::FInt, val::T = 0.0) where {T<:Number}

Set the EBCs (essential boundary conditions).

`fenids`         - array of N node identifiers
`is_fixed` = scaler Boolean: are the degrees of freedom being fixed (true)
             or released (false),
`comp` = integer, which  degree of freedom (component),
`val` = scalar of type T

Note:  Any call to `setebc!()` potentially changes the current assignment
which degrees of freedom are free and which are fixed and therefore is
presumed to invalidate the current degree-of-freedom numbering. In such a case
this method sets `nfreedofs = 0`; and  `dofnums=0`.
"""
function setebc!(self::AbstractField, fenids::FIntVec, is_fixed::Bool, comp::FInt, val::T = 0.0) where {T<:Number}
    @assert (comp >= 1 && comp <= size(self.values,2)) "Requested  nonexistent  degree of freedom"
    @assert maximum(fenids) <= size(self.values,1) "Requested nonexistent node"
    @assert minimum(fenids) >= 1 "Requested nonexistent node"
    for  j = 1:length(fenids)
        _setebc!(self, fenids[j], is_fixed, comp, val)
    end
    self.nfreedofs = 0
    fill!(self.dofnums, 0)
    return  self
end

"""
    setebc!(self::AbstractField, fenids::FIntVec, comp::FInt,
      val::FVec{T}) where {T<:Number}

Set the EBCs (essential boundary conditions).

`fenids` = array of N node identifiers
`comp` = integer, which  degree of freedom (component),
`val` = array of N values of type `T`

Note:  Any call to `setebc!()` potentially changes the current assignment
which degrees of freedom are free and which are fixed and therefore is
presumed to invalidate the current degree-of-freedom numbering. In such a case
this method sets `nfreedofs = 0`; and  `dofnums=0`.
"""
function setebc!(self::AbstractField, fenids::FIntVec, comp::FInt, val::FVec{T}) where {T<:Number}
    return setebc!(self, fenids, true, comp, val)
end

"""
    setebc!(self::AbstractField, fenids::FIntVec, comp::FInt, val::T=0.0) where {T<:Number}

Set the EBCs (essential boundary conditions).

`fenids` = array of N node identifiers
`comp` = integer, which  degree of freedom (component),
`val` = scalar of type `T`

Note:  Any call to `setebc!()` potentially changes the current assignment
which degrees of freedom are free and which are fixed and therefore is
presumed to invalidate the current degree-of-freedom numbering. In such a case
this method sets `nfreedofs = 0`; and  `dofnums=0`.
"""
function setebc!(self::AbstractField, fenids::FIntVec, comp::FInt, val::T=0.0) where {T<:Number}
    return setebc!(self, fenids, true, comp, val)
end

"""
    setebc!(self::AbstractField, fenids::FIntVec, is_fixed::Bool, comp::FIntVec, val::T=0.0) where {T<:Number}

Set the EBCs (essential boundary conditions).

`fenids` = array of N node identifiers
`comp` = integer vector, which degree of freedom (component),
`val` = scalar of type `T`, default is `0.0`

Note:  Any call to `setebc!()` potentially changes the current assignment which
degrees of freedom are free and which are fixed and therefore is presumed to
invalidate the current degree-of-freedom numbering. In such a case this method
sets `nfreedofs = 0`; and  `dofnums=0`.
"""
function setebc!(self::AbstractField, fenids::FIntVec, is_fixed::Bool, comp::FIntVec, val::T=0.0) where {T<:Number}
	for j in comp
		setebc!(self, fenids, is_fixed, j, val)
	end
    return self
end

"""
    setebc!(self::AbstractField, fenids::FIntVec)

Set the EBCs (essential boundary conditions).

Suppress all degrees of freedom at the given nodes.

`fenids`         - array of N node identifiers

Note:  Any call to `setebc!()` potentially changes the current assignment
which degrees of freedom are free and which are fixed and therefore is
presumed to invalidate the current degree-of-freedom numbering. In such a case
this method sets `nfreedofs = 0`; and  `dofnums=0`.
"""
function setebc!(self::AbstractField, fenids::FIntVec)
    zer = zero(eltype(self.fixed_values[1]))
    for comp in 1:size(self.values, 2)
        setebc!(self, fenids, true, comp, zer)
    end
    return self
end

"""
    setebc!(self::AbstractField, fenid::FInt)

Set the EBCs (essential boundary conditions).

Suppress all degrees of freedom at the given node.

`fenid`         - One integer as a node identifier

Note:  Any call to setebc!() potentially changes the current assignment
which degrees of freedom are free and which are fixed
and therefore is presumed to invalidate the
current degree-of-freedom numbering. In such a case this method sets
`nfreedofs = 0`; and  `dofnums=0`.
"""
function setebc!(self::AbstractField, fenid::FInt)
    return setebc!(self, [fenid])
end

"""
    setebc!(self::AbstractField)

Set the EBCs (essential boundary conditions).

All essential boundary conditions are CLEARED.

Note:  Any call to setebc!() potentially changes the current assignment
which degrees of freedom are free and which are fixed
and therefore is presumed to invalidate the
current degree-of-freedom numbering. In such a case this method sets
`nfreedofs = 0`; and  `dofnums=0`.
"""
function setebc!(self::AbstractField)
    self.nfreedofs = 0
    fill!(self.dofnums, 0)
    fill!(self.is_fixed, false)
    fill!(self.fixed_values, zero(eltype(self.fixed_values[1])))
    return  self
end

"""
    applyebc!(self::AbstractField)

Apply EBCs (essential boundary conditions).
"""
function applyebc!(self::AbstractField)
    nents,dim = size(self.values);
    for i = 1:nents
        for j = 1:dim
            if self.is_fixed[i,j]
                self.values[i,j] = self.fixed_values[i,j];
            end
        end
    end
    return  self
end

"""
    scattersysvec!(self::AbstractField, vec::FVec{T}) where {T<:Number}

Scatter values to the field from a system vector.
"""
function scattersysvec!(self::AbstractField, vec::FVec{T}) where {T<:Number}
    nents,dim = size(self.values);
    for i = 1:nents
        for j = 1:dim
            dn = self.dofnums[i,j];
            if (dn > 0) && (dn <= self.nfreedofs)
                self.values[i,j] = vec[dn];
            end
        end
    end
    return  self
end

"""
    incrscattersysvec!(self::AbstractField, vec::FVec{T}) where {T<:Number}

Increment values of the field by scattering a system vector.
"""
function incrscattersysvec!(self::AbstractField, vec::FVec{T}) where {T<:Number}
    nents,dim = size(self.values);
    for i = 1:nents
        for j = 1:dim
            dn = self.dofnums[i,j];
            if (dn > 0) && (dn <= self.nfreedofs)
                self.values[i,j] += vec[dn];
            end
        end
    end
    return  self
end

"""
    prescribeddofs(uebc, u)

Find which degrees of freedom are prescribed.
`uebc` = field which defines the constraints (is the dof fixed and to which value),
`u` = field which does not have the constraints applied, and serves as the source of equation numbers,
`uebc` and `u` may be one and the same field.

"""
function prescribeddofs(uebc::AbstractField, u::AbstractField)
    dofnums = FInt[]
	prescribedvalues = eltype(uebc.values)[]
	nents, dim = size(uebc.values)
	@assert size(uebc.values) == size(u.values)
	for i in 1:nents
		for j in 1:dim
			if uebc.is_fixed[i,j]
			    dn = u.dofnums[i,j];
                push!(prescribedvalues, uebc.fixed_values[i,j]);
			    push!(dofnums, dn)
			end
		end
	end
	return dofnums, prescribedvalues
end

end


text:
function get_page(opts::Dict=Dict())
    Juno.isactive() ? Juno.Atom.blinkplot() : Window(opts).content
end

Juno.render(::Juno.PlotPane, n::Union{Node, Scope, AbstractWidget}) =
    (body!(get_page(), n); nothing)

Juno.render(i::Juno.Editor, n::Union{Node, Scope, AbstractWidget}) =
    Juno.render(i, Text("$(n.instanceof) Node with $(n._descendants_count) descendent(s)"))

function WebIO.register_renderable(T::Type, ::Val{:atom})
    Juno.media(T, Juno.Media.Graphical)
    eval(:(Juno.Media.render(::Juno.PlotPane, x::$T) =
          (body!(get_page(), WebIO.render(x)); nothing)))
    if hasmethod(WebIO.render_inline, (T,))
        eval(:(Juno.render(i::Juno.Editor, x::$T) =
               Juno.render(i, WebIO.render_inline(x))))
    end
end

function WebIO.setup_provider(::Val{:atom})
    Juno.media(Node, Juno.Media.Graphical)
    Juno.media(Scope, Juno.Media.Graphical)
end
WebIO.setup(:atom)


text:
using MPIMeasurements

scanner = MPIScanner("MPS.toml")
daq = getDAQ(scanner)

params = toDict(daq.params)

params["studyName"]="TestTobi"
params["studyDescription"]="A very cool measurement"
params["scannerOperator"]="Tobi"
params["dfStrength"]=[20e-3]
params["acqNumAverages"]=1000

measurementCont(daq, params, controlPhase=true)


text:
# This file is part of Kpax3. License is MIT.

abstract type PriorRowPartition end

"""
# Ewens-Pitman distribution

## Description

The Ewens-Pitman distribution is a discrete probability distribution on the
partitions of the set N = {1, ..., n} (n = 1, 2, ...). It is a two parameters
generalization of the Ewens sampling formula. The latter is also known as the
Chinese Restaurant Process or as the marginal distribution of a Dirichlet
Process.

## Fields

* `α` Real number (see details)
* `θ` Real number (see details)

## Details

The two parameters must satisfy either
* α < 0 and θ = -Lα for some L ∈ {1, 2, ...}
* 0 ≤ α < 1 and θ > − α

The special case α = 0 and θ > 0 corresponds to the Ewens sampling formula.

## References

Aldous, D. J. (1985) Exchangeability and related topics. In *École d'Été de
Probabilités de Saint-Flour XIII — 1983*. Lecture Notes in Mathematics
**1117**, 1-198. Springer Berlin Heidelberg.
<http://dx.doi.org/10.1007/BFb0099421>.

Gnedin, A. and Pitman, J. (2006) Exchangeable Gibbs partitions and Stirling
triangles. *Journal of Mathematical Sciences* **138**(3), 5674-5685.
<http://dx.doi.org/10.1007/s10958-006-0335-z>.

Kerov, S. (2006) Coherent random allocations, and the Ewens-Pitman formula.
*Journal of Mathematical Sciences*, **138**(3).
<http://dx.doi.org/10.1007/s10958-006-0338-9>.

Pitman, J. (1995) Exchangeable and Partially Exchangeable Random Partitions.
*Probability Theory and Related Fields* **102**(2), 145-158.
<http://dx.doi.org/10.1007%2FBF01213386>.

Pitman, J. (2006) Combinatorial Stochastic Processes. In *Ecole d’Eté de
Probabilités de Saint-Flour XXXII – 2002*. Lecture Notes in Mathematics
**1875**. Springer Berlin Heidelberg. <http://dx.doi.org/10.1007/b11601500>.
"""
abstract type EwensPitman <: PriorRowPartition end

"""
# Ewens-Pitman distribution

## Description

Ewens-Pitman distribution with 0 < α < 1,  θ > -α and θ ≠ 0.

## Fields

* `α` Real number greater than zero and lesser than one
* `θ` Real number greater than `-α` but different from zero
"""
struct EwensPitmanPAUT <: EwensPitman
  α::Float64
  θ::Float64
end

"""
# Ewens-Pitman distribution

## Description

Ewens-Pitman distribution with 0 < α < 1 and θ = 0.

## Fields

* `α` Real number greater than zero and lesser than one
"""
struct EwensPitmanPAZT <: EwensPitman
  α::Float64
end

"""
# Ewens-Pitman distribution

## Description

Ewens-Pitman distribution with α = 0 and θ > 0. This is equivalent to the Ewens
sampling formula.

## Fields

* `θ` Real number greater than zero
"""
struct EwensPitmanZAPT <: EwensPitman
  θ::Float64
end

"""
# Ewens-Pitman distribution

## Description

Ewens-Pitman distribution with α < 0 and θ > 0. θ = -Lα for some
L ∈ {1, 2, ...}.

## Fields

* `α` Real number lesser than zero
* `L` Integer number greater than zero
"""
struct EwensPitmanNAPT <: EwensPitman
  α::Float64
  L::Int
end

# TODO: Julia v0.4
# we can't document this constructor: the type has already been documented

#=
# Constructor of an object of (super)type EwensPitman

## Description

Create an appropriate EwensPitman type according to the arguments' values.

## Usage

EwensPitman(α, θ)

## Arguments

* `α` Real number
* `θ` Real number
=#
function EwensPitman(α::Real,
                     θ::Float64)
  α = float(α)

  if α == zero(Float64)
    if θ > zero(Float64)
      EwensPitmanZAPT(θ)
    else
      throw(KDomainError(string("When argument α is zero, argument θ must be ",
                                "positive.")))
    end
  elseif zero(Float64) < α < one(Float64)
    if θ == zero(Float64)
      EwensPitmanPAZT(α)
    elseif θ > -α
      EwensPitmanPAUT(α, θ)
    else
      throw(KDomainError(string("When argument α is non-negative, argument θ ",
                                "must be greater than -α.")))
    end
  elseif α < zero(Float64)
    throw(KDomainError(string("When argument α is negative, provide an ",
                              "integer L to define parameter θ.")))
  else
    throw(KDomainError("Argument α cannot be greater than or equal to 1."))
  end
end

function EwensPitman(α::Real,
                     L::Int)
  α = float(α)

  if α < zero(Float64)
    if L > zero(Int)
      EwensPitmanNAPT(α, L)
    else
      throw(KDomainError("Argument L must be positive."))
    end
  else
    throw(KDomainError(string("When argument α is non-negative, provide a ",
                              "real value for parameter θ.")))
  end
end


text:
function indexfilter(dbdata,
                     filter_query;
                     id_key=DEFAULT_DB_ID_KEY,
                     sort_keys=nothing,
                     sort_reverse=false,
                     exclude=nothing)::Vector{Int}
    # Checks
    cols = colnames(dbdata)
    if !in(id_key, cols)
        @warn "Could not find id_key='$id_key'. Returning empty id vector..."
        return Int[]
    end

    # Function with methods for constructing data filtering functions (x is the current value for a column)
    __filter_function(val) = x -> x == val              # filter by equality to a value

    __filter_function(val::Tuple) = x -> x in val       # filter by being present in a set

    __filter_function(val::NTuple{N,T}) where {N, T<:AbstractString} =
        x -> any(v -> occursin(v, x), val)              # filter by having at least one value from a string set

    __filter_function(val::AbstractVector) =
        x -> x >= val[1] && x <= val[2]                 # filter by belonging to an interval

    __safe_filter_function(key, val::AbstractVector) = try
        @assert length(val) >= 2
        __filter_function(val)
    catch
        @warn "Failed to generate filter function for $key:$val"
        x -> true  # pass-through
    end

    __safe_filter_function(key, val) = __filter_function(val)

    # Find sorting, filtering columns i.e. keys
    safe_filter_keys = ()
    if !isempty(filter_query)
        safe_filter_keys = Tuple(key for key in keys(filter_query) if key in cols)
    end

    safe_sort_keys = ()
    if sort_keys !== nothing
        safe_sort_keys = Tuple(key for key in sort_keys if key in cols)
    end

    # Sort
    if !isempty(safe_sort_keys)
        sort_select_keys = safe_filter_keys
        !in(id_key, sort_select_keys) && (sort_select_keys= (id_key, sort_select_keys...))
        dbdata = sort(dbdata, safe_sort_keys; select=sort_select_keys, rev=sort_reverse)
    end

    # Filter
    if !isempty(safe_filter_keys)
        selectors = Tuple(key => __safe_filter_function(key, filter_query[key])
                          for key in safe_filter_keys)
        dbdata = filter(all, dbdata; select=selectors)
    end

    # Return ids (use row iterator as it contains all fields)
    return setdiff(rows(dbdata, id_key), [exclude])
end

sort(t::NDSparse, by...; select=colnames(t), kwargs...) =
    sort(table(nds), by...; select=select, kwargs...)


text:
#*******************************************************/
#* Copyright(c) 2018 by Artelys                        */
#* This source code is subject to the terms of the     */
#* MIT Expat License (see LICENSE.md)                  */
#*******************************************************/

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#  This example demonstrates how to use Knitro to solve the following
#  simple quadratic programming problem (QP).
#
#     minimize     0.5*(x0^2+x1^2+x2^2) + 11*x0 + x2
#     subject to   -6*x2  <= 5
#                 0 <= x0
#                 0 <= x1
#                -3 <= x2 <= 2
#
#  The optimal solution is:
#     obj=-0.4861   x=[0,0,-5/6]
#
#  The purpose is to illustrate how to invoke Knitro using the C
#  language API.
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

using KNITRO
using Test

function example_qp1(; verbose=true)
    # Used to specify whether linear and quadratic objective
    # terms are loaded separately or together in this example.
    bSeparate = false

    # Create a new Knitro solver instance.
    kc = KNITRO.KN_new()

    # Illustrate how to override default options by reading from
    # the knitro.opt file.
    options = joinpath(dirname(@__FILE__), "..", "examples", "knitro.opt")
    KNITRO.KN_load_param_file(kc, options)
    kn_outlev = verbose ? KNITRO.KN_OUTLEV_ITER : KNITRO.KN_OUTLEV_NONE
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_OUTLEV, kn_outlev)

    # Initialize Knitro with the problem definition.

    # Add the variables and set their bounds.
    # Note: unset bounds assumed to be infinite.
    KNITRO.KN_add_vars(kc, 3)
    KNITRO.KN_set_var_lobnds(kc, [0.0, 0.0, -3.0])
    KNITRO.KN_set_var_upbnds(kc, 2, 2.0)

    # Add the constraint and set the bound and coefficient.
    KNITRO.KN_add_cons(kc, 1)
    KNITRO.KN_set_con_upbnd(kc, 0, 5.0)
    KNITRO.KN_add_con_linear_struct(kc, 0, 2, -6.0)

    # Set the coefficients for the objective -
    # can either set linear and quadratic objective structure
    # separately or all at once.  We show both cases below.
    # Change the value of "bSeparate" to try both cases.

    if (bSeparate)
        # Set linear and quadratic objective structure separately.
        # First set linear objective structure.
        lobjIndexVars = Int32[0, 2]
        lobjCoefs = [11.0, 1.0]
        KNITRO.KN_add_obj_linear_struct(kc, lobjIndexVars, lobjCoefs)
        # Now set quadratic objective structure.
        qobjIndexVars1 = Int32[0, 1, 2]
        qobjIndexVars2 = Int32[0, 1, 2]
        qobjCoefs = [0.5, 0.5, 0.5]
        KNITRO.KN_add_obj_quadratic_struct(kc, qobjIndexVars1, qobjIndexVars2, qobjCoefs)
    else
        # Example of how to set linear and quadratic objective
        # structure at once. Setting the 2nd variable index in a
        # quadratic term to be negative, treats it as a linear term.
        indexVars1 = Int32[0, 1, 2, 0, 2]
        indexVars2 = Int32[0, 1, 2, -1, -1]  # -1 for linear coefficients
        objCoefs = [0.5, 0.5, 0.5, 11.0, 1.0]
        KNITRO.KN_add_obj_quadratic_struct(kc, indexVars1, indexVars2, objCoefs)
    end

    # Set minimize or maximize (if not set, assumed minimize)
    KNITRO.KN_set_obj_goal(kc, KNITRO.KN_OBJGOAL_MINIMIZE)

    # Enable iteration output and crossover procedure to try to
    # get more solution precision
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_BAR_MAXCROSSIT, 5)

    # Solve the problem.
    #
    # Return status codes are defined in "kn_defines.jl" and described
    # in the Knitro manual.
    nStatus = KNITRO.KN_solve(kc)
    nStatus, objSol, x, lambda_ =  KNITRO.KN_get_solution(kc)

    # An example of obtaining solution information.
    if verbose
        println("Knitro converged with final status = ", nStatus)
        println("  optimal objective value  = ", objSol)
        println("  optimal primal values x  = ",  x)
        println("  feasibility violation    = ", KNITRO.KN_get_abs_feas_error(kc))
        println("  KKT optimality violation = ", KNITRO.KN_get_abs_opt_error(kc))
    end
    # Delete the Knitro solver instance.
    KNITRO.KN_free(kc)

    @testset "Example QP1" begin
        @test nStatus == 0
        @test objSol ≈ -0.4861 atol=1e-4
        @test x ≈ [0., 0., -5/6] atol=1e-5
    end
end

example_qp1(; verbose=isdefined(Main, :KN_VERBOSE) ? KN_VERBOSE : true)


text:
include("BLTandPantelides/BLTandPantelidesUtilities.jl")
include("BLTandPantelides/BLTandPantelides.jl")

include("Utilities.jl")

include("ExactlyRemoveSingularities.jl")
include("StateSelection.jl")

include("DAEquations/SymbolicTransform.jl")
include("DAEquations/BasicStructuralTransform.jl")

include("StructuralTransform.jl")

include("ModelElaboration.jl")


text:
using Manifolds, ManifoldsBase, Documenter

makedocs(
    # for development, we disable prettyurls
    format = Documenter.HTML(prettyurls = false),
    modules = [Manifolds, ManifoldsBase],
    sitename = "Manifolds",
    pages = [
        "Home" => "index.md",
        "Manifolds" => [
            "Basic manifolds" => [
                "Euclidean" => "manifolds/euclidean.md",
                "Cholesky Space" => "manifolds/choleskyspace.md",
                "Rotations" => "manifolds/rotations.md",
                "Sphere" => "manifolds/sphere.md",
                "Symmetric Positive Definite" => "manifolds/symmetricpositivedefinite.md"
            ],
            "Combined manifolds" => [
                "Power manifold" => "manifolds/power.md",
                "Product manifold" => "manifolds/product.md",
                "Vector bundle" => "manifolds/vector_bundle.md"
            ],
            "Manifold decorators" => [
                "Array manifold" => "manifolds/array.md",
                "Metric manifold" => "manifolds/metric.md"
            ]
        ],
        "Distributions" => "distributions.md",
        "Library" => [
            "Public" => "lib/public.md",
            "Internals" => "lib/internals.md"
        ]
    ]
)

deploydocs(
    repo = "github.com/JuliaNLSolvers/Manifolds.jl.git",
)


text:
# Functions for the generation of a dipole field.
#
# Modified from the original Python version by
# Los Alamos National Security, LLC. 2010

# mean value of B at the magnetic equator on the Earth's surface, [nT]
const B₀ = 31200.

"""
	 b_mag(x, y)

Return the strength of the dipole magnetic field in nT for a position
`x`, `y`in units of planetary radius.
"""
function b_mag(x, y)
   r = @. √(x^2 + y^2)
   cosy = y./r
   B = @. B₀ * sqrt(1+3*cosy^2)/r^3
end

"""
	 b_mag(x, y, z)

Return the strength of the dipole magnetic field in nT for a position
`x`, `y`, `z` in units of planetary radius.
"""
function b_mag(x, y, z)
   r = @. √(x^2 + y^2 + z^2)
   cosθ = z ./ r
   B = @. B₀ * sqrt(1+3*cosθ^2)/r^3
end

"""
    b_hat(x, y)

Return two arrays, `bx` and `by`, corresponding to the direction of a dipole
field at (`x`, `y`). The grid is organized in meshgrid (default) or ndgrid
format.
"""
function b_hat(x, y; gridType="meshgrid")
   if gridType == "meshgrid"
      xgrid = [i for _ in y, i in x]
      ygrid = [j for j in y, _ in x]
   else
      xgrid = [i for i in x, _ in y]
      ygrid = [j for _ in x, j in y]
   end

   r = @. sqrt(xgrid^2 + ygrid^2)
   cosy = ygrid./r
   sinx = xgrid./r

   denom = @. sqrt(1.0 + 3.0*cosy^2)

   br = @. 2.0 * cosy / denom
   bθ =          sinx ./ denom

   bx = @. br*sinx + bθ*cosy
   by = @. br*cosy - bθ*sinx

   return bx, by
end

"""
    b_hat(x, y, z)

Return the direction of a 3D dipole field at given locations.
The grid is organized in meshgrid or ndgrid (default) format.
"""
function b_hat(x, y, z; gridType="ndgrid")
   if gridType == "meshgrid"
      xgrid = [i for _ in y, i in x, _ in z]
      ygrid = [j for j in y, _ in x, _ in z]
      zgrid = [k for _ in y, _ in x, k in z]
   else
      xgrid = [i for i in x, _ in y, _ in z]
      ygrid = [j for _ in x, j in y, _ in z]
      zgrid = [k for _ in x, _ in y, k in z]
   end

   r = @. sqrt(xgrid^2 + ygrid^2 + zgrid^2)

   bx = @. 3*xgrid*zgrid/r^5
   by = @. 3*ygrid*zgrid/r^5
   bz = @. (3*zgrid^2 - r^2) / r^5

   return bx, by, bz
end

"""
	 b_line(x, y; npoint=30)

For a starting X, Y point return x and y vectors that trace the dipole field
line that passes through the given point.
"""
function b_line(x, y; npoints=30)
   r = sqrt(x^2 + y^2)
   try
      theta = atan(x/y)
   catch
      @warn "ZeroDivisionError"
      theta = pi/2.0
   end

   R = r/(sin(theta)^2)

   if x < 0
      theta = π:π/npoints:2.0π
   else
      theta = 0:π/npoints:π
   end

   r_vec = @. R * sin(theta)^2
   x_out = @. r_vec * sin(theta)
   y_out = @. r_vec * cos(theta)

   return x_out, y_out
end

"A quick test of the dipole field functions."
function test_dipole_2d()

   x = -100.0:5.0:101.0
   y = -100.0:5.0:101.0

   x_vec, y_vec = b_hat(x,y)

   fig = plt.figure(figsize=(10,8))
   ax1 = plt.subplot(111)

   ax1.quiver(x, y, x_vec, y_vec)

   for i in -120:10:121
      x,y = b_line(float(i), 0.0, npoints=100)
      ax1.plot(x, y, "b")
   end
   for theta in π/2.0:π/100.0:3.0π/2.0
      x = sin(theta)
      y = cos(theta)
      x, y = b_line(x, y, npoints=100)
      ax1.plot(x, y, "r")
   end
   ax1.set_xlim([-100, 100])
   ax1.set_ylim([-100, 100])
   plt.title("Unit vectors for an arbitrary dipole field")

   return true
end


text:
export Conv_cds
export MVP_cds
export extended_obj
export extended_obj_grad!
export extended_conv_op
export var_proj_optim
export var_proj_op

mutable struct Conv_cds
    diagonals
    offsets
    mask
end

(C::Conv_cds)(x::Vector{TF}) where TF = σ.(MVP_cds(C, x))
(C::Conv_cds)(w::Array{TF, 2}) where TF = Conv_cds(w, C.offsets, C.mask)

function MVP_cds(
    R::Array{TF, 2}, offset::Vector{TI}, x::Vector{TF}
) where {TF<:Real, TI<:Integer}

    return sum(R[:, i] .* circshift(x, -offset[i]) for i=1:length(offset))
end

MVP_cds(C::Conv_cds, x::Vector{TF}) where TF = MVP_cds(C.diagonals, C.offsets, x)

function extended_obj(
    x::Vector{TF}, y::Vector{TF}, C::Conv_cds, w_ex::Array{TF, 2}, w::Array{TF,2}, λ
) where {TF<:Real}

    return 5f-1*norm(y - C(w_ex)(x))^2 + 5f-1*λ^2*norm(C.mask.*(w_ex .- vec(w)'))^2
end

function extended_obj(
    x::Vector{TF}, y::Vector{TF}, C::Chain, w_ex::Array{TF, 3}, w::Array{TF,3}, λ
) where {TF<:Real}

    nlayers = length(C)
    Cloc = Chain([C[i](w_ex[:, :, i]) for i=1:nlayers]...)

    obj = 5f-1*norm(y - Cloc(x))^2
    obj += 5f-1*λ^2*sum(norm(C[i].mask.*(w_ex[:, :, i].- vec(w[:, :, i])'))^2 for i=1:nlayers)
    return obj
end

function extended_obj_grad!(grad, x, y, C, w_ex, w, λ)
    grad .= C[1].mask .* gradient(w_ex->extended_obj(x, y, C, w_ex, w, λ), w_ex)[1]
end

function var_proj_optim(x, y, C, w_ex_0, w, λ, opt_var_proj)
    # Do outer loop
    ext_obj(w_ex) = extended_obj(x, y, C, w_ex, w, λ)
    ext_obj_grad!(grad, w_ex) = extended_obj_grad!(grad, x, y, C, w_ex, w, λ)
    opt_ex = optimize(ext_obj, ext_obj_grad!, w_ex_0, LBFGS(m=100), opt_var_proj)

    # Minimizer for extended weights w_ex
    w_ex_est= Optim.minimizer(opt_ex)

    func_val = ext_obj(w_ex_est)

    return func_val, w_ex_est
end

function var_proj_op(func_val, dw, w, x, y, C, w_ex_0, λ, opt_var_proj)
    # Do outer loop
    ext_obj(w_ex) = extended_obj(x, y, C, w_ex, w, λ)
    ext_obj_grad!(grad, w_ex) = extended_obj_grad!(grad, x, y, C, w_ex, w, λ)
    opt_ex = optimize(ext_obj, ext_obj_grad!, w_ex_0, LBFGS(m=100), opt_var_proj)

    # Minimizer for extended weights w_ex
    w_ex_est= Optim.minimizer(opt_ex)

    func_val = ext_obj(w_ex_est)

    dw .= gradient(w0->extended_obj(x, y, C, w_ex_est, w0, λ), w)[1]
    func_val = ext_obj(w_ex_est)
    func_val
end

function solve_reduced_obj(proj_sol, C, w)
    w_ = zero(w)
    for j = 1:size(w, 3)
        sol = sum(C[j].mask.*proj_sol[:, :, j], dims=1)./sum(C[j].mask, dims=1)
        sol = reshape(sol, size(w)[1:2])
        w_[:, :, j] .= deepcopy(sol)
    end
    return w_
end

function extended_conv_op(nx, ny; nc=1, batchsize=1, k=3, stride=1, pad=1, nlayers=4, w=nothing)

    W = ones(Float32, k, k, 1, 1)
    bias = zeros(Float32, 1)

    G = Conv(W, bias, identity; stride=1, pad=1)
    Gadj = ConvTranspose(W, bias, identity; stride=1, pad=1)

    N = nx * ny * nc * batchsize
    K = joLinearFunctionFwd_T(
        N, N,
        x -> vec(G(reshape(x, nx, ny, nc, batchsize))),
        y -> vec(Gadj(reshape(y, nx, ny, nc, batchsize))),
        Float32,Float32, name="Conv"
    )

    # @assert islinear(K)[1]
    # @assert islinear(transpose(K))[1]
    # @assert isadjoint(K)[1]

    K = K*jo_eye(Float32, nx*ny)
    K = sparse(K)

    R0 , offset = mat2CDS(K)
    I = findall(x->x==0, R0)

    M = ones(Float32, size(R0))
    M[I] .= 0f0

    if w == nothing
        C = Chain([Conv_cds(R0.*randn(Float32, 1, k*k), offset, M) for i=1:nlayers]...)
    else
        C = Chain([Conv_cds(R0.*reshape(w[:, :, i], 1, k*k), offset, M) for i=1:nlayers]...)
    end

    return C

end


text:
#=
Copyright 2021 BlackRock, Inc.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

@testset "Settings" begin
    # should not error
    Settings(; ρ=[1], σ=[1])

    # proper consntruction
    @test_throws AssertionError Settings(; ρ=[1], σ=[1], α=3)
    @test_throws AssertionError Settings(; ρ=[-1], σ=[1])
    @test_throws AssertionError Settings(; ρ=[1], σ=[-1])
    @test_throws AssertionError Settings(; ρ=[1], σ=[1], max_iters=-1)
    @test_throws AssertionError Settings(; ρ=[1], σ=[1], obj_tol=-1)
    @test_throws AssertionError Settings(; ρ=[1], σ=[1], res_tol=-1)
    @test_throws AssertionError Settings(; ρ=[1], σ=[1], non_improvement_iters=-1)

    s1 = Settings(3, 3)
    s2 = Settings(; ρ=[1, 1, 1], σ=[1, 1, 1])

    @test s1.ρ == s2.ρ
    @test s1.σ == s2.σ
end

@testset "Params" begin
    # should not error
    P = [1.0 0.0; 0.0 2.0]
    q = [1.0, 2.0]
    A = [1.0 2.0; 0.0 1.0; 1.0 2.0]
    b = [1.0, 2.0, 3.0]
    g = repeat([zero(PiecewiseQuadratic)], 2)
    AdmmParams(P, q, A, b, g)

    @test_throws AssertionError AdmmParams([1.0 0.0 3.0; 0.0 2.0 2.0], q, A, b, g)
    @test_throws AssertionError AdmmParams(P, [1.0, 2.0, 3.0], A, b, g)
    @test_throws AssertionError AdmmParams(P, q, A, [1.0, 2.0], g)
    @test_throws AssertionError AdmmParams(P, q, A, b,
                                           repeat([zero(PiecewiseQuadratic)], 3))
end

@testset "Vars" begin
    # should not error
    Vars(5, 1)
    Vars(zeros(5), zeros(1))
    Vars(zeros(5), zeros(1), zeros(5), zeros(1), zeros(5), zeros(1))

    @test_throws AssertionError Vars(zeros(5), zeros(2), zeros(5), zeros(1), zeros(5),
                                     zeros(1))
    @test_throws AssertionError Vars(zeros(5), zeros(1), zeros(2), zeros(1), zeros(5),
                                     zeros(1))
    @test_throws AssertionError Vars(zeros(5), zeros(1), zeros(5), zeros(2), zeros(5),
                                     zeros(1))
    @test_throws AssertionError Vars(zeros(5), zeros(1), zeros(5), zeros(1), zeros(2),
                                     zeros(1))
    @test_throws AssertionError Vars(zeros(5), zeros(1), zeros(5), zeros(1), zeros(5),
                                     zeros(2))
end

@testset "Validity" begin
    P = [1.0 0.0; 0.0 2.0]
    q = [1.0, 2.0]
    A = [1.0 2.0; 0.0 1.0; 1.0 2.0]
    b = [1.0, 2.0, 3.0]
    g = repeat([zero(PiecewiseQuadratic)], 2)
    params = AdmmParams(P, q, A, b, g)

    vars = Vars(2, 3)

    settings = Settings(; ρ=[1.0, 2.0, 3.0], σ=[1.0, 2.0])

    # should not error
    assert_valid(params, settings, vars)
end


text:
@testset "List of States" begin

    s0 = GenericState(zeros(50))
    s1 = GenericState(ones(10))
    s2 = GenericState(NaN*ones(10), current_time = 1.0, current_score = 0.0)

    @test typeof(ListofStates(s0)) <: AbstractListofStates
    @test typeof(ListofStates(-1, Val{GenericState}())) <: AbstractListofStates
    @test typeof(ListofStates(1, Val{GenericState}())) <: AbstractListofStates
    #@test typeof(ListofStates(-1, 3, [])) <: AbstractListofStates
    #@test typeof(ListofStates(-1, [])) <: AbstractListofStates

    stest = ListofStates(s0, max_vector_size = 2, pnorm = Inf)

    @test state_type(stest) == GenericState{Float64,Array{Float64,1}}

    add_to_list!(stest, s1, max_vector_size = 2, pnorm = Inf)
    add_to_list!(stest, s2, max_vector_size = 2, pnorm = Inf)

    @test length(stest) == 3

    stest2 = ListofStates(s0, n = 2, max_vector_size = 2, pnorm = Inf)

    add_to_list!(stest2, s1, max_vector_size = 2, pnorm = Inf)
    add_to_list!(stest2, s2, max_vector_size = 2, pnorm = Inf)

    @test length(stest2) == 2

    df1 = print(stest, verbose = false)

    df2 = print(stest2, verbose = false)

    df3 = print(stest2, verbose = false, print_sym = [:x])

    @test typeof(df2) <: DataFrame

    stest3 = ListofStates(-1, 3, [(s0, VoidListofStates()), (s1, VoidListofStates()), (s2, VoidListofStates())])

    @test stest3[2,1] == s1

    stest4 = ListofStates(-1, [(s0, VoidListofStates()), (s1, VoidListofStates()), (s2, VoidListofStates())])

    @test length(stest4) == 3

    #nested lists

    stest5 = ListofStates(-1, [(s0, stest3)])

    df5 = print(stest5[1,2], verbose = false)

    stest7 = ListofStates(-1, [s0, s1, s2])

    @test length(stest7) == 3

end


text:
# NOTE: Most tests related to biological symbols are located in BioSymbols.jl.
@testset "Symbols" begin
    @testset "DNA" begin
        @test DNA_A === BioSymbols.DNA_A
        @test ACGT  === BioSymbols.ACGT
        @test ACGTN === BioSymbols.ACGTN
        @test typeof(DNA_A) === BioSymbols.DNA
    end

    @testset "RNA" begin
        @test RNA_A === BioSymbols.RNA_A
        @test ACGU  === BioSymbols.ACGU
        @test ACGUN === BioSymbols.ACGUN
        @test typeof(RNA_A) === BioSymbols.RNA
    end

    @testset "AminoAcid" begin
        @test AA_A === BioSymbols.AA_A
        @test typeof(AA_A) === BioSymbols.AminoAcid
    end

    @testset "Predicate functions" begin
        @test iscompatible(DNA_A, DNA_N)
        @test isambiguous(DNA_N)
        @test iscertain(DNA_A)
        @test isgap(DNA_Gap)
        @test ispurine(DNA_A)
        @test ispyrimidine(DNA_C)
        @test isGC(DNA_G)
    end

    @testset "Misc. functions" begin
        @test length(BioSymbols.alphabet(DNA)) == 16
        @test BioSymbols.gap(DNA) === DNA_Gap
        @test BioSymbols.complement(DNA_A) === DNA_T
    end

    @testset "Encoder" begin
        encode = BioSequences.encode
        EncodeError = BioSequences.EncodeError

        @testset "DNA" begin
            # 2 bits
            @test encode(DNAAlphabet{2}, DNA_A) === 0x00
            @test encode(DNAAlphabet{2}, DNA_C) === 0x01
            @test encode(DNAAlphabet{2}, DNA_G) === 0x02
            @test encode(DNAAlphabet{2}, DNA_T) === 0x03
            @test_throws EncodeError encode(DNAAlphabet{2}, DNA_M)
            @test_throws EncodeError encode(DNAAlphabet{2}, DNA_N)
            @test_throws EncodeError encode(DNAAlphabet{2}, DNA_Gap)

            # 4 bits
            for nt in BioSymbols.alphabet(DNA)
                @test encode(DNAAlphabet{4}, nt) === reinterpret(UInt8, nt)
            end
            @test_throws EncodeError encode(DNAAlphabet{4}, reinterpret(DNA, 0b10000))
        end

        @testset "RNA" begin
            # 2 bits
            @test encode(RNAAlphabet{2}, RNA_A) === 0x00
            @test encode(RNAAlphabet{2}, RNA_C) === 0x01
            @test encode(RNAAlphabet{2}, RNA_G) === 0x02
            @test encode(RNAAlphabet{2}, RNA_U) === 0x03
            @test_throws EncodeError encode(RNAAlphabet{2}, RNA_M)
            @test_throws EncodeError encode(RNAAlphabet{2}, RNA_N)
            @test_throws EncodeError encode(RNAAlphabet{2}, RNA_Gap)

            # 4 bits
            for nt in BioSymbols.alphabet(RNA)
                @test encode(RNAAlphabet{4}, nt) === reinterpret(UInt8, nt)
            end
            @test_throws EncodeError encode(RNAAlphabet{4}, reinterpret(RNA, 0b10000))
        end

        @testset "AminoAcid" begin
            @test encode(AminoAcidAlphabet, AA_A) === 0x00
            for aa in BioSymbols.alphabet(AminoAcid)
                @test encode(AminoAcidAlphabet, aa) === convert(UInt8, aa)
            end
            @test_throws BioSequences.EncodeError encode(AminoAcidAlphabet, BioSymbols.AA_INVALID)
        end
    end

    @testset "Decoder" begin
        decode = BioSequences.decode
        DecodeError = BioSequences.DecodeError

        @testset "DNA" begin
            # 2 bits
            @test decode(DNAAlphabet{2}, 0x00) === DNA_A
            @test decode(DNAAlphabet{2}, 0x01) === DNA_C
            @test decode(DNAAlphabet{2}, 0x02) === DNA_G
            @test decode(DNAAlphabet{2}, 0x03) === DNA_T
            @test_throws DecodeError decode(DNAAlphabet{2}, 0x04)
            @test_throws DecodeError decode(DNAAlphabet{2}, 0x0e)

            # 4 bits
            for x in 0b0000:0b1111
                @test decode(DNAAlphabet{4}, x) === reinterpret(DNA, x)
            end
            @test_throws DecodeError decode(DNAAlphabet{4}, 0b10000)
        end

        @testset "RNA" begin
            # 2 bits
            @test decode(RNAAlphabet{2}, 0x00) === RNA_A
            @test decode(RNAAlphabet{2}, 0x01) === RNA_C
            @test decode(RNAAlphabet{2}, 0x02) === RNA_G
            @test decode(RNAAlphabet{2}, 0x03) === RNA_U
            @test_throws DecodeError decode(RNAAlphabet{2}, 0x04)
            @test_throws DecodeError decode(RNAAlphabet{2}, 0x0e)

            # 4 bits
            for x in 0b0000:0b1111
                @test decode(RNAAlphabet{4}, x) === reinterpret(RNA, x)
            end
            @test_throws DecodeError decode(RNAAlphabet{4}, 0b10000)
        end

        @testset "AminoAcid" begin
            @test decode(AminoAcidAlphabet, 0x00) === AA_A
            for x in 0x00:0x1b
                @test decode(AminoAcidAlphabet, x) === convert(AminoAcid, x)
            end
            @test_throws BioSequences.DecodeError decode(AminoAcidAlphabet, 0x1c)
        end
    end
end


text:
function simple_volatility_given_returns(returns::Array{R,1}) where R<:Real
    sqrt(mean(returns .^ 2))
end

"""
Calculates volatility with the simple method with a specified time grid.
"""
function simple_volatility_with_grid(ts::SortedDataFrame, assets::Vector{Symbol}, time_grid)
    voldict = Dict{Symbol, eltype(ts.df[:,ts.value])}(assets .=> repeat([0.0], length(assets)))
    for a in assets
        at_times = time_grid[a]
        dd_compiled = latest_value(ts, at_times; assets = Array{Symbol}([a]))
        returns = get_returns(dd_compiled)
        vol = simple_volatility_given_returns(Array(returns[:,a]))
        av_tick_duration = (maximum(at_times) - minimum(at_times)) / (length(at_times) - 1)
        voldict[a] = vol / sqrt(av_tick_duration)
    end
    return voldict
end

"""
    simple_volatility(ts::SortedDataFrame, assets::Vector{Symbol} = get_assets(ts);
                      time_grid::Union{Missing,Dict} = missing , fixed_spacing::Union{Missing,Dict,<:Real} = missing,
                      use_all_obs::Bool = false, rough_guess_number_of_intervals::Integer = 5)

Calculates volatility with the simple method.
### Inputs
* `ts` - The tick data.
* `assets` - The assets you want to estimate volatilities for.
* `time_grid` - The grid with which to calculate returns. If missing one is generated with a fixed spacing (if that is provided) or a default spacing.
* `fixed_spacing` - A spacing used to calculate a time grid. Not used if a `time_grid` is input or if `use_all_obs = true`.
* `use_all_obs` - Use all observations to estimate volatilities. Not used if a `time_grid` is provided.
* `rough_guess_number_of_intervals` - A rough number of intervals to calculate a default spacing. Not used if a `time_grid` or `fixed_spacing` is provided or if `use_all_obs = true`.
### Returns
* A `Dict` with an estimated volatility for each asset.
"""
function simple_volatility(ts::SortedDataFrame, assets::Vector{Symbol} = get_assets(ts);
                           time_grid::Union{Missing,Dict} = missing , fixed_spacing::Union{Missing,Dict,<:Real} = missing,
                           use_all_obs::Bool = false, rough_guess_number_of_intervals::Integer = 5)
    if ismissing(time_grid)
        time_grid = Dict{Symbol,Vector{eltype(ts.df[:,ts.time])}}()
        if use_all_obs
            for a in assets
                time_grid[a] = ts.df[ts.groupingrows[a],ts.time]
            end
        elseif !ismissing(fixed_spacing)
            if isa(fixed_spacing, Dict)
                for a in assets
                    time_grid[a] = collect(minimum(ts.df[:,ts.time]):fixed_spacing[a]:maximum(ts.df[:,ts.time]))
                end
            else
                for a in assets
                    time_grid[a] = collect(minimum(ts.df[:,ts.time]):fixed_spacing:maximum(ts.df[:,ts.time]))
                end
            end
        else
            n_grid = default_spacing(ts; rough_guess_number_of_intervals = rough_guess_number_of_intervals)
            for a in assets
                time_grid[a] = collect(minimum(ts.df[:,ts.time]):n_grid[a]:maximum(ts.df[:,ts.time]))
            end
        end
    end
    return simple_volatility_with_grid(ts, assets, time_grid)
end

"""
    default_spacing(ts::SortedDataFrame; rough_guess_number_of_intervals::Integer = 5,
                    T = duration(ts; in_dates_period = false))

Calculates a default spacing between returns to use.
This comes from the equation at section 1.2.3 of Zhang, Mykland, Ait-Sahalia 2005.
### Inputs
* `ts` - The tick data.
* `rough_guess_number_of_intervals` - A rough estimate of how many intervals to split the tick data into. This is used in a first pass to estimate the optimal interval spacing.
* `T` - The duration of the tick data.
### Returns
* A scalar representing the optimal interval spacing.
"""
function default_spacing(ts::SortedDataFrame; rough_guess_number_of_intervals::Integer = 5,
                         T::Real = duration(ts; in_dates_period = false))
    rough_vol_guess, rough_micro_guess  = two_scales_volatility(ts; num_grids = rough_guess_number_of_intervals)
    n_guess = Dict{Symbol,eltype(ts.df[:,ts.value])}()
    for a in keys(rough_vol_guess)
        n_guess[a] = ( (T/(4* rough_micro_guess[a]^2)) * T * rough_vol_guess[a]^4 )^(1/3)
        if (rough_micro_guess[a] < 0.0000001)
            n_guess[a] = duration(ts; in_dates_period = false)/20
        end
    end
    return n_guess
end


text:
function _reset_relaxation!(m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}

    m._working_problem._relaxed_evaluator.is_first_eval = true
    fill!(m._working_problem._relaxed_evaluator.subexpressions_eval, false)

    m._new_eval_constraint = true

    delete_nl_constraints!(m)
    delete_objective_cuts!(m)

    return nothing
end

"""
$(SIGNATURES)

Runs interval, linear, quadratic contractor methods followed by obbt and a
constraint programming walk up to tolerances specified in
`EAGO.GlobalOptimizer` object.
"""
function preprocess!(t::ExtensionType, m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}

    _reset_relaxation!(m)

    wp = m._working_problem

    # Sets initial feasibility
    feasible_flag = true
    m._obbt_performed_flag = false

    # compute initial volume
    m._initial_volume = prod(i -> _mid(BranchVar, m, i), 1:N)

    if m.fbbt_lp_depth >= m._iteration_count
        for i = 1:m.fbbt_lp_repetitions
            if feasible_flag
                for saf_leq in wp._saf_leq
                    !feasible_flag && break
                    feasible_flag &= fbbt!(m, saf_leq)
                end
                !feasible_flag && break

                for saf_eq in  wp._saf_eq
                    !feasible_flag && break
                    feasible_flag &= fbbt!(m, saf_eq)
                end
                !feasible_flag && break
            end
        end
    end

    # done after cp to prevent using cp specific flags in cut generation
    set_first_relax_point!(m)

    cp_walk_count = 0
    perform_cp_walk_flag = feasible_flag
    perform_cp_walk_flag &= (m.cp_depth >= m._iteration_count)
    perform_cp_walk_flag &= (cp_walk_count < m.cp_repetitions)
    while perform_cp_walk_flag
        feasible_flag &= set_constraint_propagation_fbbt!(m)
        !feasible_flag && break
        cp_walk_count += 1
        perform_cp_walk_flag = (cp_walk_count < m.cp_repetitions)
    end

    obbt_count = 0
    perform_obbt_flag = feasible_flag
    perform_obbt_flag &= (m.obbt_depth >= m._iteration_count)
    perform_obbt_flag &= (obbt_count < m.obbt_repetitions)

    while perform_obbt_flag
        feasible_flag &= obbt!(m)
        m._obbt_performed_flag = true
        !feasible_flag && break
        obbt_count += 1
        perform_obbt_flag = (obbt_count < m.obbt_repetitions)
    end

    m._final_volume = prod(i -> _mid(BranchVar, m, i), 1:N)
    m._preprocess_feasibility = feasible_flag

    return nothing
end


text:
import SparsityTracing; const AD = SparsityTracing # short name for convenience

using SparseArrays
using LinearAlgebra  # for I

using Test, Documenter
# using BenchmarkTools

@testset "SparsityTracing" begin

    doctest(SparsityTracing; manual=false)

    @testset "createviewassign" begin

        x = [24.0, 42.0]

        x_ad = AD.create_advec(x)

        x_jac = AD.jacobian(x_ad, length(x))
        @test x_jac == Matrix(I, 2, 2)

        a_ad = x_ad[1]
        b_ad = x_ad[2]

        @test a_ad.val   == 24.0
        @test AD.deriv(a_ad, 2) == [1, 0]
        @test b_ad.val   == 42.0
        @test AD.deriv(b_ad, 2) == [0, 1]

        ax_ad = view(x_ad, 1)
        bx_ad = view(x_ad, 2)

        @test ax_ad[].val   == 24.0
        @test AD.deriv(ax_ad[], 2) == [1, 0]
        @test bx_ad[].val   == 42.0
        @test AD.deriv(bx_ad[], 2) == [0, 1]

        y_ad = similar(x_ad)
        ay_ad = view(y_ad, 1)
        by_ad = view(y_ad, 2)

        ay_ad[] = ax_ad[]
        by_ad[] = ax_ad[]
        @test ay_ad[].val   == 24.0
        @test AD.deriv(ay_ad[], 2) == [1, 0]
        @test by_ad[].val   == 24.0
        @test AD.deriv(by_ad[], 2) == [1, 0]

        by_ad[] = 25
        @test by_ad[].val   == 25.0
        @test AD.deriv(by_ad[], 2) == [0, 0]
        @test ax_ad[].val   == 24.0
        @test AD.deriv(ax_ad[], 2) == [1, 0]

    end

    @testset "arithmetic" begin

        x = [24.0, 42.0]

        x_ad = AD.create_advec(x)

        a_ad = x_ad[1]
        b_ad = x_ad[2]

        # addition
        c_ad = a_ad + b_ad
        @test c_ad.val   == 24.0+42.0
        @test AD.deriv(c_ad, 2) == [1, 1]

        c_ad = a_ad + 2
        @test c_ad.val   == 24.0+2
        @test AD.deriv(c_ad, 2) == [1, 0]

        # check handling of ad variables with no derivatives
        x_no_ad = copy(x_ad)
        for i = 1:length(x_no_ad)
            x_no_ad[i] = AD.value(x_no_ad[i])
        end
        @test AD.deriv(x_no_ad[1], 2) == [0, 0]
        @test AD.deriv(x_no_ad[2], 2) == [0, 0]
        a_no_ad = x_no_ad[1]
        b_no_ad = x_no_ad[2]
        c_ad = a_no_ad + b_ad
        @test c_ad.val   == 24.0+42.0
        @test AD.deriv(c_ad, 2) == [0, 1]  # no deriv wrt x_ad[1]
        c_ad = a_ad + b_no_ad
        @test c_ad.val   == 24.0+42.0
        @test AD.deriv(c_ad, 2) == [1, 0]  # no deriv wrt x_ad[2]
        c_ad = a_no_ad + b_no_ad
        @test c_ad.val   == 24.0+42.0
        @test AD.deriv(c_ad, 2) == [0, 0]  # no deriv

        # unary minus
        c_ad = -a_ad
        @test c_ad.val   == -24.0
        @test AD.deriv(c_ad, 2) == -[1, 0]

        # subtraction
        c_ad = a_ad - b_ad
        @test c_ad.val   == 24.0-42.0
        @test AD.deriv(c_ad, 2) == [1, -1]

        c_ad = a_ad - 2
        @test c_ad.val   == 24.0-2
        @test AD.deriv(c_ad, 2) == [1, 0]

        c_ad = 2 - a_ad
        @test c_ad.val   == 2-24.0
        @test AD.deriv(c_ad, 2) == -[1, 0]

        #multiplication
        c_ad = a_ad*b_ad
        @test c_ad.val   == 24.0*42.0
        @test AD.deriv(c_ad, 2) == [42, 24]

        c_ad = a_ad*2
        @test c_ad.val   == 24.0*2
        @test AD.deriv(c_ad, 2) == [1, 0]*2

        # division
        c_ad = a_ad/b_ad
        @test c_ad.val   == 24.0/42.0
        # println(c_ad)
        @test AD.deriv(c_ad, 2) == [1/42, -24/42^2]

        c_ad = a_ad/2
        @test c_ad.val   == 24.0/2
        @test AD.deriv(c_ad, 2) == [1, 0]/2

        c_ad = 2/a_ad
        @test c_ad.val   == 2/24.0
        @test AD.deriv(c_ad, 2) == [-2/24^2, 0]

        # power
        c_ad = a_ad^2
        @test c_ad.val   == 24.0^2
        @test AD.deriv(c_ad, 2) == [1*24*2, 0]

        # max
        x_max = max.(x_ad, 30)
        @test x_max[1].val == 30
        @test x_max[2].val == 42.0
        @test AD.deriv(x_max[1], 2) == [0, 0]
        @test AD.deriv(x_max[2], 2) == [0, 1]
    end

    @testset "matmul" begin
        B = SparseArrays.sprand(10,10,0.5)
        C = AD.create_advec(ones(10))

        # A = AD.create_advec(zeros(10))
        A = B*C
        Ajac = AD.jacobian(A, length(C))

        @views A_view = B*C
        Ajac_view = AD.jacobian(A_view, length(C))
        @test Ajac_view == Ajac

        A_range = AD.create_advec(zeros(10))
        A_range .= 0.0

        A_range[2:5]  = B[2:5, :]*C
        A_range_jac = AD.jacobian(A_range, length(A_range))

        A_range .= 0.0
        @views A_range[2:5]  = B[2:5, :]*C  # fails - fills zeros in derivatives? (-> @views is removing sparsity ??)
        A_range_views_jac = AD.jacobian(A_range, length(A_range))
        @test A_range_views_jac == A_range_jac
        @test_broken nnz(A_range_views_jac) == nnz(A_range_jac) # fails -  A_range_view_jac has zero elements stored

        A_range .= 0.0
        mul!(view(A_range, 2:5), B[2:5, :], C)
        A_range_mul_jac = AD.jacobian(A_range, length(A_range))
        @test A_range_mul_jac == A_range_jac
        @test nnz(A_range_mul_jac) == nnz(A_range_jac)

        A_range .= 0.0
        @views mul!(A_range[2:5], B[2:5, :], C) # fails - fills zero in derivatives?
        A_range_mul_views_jac = AD.jacobian(A_range, length(A_range))
        @test A_range_mul_views_jac == A_range_jac
        @test_broken nnz(A_range_mul_views_jac) == nnz(A_range_jac) # fails -  A_range_mul_views_jac has zero elements stored

        # Looks like this is a symptom of a more fundamental problem with (lack of) views into sparse matrices,
        # and reverting to slow dense path:
        # https://stackoverflow.com/questions/58699267/julia-view-of-sparse-matrix
        # https://github.com/JuliaLang/julia/issues/21796

    end

    @testset "rober" begin
        function rober(du,u,p)
            y₁,y₂,y₃ = u
            k₁,k₂,k₃ = p
            du[1] = -k₁*y₁+k₃*y₂*y₃
            du[2] =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
            du[3] =  k₂*y₂^2
            nothing
        end

        function rober_jac(J,u,p)
            y₁,y₂,y₃ = u
            k₁,k₂,k₃ = p
            J[1,1] = k₁ * -1
            J[2,1] = k₁
            J[3,1] = 0
            J[1,2] = y₃ * k₃
            J[2,2] = y₂ * k₂ * -2 + y₃ * k₃ * -1
            J[3,2] = y₂ * 2 * k₂
            J[1,3] = k₃ * y₂
            J[2,3] = k₃ * y₂ * -1
            J[3,3] = 0
            nothing
        end

        p = (0.04,3e7,1e4)

        u = [1.0, 2.0, 3.0]
        Janalytic = zeros(3, 3)
        rober_jac(Janalytic, u, p)

        @info "rober J analytic\n $Janalytic"

        u_ad = AD.create_advec(u)
        du_ad = similar(u_ad)
        rober(du_ad, u_ad, p)

        Jad = AD.jacobian(du_ad, length(du_ad))
        @info "rober J AD\n $Jad"

        @test Janalytic == Jad
    end

    return nothing
end


text:
# take unaligned seed -> make a msa

function build_model(fileseed::String, filefull::String, ctype::Symbol, L::Int64;
		    filename_ins::String="LambdaOpen_LambdaExt.dat",
		    filename_par::String="Parameters_PlmDCA.dat",
		    filename_gap::String="Gap_Ext_Int.dat",
		    Mtest::Int64=0,
		    verbose::Bool=true)

	if ctype != :amino && ctype != :nbase
		error("Wrong second argument: choose between :amino and :nbase")
		return
	end
	println("### Reading seed ###")
	if ctype == :nbase
		print_pos = true
	else
		print_pos = false
	end
	seed = DCAlign.readfull(fileseed, ctype=ctype, pos = true)
	println("### Inferring insertions penalties ###")
	l_o, l_e = infer_ins_pen(seed, L)

	println("### Inferring a Potts model using PlmDCA ###")
	aligntmp = filter_insertions(seed)
	println("Temporary FASTA alignment in ", aligntmp)
	PlmData = PlmDCA.plmdca(aligntmp, theta=0.20)

	println("### Finding gap penalties ###")
	println("WARNING: Reasonable values are obtained when using many (> 500) sequences")

	full = DCAlign.readfull(filefull, ctype=ctype, pos = print_pos)
	if Mtest == 0
		Mtest = length(seed)
		println("Using all seed sequences to get the gap penalties")
	else
		println("Using ", Mtest, " out of ", length(seed), " to get the gap penalties")
	end
	fulltmp = extract_full_seq(full, seed, Mtest, ctype=ctype)
	println("Temporary full length sequences in ", fulltmp)

	if ctype == :amino
		q = 21
	elseif ctype == :nbase
		q = 5
	end

	print_results(filename_ins, l_o, l_e, filename_par, PlmData, ctype, L, filename_gap, 0.0, 0.0) # tmp
	mu = 0.00:0.50:4.00
	muint = 0.00:0.50:4.00
	d = zeros(length(mu),length(muint))
	aseed = DCAlign.readfull(aligntmp, ctype=ctype, pos = true)
	for a in 1:length(mu)
		for b in 1:length(muint)
			println("#### Aligning ", Mtest, " sequences using (μext, μint) = (", mu[a], ", ", muint[b],")")
			filename_out = tempname()
			filename_flag = tempname()
			DCAlign.align_all(q, L, filename_par, fulltmp, filename_ins, mu[a], muint[b]; typel=:plm, filename_flag=filename_flag, filename_align=aligntmp, filename_ins=fileseed, filename_out=filename_out, verbose = verbose, maxiter = 300)
			tmpseed = DCAlign.readfull(filename_out, ctype=ctype, pos = true)
			d[a,b] = compute_average_dist(aseed, tmpseed)
			rm(filename_out)
			rm(filename_flag)
		end
	end
	#println(d)
	aux = argmin(d)
	muext_best = mu[aux[1]]
	muint_best = muint[aux[2]]
	print_results(filename_ins, l_o, l_e, filename_par, PlmData, ctype, L, filename_gap, muext_best, muint_best)
	rm(fulltmp)
	rm(aligntmp)
	println("Done!")

end


text:
using ParallelAccelerator

@acc function add()
  # w = 5
  # h = 5
  # A = fill(1, w, h)
  # B = fill(1, w, h)
  # C = A + B
  a = 1 + 2
  b = 3 + 4 +5
  return 1 + 2 + a + b
end

function main()
    tic()
    D = add()
    println("SELFTIMED ", toq())

    println(D)

end

main()

max_stars_repo_name: malekbr/julia-to-tiramisu

max_issues_repo_name: malekbr/julia-to-tiramisu

max_forks_repo_name: malekbr/julia-to-tiramisu


text:
@testset "Grid Search" begin
    @test Optim.grid_search(x -> (1.0 - x)^2, [-1:0.1:1.0;]) == 1
end
# Cartesian product over 2 dimensions.
#grid_search(x -> (1.0 - x[1])^2 + (2.0 - x[2])^2, product[-1:0.1:1.0, -1:0.1:1.0])