text:
Base.summary(r::OptimizationResults) = summary(r.method) # might want to do more here than just return summary of the method used
minimizer(r::OptimizationResults) = r.minimizer
minimum(r::OptimizationResults) = r.minimum
iterations(r::OptimizationResults) = r.iterations
iteration_limit_reached(r::OptimizationResults) = r.iteration_converged
trace(r::OptimizationResults) = length(r.trace) > 0 ? r.trace : error("No trace in optimization results. To get a trace, run optimize() with store_trace = true.")

function x_trace(r::UnivariateOptimizationResults)
    tr = trace(r)
    !haskey(tr[1].metadata, "minimizer") && error("Trace does not contain x. To get a trace of x, run optimize() with extended_trace = true")
    [ state.metadata["minimizer"] for state in tr ]
end
function x_lower_trace(r::UnivariateOptimizationResults)
    tr = trace(r)
    !haskey(tr[1].metadata, "x_lower") && error("Trace does not contain x. To get a trace of x, run optimize() with extended_trace = true")
    [ state.metadata["x_lower"] for state in tr ]
end
x_lower_trace(r::MultivariateOptimizationResults) = error("x_lower_trace is not implemented for $(summary(r)).")
function x_upper_trace(r::UnivariateOptimizationResults)
    tr = trace(r)
    !haskey(tr[1].metadata, "x_upper") && error("Trace does not contain x. To get a trace of x, run optimize() with extended_trace = true")
    [ state.metadata["x_upper"] for state in tr ]
end
x_upper_trace(r::MultivariateOptimizationResults) = error("x_upper_trace is not implemented for $(summary(r)).")

function x_trace(r::MultivariateOptimizationResults)
    tr = trace(r)
    if isa(r.method, NelderMead)
        throw(ArgumentError("Nelder Mead does not operate with a single x. Please use either centroid_trace(...) or simplex_trace(...) to extract the relevant points from the trace."))
    end
    !haskey(tr[1].metadata, "x") && error("Trace does not contain x. To get a trace of x, run optimize() with extended_trace = true")
    [ state.metadata["x"] for state in tr ]
end

function centroid_trace(r::MultivariateOptimizationResults)
    tr = trace(r)
    if !isa(r.method, NelderMead)
        throw(ArgumentError("There is no centroid involved in optimization using $(r.method). Please use x_trace(...) to grab the points from the trace."))
    end
    !haskey(tr[1].metadata, "centroid") && error("Trace does not contain centroid. To get a trace of the centroid, run optimize() with extended_trace = true")
    [ state.metadata["centroid"] for state in tr ]
end
function simplex_trace(r::MultivariateOptimizationResults)
    tr = trace(r)
    if !isa(r.method, NelderMead)
        throw(ArgumentError("There is no simplex involved in optimization using $(r.method). Please use x_trace(...) to grab the points from the trace."))
    end
    !haskey(tr[1].metadata, "simplex") && error("Trace does not contain simplex. To get a trace of the simplex, run optimize() with trace_simplex = true")
    [ state.metadata["simplex"] for state in tr ]
end
function simplex_value_trace(r::MultivariateOptimizationResults)
    if !isa(r.method, NelderMead)
        throw(ArgumentError("There are no simplex values involved in optimization using $(r.method). Please use f_trace(...) to grab the objective values from the trace."))
    end
    !haskey(tr[1].metadata, "simplex_values") && error("Trace does not contain objective values at the simplex. To get a trace of the simplex values, run optimize() with trace_simplex = true")
    [ state.metadata["simplex_values"] for state in tr ]
end

f_trace(r::OptimizationResults) = [ state.value for state in trace(r) ]
g_norm_trace(r::OptimizationResults) = error("g_norm_trace is not implemented for $(summary(r)).")
g_norm_trace(r::MultivariateOptimizationResults) = [ state.g_norm for state in trace(r) ]

f_calls(r::OptimizationResults) = r.f_calls
f_calls(d) = first(d.f_calls)

g_calls(r::OptimizationResults) = error("g_calls is not implemented for $(summary(r)).")
g_calls(r::MultivariateOptimizationResults) = r.g_calls
g_calls(d::NonDifferentiable) = 0
g_calls(d) = first(d.df_calls)

h_calls(r::OptimizationResults) = error("h_calls is not implemented for $(summary(r)).")
h_calls(r::MultivariateOptimizationResults) = r.h_calls
h_calls(d::Union{NonDifferentiable, OnceDifferentiable}) = 0
h_calls(d) = first(d.h_calls)
h_calls(d::TwiceDifferentiableHV) = first(d.hv_calls)

converged(r::UnivariateOptimizationResults) = r.converged
converged(r::MultivariateOptimizationResults) = r.x_converged || r.f_converged || r.g_converged
x_converged(r::OptimizationResults) = error("x_converged is not implemented for $(summary(r)).")
x_converged(r::MultivariateOptimizationResults) = r.x_converged
f_converged(r::OptimizationResults) = error("f_converged is not implemented for $(summary(r)).")
f_converged(r::MultivariateOptimizationResults) = r.f_converged
f_increased(r::OptimizationResults) = error("f_increased is not implemented for $(summary(r)).")
f_increased(r::MultivariateOptimizationResults) = r.f_increased
g_converged(r::OptimizationResults) = error("g_converged is not implemented for $(summary(r)).")
g_converged(r::MultivariateOptimizationResults) = r.g_converged

x_abstol(r::OptimizationResults) = error("x_abstol is not implemented for $(summary(r)).")
x_reltol(r::OptimizationResults) = error("x_reltol is not implemented for $(summary(r)).")
x_tol(r::OptimizationResults) = error("x_tol is not implemented for $(summary(r)).")
x_abstol(r::MultivariateOptimizationResults) = r.x_abstol
x_reltol(r::MultivariateOptimizationResults) = r.x_reltol
x_tol(r::MultivariateOptimizationResults) = r.x_abstol
x_abschange(r::MultivariateOptimizationResults) = r.x_abschange
x_relchange(r::MultivariateOptimizationResults) = r.x_relchange
f_abstol(r::OptimizationResults) = error("f_abstol is not implemented for $(summary(r)).")
f_reltol(r::OptimizationResults) = error("f_reltol is not implemented for $(summary(r)).")
f_tol(r::OptimizationResults) = error("f_tol is not implemented for $(summary(r)).")
f_tol(r::MultivariateOptimizationResults) = r.f_reltol
f_abstol(r::MultivariateOptimizationResults) = r.f_abstol
f_reltol(r::MultivariateOptimizationResults) = r.f_reltol
f_abschange(r::MultivariateOptimizationResults) = r.f_abschange
f_relchange(r::MultivariateOptimizationResults) = r.f_relchange

g_tol(r::OptimizationResults) = error("g_tol is not implemented for $(summary(r)).")
g_tol(r::MultivariateOptimizationResults) = r.g_abstol
g_residual(r::MultivariateOptimizationResults) = r.g_residual

initial_state(r::OptimizationResults) = error("initial_state is not implemented for $(summary(r)).")
initial_state(r::MultivariateOptimizationResults) = r.initial_x

lower_bound(r::OptimizationResults) = error("lower_bound is not implemented for $(summary(r)).")
lower_bound(r::UnivariateOptimizationResults) = r.initial_lower
upper_bound(r::OptimizationResults) = error("upper_bound is not implemented for $(summary(r)).")
upper_bound(r::UnivariateOptimizationResults) = r.initial_upper

rel_tol(r::OptimizationResults) = error("rel_tol is not implemented for $(summary(r)).")
rel_tol(r::UnivariateOptimizationResults) = r.rel_tol
abs_tol(r::OptimizationResults) = error("abs_tol is not implemented for $(summary(r)).")
abs_tol(r::UnivariateOptimizationResults) = r.abs_tol

time_limit(r::MultivariateOptimizationResults) = r.time_limit
time_run(  r::MultivariateOptimizationResults) = r.time_run

time_limit(r::OptimizationResults) = error("time_limit is not implemented for $(summary(r)).")
time_run(  r::OptimizationResults) = error("time_run is not implemented for $(summary(r)).")


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

JLLWrappers.@generate_main_file_header("oneAPI_Level_Zero_Loader")
JLLWrappers.@generate_main_file("oneAPI_Level_Zero_Loader", UUID("13eca655-d68d-5b81-8367-6d99d727ab01"))
end  # module oneAPI_Level_Zero_Loader_jll

ext: jl


text:
module Monads

import Base.==

# types
export Monad, Identity, MList, Maybe, State, MPromise
# combinators
export mreturn, join, fmap, mbind, mcomp, mthen, (>>)
# utilities
export liftM
# do syntax
export @mdo
# MonadPlus
export MonadPlus, mzero, mplus, guard
# State
export runState, put, get, evalState, execState

abstract type Monad end
abstract type MonadPlus <: Monad end

## Buy two monad combinators, get the third free!
mreturn(::Type{M}, val) where {M<:Monad} = M(val)
join(m::Monad) = mbind(identity, m)
fmap(f::Function, m::M) where {M<:Monad} = mbind(m) do x
    mreturn(M, f(x))
end
mbind(f::Function, m::Monad) = join(fmap(f, m))

## Extra combinators
mcomp(g::Function, f::Function) = x -> mbind(g, f(x))
mthen(k::Monad, m::Monad) = mbind(_ -> k, m)
(>>)(m::Monad, k::Monad) = mthen(k, m)
(==)(m::Monad, k::Monad) = m.value == k.value

## A MonadPlus function
guard(::Type{M}, c::Bool) where {M<:MonadPlus} = c ? mreturn(M, nothing) : mzero(M)

## Friendly monad blocks
macro mdo(mtype, body)
    esc(mdo_desugar(mdo_patch(mtype, body)))
end

## patch up functions to insert the right monad
mdo_patch(mtype, expr) = expr
function mdo_patch(mtype, expr::Expr)
    expr.args = map(arg->mdo_patch(mtype, arg), expr.args)
    if expr.head == :return
        expr.head = :call
        insert!(expr.args, 1, :mreturn)
    end
    if expr.head == :call && any(expr.args[1] .== [:mreturn, :mzero, :guard, :liftM])
        insert!(expr.args, 2, mtype)
    end
    expr
end

## desugaring mdo syntax is a right fold
mdo_desugar(exprIn) = foldr(mdo_desugar_helper, exprIn.args)
mdo_desugar_helper(rest, expr) = rest
function mdo_desugar_helper(expr::Expr, rest)
    if (expr.head == :call
      && expr.args[1] == Symbol("<")
        && expr.args[3].args[1] == Symbol("-"))
        # replace "<-" with monadic binding
        newExpr = expr.args[3].args[2:end];
        quote
            mbind($(newExpr...)) do $(expr.args[2])
                $rest
            end
        end
    elseif expr.head == :(=)
        # replace assignment with let binding
        quote
            let
                $expr
                $rest
            end
        end
    elseif expr.head == :line
        rest
    elseif rest == :()
        expr
    else
        # replace with sequencing
        :(mthen($rest, $expr))
    end
end

## Function lifting
liftM(::Type{M}, f::Function) where {M<:Monad} = m1 -> @mdo M begin
    x1 <- m1
    return f(x1)
end

## Starting slow: Identity
struct Identity{T} <: Monad
    value :: T
end

mbind(f::Function, m::Identity) = f(m.value)

## List
struct MList <: MonadPlus
    value :: Vector

    MList(x::Array) = new(vec(x))
    MList(x) = new([x])
end

function join(m::MList)
    if !isempty(m.value)
        val = nothing
        for arr in m.value[1:end]
            if !isempty(arr.value)
                if val === nothing
                    val = arr.value
                else
                    append!(val, arr.value)
                end
            end
        end
        if val === nothing
            mzero(MList)
        else
            mreturn(MList, val)
        end
    else
        mzero(MList)
    end
end
fmap(f::Function, m::MList) = MList(map(f, m.value))

mbind(f::Function, m::MList) = join(fmap(f, m))

# It's also a MonadPlus
mzero(::Type{MList}) = MList([])
mplus(m1::MList, m2::MList) = join(MList([m1, m2]))

struct Nothing end

## Maybe
struct Maybe{T} <: Monad
    value :: Union{T, Nothing}
end

mbind(f::Function, m::Maybe) = isa(m.value, Nothing) ? Maybe(nothing) : f(m.value)

## State
struct State <: Monad
    runState :: Function # s -> (a, s)
end
state(f) = State(f)

runState(s::State) = s.runState
runState(s::State, st) = s.runState(st)

function mbind(f::Function, s::State)
      state(st -> begin
          (x, stp) = runState(s, st)
          runState(f(x), stp)
            end
            )
end
mreturn(::Type{State}, x) = state(st -> (x, st))

put(newState) = state(_ -> (nothing, newState))
get() = state(st -> (st, st))

evalState(s::State, st) = runState(s, st)[1]
execState(s::State, st) = runState(s, st)[2]

end


text:
import RestrictedBoltzmannMachines as RBMs
import ConvolutionalRBMs as ConvRBMs

using Test: @test, @testset, @inferred, @test_throws
using Random: bitrand, randn!
using LinearAlgebra: dot
using LogExpFunctions: logsumexp, logaddexp
using RestrictedBoltzmannMachines: visible, hidden, weights
using RestrictedBoltzmannMachines: energy, free_energy, interaction_energy
using RestrictedBoltzmannMachines: inputs_v_to_h, inputs_h_to_v
using RestrictedBoltzmannMachines: sample_v_from_v, sample_h_from_h, sample_v_from_h, sample_h_from_v
using ConvolutionalRBMs: vsizes, hsizes, BinaryConvRBM, output_size, out2in, hankel

@testset "free energy" begin
    rbm = BinaryConvRBM(randn(1), randn(1), randn(1,1,1); pool=false)
    β = rand()
    v = bitrand(1,2)
    hs = [[h1 h2] for h1=0:1, h2=0:1]
    -β * free_energy(rbm, v; β) ≈ logsumexp(-β * energy(rbm, v, h) for h in hs)

    rbm = BinaryConvRBM(randn(1), randn(1), randn(1,1,1); pool=true)
    β = rand()
    v = bitrand(1,2)
    hs = [[0 1], [1 0]]
    K = prod(output_size(rbm,v))
    M = prod(size(hidden(rbm)))
    @test exp(-β * free_energy(rbm, v; β)) ≈ sum(exp(-β * energy(rbm, v, h)) for h in hs) + K^M * exp(-β * energy(rbm, v, [0 0]))
end


text:
using MeasureTheory, Static, BenchmarkTools, LinearAlgebra, StaticArrays
using MultivariateMeasures
x = randn(200,200);
Σ = x' * x;

using PositiveFactorizations
CΣ = cholesky(Positive, Σ)

function benchmark()
    σtimes = Float64[]
    ωtimes = Float64[]
    y = randn(200);
    for k in 1:100
        U = UpperTriangular(SizedMatrix{k,k}(CΣ.U.data[1:k,1:k]));
        C = Cholesky(U, 'U', 0);
        yk = y[1:k]
        push!(σtimes, @belapsed logdensity(MvNormal(σ=$C), $yk))
        push!(ωtimes, @belapsed logdensity(MvNormal(ω=$C), $yk))
    end

    (σtimes, ωtimes)
end

(σtimes, ωtimes) = benchmark()

# using PDMats
# C = cholesky(Σ[1:20,1:20])
# Dists.MvNormal(PDMat(C))

# @btime Dists.logpdf(Dists.MvNormalCanon(PDMat($C)), $y)
# @btime Dists.logpdf(Dists.MvNormal(PDMat($C)), $y)

using Plots

plot(σtimes[1:end] .* 1e6; dpi=300, legend=150, label="σ parameterization")
plot!(ωtimes[1:end] .* 1e6, label="ω parameterization")
xlabel!("Dimensions")
ylabel!("Time (μs)")
title!("MeasureTheory.jl Log-density for MvNormal")
savefig("mvnormal.svg")
savefig("mvnormal.png")

# C_unsized = Cholesky(C.UL.data.data, 'U', )

# using PDMats
# PD = PDMat(C_unsized)

# julia> @btime Dists.logpdf(Dists.MvNormalCanon($PD), $y)
#   4.944 μs (6 allocations: 1.97 KiB)
# -124.71152557023147

# julia> @btime Dists.logpdf(Dists.MvNormalCanon(PDMat($C_unsized)), $y)
#   16.552 μs (10 allocations: 27.25 KiB)
# -124.71152557023147

# julia> @btime logdensity(MvNormal(ω=$C), $y)
#   214.948 ns (0 allocations: 0 bytes)
# -124.53160601776467

# L = inv(C.L)
# Σ = L' * L

# @btime Dists.logpdf(Dists.MvNormalCanon($PD), $y)
# @btime Dists.logpdf(Dists.MvNormalCanon(PDMat($C_unsized)), $y)

Dists.logpdf(Dists.MvNormal(Σ), y)
logdensity(MvNormal(ω=C), y)

# n = 40;
# t = CorrCholesky(StaticInt(n));
# C = transform(t, randn(dimension(t)));
@btime MeasureTheory.logdet_pos($C)

# @btime logdensity(MvNormal(ω=$C), $(randn(n)))


text:
export iauApco
"""
For a terrestrial observer, prepare star-independent astrometry
parameters for transformations between ICRS and observed
coordinates.  The caller supplies the Earth ephemeris, the Earth
rotation information and the refraction constants as well as the
site coordinates.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   date1  double       TDB as a 2-part...
   date2  double       ...Julian Date (Note 1)
   ebpv   double[2][3] Earth barycentric PV (au, au/day, Note 2)
   ehp    double[3]    Earth heliocentric P (au, Note 2)
   x,y    double       CIP X,Y (components of unit vector)
   s      double       the CIO locator s (radians)
   theta  double       Earth rotation angle (radians)
   elong  double       longitude (radians, east +ve, Note 3)
   phi    double       latitude (geodetic, radians, Note 3)
   hm     double       height above ellipsoid (m, geodetic, Note 3)
   xp,yp  double       polar motion coordinates (radians, Note 4)
   sp     double       the TIO locator s' (radians, Note 4)
   refa   double       refraction constant A (radians, Note 5)
   refb   double       refraction constant B (radians, Note 5)

Returned:
   astrom iauASTROM*   star-independent astrometry parameters:
      pmt    double       PM time interval (SSB, Julian years)
      eb     double[3]    SSB to observer (vector, au)
      eh     double[3]    Sun to observer (unit vector)
      em     double       distance from Sun to observer (au)
      v      double[3]    barycentric observer velocity (vector, c)
      bm1    double       sqrt(1-|v|^2): reciprocal of Lorenz factor
      bpn    double[3][3] bias-precession-nutation matrix
      along  double       longitude + s' (radians)
      xpl    double       polar motion xp wrt local meridian (radians)
      ypl    double       polar motion yp wrt local meridian (radians)
      sphi   double       sine of geodetic latitude
      cphi   double       cosine of geodetic latitude
      diurab double       magnitude of diurnal aberration vector
      eral   double       "local" Earth rotation angle (radians)
      refa   double       refraction constant A (radians)
      refb   double       refraction constant B (radians)

Notes:

   1. The TDB date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TDB)=2450123.7 could be expressed in any of these ways, among
   others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in cases
   where the loss of several decimal digits of resolution is
   acceptable.  The J2000 method is best matched to the way the
   argument is handled internally and will deliver the optimum
   resolution.  The MJD method and the date & time methods are both
   good compromises between resolution and convenience.  For most
   applications of this function the choice will not be at all
   critical.

   TT can be used instead of TDB without any significant impact on
   accuracy.

   2. The vectors eb, eh, and all the astrom vectors, are with respect
   to BCRS axes.

   3. The geographical coordinates are with respect to the WGS84
   reference ellipsoid.  TAKE CARE WITH THE LONGITUDE SIGN
   CONVENTION:  the longitude required by the present function is
   right-handed, i.e. east-positive, in accordance with geographical
   convention.

   4. xp and yp are the coordinates (in radians) of the Celestial
   Intermediate Pole with respect to the International Terrestrial
   Reference System (see IERS Conventions), measured along the
   meridians 0 and 90 deg west respectively.  sp is the TIO locator
   s', in radians, which positions the Terrestrial Intermediate
   Origin on the equator.  For many applications, xp, yp and
   (especially) sp can be set to zero.

   Internally, the polar motion is stored in a form rotated onto the
   local meridian.

   5. The refraction constants refa and refb are for use in a
   dZ = A*tan(Z)+B*tan^3(Z) model, where Z is the observed
   (i.e. refracted) zenith distance and dZ is the amount of
   refraction.

   6. It is advisable to take great care with units, as even unlikely
   values of the input parameters are accepted and processed in
   accordance with the models used.

   7. In cases where the caller does not wish to provide the Earth
   Ephemeris, the Earth rotation information and refraction
   constants, the function iauApco13 can be used instead of the
   present function.  This starts from UTC and weather readings etc.
   and computes suitable values using other SOFA functions.

   8. This is one of several functions that inserts into the astrom
   structure star-independent parameters needed for the chain of
   astrometric transformations ICRS <-> GCRS <-> CIRS <-> observed.

   The various functions support different classes of observer and
   portions of the transformation chain:

         functions         observer        transformation

      iauApcg iauApcg13    geocentric      ICRS <-> GCRS
      iauApci iauApci13    terrestrial     ICRS <-> CIRS
      iauApco iauApco13    terrestrial     ICRS <-> observed
      iauApcs iauApcs13    space           ICRS <-> GCRS
      iauAper iauAper13    terrestrial     update Earth rotation
      iauApio iauApio13    terrestrial     CIRS <-> observed

   Those with names ending in "13" use contemporary SOFA models to
   compute the various ephemerides.  The others accept ephemerides
   supplied by the caller.

   The transformation from ICRS to GCRS covers space motion,
   parallax, light deflection, and aberration.  From GCRS to CIRS
   comprises frame bias and precession-nutation.  From CIRS to
   observed takes account of Earth rotation, polar motion, diurnal
   aberration and parallax (unless subsumed into the ICRS <-> GCRS
   transformation), and atmospheric refraction.

   9. The context structure astrom produced by this function is used by
   iauAtioq, iauAtoiq, iauAtciq* and iauAticq*.

Called:
   iauAper      astrometry parameters: update ERA
   iauC2ixys    celestial-to-intermediate matrix, given X,Y and s
   iauPvtob     position/velocity of terrestrial station
   iauTrxpv     product of transpose of r-matrix and pv-vector
   iauApcs      astrometry parameters, ICRS-GCRS, space observer
   iauCr        copy r-matrix

This revision:   2013 October 9

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
function iauApco(date1::Real, date2::Real,
                 ebpv::AbstractMatrix{<:Real}, ehp::AbstractVector{<:Real},
                 x::Real, y::Real, s::Real, theta::Real,
                 elong::Real, phi::Real, hm::Real,
                 xp::Real, yp::Real, sp::Real,
                 refa::Real, refb::Real)

   # Allocate return value
   ref_astrom = Ref{iauASTROM}(iauASTROM())

   # Transpose to map Julia (FORTRAN) -> C style memory allocation
   ebpv = Array{Float64, 2}(ebpv') # Transpose input up front

   ccall((:iauApco, libsofa_c), Cvoid,
            (Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble},
            Cdouble, Cdouble, Cdouble,
            Cdouble, Cdouble, Cdouble,
            Cdouble, Cdouble, Cdouble,
            Cdouble, Cdouble, Cdouble,
            Ref{iauASTROM}),
            convert(Float64, date1),
            convert(Float64, date2),
            pointer(ebpv),
            pointer(ehp),
            convert(Float64, x),
            convert(Float64, y),
            convert(Float64, s),
            convert(Float64, theta),
            convert(Float64, elong),
            convert(Float64, phi),
            convert(Float64, hm),
            convert(Float64, xp),
            convert(Float64, yp),
            convert(Float64, sp),
            convert(Float64, refa),
            convert(Float64, refb),
            ref_astrom)

   return ref_astrom[]
end


text:
# ====================================================================
#                           USER INTERFACE
# ====================================================================

# ____________________________________________________________________
# Wrapper for user interface
#
struct UI{T}
    _w::T
end
wrappee(v::UI{T}) where T = getfield(v, :_w)

# ____________________________________________________________________
# Model
#
Model() = UI{Model}(Model(nothing))
Model(args::Vararg{Pair, N}) where N =
    addcomp!(Model(), args...)

propertynames(w::UI{Model}) = collect(keys(wrappee(w).comp))
getproperty(w::UI{Model}, s::Symbol) = UI{WComponent}(get(wrappee(w).comp, s, nothing))
getindex(w::UI{Model}, s::Symbol) = getproperty(w, s)

code(w::UI{Model}, id::Int=1) = wrappee(w).instruments[id].code

function (w::UI{Model})(id::Int, expr::Symbol)
    model = wrappee(w)
    c = model.instruments
    @assert length(c) >= 1 "No domain in model"
    @assert 1 <= id <= length(c) "Invalid index (allowed range: 1 : " * string(length(c)) * ")"
    instr = c[id]

    if expr == :domain
        ret = instr.domain
        (ndims(ret) == 1)  &&  (return ret[1])
        return ret
    end
    i = findall(expr .== instr.compnames)
    if length(i) == 1
        return instr.compevals[i[1]].result
    end
    i = findall(expr .== instr.exprnames)
    if length(i) == 1
        return instr.exprevals[i[1]]
    end
    error("No expression $expr available for instrument $id")
end

(w::UI{Model})(expr::Symbol) = w(1, expr)

function (w::UI{Model})(id::Int=1)
    model = wrappee(w)
    c = model.instruments[id]
    i = findall(c.exprcmp)
    @assert length(i) == 1 "No unambiguous expression for instrument $id"
    return c.exprevals[i[1]]
end

getparamvalues(w::UI{Model}) = getparamvalues(wrappee(w))
setparamvalues!(w::UI{Model}, pval::Vector{Float64}) = setparamvalues!(wrappee(w), pval)

function resetcounters!(w::UI{Model})
    model = wrappee(w)
    for instr in model.instruments
        instr.counter = 0
        for ce in instr.compevals
            ce.counter = 0
        end
    end
    return w
end

# ____________________________________________________________________
# Components
#
propertynames(w::UI{WComponent}) = [:fixed, collect(fieldnames(typeof(wrappee(w).comp)))...]
componenttype(w::UI{WComponent}) = typeof(wrappee(w).comp)

function getproperty(w::UI{WComponent}, s::Symbol)
    (s == :fixed)  &&  (return wrappee(w).fixed)
    return getfield(wrappee(w).comp, s)
end

function setproperty!(w::UI{WComponent}, s::Symbol, value)
    @assert s == :fixed "Can't assign a value to $s"
    wrappee(w).fixed = value
    return value
end

function addcomp!(w::UI{Model}, args::Vararg{Pair, N}) where N
    model = wrappee(w)

    for c in args
        @assert typeof(c[1]) == Symbol
        @assert typeof(c[2]) <: AbstractComponent
        @assert !(c[1] in keys(model.comp)) "Component $(c[1]) already exists"
    end

    for c in args
        newcomp = deepcopy(c[2])
        for pname in fieldnames(typeof(newcomp))
            par = getfield(newcomp, pname)
            if typeof(par) == Parameter
                par._private.model = model
                par._private.cname = c[1]
                par._private.pname = pname
                par._private.index = 0
            elseif typeof(par) == Vector{Parameter}
                for i in 1:length(par)
                    par[i]._private.model = model
                    par[i]._private.cname = c[1]
                    par[i]._private.pname = pname
                    par[i]._private.index = i
                end
            end
        end
        model.comp[c[1]] = WComponent(c[1], newcomp, false)
    end
    return w
end

# ____________________________________________________________________
# Parameter
#
function propertynames(p::Parameter)
    out = collect(fieldnames(Parameter))
    out = out[findall((out .!= :_private)  &&
                      (out .!= :cfixed))]
end

function setproperty!(p::Parameter, s::Symbol, value)
    bkp = p.expr
    setfield!(p, s, convert(typeof(getfield(p, s)), value))
    if (s == :expr)  &&  (p._private.model != nothing)
        for instr in p._private.model.instruments
            (p._private.cname in instr.compnames)  &&  (instr.compile = true)
        end
    end
    return value
end

# ____________________________________________________________________
# Instruments/domains
#
function add_dom!(w::UI{Model}, dom::AbstractDomain)
    model = wrappee(w)
    push!(model.instruments, Instrument(flatten(dom)))
    return length(model.instruments)
end

function add_dom!(w::UI{Model}, args...)
    model = wrappee(w)
    push!(model.instruments, Instrument(Domain(args...)))
    return length(model.instruments)
end

rm_dom!(w::UI{Model}, id::Int) = deleteat!(wrappee(w).instruments, id)

dom_count(w::UI{Model}) = length(wrappee(w).instruments)

# ____________________________________________________________________
# Expressions
#
addexpr!(w::UI{Model}, args...; kw...) = addexpr!(wrappee(w), args...; kw...)

replaceexpr!(w::UI{Model}, label::Symbol, expr::Expr) = replaceexpr!(w::UI{Model}, 1, label, expr)

function replaceexpr!(w::UI{Model}, expr::Expr)
    model = wrappee(w)
    i = findall(model.instruments[1].exprcmp)
    @assert length(i) == 1
    label = model.instruments[1].exprnames[i[1]]
    return replaceexpr!(w, 1, label, expr)
end

function replaceexpr!(w::UI{Model}, id::Int, label::Symbol, expr::Expr)
    model = wrappee(w)
    ii = findall(label .== model.instruments[id].exprnames)
    @assert length(ii) == 1 "No expression labelled $label in domain $id"
    model.instruments[id].exprs[ii[1]] = expr
    model.instruments[id].compile = true
    return model
end

setflag!(w::UI{Model}, label::Symbol, flag::Bool) = setflag!(w, 1, label, flag)
function setflag!(w::UI{Model}, id::Int, label::Symbol, flag::Bool)
    model = wrappee(w)
    ii = findall(label .== model.instruments[id].exprnames)
    @assert length(ii) == 1 "No expression labelled $label in domain $id"
    model.instruments[id].exprcmp[ii[1]] = flag
    model.instruments[id].compile = true
end

# ____________________________________________________________________
# Evaluate & fit
#
evaluate!(w::UI{Model}) = _evaluate!(wrappee(w))

fit!(w::UI{Model}, data::AbstractData; kw...) = _fit!(wrappee(w), [data]; kw...)
fit!(w::UI{Model}, data::Vector{T}; kw...) where T<:AbstractMeasures = _fit!(wrappee(w), data; kw...)
fit(w::UI{Model}, data::AbstractData; kw...) = _fit(wrappee(w), [data]; kw...)
fit(w::UI{Model}, data::Vector{T}; kw...) where T<:AbstractMeasures = _fit(wrappee(w), data; kw...)

function propertynames(w::UI{FitResult})
    res = wrappee(w)
    out = collect(fieldnames(typeof(res)))
    out = out[findall((out .!= :bestfit)  .&
                      (out .!= :fitter))]
    out = [out; collect(keys(res.bestfit))]
    return out
end

function getproperty(w::UI{FitResult}, s::Symbol)
    res = wrappee(w)
    if s in fieldnames(typeof(res))
        return getfield(res, s)
    end
    return UI{FitComp}(res.bestfit[s])
end

propertynames(w::UI{FitComp}) = collect(keys(wrappee(w).params))

getproperty(w::UI{FitComp}, s::Symbol) = wrappee(w).params[s]

# ____________________________________________________________________
# Miscellaneous
#
function test_component(comp::AbstractComponent, args...; iter=10)
    model = Model(:test => comp)
    add_dom!(model, args...)
    addexpr!(model, :test)

    println("Warmup:")
    @time result = evaluate!(model)

    if iter > 0
        println("Evaluating component $iter times:")
        @time begin
            for i in 1:iter
                result = evaluate!(model)
                wrappee(model).instruments[1].compevals[1].lastParams[1] = NaN  # Force re-calculation
            end
        end
    end
    return model
end
test_component(comp::AbstractComponent, domain::AbstractCartesianDomain, iter=1) =
    test_component(comp, flatten(domain); iter=iter)

probe(lparams::Parameter, data::AbstractMeasures; kw...) = probe([lparams], [data]; kw...)
probe(lparams::Vector{Parameter}, data::AbstractMeasures; kw...) = probe(lparams, [data]; kw...)
probe(lparams::Parameter, data::Vector{T}; kw...) where T <: AbstractMeasures = probe([lparams], data; kw...)

function probe(lparams::Vector{Parameter}, data::Vector{T}; delta=3., kw...) where T <: AbstractMeasures
    (length(lparams) == 0)  &&  (return nothing)
    (length(data) == 0)  &&  (return nothing)
    for par in lparams
        @assert par._private.model == lparams[1]._private.model
    end
    model = lparams[1]._private.model
    @assert model != nothing
    @assert length(model.instruments) >= 1
    params = collect(values(getparams(model)))
    pvalues0 = getfield.(params, :val)
    pvalues = deepcopy(pvalues0)

    data1d = data1D(model, data)
    ipar = Vector{Int}()
    range = Vector{Float64}()
    _evaluate!(model, pvalues0)
    cost0 = sum(abs2, residuals1d(model, data1d))
    for par in lparams
        for ii in 1:length(params)
            if isequal(params[ii], par)
                push!(ipar, ii)
                @assert isfinite(par._private.fitunc) "The model has not been fit: can't guess step for parameter"
                pvalues = deepcopy(pvalues0)
                lo = pvalues0[ii] - par._private.fitunc/5
                while true
                    if lo < par.low
                        lo = par.low
                        break
                    end
                    pvalues[ii] = lo
                    _evaluate!(model, pvalues)
                    cost = sum(abs2, residuals1d(model, data1d))
                    (abs(cost - cost0) >= delta)  &&  break
                    lo -= par._private.fitunc/5
                end
                hi = pvalues0[ii] + par._private.fitunc/5
                while true
                    if lo > par.high
                        hi = par.high
                        break
                    end
                    pvalues[ii] = hi
                    _evaluate!(model, pvalues)
                    cost = sum(abs2, residuals1d(model, data1d))
                    (abs(cost - cost0) >= delta)  &&  break
                    hi += par._private.fitunc/5
                end
                append!(range, [lo, hi])
            end
        end
    end
    @assert length(ipar) == length(lparams)
    _evaluate!(model, pvalues0)
    rr = collect(reshape(range, 2, div(length(range),2))')
    return probe(lparams, data, rr; kw...)
end

probe(lparams::Parameter, data::AbstractMeasures, rr::Matrix{Float64}; kw...) = probe([lparams], [data], rr; kw...)
probe(lparams::Vector{Parameter}, data::AbstractMeasures, rr::Matrix{Float64}; kw...) = probe(lparams, [data], rr; kw...)
probe(lparams::Parameter, data::Vector{T}, rr::Matrix{Float64}; kw...) where T <: AbstractMeasures = probe([lparams], data, rr; kw...)

function probe(lparams::Vector{Parameter}, data::Vector{T}, rr::Matrix{Float64}; nstep=11) where {T<:AbstractMeasures}
    (length(lparams) == 0)  &&  (return nothing)
    (length(data) == 0)  &&  (return nothing)
    for par in lparams
        @assert par._private.model == lparams[1]._private.model
    end
    model = lparams[1]._private.model
    @assert model != nothing
    @assert length(model.instruments) >= 1
    params = collect(values(getparams(model)))
    pvalues0 = getfield.(params, :val)
    pvalues = deepcopy(pvalues0)
    @assert size(rr)[1] == length(lparams)

    data1d = data1D(model, data)
    _evaluate!(model, pvalues0)
    cost0 = sum(abs2, residuals1d(model, data1d))

    ranges = Vector{AbstractArray}()
    lnstep = nstep
    if (length(nstep) == 1)  &&  (length(lparams) > 1)
        lnstep = fill(nstep[1], length(lparams))
    end
    @assert length(lnstep) == length(lparams)
    ipar = Vector{Int}()
    for par in lparams
        for ii in 1:length(params)
            if params[ii] == par
                push!(ipar, ii)
                push!(ranges, range(rr[length(ipar), 1], stop=rr[length(ipar), 2],
                                    length=lnstep[length(ipar)]))
            end
        end
    end
    @assert length(ipar) == length(lparams)
    if length(lparams) > 1
        cd = CartesianDomain(ranges...)
    else
        cd = Domain(ranges...)
    end
    dom = flatten(cd)

    out = Matrix{Float64}(undef, length(dom), length(lparams)+1)
    for ii in 1:length(dom)
        if length(lparams) > 1
            pvalues[ipar] .= dom.axis[:,ii]
            out[ii, 1:end-1] .= dom.axis[:,ii]
        else
            pvalues[ipar] .= dom.axis[ii]
            out[ii, 1] = dom.axis[ii]
        end
        _evaluate!(model, pvalues)
        out[ii, end] = sum(abs2, residuals1d(model, data1d)) .- cost0
    end
    _evaluate!(model, pvalues0)
    return out
end


text:
using Qsosed, Test

@testset "Test warm component" begin

    bh = BlackHole(1e8 * M_SUN, 0.5, 0.0)
    corona =
        Corona(bh, hard_xray_fraction = 0.02, electron_energy = 100, reprocessing = true)
    warm = Warm(corona, electron_energy = 0.2, photon_index = 2.5)
    @test compute_warm_photon_index(warm) == 2.5
    @test warm.electron_energy == 0.2
    @test warm.radius ≈ 17.801792397994166 rtol = 1e-2

    temp_warm = compute_disk_temperature(bh, warm.radius)
    @test temp_warm ≈ 87745.726703792374 rtol = 1e-1
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

abstract AbstractSparseArray{Tv,Ti,N} <: AbstractArray{Tv,N}

typealias AbstractSparseVector{Tv,Ti} AbstractSparseArray{Tv,Ti,1}
typealias AbstractSparseMatrix{Tv,Ti} AbstractSparseArray{Tv,Ti,2}

"""
    issparse(S)

Returns `true` if `S` is sparse, and `false` otherwise.
"""
issparse(A::AbstractArray) = false
issparse(S::AbstractSparseArray) = true

issparse{T, A<:AbstractSparseMatrix}(S::Symmetric{T, A}) = true
issparse{T, A<:AbstractSparseMatrix}(S::Hermitian{T, A}) = true
issparse{T, A<:AbstractSparseMatrix}(S::LowerTriangular{T, A}) = true
issparse{T, A<:AbstractSparseMatrix}(S::LinAlg.UnitLowerTriangular{T, A}) = true
issparse{T, A<:AbstractSparseMatrix}(S::UpperTriangular{T, A}) = true
issparse{T, A<:AbstractSparseMatrix}(S::LinAlg.UnitUpperTriangular{T, A}) = true

indtype{Tv,Ti}(S::AbstractSparseArray{Tv,Ti}) = Ti

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
using Plots, Distributions, LinearAlgebra, ForwardDiff

n = 100
m = 50

Q = Diagonal(ones(n))

A = rand(m, n)
b = A * randn(n)
z_sol = 1.0 * [2 * Q A'; A zeros(m, m)] \ [zeros(n); b]
# @show norm(lagrangian_gradient(z_sol))
# @show lagrangian(z_sol)

function obj(x)
	x' * Q * x
end

function gradient(x)
	2.0 * Q * x
end

function lagrangian(z)
	x = view(z, 1:n)
	y = view(z, n .+ (1:m))

	obj(x) + y' * (A * x - b)
end

function lagrangian_gradient(z)
	x = view(z, 1:n)
	y = view(z, n .+ (1:m))

	[gradient(x) + A' * y; A * x - b]
end

dir = Diagonal([ones(n); -1.0 * ones(m)])

# adam
function adam(obj, gradient, x;
		max_iter = 1000)

	# history
	J_hist = [obj(x)]

	# parameters
	α = 1.0e-4
	β1 = 0.9
	β2 = 0.999
	ϵ = 10.0^(-8.0)
	m = zero(x)
	v = zero(x)
	m̂ = zero(x)
	v̂ = zero(x)

	# iterate
	for i = 1:max_iter
		g = gradient(x)
		# m .= β1 .* m + (1.0 - β1) .* g
		# v .= β2 .* v + (1.0 - β2) .* (g.^2.0)
		# m̂ .= m ./ (1.0 - β1^i)
		# v̂ .= v ./ (1.0 - β2^i)
		# x .-= α * m̂ ./ (sqrt.(v̂) .+ ϵ)
		x .-= α * dir * g
		push!(J_hist, obj(x))
		println("norm $(norm(gradient(x)))")
		i % 1000 == 0 && println("iter: $i")
	end

	return x, J_hist
end

x_sol, J_hist = adam(lagrangian, lagrangian_gradient, ones(n + m),
	max_iter = 100000)

plt = plot(J_hist,
	xlabel = "iteration", ylabel = "obj.", title = "sgd",
	# label = "α = $α (adam)",
	color = :red, width = 2.0)

display(plt)


text:
using Plots, JLD, ProgressBars, Statistics
Directions=[[0,1],[0,-1],[1,0],[-1,0]]
function f(N, Network_,Pos)
    if N==0
        return 1
    end
    Network_[Pos...]=true
    c=0
    for path in Directions
        NextPos=Pos .+ path
        if !Network_[NextPos...]
            c+=f(N-1, Network_, NextPos)
        end
    end
    Network_[Pos...]=false
    return c
end
Totdata=[]
N=16
for n in ProgressBar(0:N)
    L=2n+1
    Network_=falses(L,L)
    FirstPos=[n+1, n+1]
    push!(Totdata, f(n, Network_,FirstPos))
end
plot(0:N,Totdata, legend=false, dpi=400)
scatter!(0:N,Totdata, legend=false, xlabel="N",ylabel="Number of Paths",
    title="Number of Paths versus N for Self-Avoiding Walker", titlefontsize=10, dpi=400)
save("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter5\\4.7\\SAW.jld", "Data", Totdata)
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter5\\4.7\\Figs\\NumberOfPaths_N.png")
plot(0:N, Totdata./ 4 .^(0:N), legend=false, dpi=400)
scatter!(0:N, Totdata./ 4 .^(0:N), legend=false, xlabel="N",ylabel="Ratio of Number of Paths",
    title="Ratio of Number of Paths of SAW and RW ", dpi=400)
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter5\\4.7\\Figs\\Ratio_N.png")


text:
"""
    estimateKernelSize{T <: Real}(X::Array{T}; sampleSize = 100)

    Estimate kernel size by taking median distance between points.
"""
function estimateKernelSize(X::Array{T}; sampleSize = 100) where T <: Real

	M = length(X)

	# set kernel size to median distance between points
	if M > sampleSize
		Xmed = reshape(sample(X, sampleSize, replace=false), (sampleSize, 1))
		S = sampleSize
	else
		Xmed = reshape(X, (M, 1))
		S = M
	end

	dists = pairwise(SqEuclidean(), Xmed, Xmed, dims = 1)
    sig = sqrt(0.5 * median(dists))

    return sig
end

"""
    rbfDotProduct(X::Array, X::Array, kernelSize::Float64)

"""
function rbfDotProduct(X::Array, Y::Array, kernelSize)

	G = sum((X.*Y), dims = 2)

	Q = repeat(G, 1, length(Y))
	R = repeat(G', length(X), 1)

	H = Q + R - 2*X*Y'

	return exp.(-H/2/kernelSize^2)

end

eye(n) = Matrix(I, n, n)


text:
using DataStructure101
const DS = DataStructure101

export find_entrances, topological_sort

function calc_in_degree(g::DS.AbstractGraph)
    in_degrees = zeros(DS.nv(g))

    for v = 1:DS.nv(g)
        for n in DS.neighbor(g, v)
            in_degrees[n] += 1
        end # for
    end # for

    return in_degrees
end

function find_entrances(g::DS.AbstractGraph)
    entrances = Int[]

    in_degrees = calc_in_degree(g)
    for (v, in_degree) in enumerate(in_degrees)
        in_degree == 0 && push!(entrances, v)
    end # for

    return entrances
end # function

function topological_sort(g::DS.AbstractGraph)
    sorted = Int[]
    entrances = find_entrances(g)

    while length(entrances) > 0
        v = popfirst!(entrances)
        push!(sorted, v)

        ns = DS.neighbor(g, v)
        for n in ns
            if !(n in sorted || n in entrances) push!(entrances, n) end
        end
    end # while

    return sorted
end


text:
# Test features specific to PermMat

m = PermMat([2, 1, 4, 5, 3])
m3 = m * 3
@test size(m3) == (5,5)
@test sum(m3) == 15

p = list(m)
@test p * 3 == 4
@test m * "cat" == p * "cat" == "act"
mf = full(m)
@test p * mf == m * mf

@test p[1] == 2
@test m[1] == 0

@test length(m) == 25
@test size(m) == (5,5)
# Can't do this with the other permutation types
@test repmat(m,3,4) == repmat(mf,3,4)
@test m + eye(Int,5) == mf + eye(Int,5)
v = [1,2,3,4,5]
@test broadcast(+,v,m) == broadcast(+,v,mf)
@test broadcast(-,v,m) == broadcast(-,v,mf)
@test m + m == mf + mf
e1 = [1,2,3]
@test m[e1] == Int[0,1,0]
@test m[1:3] == Int[0,1,0]
v1 = [.2,.3,.1,.7,.9]
@test m * v1 == mf * v1

@test length([i for i in m]) == 25
@test length([i for i in p]) == 5

# wikipedia example
p1 = PermList([1,4,2,5,3])
m1 = matrix(p1)
@test m1 * v == p1.data
@test full(m1) * v == p1.data
@test sparse(m1) * v == p1.data


text:
abstract type AbstractSecondaryEdgeCutModel end

@inline (mdl::AbstractSecondaryEdgeCutModel)(n_vertices, n_self_pins, n_cut_pins, k) = mdl(n_vertices, n_self_pins, n_cut_pins)

struct AffineSecondaryEdgeCutModel{Tv} <: AbstractSecondaryEdgeCutModel
    α::Tv
    β_vertex::Tv
    β_self_pin::Tv
    β_cut_pin::Tv
end

function AffineSecondaryEdgeCutModel(; α = false, β_vertex = false, β_self_pin = false, β_cut_pin = false)
    AffineSecondaryEdgeCutModel(promote(α, β_vertex, β_self_pin, β_cut_pin)...)
end

@inline cost_type(::Type{AffineSecondaryEdgeCutModel{Tv}}) where {Tv} = Tv

(mdl::AffineSecondaryEdgeCutModel)(n_vertices, n_self_pins, n_cut_pins, k) = mdl.α + n_vertices * mdl.β_vertex + n_self_pins * mdl.β_self_pin + n_cut_pins * mdl.β_cut_pin

function bound_stripe(A::SparseMatrixCSC, K, Π, mdl::AffineSecondaryEdgeCutModel)
    @assert mdl.β_vertex >= 0
    @assert mdl.β_self_pin >= 0
    @assert mdl.β_cut_pin >= 0
    adj_A = adjointpattern(A)
    c_hi = bottleneck_value(adj_A, Π, AffineWorkModel(mdl.α, mdl.β_vertex, mdl.β_cut_pin))
    c_lo = bottleneck_value(adj_A, Π, AffineWorkModel(mdl.α, mdl.β_vertex, mdl.β_self_pin))
    return (c_lo, c_hi)
end

struct SecondaryEdgeCutOracle{Ti, Lcv, Mdl} <: AbstractOracleCost{Mdl}
    Π::SplitPartition{Ti}
    pos::Vector{Ti}
    πos::Vector{Ti}
    lcv::Lcv
    mdl::Mdl
end

oracle_model(ocl::SecondaryEdgeCutOracle) = ocl.mdl

function bound_stripe(A::SparseMatrixCSC, K, Π::SplitPartition, ocl::SecondaryEdgeCutOracle{<:Any, <:Any, <:AffineSecondaryEdgeCutModel})
    @inbounds begin
        c_lo = 0
        c_hi = 0
        (m, n) = size(A)
        mdl = oracle_model(ocl)
        @assert mdl.β_vertex >= 0
        @assert mdl.β_self_pin >= 0
        @assert mdl.β_cut_pin >= 0
        for k = 1:K
            n_vertices = Π.spl[k + 1] - Π.spl[k]
            n_pins = ocl.pos[Π.spl[k + 1]] - ocl.pos[Π.spl[k]]
            c_lo = max(c_lo, mdl.α + n_vertices * mdl.β_vertex + n_pins * min(mdl.β_self_pin, mdl.β_cut_pin))
            c_hi = max(c_hi, mdl.α + n_vertices * mdl.β_vertex + n_pins * max(mdl.β_self_pin, mdl.β_cut_pin))
        end
        return (c_lo, c_hi)
    end
end

function oracle_stripe(hint::AbstractHint, mdl::AbstractSecondaryEdgeCutModel, A::SparseMatrixCSC, Π::SplitPartition; adj_pos=nothing, adj_A=nothing, pin=nothing, kwargs...)
    @inbounds begin
        m, n = size(A)

        if adj_A === nothing
            adj_A = adjointpattern(A)
        end

        if adj_pos === nothing
            adj_pos = adj_A.colptr
        end

        (n, K, n′, πos, prm), Ap = partwise(A, convert(MapPartition, Π))
        lcv = partwisecount!(hint, n, K, n′, πos, prm, pincount(hint, Ap); kwargs...)

        return SecondaryEdgeCutOracle(Π, adj_pos, πos, lcv, mdl)
    end
end

@inline function (cst::SecondaryEdgeCutOracle{Ti, Mdl})(i::Ti, i′::Ti, k) where {Ti, Mdl}
    @inbounds begin
        w = cst.pos[cst.Π.spl[k + 1]] - cst.pos[cst.Π.spl[k]]
        d = cst.πos[k + 1] - cst.πos[k]
        l = cst.lcv(i, i′, k)
        return cst.mdl(cst.Π.spl[k + 1] - cst.Π.spl[k], l, w - l, k)
    end
end

@inline function (stp::Step{Ocl})(_j, _j′, _k) where {Ti, Mdl, Ocl <: SecondaryEdgeCutOracle{Ti, Mdl}}
    @inbounds begin
        cst = stp.ocl
        j = destep(_j)
        j′ = destep(_j′)
        k = destep(_k)
        w = cst.pos[cst.Π.spl[k + 1]] - cst.pos[cst.Π.spl[k]]
        d = cst.πos[k + 1] - cst.πos[k]
        l = Step(cst.lcv)(_j, _j′, _k)
        return cst.mdl(cst.Π.spl[k + 1] - cst.Π.spl[k], l, w - l, k)
    end
end


text:
"""
  edge

Edge of bipartite graph with star expansion of a hypergraph.

**Arguments**

* from : Vertex number that the edge connects to.  It is a number representing a vertex of the hypergraph.
* to : Vertex number that the edge connects to. It is a number representing a hyperedge of the hypergraph.
* weight : Edge weight.
* id : The number of edges that this structure represents.
"""
mutable struct edge
  from::Int
  to::Int
  weight::Float64
  id::Int
end

"""
  edge_comp(a::edge, b::edge)

Edge comparison function.
Sort in descending order of weight.
If the weights are the same, they are compared according to the IDs of the edges.
"""
function edge_comp(a::edge, b::edge)
  if a.weight == b.weight
    return a.id < b.id
  end
  return a.weight < b.weight
end

"""
partition(uf::UnionFind{Int},
          to::Int=length(uf.parent),
          from::Int=1,
          )::Vector{Set{Int}}

Compute the clustering result from the disjoint set.

**Arguments**

* uf : The disjoint set obtained by clustering.
* to : Disjoint sets up to `to` are used. For example, if you want to get a vertex-only clustering result,
specify the number of vertices.
* from : Use a disjoint set from `from`.
"""
function partition(uf::UnionFind{Int64},
                   to::Int=length(uf.parent),
                   from::Int=1
                  )::Vector{Set{Int}}
  d = Dict()
  for (i, elem) in (enumerate(uf.parent[from:to]))
    cluster_num = root(uf, from-1+i)
    d[cluster_num] = get(d, cluster_num, Set())
    push!(d[cluster_num], from-1+i)
  end
  d = [(length(v), Set{Int64}(v)) for (k, v) in d]
  sort!(d, by=x->x[1], rev=true)
  d = Vector([s for (l, s) in d])
  return d
end

"""
  my_mod(h::Hypergraph, part::Array{Set{any}})

Computes the modularity when `h` is divided by `part`.
"""
function my_mod(h::Hypergraph, part)
  ml = 0
  mr = 0
  v = 0
  for i in 1:nhv(h)
    v += length(gethyperedges(h, i))
  end
  dict = Dict([])
  for i in 1:nhe(h)
    dict[length(getvertices(h, i))] = get(dict, length(getvertices(h, i)), 0) + 1
  end

  # vav = 0
  va = Dict([i => 0.0 for i in 1:length(part)])
  for (i, cluster) in enumerate(part)
    for node in cluster
      va[i] += length(gethyperedges(h, node))
      # vav += length(gethyperedges(h, node))
    end
    va[i] /= v
  end

  for (key, val) in dict
    tmp = 0
    for (key2, val2) in va
      tmp += (val2 ^ key)
    end
    tmp *= val
    mr += tmp
  end

  n2c = [-1 for i in 1:nhv(h)]
  for (num, cluster) in enumerate(part)
    for node in cluster
      n2c[node] = num
    end
  end

  for he_i in 1:nhe(h)
    he = [k for (k, v) in getvertices(h, he_i)]
    if length(he) == 0
      error("There are hypernodes that are not included in the hyperedge.")
      exit(1)
    end
    flag = true
    fnode_cluster = n2c[he[1]]
    for node in he
      if n2c[node] != fnode_cluster
        flag = false
        break
      end
    end
    if flag ml += 1 end
    # if issubset(he, cluster) ml += 1 end
  end

  (ml - mr) / nhe(h)
end

"""
  star_expansion(h::Hypergraph, weighted_f::Function=tfidf, param=Dict())

Star expandion `h` and construct a bipartite graph. Each edge is weighted by `weighted_f`.
"""
function star_expansion(h::Hypergraph, weighted_f=tfidf, params=Dict())::Array{edge}
  edges::Array{edge} = []
  params["dl_ave"] = 0

  for he in 1:nhe(h)
    params["dl_ave"] += length(getvertices(h, he))
  end
  params["dl_ave"] /= nhe(h)

  id = 1
  @showprogress 1 "computing..." for node in (1:nhv(h))
    hes = gethyperedges(h, node)
    hes = [key for (key, val) in hes]
    for he in hes
      w = weighted_f(h, he, node, params)
      # edgeはnode, edge, tfidfの構造体
      push!(edges, edge(node, nhv(h)+he, w, id))
      id += 1
    end
  end
  # edges = [i for i in edges]
  sort!(edges, rev=true, lt=edge_comp)
  return edges
end

# 何番目にnoise nodeが追加されたか
function noise_order(edges, noise_indexes)
  orders = []
  for (i, e) in enumerate(edges)
    node = e.from
    he = e.to
    if noise_indexes[node][he] push!(orders, i) end
  end

  return orders
end

function plot_incidence(h::Hypergraph, name="", weighted_f=tfidf, params=Dict())
  pyplot()
  bg = star_expansion(h, weighted_f, params)
  plot_arr::Array{Tuple{Int64, Int64}} = [(i.to-nhv(h), i.from) for i in bg]
  maxw = maximum([i.weight for i in bg])
  minw = minimum([i.weight for i in bg])
  f(w) = 1-sqrt((w-minw) / (maxw-minw))
  color = [RGB(f(i.weight), f(i.weight), f(i.weight)) for i in bg]
  # color = [RGB(0, f(i.weight), 1.) for i in bg]
  return Plots.scatter(
                plot_arr,
                markersize=6,
                markerstrokewidth=0,
                color=color,
                xlabel="Hyperedges",
                ylabel="Node",
                label="",
                title=name,
               )
end

function plot_incidence2(h::Hypergraph, arr)
  # pyplot()
  # gr()
  bg = arr
  # bg = build_bg(h, weighted_f, params)
  plot_arr::Array{Tuple{Int64, Int64}} = [(i.to-nhv(h), i.from) for i in bg]
  maxw = length(bg)
  minw = 1
  f(w) = sqrt((w-minw) / (maxw-minw))
  color = [RGB(f(i), f(i), f(i)) for (i, e) in enumerate(bg)]
  # println(color)
  # color = [RGB(0, f(i.weight), 1.) for i in bg]
  return Plots.scatter(
                plot_arr,
                markersize=6,
                markerstrokewidth=0,
                color=color,
                xlabel="hyperedges",
                ylabel="node",
                label="",
               )
end

function curve(h::Hypergraph, f1, f2, similar_f=f1_score, params1=Dict(), params2=Dict())
  bg1 = star_expansion(h, f1, params1)
  bg2 = star_expansion(h, f2, params2)
  # bg1 = map(i->i.id, bg1[1:Int(floor(length(bg1)/100)):end])
  # bg2 = map(i->i.id, bg2[1:Int(floor(length(bg1)/100)):end])
  bg1 = map(i->i.id, bg1)
  bg2 = map(i->i.id, bg2)
  n = length(bg1)

  # arr1 = [0 for i in 1:n]
  # println(length(arr1[bg1]))
  res = []
  @showprogress 1 "computing..."  for i in 1:n
    arr1 = [0 for i in 1:n]
    arr2 = [0 for i in 1:n]
    arr1[bg1[1:i]] .= 1
    arr2[bg2[1:i]] .= 1
    # arr2 = bg2[1:i]
    push!(res, similar_f(arr1, arr2))
  end

  return res
end

function disp_basic(h::Hypergraph)
  node_degree_dist = [0 for i in 1:nhe(h)]
  he_degree_dist = [0 for i in 1:nhv(h)]
  for node in 1:nhv(h)
    size = length(gethyperedges(h, node))
    node_degree_dist[size] += 1
  end
  for he in 1:nhe(h)
    size = length(getvertices(h, he))
    he_degree_dist[size] += 1
  end
  # println(node_degree_dist)
  # println(he_degree_dist)
  p = Plots.bar(node_degree_dist, show=true)
  display(p)
  p = Plots.bar(he_degree_dist, show=true)
  display(p)
end

"""
  epart2cluster(h::Hypergraph, epart::Array{Set{Int}})

"""
function epart2cluster(h::Hypergraph, epart)
  part::Array{Set{Int64}} = []
  @showprogress for ecluster in epart
    p = Set([])
    for he in ecluster
      p = union(p, keys(getvertices(h, he-nhv(h))))
    end
    push!(part, p)
  end
  return part
end

function epart2probdist(h::Hypergraph, uf, p; weighted_f=tfidf, params=Dict())
  cluster_set = Set()
  for i in nhv(h)+1:nhv(h)+nhe(h) push!(cluster_set, root(uf, i)) end
  cluster_nums = Dict([val => i for (i, val) in enumerate(cluster_set)])
  entropies = Dict([i => [.0 for j in 1:length(cluster_nums)] for i in 1:nhv(h)])
  bg = star_expansion(h, weighted_f)
  @showprogress 1 "computing.." for e in bg
    cluster_num = cluster_nums[root(uf, e.to)]
    entropies[e.from][cluster_num] += e.weight
  end
  for e in bg
    entropies[e.from] ./= sum(entropies[e.from])
  end
  return entropies
end

function avg_entropy(prob_dist)
  avg = 0
  for prob in prob_dist
    if prob == 0 continue end
    avg += prob * log(length(prob_dist), prob)
  end
  avg *= -1
  if avg == -0.0 return 0
  else return avg end
end

function prob_dist2part(prob_dists, threshold=0.5)
  part::Array{Set{Int64}} = [Set() for i in 1:length(prob_dists[1])]
  @showprogress 1 "computing..." for (node, prob_dist) in prob_dists
    for (i, prob) in enumerate(prob_dist)
      if prob > threshold push!(part[i], node) end
    end
  end
  return part
end

# function cluster_visualization(h::Hypergraph, p)
#   colors = ["red" "blue" "yellow" "green" "brown" "cyan" "magenta" "pink" "white" "gray" ]
#   ec = ["gray" for i in nhe(h)]
#   c = ["" for i in 1:nhv(h)]
#   sort!(p, rev=true, by=i->length(i))
#
#   for (i, cluster) in enumerate(p)
#     for node in cluster
#       c[node] = length(cluster) == 1 ? "black" : colors[i]
#     end
#   end
#
#   SimpleHypergraphs.draw(h,
#                          HyperNetX,
#                          nodes_kwargs=Dict(["facecolors"=>c]),
#                          edges_kwargs=Dict(["edgecolors"=>ec]),
#                          layout_kwargs=Dict(["seed"=>0]),
#                          with_node_labels=false,
#                          with_edge_labels=false,
#                         )
#
#   return c
# end

function disp_cluster_bias(p)
  Plots.pie(length.(p), label="")
end

function h2correlation(h::Hypergraph, f1, f2)
  okapi_e = star_expansion(h, f1)
  tfidf_e = star_expansion(h, f2)

  rank1 = Dict([i => Dict([j => 0 for j in 1:nhe(h)]) for i in 1:nhv(h)])
  rank2 = Dict([i => Dict([j => 0 for j in 1:nhe(h)]) for i in 1:nhv(h)])
  for (i, e) in enumerate(okapi_e)
    rank1[e.from][e.to-nhv(h)] = i
  end

  for (i, e) in enumerate(tfidf_e)
    rank2[e.from][e.to-nhv(h)] = i
  end

  arr1 = Array{Int64}([])
  arr2 = Array{Int64}([])
  @showprogress 1 "computing..." for node in 1:nhv(h)
    for he in 1:nhe(h)
      if rank1[node][he] != 0 push!(arr1, rank1[node][he]::Int64) end
      if rank2[node][he] != 0 push!(arr2, rank2[node][he]::Int64) end
    end
  end

  return (arr1, arr2)
end

# Or later, feature functoins
function gini(h::Hypergraph, cluster_dict)
  k = 0
  sum = 0
  for (key, val) in cluster_dict
    k += val
    sum += key * val
  end
  if k == 1 return (1, 1) end

  gini = 0

  for (i, val1) in cluster_dict
    for (j, val2) in cluster_dict
      gini += abs(i - j) * (val1 * val2)
    end
  end
  gini /= 2

  gini = (gini / binomial(k, 2)) / (sum / k)
  gini /= 2

  return (gini, k)
end

function gini2(h::Hypergraph, cluster_dict)
  k = 0
  sum = 0

  for (key, val) in cluster_dict
    if(key == 1) continue end
    k += val
    sum += key * val
  end
  if k <= 1 return (1, 1) end

  gini = 0

  for (i, val1) in cluster_dict
    for (j, val2) in cluster_dict
      if(i == 1 || j == 1) continue end
      gini += abs(i - j) * (val1 * val2)
    end
  end
  gini /= 2

  gini = (gini / binomial(k, 2)) / (sum / k)
  gini /= 2

  return (gini, k)
end
# function testmod(h::Hypergraph, method::CFModularityCNMLike, fx)
#     ha = HypergraphAggs(h)
#     best_modularity = 0
#     comms = [Set(i) for i in 1:nhv(h)]
#     mod_history = Vector{Float64}(undef, method.reps)
#     for rep in 1:method.reps
#         he = rand(1:nhe(h))
#         vers = collect(keys(getvertices(h, he)))
#         if length(vers) == 0
#             continue
#         end;
#         c = deepcopy(comms)
#         i0 = fx(c, vers)
#         max_i = length(c)
#         i_cur = i0
#         while i_cur < max_i
#             i_cur += 1
#             if length(intersect(c[i_cur],vers)) > 0
#                 union!(c[i0], c[i_cur])
#                 c[i_cur]=c[max_i]
#                 max_i += -1
#             end
#                 println(c)
#         end
#         resize!(c,max_i)
#         m = my_mod(h, c)
#         if m > best_modularity
#             best_modularity = m
#             comms = c
#         end
#         mod_history[rep] = best_modularity
#     end
#     return (bm=best_modularity, bp=comms, mod_history=mod_history)
# end


text:
module NLSat

export NLSolver, Z3Solver, SMTSolverZ3, SMTSolverYices, SMTSolver, CFiniteSolver
export NLStatus, NLModel
export add_vars!, add!, add_soft!, solve

using DataStructures
using DelimitedFiles
using Distributed
using MacroTools: walk, postwalk, @capture, @match, replace
using Dates
using Z3

const NLModel = Dict{Symbol,<:Number}

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

abstract type AlgebraicNumber end
export AlgebraicNumber

function AlgebraicNumber(x::Expr, n::Int)
    @warn "Algebraic numbers not yet supported, got $((x, n)), returning $(NaN)"
    NaN
end

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

include("../utils.jl")
include("clauseset.jl")
include("smtlib.jl")
include("lisp.jl")
include("z3utils.jl")

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

@enum NLStatus sat unsat unknown timeout

abstract type NLSolver end

function add_vars!(s::NLSolver, d::Dict{Symbol,Type}) end
function add!(s::NLSolver, cs::ClauseSet) end
function add_soft!(s::NLSolver, cs::ClauseSet) end
function solve(s::NLSolver; timeout::Int = -1) end

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

# include("mathematica.jl")
include("z3solver.jl")
include("smtsolver.jl")
include("cfinitesolver.jl")

end # module


text:
# ===================================
# Homogeneous case without invariant
# ===================================

# note that each of the computations in the loop over directions is independent
function reach_homog_LGG09!(F::Vector{RT},
                            dirs::TN,
                            Ω₀::LazySet,
                            Φ::AbstractMatrix,
                            NSTEPS::Integer,
                            δ::N,
                            X::Universe, # no invariant
                            Δt0::TimeInterval,
                            cache,
                            threaded) where {N, VN, TN, SN, RT<:TemplateReachSet{N, VN, TN, SN}}

    # transpose coefficients matrix
    Φᵀ = copy(transpose(Φ))

    # preallocate output sequence
    ρℓ = Matrix{N}(undef, length(dirs), NSTEPS)

    _reach_homog_dir_LGG09!(ρℓ, Ω₀, Φᵀ, dirs, NSTEPS, cache, threaded)

    # fill template reach-set sequence
    Δt = (zero(N) .. δ) + Δt0
    @inbounds for k in 1:NSTEPS
        F[k] = TemplateReachSet(dirs, view(ρℓ, :, k), Δt)
        Δt += δ
    end

    return ρℓ
end

function _reach_homog_dir_LGG09!(ρℓ, Ω₀, Φᵀ, dirs, NSTEPS, cache, threaded::Val{true})
    ℓ = _collect(dirs)
    Threads.@threads for j in 1:length(ℓ)
        reach_homog_dir_LGG09!(ρℓ, j, Ω₀, Φᵀ, ℓ[j], NSTEPS, cache)
    end
end

function _reach_homog_dir_LGG09!(ρℓ, Ω₀, Φᵀ, dirs, NSTEPS, cache, threaded::Val{false})
    #for each direction, compute NSTEPS iterations
    @inbounds for (j, ℓ) in enumerate(dirs)
        reach_homog_dir_LGG09!(ρℓ, j, Ω₀, Φᵀ, ℓ, NSTEPS, cache)
    end
end

function reach_homog_dir_LGG09!(ρvec_ℓ::AbstractMatrix{N}, j, Ω₀, Φᵀ, ℓ, NSTEPS, cache::Val{true}) where {N}
    rᵢ = copy(ℓ)
    rᵢ₊₁ = similar(rᵢ)

    @inbounds for i in 1:NSTEPS
        ρvec_ℓ[j, i] = ρ(rᵢ, Ω₀)

        # update cache for the next iteration
        mul!(rᵢ₊₁, Φᵀ, rᵢ)
        copy!(rᵢ, rᵢ₊₁)   # rᵢ .= rᵢ₊₁
    end
    return ρvec_ℓ
end

function reach_homog_dir_LGG09!(ρvec_ℓ::AbstractMatrix{N}, j, Ω₀, Φᵀ, ℓ, NSTEPS, cache::Val{false}) where {N}
    rᵢ = copy(ℓ)

    @inbounds for i in 1:NSTEPS
        ρvec_ℓ[j, i] = ρ(rᵢ, Ω₀)

        # update cache for the next iteration
        rᵢ = Φᵀ * rᵢ
    end
    return ρvec_ℓ
end

# TODO: needs specialization for static vector / static matrix ?
# compute NSTEPS iterations support function along direction ℓ
# ``ρ(ℓ, Ω₀), ρ(ℓ, Φ * Ω₀), ρ(ℓ, Φ^2 * Ω₀), ..., ρ(ℓ, Φ^NSTEPS * Ω₀)``
function reach_homog_dir_LGG09!(ρvec_ℓ::AbstractVector{N}, Ω₀, Φᵀ, ℓ, NSTEPS, cache::Val{true}) where {N}
    rᵢ = copy(ℓ)
    rᵢ₊₁ = similar(rᵢ)

    @inbounds for i in 1:NSTEPS
        ρvec_ℓ[i] = ρ(rᵢ, Ω₀)

        # update cache for the next iteration
        mul!(rᵢ₊₁, Φᵀ, rᵢ)
        copy!(rᵢ, rᵢ₊₁)   # rᵢ .= rᵢ₊₁
    end
    return ρvec_ℓ
end

# ``ρ(ℓ, Ω₀), ρ(ℓ, Φ * Ω₀), ρ(ℓ, Φ^2 * Ω₀), ..., ρ(ℓ, Φ^NSTEPS * Ω₀)``
function reach_homog_dir_LGG09!(ρvec_ℓ::AbstractVector{N}, Ω₀, Φᵀ, ℓ, NSTEPS, cache::Val{false}) where {N}
    rᵢ = copy(ℓ)
    @inbounds for i in 1:NSTEPS
        ρvec_ℓ[i] = ρ(rᵢ, Ω₀)

        # update cache for the next iteration
        rᵢ = Φᵀ * rᵢ
    end
    return ρvec_ℓ
end

# NOTES:
# - if ℓ and -ℓ are in the template, they could share (Φᵀ)^k * ℓ
# - the loop over directions can be parallelized

#=---- version using vector-of-vectors
    # preallocate output sequence
    ρℓ = [Vector{N}(undef, NSTEPS) for _ in 1:length(dirs)]

    # for each direction, compute NSTEPS iterations
    @inbounds for (i, ℓ) in enumerate(dirs)
        reach_homog_dir_LGG09!(ρℓ[i], Ω₀, Φᵀ, ℓ, NSTEPS)
    end

    # fill template reach-set sequence
    Δt = (zero(N) .. δ) + time_shift
    @inbounds for k in 1:NSTEPS
        sf = [ρℓ[i][k] for i in 1:length(dirs)]
        F[k] = TemplateReachSet(dirs, sf, Δt)
        Δt += δ
    end
=#

# ------------------------------------------------------------
# Functionality that requires ExponentialUtilities.jl
# ------------------------------------------------------------

function load_krylov_LGG09_homog()
return quote

"""
    reach_homog_krylov_LGG09!(out, Ω₀::LazySet, Aᵀδ::AbstractMatrix,
                              ℓ::AbstractVector, NSTEPS;
                              hermitian=false, m=min(30, size(Aᵀδ, 1)), tol=1e-7)

### Algorithm

We compute the sequence:

```math
    ρ(ℓ, Ω₀), ρ(ℓ, Φ Ω₀), ρ(ℓ, Φ^2 Ω₀), ρ(ℓ, Φ^3 Ω₀), ...
```

Using Krylov subspace approximations to compute the action of Φ := exp(Aδ) over
the direction ℓ.

The method is (see [1]):

```julia
out[1] <- ρ(ℓ, Ω₀)

out[2] <- ρ(ℓ, Φ Ω₀) = ρ(Φᵀ ℓ, Ω₀)

out[3] <- ρ(ℓ, Φ^2 Ω₀) = ρ((Φᵀ)^2 ℓ, Ω₀)

out[4] <- ρ(ℓ, Φ^3 Ω₀) = ρ((Φᵀ)^3 ℓ, Ω₀)
```
 and so on.

### References

[1] Reach Set Approximation through Decomposition with Low-dimensional Sets and
    High-dimensional Matrices. Sergiy Bogomolov, Marcelo Forets, Goran Frehse,
    Frédéric Viry, Andreas Podelski and Christian Schilling (2018) HSCC'18
    Proceedings of the 21st International Conference on Hybrid Systems: Computation
    and Control: 41–50.
"""
function reach_homog_krylov_LGG09!(out, Ω₀::LazySet, Aᵀδ::AbstractMatrix,
                                   ℓ::AbstractVector, NSTEPS;
                                   hermitian=false, m=min(30, size(Aᵀδ, 1)), tol=1e-7)

    # initialization of the krylov subspace
    TA, Tb = eltype(Aᵀδ), eltype(ℓ)
    T = promote_type(TA, Tb)
    Ks = KrylovSubspace{T, real(T)}(length(ℓ), m)
    arnoldi!(Ks, Aᵀδ, ℓ; m=m, ishermitian=hermitian, tol=tol)

    # rᵢ stores is the cache for each vector: (Φᵀ)^i ℓ
    rᵢ = deepcopy(ℓ)

    @inbounds for i in 1:NSTEPS
        out[i] = ρ(rᵢ, Ω₀)
        expv!(rᵢ, i*1.0, Ks)
    end
    return out
end

end end  # quote / load_krylov_LGG09_homog()

# ------------------------------------------------------------
# Methods using eigenvalues of the transition matrix
# ------------------------------------------------------------

# it is assumed that (λ, d) is an eigenvalue-eigenvector pair of the matrix Φᵀ
function reach_homog_dir_eig_LGG09!(out::AbstractVector{N}, X₀, d::AbstractVector{N}, λ::N, NSTEPS) where {N}
    if iszero(λ)
        _reach_homog_dir_eig_LGG09_zero!(out, X₀, d, NSTEPS)
    elseif λ > zero(N)
        _reach_homog_dir_eig_LGG09_positive!(out, X₀, d, λ, NSTEPS)
    else
        _reach_homog_dir_eig_LGG09_negative!(out, X₀, d, λ, NSTEPS)
    end
    return out
end

function _reach_homog_dir_eig_LGG09_zero!(out::AbstractVector{N}, X₀, d, NSTEPS) where {N}
    @inbounds begin
        out[1] = ρ(d, X₀)
        for i in 2:NSTEPS
            out[i] = zero(N)
        end
    end
end

function _reach_homog_dir_eig_LGG09_positive!(out::AbstractVector{N}, X₀, d, λ, NSTEPS) where {N}
    ρ₀ = ρ(d, X₀)
    λⁱ = one(N)
    @inbounds for i in 1:NSTEPS
        out[i] = λⁱ * ρ₀
        λⁱ = λⁱ * λ
    end
end

function _reach_homog_dir_eig_LGG09_negative!(out::AbstractVector{N}, X₀, d, λ, NSTEPS) where {N}
    ρ₀ = ρ(d, X₀)
    ρ₀₋ = ρ(-d, X₀)
    λⁱ = λ
    @inbounds begin
        out[1] = ρ₀
        for i in 2:NSTEPS
            if iseven(i)
                out[i] = -λⁱ * ρ₀₋
            else
                out[i] = λⁱ * ρ₀
            end
            λⁱ = λⁱ * λ
        end
    end
end

# Let A be an `n x n` matrix with real eigenvalues.
# Given the eigendecomposition: Aδ = P * Λ * P^{-1},
# where Λ is a diagonal matrix that contains the eigenvalues of Aδ
# and P is an invertible matrix.
# Let Q = (P^{-1})^T, and let Φ = exp(Aδ). Then it holds that:
# Φ^T = Q * exp(Λ) * Q^{-1}
# Moreover, the eigenvectors of Φ^T are the columns of Q.
#
# This function returns a matrix `ρmat` of size `2n × NSTEPS` such that
# the i-th row of `ρmat` contains all support function evaluations
# ρ(Φ^k dᵢ, Ω₀) for each k = 0, …, NSTEPS-1,
# where dᵢ = vᵢ and vᵢ is the i-th eigenvector of Φ^T if 1 ≤ i ≤ n and for i > n
# dᵢ = -vᵢ is the negative of the i-th eigenvector of Φ^T.
function reach_homog_eig_LGG09_posneg(Λ::Vector{N}, Q, Ω₀, NSTEPS) where {N}

    n = length(Λ)
    @assert n == size(Q, 1) == size(Q, 2)

    ρmat = Matrix{N}(undef, 2*n, NSTEPS)
    Q₋ = -Q
    for j in 1:n
        λj = exp(Λ[j])

        Qj₊ = view(Q, :, j)
        reach_homog_dir_eig_LGG09!(view(ρmat, j, :), Ω₀,  Qj₊, λj, NSTEPS)

        Qj₋ = view(Q₋, :, j)
        reach_homog_dir_eig_LGG09!(view(ρmat, j+n, :), Ω₀, Qj₋, λj, NSTEPS)
    end
    return ρmat
end

# same as reach_homog_eig_LGG09_posneg, but only computes along Q
function reach_homog_eig_LGG09(Λ::Vector{N}, Q, Ω₀, NSTEPS) where {N}

    n = length(Λ)
    @assert n == size(Q, 1) == size(Q, 2)

    ρmat = Matrix{N}(undef, n, NSTEPS)
    for j in 1:n
        λj = exp(Λ[j])
        Qj = view(Q, :, j)
        reach_homog_dir_eig_LGG09!(view(ρmat, j, :), Ω₀,  Qj, λj, NSTEPS)
    end
    return ρmat
end

# compute an upper bound on ρ(eᵢ, Xₖ)
# ρmat is a matrix of size 2n x NSTEPS
# P is the eigenvectors matrix
# k is the time index

# `ρmat` is a matrix of size `2n × NSTEPS` such that
# the i-th row of `ρmat` contains all support function evaluations
# ρ(Φ^k dᵢ, Ω₀) for each k = 0, …, NSTEPS-1,
# where dᵢ = vᵢ and vᵢ is the i-th eigenvector of Φ^T if 1 ≤ i ≤ n and for i > n
# dᵢ = -vᵢ is the negative of the i-th eigenvector of Φ^T.
#
# This function receives:
# - `ρmat` -- matrix of support functions
# - `i`    -- an index `i` between `1` and `n`
# - `Qinv` -- the inverse of the matrix of eigenvectors of `Φ^T`, `Qinv`
# - `k`    -- time index `k` such that 1 ≤ k ≤ size(ρmat, 2) = NSTEPS
#
# The function returns an upper bound on ρ((Φ^T)^k eᵢ, Ω₀).
function _upper_bound_eig_dir(ρmat, i, Qinv, k)
    res = zero(eltype(Qinv))
    α = view(Qinv, :, i) # i-th column of Q^{-1}
    n = size(Qinv, 1)
    for (j, αⱼ) in enumerate(α)
        if αⱼ > 0
            res += αⱼ * ρmat[j, k]
        else
            res += abs(αⱼ) * ρmat[j + n, k]
        end
    end
    return res
end

# similar to _upper_bound_eig_dir but with eᵢ and -eᵢ
function _upper_bound_eig_dir_posneg(ρmat, i, Qinv, k)
    N = eltype(Qinv)
    a = zero(N)
    b = zero(N)

    α = view(Qinv, :, i) # i-th column of Q^{-1}
    n = size(Qinv, 1)
    for (j, αⱼ) in enumerate(α)
        if αⱼ > 0
            a += αⱼ * ρmat[j, k]
            b += αⱼ * ρmat[j + n, k]
        else
            a += abs(αⱼ) * ρmat[j + n, k]
            b += abs(αⱼ) * ρmat[j, k]
        end
    end
    return a, b
end

function _upper_bound_eig(ρmat::Matrix{N}, Qinv, NSTEPS) where {N}
    n = size(Qinv, 1)
    ρmat_box = Matrix{N}(undef, 2n, NSTEPS)
    @inbounds begin
        for k in 1:NSTEPS
            for i in 1:n
                a, b = _upper_bound_eig_dir_posneg(ρmat, i, Qinv, k)
                ρmat_box[i, k] = a
                ρmat_box[i + n, k] = b
            end
        end
    end
    return ρmat_box
end

# let x' = Ax, s.t. A has real eigenvalues of size n x n
# Aδ = P Λ P^{-1}, and Q = (P^{-1})^T
# Q has the eigenvectors of Φ^T = exp(A^T δ)
# the function assumes that the matrix Q is orthogonal and returns the box overapproximation
# of the flowipe by support function evaluation matrix
function reach_homog_eig_LGG09_box(Λ::Vector{N}, Q, Ω₀, NSTEPS) where {N}
    # both M₊ and M₋ are matrices of size n x NSTEPS
    # compute support functions along +vj where vj is column of Q (eigenvector of Φ^T)
    M₊ = reach_homog_eig_LGG09(Λ, Q, Ω₀, NSTEPS)

    # compute support functions along -vj where vj is column of Q (eigenvector of Φ^T)
    M₋ = reach_homog_eig_LGG09(Λ, -Q, Ω₀, NSTEPS)

    # center
    C = Q * (M₊ - M₋)/2

    # radius
    R = abs.(Q) * (M₊ + M₋)/2

    # support function of the hyperrectangular approximations
    B₊ = R + C
    B₋ = R - C

    return B₊, B₋
end

# ===================================
# Homogeneous case with invariant
# ===================================

function reach_homog_LGG09!(F::Vector{RT},
                            dirs::TN,
                            Ω₀::LazySet{N},
                            Φ::AbstractMatrix{N},
                            NSTEPS::Integer,
                            δ::N,
                            X::LazySet,
                            Δt0::TimeInterval,
                            cache,
                            threaded) where {N, VN, TN, SN, RT<:TemplateReachSet{N, VN, TN, SN}}

    # transpose coefficients matrix
    Φᵀ = copy(transpose(Φ))

    # preallocate output sequence
    ndirs = length(dirs)
    ρmat = Matrix{N}(undef, ndirs, NSTEPS)

    rᵢ = reduce(hcat, dirs) # one column per template direction
    rᵢ₊₁ = similar(rᵢ)

    # fill template reach-set sequence
    Δt = (zero(N) .. δ) + Δt0
    k = 1
    @inbounds while k <= NSTEPS

        # loop dispatch in threaded
        _reach_homog_LGG09_invariant!(Ω₀, rᵢ, ρmat, k, ndirs, threaded)

        # update cache for the next iteration
        mul!(rᵢ₊₁, Φᵀ, rᵢ)
        copy!(rᵢ, rᵢ₊₁)

        F[k] = TemplateReachSet(dirs, view(ρmat, :, k), Δt)
        _is_intersection_empty(X, set(F[k])) && break  # TODO pass disjointness method
        Δt += δ
        k += 1
    end
    if k < NSTEPS
        resize!(F, k-1)
    end
    return ρmat
end

function _reach_homog_LGG09_invariant!(Ω₀, rᵢ, ρmat, k, ndirs, threaded::Val{false})
    for j in 1:ndirs
        d = view(rᵢ, :, j)
        ρmat[j, k] = ρ(d, Ω₀)
    end
end

function _reach_homog_LGG09_invariant!(Ω₀, rᵢ, ρmat, k, ndirs, threaded::Val{true})
    Threads.@threads for j in 1:ndirs
        d = view(rᵢ, :, j)
        ρmat[j, k] = ρ(d, Ω₀)
    end
end


text:
v = rand(10)

a, b = 4, 8
c, d = 5, 9

M = randn(a, b)
N = rand(c, d)

p = rand(1:a, 6)
q = rand(1:c, 6)

r = rand(1:b, 10)
t = rand(1:d, 10)

kronprod = N ⊗ M
ikp = kronprod[p,q,r,t]

subsystem = kron(N, M)[a * (q .- 1) .+ p, b * (t .- 1) .+ r]
@test subsystem ≈ collect(ikp)

# result shortcut
#u = genvectrick(M, N, v, p, q, r, t)
u = ikp * v
# result naive
unaive = kron(N, M)[a * (q .- 1) .+ p, b * (t .- 1) .+ r] * v
@test all(u .≈ unaive)


text:
# Product between PolyVar and Monomial -> Monomial
function Base.:(*)(x::PolyVar{true}, y::PolyVar{true})
    if x === y
        Monomial{true}([x], [2])
    else
        Monomial{true}(x > y ? [x,y] : [y,x], [1,1])
    end
end
function multiplyvar(v::Vector{PolyVar{true}}, x::PolyVar{true})
    i = findfirst(w->w <= x, v)
    if (i != nothing && i > 0) && v[i] == x
        multiplyexistingvar(v, x, i)
    else
        insertvar(v, x, (i == nothing || i == 0) ? length(v)+1 : i)
    end
end
function Base.:(*)(x::PolyVar{true}, y::Monomial{true})
    w, updatez = multiplyvar(y.vars, x)
    Monomial{true}(w, updatez(y.z))
end
Base.:(*)(y::MonomialVector{true}, x::PolyVar{true}) = x * y
function Base.:(*)(x::PolyVar{true}, y::MonomialVector{true})
    w, updatez = multiplyvar(y.vars, x)
    MonomialVector{true}(w, updatez.(y.Z))
end
function multdivmono!(output_variables::Vector{PolyVar{true}},
                      v::Vector{PolyVar{true}}, x::Monomial{true}, op)
    if v == x.vars
        updatez! = (output, z) -> @. output = op(z, x.z)
    else
        w, maps = mergevars([v, x.vars])
        n = length(v)
        resize!(output_variables, length(w))
        output_variables[maps[1]] = v[1:n]
        updatez! = (output, z) -> begin
            resize!(output, length(w))
            z[maps[1]] = z[1:n]
            I = maps[1]; i = 1; lI = length(I)
            J = maps[2]; j = 1; lJ = length(J)
            while i <= lI || j <= lJ
                if i > lI || (j <= lJ && J[j] < I[i])
                    output[J[j]] = op(0, x.z[j])
                    j += 1
                elseif j > lJ || (i <= lI && I[i] < J[j])
                    output[I[i]] = op(z[I[i]], 0)
                    i += 1
                else
                    @assert I[i] == J[j]
                    output[I[i]] = op(z[I[i]], x.z[j])
                    i += 1
                    j += 1
                end
            end
        end
    end
    return updatez!
end
function multdivmono(v::Vector{PolyVar{true}}, x::Monomial{true}, op)
    if v == x.vars
        w = copy(v)
        updatez = z -> op.(z, x.z)
    else
        w, maps = mergevars([v, x.vars])
        updatez = z -> begin
            newz = zeros(Int, length(w))
            I = maps[1]; i = 1; lI = length(I)
            J = maps[2]; j = 1; lJ = length(J)
            while i <= lI || j <= lJ
                if i > lI || (j <= lJ && J[j] < I[i])
                    newz[J[j]] = op(0, x.z[j])
                    j += 1
                elseif j > lJ || (i <= lI && I[i] < J[j])
                    newz[I[i]] = op(z[i], 0)
                    i += 1
                else
                    @assert I[i] == J[j]
                    newz[I[i]] = op(z[i], x.z[j])
                    i += 1
                    j += 1
                end
            end
            newz
        end
    end
    w, updatez
end
function MP.mapexponents_to!(output::Monomial{true}, f::Function, x::Monomial{true}, y::Monomial{true})
    updatez! = multdivmono!(output.vars, x.vars, y, f)
    updatez!(output.z, x.z)
    return x
end
function MP.mapexponents!(f::Function, x::Monomial{true}, y::Monomial{true})
    updatez! = multdivmono!(x.vars, x.vars, y, f)
    updatez!(x.z, x.z)
    return x
end
function MP.mapexponents(f::Function, x::Monomial{true}, y::Monomial{true})
    w, updatez = multdivmono(x.vars, y, f)
    Monomial{true}(w, updatez(x.z))
end
function Base.:(*)(x::Monomial{true}, y::MonomialVector{true})
    w, updatez = multdivmono(y.vars, x, +)
    MonomialVector{true}(w, updatez.(y.Z))
end
Base.:(*)(y::MonomialVector{true}, x::Monomial{true}) = x * y
Base.:(*)(x::Monomial{true}, y::PolyVar{true}) = y * x


text:
EPSVALUE = 0.1

"""
    columngeneration(patients, visits, resources, timeslots, mastercalendar, timeDelta; setuponly = false, multithreading = false)

Builds the optimization problem and solves it.
If the problem should only be set up and not solved, apply 'setuponly = true.'
The problem can be solved with multithreading by setting 'multithreading = true'

"""
function columngeneration(patients, visits, resources, timeslots, mastercalendar, timeDelta; setuponly = false, multithreading = false)
    println("Building sets")
    timeslotsNDSparse = ndsparse(timeslots)
    subproblems = Subproblems()
    sets = Sets()
    sets.Vp = buildVp(visits)
    sets.Dv = buildDv(resources,visits)
    sets.Dp = buildDp(resources,visits)
    sets.J = keys(mastercalendar)
    sets.Jd = buildJd(timeslots,mastercalendar)
    sets.I = buildI(timeslots)
    sets.Pg = buildPg(patients,visits,sets.Vp,timeDelta)
    println("Setting up sub-problems")
    setup_sub!(subproblems,patients,visits,resources,timeslotsNDSparse,mastercalendar,timeDelta,sets,1)
    println("Setting up master")
    mp = setupmaster(subproblems,patients,resources,timeslots,mastercalendar,sets)
    generateInitialColumns!(mp,sets,timeslots,timeDelta)

    if setuponly
        return mp, subproblems, sets
    end
    iteration = 0
    done = false
    while !done
        iteration +=1
        println(iteration)
        println("Solving master")
        optimize!(mp.model)
        println("Master objective value = $(JuMP.objective_value(mp.model))")
        done = true
        status = termination_status(mp.model)
        if status != MOI.TerminationStatusCode(1)
            throw("Error: Non optimal masterproblem")
        end

        #ϕ = dual.(mp.consref_offtime)
        ϕ = []
        θ = dual.(mp.consref_onepatient)
        κ = dual.(mp.convexitycons)

        println("Solving subproblems")
        if multithreading
            subthreads = []
            for sub in values(subproblems.pricingproblems)
                push!(subthreads,Threads.@spawn solveSub!(sub,sets,ϕ,θ,κ))
            end
            #end
            for thread in subthreads
                wait(thread)
            end
        else
            for sub in values(subproblems.pricingproblems)
                solveSub!(sub,sets,ϕ,θ,κ)
                # solveSub!(sub,ϕ,θ,κ)
                # done2 = addcolumntomaster!(mp,sub,iteration,EPSVALUE) #TODO add patient info
                # done &= done2
                # count = 0
                # while !done2 && count < 0#length(sub.patients)
                #     count += 1
                #     solveSub!(sub,getIndexofPositiveVariables(sub.xvars)) #TODO time this
                #     done2 = addcolumntomaster!(mp,sub,iteration,EPSVALUE)
                # end
            end
            #end

        end
        done = addcolumntomaster!(mp,subproblems,iteration,EPSVALUE)

    end

    println("LP objective value = $(JuMP.objective_value(mp.model))")
    Gurobi.setparams!(mp.env, OutputFlag=1,TimeLimit=600)
    set_integer.(mp.lambda)
    optimize!(mp.model)
    println("MIP objective value = $(JuMP.objective_value(mp.model))")

    return mp, subproblems, sets
end
"Find a visit that fits with a patient and a visittype"
function findvisit(visits,type,patient)
   result = @from i in visits begin
      @where i.req_type == type && i.patientID == patient
      @select i.intID
      @collect
   end
   if length(result)== 1
      return first(result)
   end
   if length(result) == 0
       @warn "no visits found of type $type, patient has following types: $(visits |> @filter(_.patientID ==patient) |> @map(_.req_type) |> collect )"
       return -1
   end
   @warn "Multiple visits of same type found"
   return -1
end

"Extracts the solution into an indexed table"
function extractsolution(mp,sets,timeslots,visits)
    plannedappointments = table((pattern = Int64[],patientgroupID = Int64[],patientID = Int64[],visitID = Int64[], resourceID = Int64[], dayID= Int64[], timeslotID = Int64[]);pkey=[:visitID] )
    tempPg = sets.Pg |>@map(copy(_.patients)) |> collect
    for x in getPositiveVariables(mp.lambda)
       patientgroup = first(filter(λ->normalized_coefficient(mp.convexitycons[λ],x)> 0,eachindex(mp.convexitycons)))[1]
       patient = pop!(tempPg[patientgroup])
       for y in eachindex(mp.consref_onepatient)
           if JuMP.normalized_coefficient(mp.consref_onepatient[y],x) > 0
               push!(rows(plannedappointments),(pattern=x.index.value,patientgroupID = patientgroup,patientID = patient,visitID = findvisit(visits,timeslots[y[3]].type,patient), resourceID = y[1],dayID= y[2], timeslotID = y[3]))
           end
       end
    end
    plannedappointments
end

"Used to validate that the solution contain all visits"
function allvisits(visits,plannedappointments)
    pass = true
    novisitsfound = @from i in plannedappointments begin
        @left_outer_join j in visits on i.visitID equals j.intID
        @where j.intID != i.visitID
        @select (i.pattern,i.timeslot)
        @collect
    end
    if length(novisitsfound) > 0
        for i in novisitsfound
            @warn "no corresponding visit found for pattern $(i[1]) timeslot $(i[2])"
        end
        pass = false
    end
    badvisits = @from i in visits begin
        @join j in plannedappointments on i.intID equals j.visitID into k
        @where length(k) != 1
        @select (intID = i.intID,number=length(k))
        @collect
    end
    if length(badvisits) > 0
        toomanyplanned = badvisits |> @filter(_.number > 1) |> @map(_.intID) |> collect
        if length(toomanyplanned) > 0
                @warn "To many planned appointments for visits $(toomanyplanned)"
        end
        unplannedvisits = badvisits |> @filter(_.number ==0 ) |> @map(_.intID) |> collect
        if length(unplannedvisits) > 0
            @warn "No appointments planned for vistis $(unplannedvisits)"
        end
        pass = false
    end
    if pass
        println("test for all visits passed")
    end
    return pass
end

"Used to validate that every visit is in the correct order according to timeDelta"
function correctorder(plannedappointments,timeslots,timeDelta)
    pass = true
    for patient in plannedappointments |> @groupby(_.patientID)|> @map((patientID =key(_), timeslots= map(x-> x.timeslotID,_))) |> collect
        appointmenttimeslots = timeslots |> @filter(_.timeslotID in patient.timeslots) |> collect
        for i in appointmenttimeslots, j in appointmenttimeslots
            if haskey(timeDelta,(i.type,j.type))
                if  j.dayID-i.dayID < timeDelta[i.type,j.type].delta
                    @warn "Patient $(patient.patientID) has visits of type $(i.type) and $(j.type) that should have $(timeDelta[i.type,j.type]) day between them bu have only $(j.dayID-i.dayID)"
                    pass = false
                elseif(j.dayID == i.dayID && i.endTime > j.startTime)
                    @warn "Patient $(patient.patientID) has visits of type $(i.type) and $(j.type) planned on the same day but the one that should be done first ends at $(i.endTime) and the other starts at $(j.startTime)"
                    pass = false
                end
            elseif j.dayID == i.dayID && j.startTime < i.startTime && j.endTime > i.startTime
                pass = false
            end
        end
    end
    if pass
        println("test for correct order passed")
    end
    return pass
end

function deadline()
    pass = true
    if true
        pass = false
        println("deadline not implemented")
    end
    return pass
end

"Used to validate that the objective of the masterproblem is calculated correctly"
function daysused(mp,plannedappointments)
    pass = true
    if sum(plannedappointments |> @groupby(_.patientID) |> @map(length(unique(map(x->x.dayID,_)))) |> collect) != JuMP.objective_value(mp.model)
        pass = false
        println("Mip model objective wrong")
    else
        println("Mip model objective correct")
    end
    return pass
end

"Validates solutions "
function solutionvalidator(mp,sets,timeslots,visits,timeDelta)
    plannedappointments = extractsolution(mp,sets,timeslots,visits)
    pass = true
    pass = pass && allvisits(visits,plannedappointments)
    pass = pass && correctorder(plannedappointments,timeslots,timeDelta)
    pass = pass && daysused(mp,plannedappointments)
    pass = pass && deadline()
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

@testset "OpticalSystem" begin
    @testset "Single threaded trace makes sure function executes properly" begin
        conv = Examples.doubleconvex()
        rays = Emitters.Sources.CompositeSource(translation(-unitZ3()), repeat([Emitters.Sources.Source()], 100))
        trace(conv, rays)
        @test true #just want to verify that the trace function executed properly
    end
end


text:
using Cards
using Test

@testset "char" begin
    @test Cards.char(Card(5, Suit(1))) == '🃅'
    @test Cards.char(Card(14, Suit(0))) == '🃑'
    @test Cards.char(Card(12, Suit(2))) == '🂽'
    for c in Cards.deck
        @test parse(Card, Cards.char(c)) == c
    end
end

@testset "parsing" begin
    @test parse(Card, "2H") == Card(2, Suit(2))
    @test parse(Card, "10S") == Card(10, Suit(3))
    @test parse(Card, "A♦") == Card(14, Suit(1))
end


text:
#test/runtests.jl
#import Pkg; Pkg.add("Test")
using Test, DataFrames,SalesForceBulkApi

## Pulling the data ##
session = login("test@jltest-dev-ed.com", "9d3T67hTK8DwKjApVAiwZL4nmBmPGqFpMNnK2YoRE4B7Sgf78", "45.0")
all_object_fields_return = all_object_fields(session)
all_object_fields_return[[:name, :object]]
queries = ["Select Name From Account Limit 10", "Select LastName From Contact limit 10"]
res1 = sf_bulkapi_query(session, "Select LastName From Contact limit 10")
multi_result = multiquery(session, queries)
multi_result_all = multiquery(session, queries, true)

# Testing content #
@test eltype(session["sessionId"]) == Char
@test isa(all_object_fields_return, DataFrame)
@test size(all_object_fields_return, 1) > 1
@test size(all_object_fields_return, 2) > 50
@test size(res1) == (10,1)
@test typeof(multi_result) == Dict{Any,Any}
@test multi_result[queries[2]] == res1
@test multi_result[queries[1]][1,1] == "GenePoint"
@test size(multi_result[queries[1]],1) <= size(multi_result_all[queries[1]],1)


text:
@doc raw"""
`TableauIPRK`: Tableau of an Implicit Partitioned Runge-Kutta method
```math
\begin{align*}
P_{n,i} &= \dfrac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , &
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, V_{n,j} , &
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V_{n,i} , \\
F_{n,i} &= \dfrac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , &
P_{n,i} &= p_{n} + h  \sum \limits_{i=1}^{s} \bar{a}_{ij} \, F_{n,j} , &
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{b}_{i} \, F_{n,i} ,
\end{align*}
```
usually satisfying the symplecticity conditions
```math
\begin{align*}
b_{i} \bar{a}_{ij} + b_{j} a_{ji} &= b_{i} b_{j} , &
\bar{b}_i &= b_i .
\end{align*}
```
"""
struct TableauIPRK{T} <: AbstractTableauPRK{T}
    @HeaderTableau

    q::CoefficientsRK{T}
    p::CoefficientsRK{T}

    function TableauIPRK{T}(name, o, q, p) where {T}
        @assert q.s==p.s
        new(name, o, q.s, q, p)
    end
end

function TableauIPRK(name::Symbol, order::Int, q::CoefficientsRK{T}, p::CoefficientsRK{T}) where {T}
    TableauIPRK{T}(name, order, q, p)
end

function TableauIPRK(name::Symbol, order::Int, q::CoefficientsRK{T}) where {T}
    TableauIPRK{T}(name, order, q, q)
end

# TODO function readTableauIPRKFromFile(dir::AbstractString, name::AbstractString)

"Parameters for right-hand side function of implicit partitioned Runge-Kutta methods."
mutable struct ParametersIPRK{DT, TT, ET <: PODE{DT,TT}, D, S} <: Parameters{DT,TT}
    equ::ET
    tab::TableauIPRK{TT}
    Δt::TT

    t::TT
    q::Vector{DT}
    p::Vector{DT}
end

function ParametersIPRK(equ::ET, tab::TableauIPRK{TT}, Δt::TT) where {DT, TT, ET <: PODE{DT,TT}}
    q = zeros(DT, equ.d)
    p = zeros(DT, equ.d)

    ParametersIPRK{DT, TT, ET, equ.d, tab.s}(equ, tab, Δt, 0, q, p)
end

"""
Implicit partitioned Runge-Kutta integrator cache.

### Fields

* `q̃`: initial guess of q
* `p̃`: initial guess of p
* `ṽ`: initial guess of v
* `f̃`: initial guess of f
* `s̃`: holds shift due to periodicity of solution
* `Q`: internal stages of q
* `P`: internal stages of p
* `V`: internal stages of v
* `F`: internal stages of f
* `Y`: vector field of internal stages of q
* `Z`: vector field of internal stages of p
"""
struct IntegratorCacheIPRK{ST,D,S} <: PODEIntegratorCache{ST,D}
    q̃::Vector{ST}
    p̃::Vector{ST}
    ṽ::Vector{ST}
    f̃::Vector{ST}
    s̃::Vector{ST}

    Q::Vector{Vector{ST}}
    P::Vector{Vector{ST}}
    V::Vector{Vector{ST}}
    F::Vector{Vector{ST}}
    Y::Vector{Vector{ST}}
    Z::Vector{Vector{ST}}

    function IntegratorCacheIPRK{ST,D,S}() where {ST,D,S}
        # create temporary vectors
        q̃ = zeros(ST,D)
        p̃ = zeros(ST,D)
        ṽ = zeros(ST,D)
        f̃ = zeros(ST,D)
        s̃ = zeros(ST,D)

        # create internal stage vectors
        Q = create_internal_stage_vector(ST, D, S)
        P = create_internal_stage_vector(ST, D, S)
        V = create_internal_stage_vector(ST, D, S)
        F = create_internal_stage_vector(ST, D, S)
        Y = create_internal_stage_vector(ST, D, S)
        Z = create_internal_stage_vector(ST, D, S)

        new(q̃, p̃, ṽ, f̃, s̃, Q, P, V, F, Y, Z)
    end
end

"Implicit partitioned Runge-Kutta integrator."
struct IntegratorIPRK{DT, TT, PT <: ParametersIPRK{DT,TT},
                              ST <: NonlinearSolver{DT},
                              IT <: InitialGuessPODE{DT,TT}, D, S} <: IntegratorPRK{DT,TT}
    params::PT
    solver::ST
    iguess::IT
    cache::IntegratorCacheIPRK{DT,D,S}
end

function IntegratorIPRK(equation::PODE{DT,TT,VT,FT}, tableau::TableauIPRK{TT}, Δt::TT) where {DT,TT,VT,FT}
    D = equation.d
    M = equation.n
    S = tableau.s

    # create params
    params = ParametersIPRK(equation, tableau, Δt)

    # create solver
    solver = create_nonlinear_solver(DT, 2*D*S, params)

    # create initial guess
    iguess = InitialGuessPODE(get_config(:ig_interpolation), equation, Δt)

    # create cache
    cache = IntegratorCacheIPRK{DT,D,S}()

    # create integrator
    IntegratorIPRK{DT, TT, typeof(params), typeof(solver), typeof(iguess), D, S}(params, solver, iguess, cache)
end

function update_params!(int::IntegratorIPRK, sol::AtomicSolutionPODE)
    # set time for nonlinear solver and copy previous solution
    int.params.t  = sol.t
    int.params.q .= sol.q
    int.params.p .= sol.p
end

"Compute stages of implicit partitioned Runge-Kutta methods."
function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersIPRK{DT,TT,ET,D,S}) where {ST,DT,TT,ET,D,S}
    # create caches for internal stages
    Q = create_internal_stage_vector(ST, D, S)
    P = create_internal_stage_vector(ST, D, S)
    V = create_internal_stage_vector(ST, D, S)
    F = create_internal_stage_vector(ST, D, S)
    Y = create_internal_stage_vector(ST, D, S)
    Z = create_internal_stage_vector(ST, D, S)

    # compute stages from nonlinear solver solution x
    compute_stages!(x, Q, V, Y, P, F, Z, params)

    # compute b = - [(Y-AV), (Z-AF)]
    for i in 1:S
        for k in 1:D
            b[2*(D*(i-1)+k-1)+1] = - Y[i][k]
            b[2*(D*(i-1)+k-1)+2] = - Z[i][k]
            for j in 1:S
                b[2*(D*(i-1)+k-1)+1] += params.tab.q.a[i,j] * V[j][k]
                b[2*(D*(i-1)+k-1)+2] += params.tab.p.a[i,j] * F[j][k]
            end
        end
    end
end

function compute_stages!(x::Vector{ST}, Q::Vector{Vector{ST}}, V::Vector{Vector{ST}}, Y::Vector{Vector{ST}},
                                        P::Vector{Vector{ST}}, F::Vector{Vector{ST}}, Z::Vector{Vector{ST}},
                                        params::ParametersIPRK{DT,TT,ET,D,S}) where {ST,DT,TT,ET,D,S}
    local tqᵢ::TT
    local tpᵢ::TT

    for i in 1:S
        for k in 1:D
            # copy y to Y and Z
            Y[i][k] = x[2*(D*(i-1)+k-1)+1]
            Z[i][k] = x[2*(D*(i-1)+k-1)+2]

            # compute Q and P
            Q[i][k] = params.q[k] + params.Δt * Y[i][k]
            P[i][k] = params.p[k] + params.Δt * Z[i][k]
        end

        # compute time of internal stage
        tqᵢ = params.t + params.Δt * params.tab.q.c[i]
        tpᵢ = params.t + params.Δt * params.tab.p.c[i]

        # compute v(Q,P) and f(Q,P)
        params.equ.v(tqᵢ, Q[i], P[i], V[i])
        params.equ.f(tpᵢ, Q[i], P[i], F[i])
    end
end

function initialize!(int::IntegratorIPRK, sol::AtomicSolutionPODE)
    sol.t̅ = sol.t - timestep(int)

    int.params.equ.v(sol.t, sol.q, sol.p, sol.v)
    int.params.equ.f(sol.t, sol.q, sol.p, sol.f)

    initialize!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f,
                            sol.t̅, sol.q̅, sol.p̅, sol.v̅, sol.f̅)
end

function initial_guess!(int::IntegratorIPRK, sol::AtomicSolutionPODE)
    for i in eachstage(int)
        evaluate!(int.iguess, sol.q, sol.p, sol.v, sol.f,
                              sol.q̅, sol.p̅, sol.v̅, sol.f̅,
                              int.cache.q̃, int.cache.p̃, int.cache.ṽ, int.cache.f̃,
                              tableau(int).q.c[i],
                              tableau(int).p.c[i])

        for k in eachdim(int)
            int.cache.V[i][k] = int.cache.ṽ[k]
            int.cache.F[i][k] = int.cache.f̃[k]
        end
    end
    for i in eachstage(int)
        for k in eachdim(int)
            int.solver.x[2*(ndims(int)*(i-1)+k-1)+1] = 0
            int.solver.x[2*(ndims(int)*(i-1)+k-1)+2] = 0
            for j in eachstage(int)
                int.solver.x[2*(ndims(int)*(i-1)+k-1)+1] += tableau(int).q.a[i,j] * int.cache.V[j][k]
                int.solver.x[2*(ndims(int)*(i-1)+k-1)+2] += tableau(int).p.a[i,j] * int.cache.F[j][k]
            end
        end
    end
end

"Integrate ODE with implicit partitioned Runge-Kutta integrator."
function integrate_step!(int::IntegratorIPRK{DT,TT}, sol::AtomicSolutionPODE{DT,TT}) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int, sol)

    # compute initial guess
    initial_guess!(int, sol)

    # reset cache
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector fields at internal stages
    compute_stages!(int.solver.x, int.cache.Q, int.cache.V, int.cache.Y,
                                  int.cache.P, int.cache.F, int.cache.Z, int.params)

    # compute final update
    update_solution!(sol.q, sol.q̃, int.cache.V, tableau(int).q.b, tableau(int).q.b̂, timestep(int))
    update_solution!(sol.p, sol.p̃, int.cache.F, tableau(int).p.b, tableau(int).p.b̂, timestep(int))

    # copy solution to initial guess
    update!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f)
end


text:
using BluseBeamformerRecipes
using Test

@testset "BluseBeamformerRecipes.jl" begin
    # Write your tests here.
end


text:
# Display the entire type hierarchy starting from the specified `roottype`
function subtypetree(roottype, level = 1, indent = 4)
    level == 1 && println(roottype)
    for s in subtypes(roottype)
        println(join(fill(" ", level * indent)) * string(s))
        subtypetree(s, level + 1, indent)
    end
end


text:
import LinearAlgebra: transpose
import Splittings: Geometry
import Splittings: SplineNN, interpolate!

@testset "Natural Spline" begin

    geom = Geometry( 101, 101, -5.0, -5.0, 0.1, 0.1 )
    x1, x2 = geom.x1grid, geom.x2grid
    dt = 0.5
    fc = exp.(-x1.^2) .* transpose(exp.(-x2.^2))
    spline = SplinePP( geom )
    interpolate!( spline, fc, dt, dt )
    fe = exp.(-(x1.-dt).^2) .* transpose(exp.(-(x2.-dt).^2))
    error = maximum( abs.( fe .- fc ))
    println( " error  = $error ")
    @test fc ≈ fe

    geom = Geometry( 51, 101, -5.0, -5.0, 0.2, 0.1 )
    x1, x2 = geom.x1grid, geom.x2grid
    dt = 0.5
    fc = exp.(-x1.^2) .* transpose(exp.(-x2.^2))
    spline = SplineNN( geom )
    interpolate!( spline, fc, dt, dt )
    fe = exp.(-(x1.-dt).^2) .* transpose(exp.(-(x2.-dt).^2))
    error = maximum( abs.( fe .- fc ))
    println( " error  = $error ")
    @test error < 1e-4

end


text:
# parameters that govern counterfactual

## kappa remains at 1972 level
for t=1:parameters[:T]
	parameters[:kappa_mnjt][:,:,:,t] = parameters[:kappa_mnjt][:,:,:,1]
end


text:
"""
Colored console appender

It appends log events to `STDOUT` or `STDERR` using a layout specified by the user.
"""
type ColorConsole <: Appender
    name::AbstractString
    layout::LAYOUT
    filter::FILTER
    state::LifeCycle.State

    io::IO
end
function ColorConsole(config::Dict{AbstractString,Any})
    nm = get(config,  "name", "STDOUT")
    lyt= get(config,  "layout", nothing)
    io = get(config,  "target", :STDOUT)
    flt = get(config, "filter", nothing)
    ColorConsole(nm, LAYOUT(lyt), FILTER(flt), LifeCycle.INITIALIZED, io == :STDOUT ? STDOUT : STDERR)
end
ColorConsole(;kwargs...) = map(e->(string(e[1]),e[2]), kwargs) |> Dict{AbstractString,Any} |> ColorConsole
show(io::IO, apnd::ColorConsole) = print(io, "ColorConsole($(apnd.name))")

const ColorReset = UInt8[0x1b, 0x5b, 0x30, 0x6d]
const LevelColors = Dict(
    Level.ALL   => "\e[0m",
    Level.TRACE => "\e[33m",
    Level.DEBUG => "\e[32;1m",
    Level.INFO  => "\e[36;1m",
    Level.WARN  => "\e[33;1m",
    Level.ERROR => "\e[31;1m",
    Level.FATAL => "\e[31m",
    Level.OFF   => "\e[8m"
)
name(apnd::ColorConsole) = isempty(apnd.name) ? string(typeof(apnd)) : apnd.name
layout(apnd::ColorConsole) = apnd.layout

function append!(apnd::ColorConsole, evnt::Event)
    isnull(apnd.layout) && return
    col = get(LevelColors, level(evnt), "\e[39;1m")
    write(apnd.io, col)
    write(apnd.io, serialize(apnd.layout, evnt))
    write(apnd.io, ColorReset)
end