text:
abstract type VSIEOperator <: IntegralOperator end
abstract type VolumeSurfaceOperator <: VSIEOperator end
abstract type BoundarySurfaceOperator <: VSIEOperator end
abstract type VSIEOperatorT <: VSIEOperator end
abstract type VolumeSurfaceOperatorT <: VolumeSurfaceOperator end
abstract type BoundarySurfaceOperatorT <: BoundarySurfaceOperator end

struct KernelValsVSIE{T,U,P,Q,K}
    gamma::U
    vect::P
    dist::T
    green::U
    gradgreen::Q
    tau::K
end

function kernelvals(viop::VSIEOperator, p ,q)
    Y = viop.gamma;
    r = cartesian(p)-cartesian(q)
    R = norm(r)
    yR = Y*R

    expn = exp(-yR)
    green = expn / (4*pi*R)
    gradgreen = - (Y +1/R)*green/R*r

    tau = viop.tau(cartesian(q))

    KernelValsVSIE(Y,r,R, green, gradgreen,tau)
end

struct VSIESingleLayer{T,U,P} <: VolumeSurfaceOperator
    gamma::T
    α::U
    β::U
    tau::P
end

struct VSIEBoundary{T,U,P} <: BoundarySurfaceOperator
    gamma::T
    α::U
    tau::P
end

struct VSIEDoubleLayer{T,U,P} <: VolumeSurfaceOperator
    gamma::T
    α::U
    tau::P
end

#=
struct VSIESingleLayer2{T,U,P} <: VolumeSurfaceOperator
    gamma::T
    α::U
    β::U
    tau::P
end

struct VSIEBoundaryT{T,U,P} <: BoundarySurfaceOperatorT
    gamma::T
    α::U
    tau::P
end

struct VSIEDoubleLayerT{T,U,P} <: VolumeSurfaceOperatorT
    gamma::T
    α::U
    tau::P
end
=#
scalartype(op::VSIEOperator) = typeof(op.gamma)

export VSIE

struct VSIEIntegrand{S,T,O,K,L}
    test_tetrahedron_element::S
    trial_triangle_element::T
    op::O
    test_local_space::K
    trial_local_space::L
end

#=
struct VSIEIntegrandT{S,T,O,K,L}
    test_tetrahedron_element::S
    trial_triangle_element::T
    op::O
    test_local_space::K
    trial_local_space::L
end
=#

function (igd::VSIEIntegrand)(u,v)

    #mesh points
    tgeo = neighborhood(igd.test_tetrahedron_element,v)
    bgeo = neighborhood(igd.trial_triangle_element,u)

    #kernel values
    kerneldata = kernelvals(igd.op,tgeo,bgeo)

    #values & grad/div/curl of local shape functions
    tval = igd.test_local_space(tgeo)
    bval = igd.trial_local_space(bgeo)

    #jacobian
    j = jacobian(tgeo) * jacobian(bgeo)

    integrand(igd.op, kerneldata,tval,tgeo,bval,tgeo) * j
end

#=
function (igd::VSIEIntegrandT)(u,v)

    #mesh points
    tgeo = neighborhood(igd.test_tetrahedron_element,u)
    bgeo = neighborhood(igd.trial_triangle_element,v)

    #kernel values
    kerneldata = kernelvals(igd.op,tgeo,bgeo)

    #values & grad/div/curl of local shape functions
    tval = igd.test_local_space(tgeo)
    bval = igd.trial_local_space(bgeo)

    #jacobian
    j = jacobian(tgeo) * jacobian(bgeo)

    integrand(igd.op, kerneldata,tval,tgeo,bval,tgeo) * j
end
=#

function integrand(viop::VSIESingleLayer, kerneldata, tvals, tgeo, bvals, bgeo)

    gx = @SVector[tvals[i].value for i in 1:3]
    fy = @SVector[bvals[i].value for i in 1:4]

    dgx = @SVector[tvals[i].divergence for i in 1:3]
    dfy = @SVector[bvals[i].divergence for i in 1:4]

    G = kerneldata.green

    Ty = kerneldata.tau

    α = viop.α
    β = viop.β

    @SMatrix[α * dot(gx[i],Ty*fy[j]) * G + β * (dgx[i] * Ty*dfy[j])*G for i in 1:3, j in 1:4]
end

function integrand(viop::VSIEBoundary, kerneldata, tvals, tgeo, bvals, bgeo)

    dgx = @SVector[tvals[i].divergence for i in 1:3]
    fy = @SVector[bvals[i].value for i in 1:1]

    G = kerneldata.green

    Ty = kerneldata.tau

    α = viop.α

    @SMatrix[α * Ty * dgx[i] * fy[j] * G for i in 1:3, j in 1:1]
end

#=
function integrand(viop::VSIESingleLayer2, kerneldata, tvals, tgeo, bvals, bgeo)

    gx = @SVector[tvals[i].value for i in 1:3]
    fy = @SVector[bvals[i].value for i in 1:4]

    dgx = @SVector[tvals[i].divergence for i in 1:3]
    dfy = @SVector[bvals[i].divergence for i in 1:4]

    G = kerneldata.green

    Ty = kerneldata.tau

    α = viop.α
    β = viop.β

    @SMatrix[α * dot(gx[i],Ty*fy[j]) * G - β * (dgx[i] * Ty*dfy[j])*G for i in 1:3, j in 1:4]
end

function integrand(viop::VSIEBoundaryT, kerneldata, tvals, tgeo, bvals, bgeo)

    gx = @SVector[tvals[i].value for i in 1:1]
    dfy = @SVector[bvals[i].divergence for i in 1:3]

    G = kerneldata.green

    Tx = kerneldata.tau

    α = viop.α

    @SMatrix[α * Tx * gx[i] * dfy[j] * G for i in 1:3, j in 1:1]
end
=#

function integrand(viop::VSIEDoubleLayer, kerneldata, tvals, tgeo, bvals, bgeo)
    gx = @SVector[tvals[i].value for i in 1:3]
    fy = @SVector[bvals[i].value for i in 1:4]

    gradG = kerneldata.gradgreen

    Ty = kerneldata.tau

    α = viop.α

    @SMatrix[α * dot(cross(Ty*fy[j],gx[i]),gradG) for i in 1:3, j in 1:4]
end

#=
function integrand(viop::VSIEDoubleLayerT, kerneldata, tvals, tgeo, bvals, bgeo)
    gx = @SVector[tvals[i].value for i in 1:4]
    fy = @SVector[bvals[i].value for i in 1:3]

    gradG = kerneldata.gradgreen

    Ty = kerneldata.tau

    α = viop.α

    @SMatrix[α * transpose(cross(fy[j],gx[i])) *gradG for i in 1:4, j in 1:3]
end
=#

defaultquadstrat(op::VSIEOperator, tfs, bfs) = SauterSchwab3DQStrat(3,3,3,3,3,3)

function quaddata(op::VSIEOperator,
    test_local_space::RefSpace, trial_local_space::RefSpace,
    test_charts, trial_charts, qs::SauterSchwab3DQStrat)

    #The combinations of rules (6,7) and (5,7 are) BAAAADDDD
    # they result in many near singularity evaluations with any
    # resemblence of accuracy going down the drain! Simply don't!
    # (same for (5,7) btw...).
    t_qp = quadpoints(test_local_space,  test_charts,  (qs.outer_rule))
    b_qp = quadpoints(trial_local_space, trial_charts, (qs.inner_rule))

    sing_qp = (SauterSchwab3D._legendre(qs.sauter_schwab_1D,0,1),
               SauterSchwab3D._shunnham2D(qs.sauter_schwab_2D),
               SauterSchwab3D._shunnham3D(qs.sauter_schwab_3D),
               SauterSchwab3D._shunnham4D(qs.sauter_schwab_4D),)

    return (tpoints=t_qp, bpoints=b_qp, sing_qp=sing_qp)
end

quadrule(op::VolumeSurfaceOperator, g::RefSpace, f::RefSpace, i, τ, j, σ, qd, qs) = qr_volume(op, g, f, i, τ, j, σ, qd, qs)

function qr_volume(op::VolumeSurfaceOperator, g::RefSpace, f::RefSpace, i, τ, j, σ, qd,
    qs::SauterSchwab3DQStrat)

    @assert (length(τ.vertices)==3 && length(σ.vertices)==4)  "Expected simplex wrong"

    dtol = 1.0e3 * eps(eltype(eltype(τ.vertices)))

    hits = 0
    idx_t = Int64[]
    idx_s = Int64[]
    sizehint!(idx_t,4)
    sizehint!(idx_s,4)
    dmin2 = floatmax(eltype(eltype(τ.vertices)))
    D = dimension(τ)+dimension(σ)
    for (i,t) in enumerate(τ.vertices)
        for (j,s) in enumerate(σ.vertices)
            d2 = LinearAlgebra.norm_sqr(t-s)
            dmin2 = min(dmin2, d2)
            if d2 < dtol
                push!(idx_t,j)
                push!(idx_s,i)
                hits +=1
                break
            end
        end
    end

    #singData = SauterSchwab3D.Singularity{D,hits}(idx_t, idx_s )

    hits == 3 && return SauterSchwab3D.CommonFace5D_S(SauterSchwab3D.Singularity5DFace(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3]))
    hits == 2 && return SauterSchwab3D.CommonEdge5D_S(SauterSchwab3D.Singularity5DEdge(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3]))
    hits == 1 && return SauterSchwab3D.CommonVertex5D_S(SauterSchwab3D.Singularity5DPoint(idx_t,idx_s),(qd.sing_qp[3],qd.sing_qp[2]))
    #hits == 0 && return SauterSchwab3D.PositiveDistance5D_S(SauterSchwab3D.Singularity5DPositiveDistance(),(qd.sing_qp[3],qd.sing_qp[2]))

    return DoubleQuadRule(
        qd[1][1,i],
        qd[2][1,j])

end

quadrule(op::BoundarySurfaceOperator, g::RefSpace, f::RefSpace, i, τ, j, σ, qd, qs) = qr_boundary(op, g, f, i, τ, j, σ, qd, qs)

function qr_boundary(op::BoundarySurfaceOperator, g::RefSpace, f::RefSpace, i, τ, j,  σ, qd,
    qs::SauterSchwab3DQStrat)

    dtol = 1.0e3 * eps(eltype(eltype(τ.vertices)))

    hits = 0
    idx_t = Int64[]
    idx_s = Int64[]
    sizehint!(idx_t,4)
    sizehint!(idx_s,4)
    dmin2 = floatmax(eltype(eltype(τ.vertices)))
    D = dimension(τ)+dimension(σ)
    for (i,t) in enumerate(τ.vertices)
        for (j,s) in enumerate(σ.vertices)
            d2 = LinearAlgebra.norm_sqr(t-s)
            dmin2 = min(dmin2, d2)
            if d2 < dtol
                push!(idx_t,i)
                push!(idx_s,j)
                hits +=1
                break
            end
        end
    end

    #singData = SauterSchwab3D.Singularity{D,hits}(idx_t, idx_s )

    hits == 3 && return SauterSchwab3D.CommonFace4D_S(SauterSchwab3D.Singularity4DFace(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[3]))
    hits == 2 && return SauterSchwab3D.CommonEdge4D_S(SauterSchwab3D.Singularity4DEdge(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2]))
    hits == 1 && return SauterSchwab3D.CommonVertex4D_S(SauterSchwab3D.Singularity4DPoint(idx_t,idx_s),(qd.sing_qp[2]))

    return DoubleQuadRule(
        qd[1][1,i],
        qd[2][1,j])

end


text:
__precompile__(false)

#https://github.com/shiena/ansicolor/blob/master/README.md

#ThreadMutex = Threads.Mutex()
#PrintMutex = Threads.Mutex()

#function call(f::Function, args...;mutex=ThreadMutex)
#	Threads.lock(mutex)
#	result = f(args...)
#	Threads.unlock(mutex)
#  result
#end

#print(args...) = call(Base.print, args...; mutex=PrintMutex)
#println(args...) = call(Base.println, args...; mutex=PrintMutex)
#error(args...) = call(Base.error, args...; mutex=PrintMutex)
#warn(args...) = call(Base.warn, args...; mutex=PrintMutex)
#sleep(n) = Libc.systemsleep(n)

#include("Includes.jl")
#using ThreadManager

## COMPILE C File
const compileAndLink = isdefined(@__MODULE__,:USE_COMPILE)
if compileAndLink
  include(joinpath(@__DIR__,"compileAndLink.jl"))
  compileWithGCC()
end

ext: jl

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
#!/bin/bash
#=
exec julia --project=. -e 'include(popfirst!(ARGS))' \
    "${BASH_SOURCE[0]}" "$@"
=#

using MagneticResonanceSignals
using StaticArrays

function strip_twix_adjustment(in_file, out_file)
    open(in_file, "r") do src
        open(out_file, "w") do dest
            twix_id, num_measurements, mdata = MagneticResonanceSignals.read_meas_headers(src)

            # Create temp buffer
            buf = IOBuffer()

            @assert length(mdata) == 2 "Expected to have 2 measurements for adjustment and actual scan"
            meas = mdata[2]

            meas_offset = mdata[1].meas_offset

            # Change num_measurement to 1
            write(buf, SVector{2,Int32}(twix_id, 1))
            write(buf, meas.meas_uid)
            write(buf, meas.file_id)
            write(buf, meas_offset)
            write(buf, meas.meas_length)
            write(buf, meas.patient_name)
            write(buf, meas.protocol_name)

            # Seek to selected data
            seek(src, meas.meas_offset)

            write(buf, zeros(UInt8, meas_offset - position(buf)))
            @assert position(buf) == meas_offset
            write(buf, read(src, meas.meas_length))

            seekstart(buf)
            write(dest, read(buf))
        end
    end
end

using ArgParse

argdef = ArgParseSettings(exc_handler=isinteractive() ? ArgParse.debug_handler : ArgParse.default_handler,
                          prog="Twix Adjustment Remover",
                          description="""
                            Small program to remove adjustment data on Siemens twix, which relies heavily on MagneticResonanceSignals.
                            This program will take the in_file and write the no adjustment twix in out_file.
                            """)
@add_arg_table argdef begin
    "in_file"
        help = "path to the input twix file (Siemens meas_*.dat RAID format)"
        required = true
    "out_file"
        help = "path to the output file Siemens twix without adjustment data"
        required = true
end

function process_args(argvec)
    parsed_args = parse_args(argvec, argdef)

    in_file = parsed_args["in_file"]
    out_file = parsed_args["out_file"]

    (in_file=in_file, out_file=out_file)
end

args = process_args(ARGS)
strip_twix_adjustment(args.in_file, args.out_file)

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
using DataFrames, DataFramesMeta, Statistics
import FileIO:save
using PyPlot

# Create the parameterDict with the optimal parameters,
# after having ran the gridsearch
# The parameterDict has the following form:
# parameterDict = Dict(
# 		(stochasticityvalue1, pcvalue1) => optimalvalueparameter1
# 		(stochasticityvalue2, pcvalue2) => optimalvalueparameter2
# 		(stochasticityvalue3, pcvalue3) => optimalvalueparameter3
#                 ...
#                 )

# For the Categorical task (Section 2.3.2)
function findoptimalparams_multi(df;
                        savepath = homedir(),
                        filename = "multi_varsmile")
    parameterDict = Dict{Tuple{Float64, Float64},  Float64}()

    changeprobvals = unique(df[:changeprobability])
    stochvals = unique(df[:stochasticity])
    for i in 1:length(stochvals)
        for j in 1:length(changeprobvals)
            dftrimmed = @where(df, :stochasticity .== stochvals[i],
                        :changeprobability .== changeprobvals[j])
            minerror = minimum(dftrimmed[:meanterrors])
            if isnan(minerror)
                parameterDict[(stochvals[i], changeprobvals[j])] = -1.
                @show (stochvals[i], changeprobvals[j])
            else
                optimparam = [dftrimmed[:param][i] for i in 1:length(dftrimmed[:meanterrors])
                            if dftrimmed[:meanterrors][i]==minerror][1]
                parameterDict[(stochvals[i], changeprobvals[j])] = optimparam
            end
        end
    end
    save(joinpath(savepath, "parameterDict_"*filename*".jld2"),
        Dict("parameterDict" => parameterDict))
    parameterDict
end

#[parameterDict[(i, 0.1)] for i in stochvals]
# For the Gaussian task (Section 2.3.1)
function findoptimalparams_gauss(df;
                        savepath = homedir(),
                        filename = "gaussian_pf20")
    parameterDict = Dict{Tuple{Float64, Float64},  Float64}()

    changeprobvals = unique(df[:changeprobability])
    sigmavals = unique(df[:sigma])
    for i in 1:length(sigmavals)
        for j in 1:length(changeprobvals)
            dftrimmed = @where(df, :sigma .== sigmavals[i],
                        :changeprobability .== changeprobvals[j])
            minerror = minimum(dftrimmed[:meanmerrors])
            if isnan(minerror)
                parameterDict[(sigmavals[i], changeprobvals[j])] = -1.
                @show (sigmavals[i], changeprobvals[j])
            else
                optimparam = [dftrimmed[:param][i] for i in 1:length(dftrimmed[:meanmerrors])
                            if dftrimmed[:meanmerrors][i]==minerror][1]
                parameterDict[(sigmavals[i], changeprobvals[j])] = optimparam
            end
        end
    end
    save(joinpath(savepath, "parameterDict_"*filename*".jld2"),
        Dict("parameterDict" => parameterDict))
    parameterDict
end

function exploreparameterDict(parameterDict; svalue = 0.1)

    pcs = [j  for (i,j) in keys(parameterDict) if i == svalue]
    optimalvals = [parameterDict[(i,j)] for (i,j) in keys(parameterDict) if i == svalue]

    fig = figure(); ax = gca()
    ax.scatter(pcs, optimalvals)
    ax.grid(true)
end


text:
@testset "Default" begin
    @testset "Library" begin
        function default_library_test(d)
            @test default_library(d(2), Float64) isa Polyhedra.DefaultLibrary
            @test default_library(d(1), Float64) isa IntervalLibrary
            @test default_library(d(1), Int) isa IntervalLibrary
            @test default_library(d(1), BigInt) isa IntervalLibrary
            @test default_library(d(2), AbstractFloat) isa Polyhedra.DefaultLibrary{Float64}
            @test default_library(d(1), AbstractFloat) isa IntervalLibrary{Float64}
        end
        default_library_test(N -> StaticArrays.Size((N,)))
        default_library_test(identity)
    end
    @testset "LP Solver" begin
        @test Polyhedra.default_solver(hrep([HalfSpace((@SVector [1, 2]), 3)])) === nothing
        @test Polyhedra.solver(hrep([HalfSpace((@SVector [1, 2]), 3)]), lp_solver) === lp_solver
    end
end


text:
# See also fastmath_able.jl for where rules are defined simple base functions
# that also have FastMath versions.

@scalar_rule copysign(y, x) (ifelse(signbit(x)!=signbit(y), -one(y), +one(y)), NoTangent())

@scalar_rule one(x) zero(x)
@scalar_rule zero(x) zero(x)
@scalar_rule transpose(x) true

# `adjoint`

frule((_, Δz), ::typeof(adjoint), z::Number) = (z', Δz')

function rrule(::typeof(adjoint), z::Number)
    adjoint_pullback(ΔΩ) = (NoTangent(), ΔΩ')
    return (z', adjoint_pullback)
end

# `real`

@scalar_rule real(x::Real) true

frule((_, Δz), ::typeof(real), z::Number) = (real(z), real(Δz))

function rrule(::typeof(real), z::Number)
    # add zero(z) to embed the real number in the same number type as z
    real_pullback(ΔΩ) = (NoTangent(), real(ΔΩ) + zero(z))
    return (real(z), real_pullback)
end

# Conversions to Float

@scalar_rule float(x) true
@scalar_rule Float64(x::Real) true
@scalar_rule Float32(x::Real) true
@scalar_rule AbstractFloat(x::Real) true

# `imag`

@scalar_rule imag(x::Real) ZeroTangent()

frule((_, Δz), ::typeof(imag), z::Complex) = (imag(z), imag(Δz))

function rrule(::typeof(imag), z::Complex)
    imag_pullback(ΔΩ) = (NoTangent(), real(ΔΩ) * im)
    return (imag(z), imag_pullback)
end

# `Complex`

frule((_, Δz), ::Type{T}, z::Number) where {T<:Complex} = (T(z), Complex(Δz))
function frule((_, Δx, Δy), ::Type{T}, x::Number, y::Number) where {T<:Complex}
    return (T(x, y), Complex(Δx, Δy))
end

function rrule(::Type{T}, z::Complex) where {T<:Complex}
    Complex_pullback(ΔΩ) = (NoTangent(), Complex(ΔΩ))
    return (T(z), Complex_pullback)
end
function rrule(::Type{T}, x::Real) where {T<:Complex}
    Complex_pullback(ΔΩ) = (NoTangent(), real(ΔΩ))
    return (T(x), Complex_pullback)
end
function rrule(::Type{T}, x::Number, y::Number) where {T<:Complex}
    project_x = ProjectTo(x)
    project_y = ProjectTo(y)

    function Complex_pullback(Ω̄)
        ΔΩ = unthunk(Ω̄)
        return (NoTangent(), project_x(real(ΔΩ)), project_y(imag(ΔΩ)))
    end
    return (T(x, y), Complex_pullback)
end

# `hypot`

@scalar_rule hypot(x::Real) sign(x)

function frule((_, Δz), ::typeof(hypot), z::Complex)
    Ω = hypot(z)
    ∂Ω = realdot(z, Δz) / ifelse(iszero(Ω), one(Ω), Ω)
    return Ω, ∂Ω
end

function rrule(::typeof(hypot), z::Complex)
    Ω = hypot(z)
    function hypot_pullback(ΔΩ)
        return (NoTangent(), (real(ΔΩ) / ifelse(iszero(Ω), one(Ω), Ω)) * z)
    end
    return (Ω, hypot_pullback)
end

@scalar_rule fma(x, y, z) (y, x, true)
@scalar_rule muladd(x, y, z) (y, x, true)
@scalar_rule rem2pi(x, r::RoundingMode) (true, NoTangent())
@scalar_rule(
    mod(x, y),
    @setup((u, nan) = promote(x / y, NaN16), isint = isinteger(x / y)),
    (ifelse(isint, nan, one(u)), ifelse(isint, nan, -floor(u))),
)

@scalar_rule deg2rad(x) deg2rad(one(x))
@scalar_rule rad2deg(x) rad2deg(one(x))

@scalar_rule(ldexp(x, y), (2^y, NoTangent()))

# Can't multiply though sqrt in acosh because of negative complex case for x
@scalar_rule acosh(x) inv(sqrt(x - 1) * sqrt(x + 1))
@scalar_rule acoth(x) inv(1 - x ^ 2)
@scalar_rule acsch(x) -(inv(x ^ 2 * sqrt(1 + x ^ -2)))
@scalar_rule acsch(x::Real) -(inv(abs(x) * sqrt(1 + x ^ 2)))
@scalar_rule asech(x) -(inv(x * sqrt(1 - x ^ 2)))
@scalar_rule asinh(x) inv(sqrt(x ^ 2 + 1))
@scalar_rule atanh(x) inv(1 - x ^ 2)

@scalar_rule acosd(x) -inv(deg2rad(sqrt(1 - x ^ 2)))
@scalar_rule acotd(x) -inv(deg2rad(1 + x ^ 2))
@scalar_rule acscd(x) -inv(deg2rad(x^2 * sqrt(1 - x ^ -2)))
@scalar_rule acscd(x::Real) -inv(deg2rad(abs(x) * sqrt(x ^ 2 - 1)))
@scalar_rule asecd(x) inv(deg2rad(x ^ 2 * sqrt(1 - x ^ -2)))
@scalar_rule asecd(x::Real) inv(deg2rad(abs(x) * sqrt(x ^ 2 - 1)))
@scalar_rule asind(x) inv(deg2rad(sqrt(1 - x ^ 2)))
@scalar_rule atand(x) inv(deg2rad(1 + x ^ 2))

@scalar_rule cot(x) -((1 + Ω ^ 2))
@scalar_rule coth(x) -(csch(x) ^ 2)
@scalar_rule cotd(x) -deg2rad(1 + Ω ^ 2)
@scalar_rule csc(x) -Ω * cot(x)
@scalar_rule cscd(x) -deg2rad(Ω * cotd(x))
@scalar_rule csch(x) -(coth(x)) * Ω
@scalar_rule sec(x) Ω * tan(x)
@scalar_rule secd(x) deg2rad(Ω * tand(x))
@scalar_rule sech(x) -(tanh(x)) * Ω

@scalar_rule acot(x) -(inv(1 + x ^ 2))
@scalar_rule acsc(x) -(inv(x ^ 2 * sqrt(1 - x ^ -2)))
@scalar_rule acsc(x::Real) -(inv(abs(x) * sqrt(x ^ 2 - 1)))
@scalar_rule asec(x) inv(x ^ 2 * sqrt(1 - x ^ -2))
@scalar_rule asec(x::Real) inv(abs(x) * sqrt(x ^ 2 - 1))

@scalar_rule cosd(x) -deg2rad(sind(x))
@scalar_rule cospi(x) -π * sinpi(x)
@scalar_rule sind(x) deg2rad(cosd(x))
@scalar_rule sinpi(x) π * cospi(x)
@scalar_rule tand(x) deg2rad(1 + Ω ^ 2)

@scalar_rule sinc(x) cosc(x)

# the position of the minus sign below warrants the correct type for π
if VERSION ≥ v"1.6"
    @scalar_rule sincospi(x) @setup((sinpix, cospix) = Ω) (π * cospix)  (π * (-sinpix))
end

@scalar_rule(
    clamp(x, low, high),
    @setup(
        islow = x < low,
        ishigh = high < x,
    ),
    (!(islow | ishigh), islow, ishigh),
)
@scalar_rule x \ y (-(Ω / x), one(y) / x)

function frule((_, ẏ), ::typeof(identity), x)
    return (x, ẏ)
end

function rrule(::typeof(identity), x)
    function identity_pullback(ȳ)
        return (NoTangent(), ȳ)
    end
    return (x, identity_pullback)
end

# rouding related,
# we use `zero` rather than `ZeroTangent()` for scalar, and avoids issues with map etc
@scalar_rule round(x) zero(x)
@scalar_rule floor(x) zero(x)
@scalar_rule ceil(x) zero(x)

# `literal_pow`
# This is mostly handled by AD; it's a micro-optimisation to provide a gradient for x*x*x
# Note that rules for `^` are defined in the fastmath_able.jl

function frule((_, _, Δx, _), ::typeof(Base.literal_pow), ::typeof(^), x::Real, ::Val{2})
    return x * x, 2 * x * Δx
end
function frule((_, _, Δx, _), ::typeof(Base.literal_pow), ::typeof(^), x::Real, ::Val{3})
    x2 = x * x
    return x2 * x, 3 * x2 * Δx
end

function rrule(::typeof(Base.literal_pow), ::typeof(^), x::Real, ::Val{2})
    square_pullback(dy) = (NoTangent(), NoTangent(), ProjectTo(x)(2 * x * dy), NoTangent())
    return x * x, square_pullback
end
function rrule(::typeof(Base.literal_pow), ::typeof(^), x::Real, ::Val{3})
    x2 = x * x
    cube_pullback(dy) = (NoTangent(), NoTangent(), ProjectTo(x)(3 * x2 * dy), NoTangent())
    return x2 * x, cube_pullback
end


text:
"""
    paged_mango_query(db::Database, selector::Selector;
      fields::Vector{T} where T<:AbstractString          = Vector{AbstractString}(),
      sort::Vector{Dict{T, Any}} where T<:AbstractString = Vector{Dict{AbstractString, Any}}(),
      pagesize                                           = 100)

Perform multiple HTTP requests against a JSON-type index producing the
intermediate results. This is a wrapper around `query()` using the
`skip` and `limit` parameters.

### Examples

    for page in Channel(paged_mango_query(db, q"data = ..."; pagesize=10))
      for doc in page.docs
        # ...
      end
    end
"""
function paged_mango_query(db::Database, selector::Selector;
  fields::Vector{T} where T<:AbstractString          = Vector{AbstractString}(),
  sort::Vector{Dict{T, Any}} where T<:AbstractString = Vector{Dict{AbstractString, Any}}(),
  pagesize                                           = 100)

  skip = 0
  function producer(chan::Channel)
    while true
      result = mango_query(db, selector; fields=fields, sort=sort, skip=skip, limit=pagesize)
      skip += pagesize
      if length(result.docs) == 0
        break
      end
      put!(chan, result)
    end
  end

end

max_stars_repo_name: xpqz/couchzilla

max_issues_repo_name: xpqz/couchzilla

max_forks_repo_name: xpqz/couchzilla


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function geo2eci(lat::T, long::T, alt::T, jd::T) where {T<:AbstractFloat}
    Re    = planets["earth"].eqradius / 1000
    sin_lat, cos_lat = sincos(deg2rad(lat))
    long  = deg2rad(long)
    gst   = ct2lst(zero(T), jd)
    sid_angle = gst*pi/12 # Sidereal angle.
    theta = long + sid_angle # Azimuth
    altRe = alt + Re
    r     = altRe * cos_lat
    sin_theta, cos_theta = sincos(theta)
    return r * cos_theta, r * sin_theta, altRe * sin_lat
end

"""
    geo2eci(latitude, longitude, altitude, jd) -> x, y, z

### Purpose ###

Convert geographic spherical coordinates to Earth-centered inertial coordinates.

### Explanation ###

Converts from geographic spherical coordinates (latitude, longitude, altitude)
to ECI (Earth-Centered Inertial) (x, y, z) rectangular coordinates.  Julian days
is also needed.

Geographic coordinates assume the Earth is a perfect sphere, with radius equal
to its equatorial radius.  ECI coordinates are in km from Earth center at epoch
TOD (True of Date).

### Arguments ###

* `latitude`: geographic latitude, in degrees.
* `longitude`: geographic longitude, in degrees.
* `altitude`: geographic altitude, in kilometers.
* `jd`: Julian days.

The three coordinates can be passed as a 3-tuple `(latitude, longitude,
altitude)`.  In addition, `latitude`, `longitude`, `altitude`, and `jd` can be
given as arrays of the same length.

### Output ###

The 3-tuple of ECI (x, y, z) coordinates, in kilometers.  The TOD epoch is the
supplied `jd` time.

If geographical coordinates are given as arrays, a 3-tuple of arrays of the same
length is returned.

### Example ###

Obtain the ECI coordinates of the intersection of the equator and Greenwich's
meridian on 2015-06-30T14:03:12.857

```jldoctest
julia> using AstroLib

julia> geo2eci(0, 0, 0, jdcnv("2015-06-30T14:03:12.857"))
(-4024.8671780315185, 4947.835465127513, 0.0)
```

### Notes ###

`eci2geo` converts Earth-centered inertial coordinates to geographic spherical
coordinates.

Code of this function is based on IDL Astronomy User's Library.
"""
geo2eci(lat::Real, long::Real, alt::Real, jd::Real) =
    geo2eci(promote(float(lat), float(long), float(alt), float(jd))...)

geo2eci(lla::Tuple{Real, Real, Real}, jd::Real) =
    geo2eci(lla..., jd)

function geo2eci(lat::AbstractArray{LA}, long::AbstractArray{<:Real},
                 alt::AbstractArray{<:Real},
                 jd::AbstractArray{<:Real}) where {LA<:Real}
    @assert length(lat) == length(long) == length(alt) == length(jd)
    typela = float(LA)
    x = similar(lat, typela)
    y = similar(lat, typela)
    z = similar(lat, typela)
    for i in eachindex(lat)
        x[i], y[i], z[i] =
            geo2eci(lat[i], long[i], alt[i], jd[i])
    end
    return x, y, z
end


text:
using POMDPs
using ContinuousPOMDPTreeSearchExperiments
using POMDPToolbox
using POMCP
using POMCPOW
using LightDarkPOMDPs
using ProgressMeter
using PmapProgressMeter
using ParticleFilters
using JLD

N = 100

@everywhere begin
    using LightDarkPOMDPs
    using POMDPs
    using ContinuousPOMDPTreeSearchExperiments
    using POMDPToolbox
    using POMCP
    using POMCPOW
    using ParticleFilters

    pomdp = LightDark2DTarget(term_radius=0.1,
                              init_dist=SymmetricNormal2([2.0, 2.0], 5.0),
                              discount=0.95)
    p_rng = MersenneTwister(4)
    solvers = Dict{String, Union{Solver,Policy}}(
        "heuristic" => begin
            homing_policy = SimpleFeedback(gain=1.0, max_radius=10.0)
            exploring_up = ObsAdaptiveParticleFilter(pomdp, LowVarianceResampler(1000), 0.05, p_rng)
            homing_up = ObsAdaptiveParticleFilter(pomdp, LowVarianceResampler(1000), 0.0, p_rng)
            InfoGatherHeur(5.0, 0.03, homing_policy, exploring_up, homing_up)
        end,

        "vanilla_pomcp" => begin
            rng3 = copy(p_rng)
            action_gen = AdaptiveRadiusRandom(max_radius=6.0, to_zero_first=true, to_light_second=true, rng=rng3)
            rollout_policy = RadiusRandom(radius=10.0, rng=rng3)
            tree_queries = 50_000
            solver = POMCPDPWSolver(next_action=action_gen,
                                    tree_queries=tree_queries,
                                    c=5.0,
                                    max_depth=40,
                                    k_action=10.0,
                                    alpha_action=1/8,
                                    k_observation=4.0,
                                    alpha_observation=1/8,
                                    estimate_value=RolloutEstimator(rollout_policy),
                                    default_action=(b, ex)->Vec2(0.1, 0.1),
                                    rng=rng3
                                   )
        end,

        "modified_pomcp" => begin
            rng3 = copy(p_rng)
            action_gen = AdaptiveRadiusRandom(max_radius=6.0, to_zero_first=true, to_light_second=true, rng=rng3)
            rollout_policy = SimpleFeedback(max_radius=10.0)
            tree_queries = 50_000
            node_updater = ObsAdaptiveParticleFilter(pomdp, LowVarianceResampler(100), 0.05, rng3)
            # node_updater = SimpleParticleFilter(pomdp, LowVarianceResampler(100), rng=rng3)
            solver = POMCPDPWSolver(next_action=action_gen,
                                    tree_queries=tree_queries,
                                    c=5.0,
                                    max_depth=40,
                                    k_action=10.0,
                                    alpha_action=1/8,
                                    k_observation=4.0,
                                    alpha_observation=1/8,
                                    estimate_value=PORolloutEstimator(rollout_policy, node_updater),
                                    node_belief_updater=node_updater,
                                    default_action=(b, ex)->Vec2(0.1, 0.1),
                                    rng=rng3
                                   )
        end,

        "pomcpow" => begin
            rng3 = copy(p_rng)
            action_gen = AdaptiveRadiusRandom(max_radius=6.0, to_zero_first=true, to_light_second=true, rng=rng3)
            tree_queries = 1_000_000
            node_updater = ObsAdaptiveParticleFilter(pomdp, LowVarianceResampler(100), 0.05, rng3)
            solver = POMCPOWSolver(next_action=action_gen,
                                    tree_queries=tree_queries,
                                    criterion=MaxUCB(5.0),
                                    final_criterion=MaxTries(),
                                    max_depth=40,
                                    k_action=10.0,
                                    alpha_action=1/8,
                                    k_observation=4.0,
                                    alpha_observation=1/8,
                                    estimate_value=OneStepValue(),
                                    node_belief_updater=node_updater,
                                    rng=rng3
                                   )
        end,

        "greedy" => SimpleFeedback(gain=1.0, max_radius=10.0)
    )

    up_rng = MersenneTwister(5)
    standard_up = ObsAdaptiveParticleFilter(pomdp, LowVarianceResampler(100_000), 0.05, up_rng)
    updaters = Dict{String, Updater}(
        "vanilla_pomcp" => standard_up,
        "modified_pomcp" => standard_up,
        "greedy" => standard_up,
        "pomcpow" => standard_up,
        "heuristic" => updater(solvers["heuristic"])
    )
end

solver_keys = keys(solvers)
rewards = Dict{String, AbstractVector{Float64}}()

for (j, sk) in enumerate(solver_keys)
    s_rewards = SharedArray(Float64, N)
    prog = Progress(N, desc="$sk ($j of $(length(solver_keys)))...")
    pmap(prog, 1:N) do i
        sim_rng = MersenneTwister(i)
        if isa(solvers[sk], Policy)
            policy = deepcopy(solvers[sk])
        else
            policy = solve(solvers[sk], pomdp)
        end
        up = deepcopy(updaters[sk])
        sim = RolloutSimulator(max_steps=40, rng=sim_rng)
        s_rewards[i] = simulate(sim, pomdp, policy, up)
    end
    rewards[sk] = sdata(s_rewards)
end

for k in solver_keys
    println("$k mean: $(mean(rewards[k])) sem: $(std(rewards[k])/sqrt(N))")
end

filename = Pkg.dir("ContinuousPOMDPTreeSearchExperiments", "data", "compare_$(Dates.format(now(), "E_d_u_HH_MM")).jld")
# @save(filename, pomdp, solver_keys, solvers, rewards, updaters)


text:
@doc raw"""
    Ray2

An object `r` of the data type [`Ray2`](@ref) is a directed straight ray in the
two-dimensional Euclidean plane ``Ε^2``.

It starts in a point called the *source* of `r` and goes to infinity.
"""
Ray2

"""
    Ray2(p::Point2, q::Point2)

Introduces a ray `r` with source `p` and passing through point `q`.
"""
Ray2(p::Point2, q::Point2)

"""
    Ray2(p::Point2, d::Direction2)

Introduces a ray `r` starting at source `p` with direction `d`.
"""
Ray2(p::Point2, d::Direction2)

"""
    Ray2(p::Point2, v::Vector2)

Introduces a ray `r` starting at source `p` with the direction of `v`.
"""
Ray2(p::Point2, v::Vector2)

"""
    Ray2(p::Point2, l::Line2)

Introduces a ray `r` starting at source `p` with the same direction as `l`.
"""
Ray2(p::Point2, l::Line2)

"""
    ==(r::Ray2, h::Ray2)

Test for equality: two rays are equal, iff they have the same source and the
same direction.
"""
==(r::Ray2, h::Ray2)

"""
    source(r::Ray2)

Returns the source of `r`.
"""
source(r::Ray2)

@doc raw"""
    point(r::Ray2, i::Real)

Returns a point on `r`.

!!! info "Precondition"

    ``i ≥ 0``.
"""
point(r::Ray2, i::Real)

@cxxdereference point(r::Ray2, i::Real) = point(r, convert(FT, i))

"""
    direction(r::Ray2)

Returns the direction of `r`.
"""
direction(r::Ray2)

"""
    to_vector(r::Ray2)

Returns a vector giving the direction of `r`.
"""
to_vector(r::Ray2)

"""
    supporting_line(r::Ray2)

Returns the line supporting `r` which has the same direction.
"""
supporting_line(r::Ray2)

"""
    opposite(r::Ray2)

Returns the ray with the same source and opposite direction.
"""
opposite(r::Ray2)

"""
    is_degenerate(r::Ray2)

Ray `r` is degenerate, if the source and the second defining point fall
together (that is if the direction is degenerate).
"""
is_degenerate(r::Ray2)

"""
    is_horizontal(r::Ray2)

Ray `r` is horizontal, if coefficient `a` of the supporting line's equation is
zero.
"""
is_horizontal(r::Ray2)

"""
    is_vertical(r::Ray2)

Ray `r` is vertical, if coefficient `b` of the supporting line's equation is
zero.
"""
is_vertical(r::Ray2)

"""
    has_on(r::Ray2, p::Point2)

A point is on `r`, iff it is equal to the source of `r`, or if it is in the
interior of `r`.
"""
has_on(r::Ray2, p::Point2)

"""
    collinear_has_on(r::Ray2, p::Point2)

Checks if point `p` is on `r`.

!!! info "Precondition"

    `p` is on the supporting line of `r`.
"""
collinear_has_on(r::Ray2, p::Point2)

"""
    transform(r::Ray2, t::AffTransformation2)

Returns the ray obtained by applying `t` on the source and on the direction of
`r`.
"""
transform(r::Ray2, t::AffTransformation2)


text:
module TestCaret
using TSML
using TSMLextra
using Test

const IRIS = getiris()
const X = IRIS[:,1:4] |> Matrix
const Y = IRIS[:,5] |> Vector

const XX = IRIS[:,1:1] |> Matrix
const YY = IRIS[:,4] |> Vector

#const learners=["rf","treebag","svmRadialWeights","pls","svmLinear","bagFDA","rpart"]
const learners=["rf"]

function test_caret_fit(learner::String)
    #crt = CaretLearner(Dict(:learner=>learner,:fitControl=>"trainControl(method='cv')"))
    @info learner
    crt = CaretLearner(Dict(:learner=>learner))
    fit!(crt,X,Y)
    @test crt.model != nothing
end

function test_caret_transform(learner::String)
    #crt = CaretLearner(Dict(:learner=>learner,:fitControl=>"trainControl(method='cv')"))
    @info learner
    crt = CaretLearner(Dict(:learner=>learner))
    fit!(crt,X,Y)
    @test sum(transform!(crt,X) .== Y)/length(Y) > 0.10
end

if !Sys.islinux() # bug in travis using old linux and R versions
    @testset "caret training classifiers" begin
	for lrn in learners
	    test_caret_fit(lrn)
	end
    end

    @testset "caret prediction classifiers" begin
	for lrn in learners
	    test_caret_transform(lrn)
	end
    end
end

end


text:
# This file was generated, do not modify it. # hide
figure(figsize=(12, 6))
stem(res)
xlim(0, length(res))

savefig(joinpath(@OUTPUT, "residuals.png")) # hide


text:
# TODO
# - transparency
# - fancy sources
# - text layout attributes
# - ellipses, bezier curves or similar
# - actions that take lists of points (or generators?)
# - other shapes
# - align consistent
# - asymmetric scaling of text?
# - rotation only as a keyword arg for shapes?
# - axes shifting separate from scale / shift?

module Drawing

import Cairo; const X = Cairo
import Graphics; const G = Graphics
import Colors; const C = Colors
import Tk; const T = Tk

export DrawingError,
       current_context,
       with, cairo, draw, paint,
       cm, mm, inch, pts, rad, deg,
       PNG, PDF, TK, Paper, Axes, Pen, Ink, Scale, Translate, Font,
       move, line, circle, text, rectangle, square,
       print_fonts

# add additional calls to cairo
include("cairo.jl")

# --- utilities

type DrawingError <: Exception
    msg::AbstractString
end
Base.showerror(io::IO, e::DrawingError) = print(io, e.msg)

# allow sizes to be given as 30cm etc (convert to mm)
const cm = 10
const mm = 1
const inch = 25.4
const pts = inch/72

# similarly for angles (convert to degrees)
const rad = 180/pi
const deg = 1

int(x) = round(Int, x)

# lookup value from table, using lower case string as key
function lookup(name, table, value::AbstractString)
    v = lowercase(value)
    if haskey(table, v)
        table[v]
    else
        throw(DrawingError("Unknown $(name): $(value)"))
    end
end

parse_color(c::AbstractString) = parse(C.Colorant, c)
parse_color(c::C.Colorant) = c

const RED = parse_color("red")
const GREEN = parse_color("green")
const BLUE = parse_color("blue")
const WHITE = parse_color("white")
const BLACK = parse_color("black")

# --- global / thread-local state

# this is tracked via the stage variable
typealias Context Nullable{X.CairoContext}

const SCOPE_NONE = 0
const SCOPE_INNER = 1
const SCOPE_OUTER_DRAW = 2
const SCOPE_OUTER_PAINT = 3
const SCOPE_OUTER_OTHER = 3

const STAGE_NONE = 0     # context uninitialized
const STAGE_OUTPUT = 1   # create context (and on exit save file)
const STAGE_PAPER = 2    # background colour
const STAGE_AXES = 3     # axis scaling
const STAGE_DRAW = 4     # general drawing

type FontState
    fd::FontDescription
    align::Int
    size::Float64        # not applied directly as transform dependent
    FontState() = new("sans", 1, -1)
    FontState(f::FontState) = new(copy(f.fd), f.align, f.size)
end

# this needs to be in a thread local variable once such things exist in julia
type ThreadContext
    context::Context
    font::Vector{FontState}  # separate from context
    stage::Int           # used to track context creation
    scope::Vector{Int}   # used to track scoping rules
    ThreadContext() = new(Context(), FontState[FontState()], STAGE_NONE, Int[SCOPE_NONE])
end

const thread_context = ThreadContext()
current_context() = get(thread_context.context)

# --- scoping infrastructure

# before and after are functions of (ThreadContext, Attribute) -> nothing
type Attribute
    name::AbstractString
    stage::Int
    before::Vector{Function}
    after::Vector{Function}
    Attribute(name, stage, before, after) = new(name, stage, before, after)
end

const NO_ACTIONS = Function[]
ctx(f) = (c, a) -> f(get(c.context))

function save_once(ctx, saved)
    # save context if it exists
    if !saved && ctx.stage > STAGE_NONE
        X.save(get(ctx.context))
        true
    else
        saved
    end
end

function add_defaults(stage, attributes)
    i = 1
    while i <= length(attributes) && stage <= attributes[end].stage
        if attributes[i].stage == stage
            stage += 1
            i += 1
        elseif attributes[i].stage > stage
            splice!(attributes, i:i-1, [DEFAULT_ATTRIBUTES[stage]])
            stage += 1
        else
            i += 1
        end
    end
    while stage < STAGE_DRAW
        push!(attributes, DEFAULT_ATTRIBUTES[stage])
        stage += 1
    end
end

function make_scope(verify, before, after)

    function scope(f, attributes::Attribute...)

        c = thread_context
        saved = false
        pushed = false

        try

            saved = save_once(c, saved)

            # stable sort, respecting user where possible
            a = sort!([attributes...], by=a -> a.stage, alg=InsertionSort)
            add_defaults(c.stage+1, a)

            # check for conflicts and push scope
            scope = verify(c, a)
            push!(c.scope, scope)
            push!(c.font, FontState(c.font[end]))
            pushed = true

            for attribute in a

                if min(c.stage+1, STAGE_DRAW) != attribute.stage
                    throw(DrawingError("Initialization attributes must be in outermost scope."))
                end

                for x in attribute.before
                    x(c, attribute)
                end
                c.stage = attribute.stage
                saved = save_once(c, saved)
            end

            args = before(get(c.context))

            args == nothing ? f() : f(args...)

            after(get(c.context))

            for attribute in reverse(a)
                for x in attribute.after
                    x(c, attribute)
                end
            end

        finally

            if pushed
                pop!(c.scope)
                pop!(c.font)
            end

            if c.scope[end] == SCOPE_NONE
                c.context = Context()
                c.stage = STAGE_NONE
            else
                saved && X.restore(get(c.context))
            end

        end
    end
end

# --- actual scopes

function preserve_current_point(action)
    function (c)
        has_point = has_current_point(c)
        x, y = get_current_point(c)
        action(c)
        if has_point
            X.move_to(c, x, y)
        end
    end
end

stroke = preserve_current_point(X.stroke)
fill = preserve_current_point(X.fill)

function make_verify(scope)
    function verify(c, a)
        if c.scope[end] > SCOPE_INNER
            throw(DrawingError("Cannot nest a scope inside an action scope"))
        end
        scope
    end
end

INACTIVE = c -> nothing
with = make_scope(make_verify(SCOPE_INNER), INACTIVE, INACTIVE)
cairo = make_scope(make_verify(SCOPE_OUTER_OTHER), c -> [c], INACTIVE)
draw = make_scope(make_verify(SCOPE_OUTER_DRAW), INACTIVE, stroke)
paint = make_scope(make_verify(SCOPE_OUTER_PAINT), INACTIVE, fill)

# --- output (file, display, etc)

# these should be as portrait, in mm (x, y)
const PAPER_SIZES = Dict("a0" => [841, 1189],
                         "a1" => [594, 841],
                         "a2" => [420, 594],
                         "a3" => [297, 420],
                         "a4" => [210, 297],
                         "a5" => [148, 210],
                         "a6" => [105, 148],
                         "a7" => [74, 105],
                         "a8" => [52, 74],
                         "a9" => [37, 52],
                         "a10" => [26, 37],
                         "letter" => [216, 279],
                         "legal" => [216, 356],
                         "junior" => [127, 203],
                         "ledger" => [279, 432])
parse_paper_size(size::AbstractString) = lookup("paper size", PAPER_SIZES, size)

const LANDSCAPE = 1
const PORTRAIT = 2
const ORIENTATIONS = Dict("landscape" => LANDSCAPE,
                          "porttrait" => PORTRAIT)
parse_orientation(orientation::AbstractString) = lookup("orientation", ORIENTATIONS, orientation)
parse_orientation(orientation::Integer) = orientation

# set these at the start so thing look pretty
function with_defaults(before)
    vcat(before,
         Pen(join="round", cap="round").before,
         Font("sans").before)
end

# width and heigh are in mm
function PDF(path, width_mm, height_mm)
    Attribute("PDF", STAGE_OUTPUT,
              with_defaults([(c, a) -> c.context = X.CairoContext(X.CairoPDFSurface(path, width_mm/pts, height_mm/pts))]),
              [ctx(c -> X.finish(c.surface)),
               ctx(c -> X.destroy(c))])
end

function PDF(path; size="a4", orientation=LANDSCAPE)
    w, h = parse_paper_size(size)
    w, h = parse_orientation(orientation) == LANDSCAPE ? (h, w) : (w, h)
    PDF(path, w, h)
end

function PNG(path, width_px::Int, height_px::Int)
    Attribute("PNG", STAGE_OUTPUT,
              with_defaults([(c, a) -> c.context = X.CairoContext(X.CairoRGBSurface(width_px, height_px))]),
              [ctx(c -> X.write_to_png(c.surface, path)),
               (c,a) -> println("wrote $(path)"),
               ctx(c -> X.destroy(c))])
end

function press_return()
    println("Press RETURN to close window")
    readline(STDIN)
end

function TK(width_px, height_px; name="Drawing", destroy=true)
    window, canvas = nothing, nothing
    function create(c, a)
        window = Tk.Toplevel(name, width_px, height_px)
        canvas = Tk.Canvas(window)
        Tk.pack(canvas, expand=true, fill="both")
        c.context = Tk.getgc(canvas)
    end
    function close(c, a)
        Tk.reveal(canvas)
        Tk.update()
        if destroy
            press_return()
            Tk.destroy(window)
        end
    end
    Attribute("TK", STAGE_OUTPUT,
              with_defaults([create]),
              [close])
end

# --- paper (background colour)

function Paper(background="white")
    bg = parse_color(background)
    Attribute("Paper", STAGE_PAPER,
              [ctx(c -> set_background(c, bg))],
              NO_ACTIONS)
end

function set_background(c, bg)
    X.save(c)
    X.rectangle(c, 0, 0, X.width(c.surface), X.height(c.surface))
    X.set_source(c, bg)
    X.fill(c)
    X.restore(c)
end

# --- axes (coordinate scaling)

function Axes(; scale=1, border=0.1, centred=true, negative=false)
    Attribute("Axes", STAGE_AXES,
              [ctx(c -> set_coords(c, scale, border, centred, negative)),
               # start at origin
               ctx(c -> X.move_to(c, 0, 0))],
              NO_ACTIONS)
end

function set_coords(c, scale, border, centred, negative)

    if isa(scale, Number)
        scale = (scale, scale)
    end

    nx, ny = X.width(c.surface), X.height(c.surface)
    portrait = nx < ny
    sm, lg = portrait ? (nx, ny) : (ny, nx)
    b = border * sm
    u = sm - 2 * b

    # used (non-border) lengths
    ux = nx - 2 * b
    uy = ny - 2 * b
    sc = max(scale[1] / ux, scale[2] / uy)

    if centred
        if negative
            G.set_coords(c, 0, 0, nx, ny, -sc*nx, sc*nx, sc*ny, -sc*ny)
        else
            dx, dy = scale[1]-sc*ux, scale[2]-sc*uy
            G.set_coords(c, 0, 0, nx, ny, -sc*b+dx/2, sc*(nx-b)+dx/2,
                         sc*(ny-b)+dy/2, -sc*b+dy/2)
        end
    else
        if negative
            dx, dy = scale[1]-sc*ux, scale[2]-sc*uy
            G.set_coords(c, 0, 0, nx, ny, -sc*nx-dx, sc*nx-dx,
                         sc*ny-dy, -sc*ny-dy)
        else
            G.set_coords(c, 0, 0, nx, ny, -sc*b, sc*(nx-b), sc*(ny-b), -sc*b)
        end
    end
end

# --- source/stroke attributes

function Ink(foreground)
    f = parse_color(foreground)
    Attribute("Ink", STAGE_DRAW, [ctx(c -> X.set_source(c, f))], NO_ACTIONS)
end

const LINE_CAPS = Dict("butt" => X.CAIRO_LINE_CAP_BUTT,
                       "round" => X.CAIRO_LINE_CAP_ROUND,
                       "square" => X.CAIRO_LINE_CAP_SQUARE)
parse_line_cap(cap::AbstractString) = lookup("line cap", LINE_CAPS, cap)
parse_Line_cap(cap::Integer) = cap

const LINE_JOINS = Dict("mitre" => X.CAIRO_LINE_JOIN_MITER,
                        "miter" => X.CAIRO_LINE_JOIN_MITER,
                        "round" => X.CAIRO_LINE_JOIN_ROUND,
                        "bevel" => X.CAIRO_LINE_JOIN_BEVEL)
parse_line_join(join::AbstractString) = lookup("line join", LINE_JOINS, join)
parse_line_join(join::Integer) = join

function Pen(width; cap=nothing, join=nothing)
    Attribute("Pen", STAGE_DRAW,
              vcat(width >= 0 ? [ctx(c -> set_width(c, width))] : [],
                   cap != nothing ? [ctx(c -> X.set_line_cap(c, parse_line_cap(cap)))] : [],
                   join != nothing ? [ctx(c -> X.set_line_join(c, parse_line_join(join)))] : []),
              NO_ACTIONS)
end

Pen(; cap=X.CAIRO_LINE_CAP_ROUND, join=X.CAIRO_LINE_JOIN_ROUND) = Pen(-1; cap=cap, join=join)

function set_width(c, width)
    # width is in user coords, so scale to device coords
    v = [width, width]
    X.user_to_device_distance!(c, v)
    v = map(abs, v)
    X.set_line_width(c, sum(v)/2)
end

# --- transforms

# this is much, much more limited that raw cairo.  and may be reduced further,
# so that all you can do is specify new axes on some sub-region.

# instead, individual commands support rotation and scaling in varius ways.

Scale(factor) = Attribute("Scale", STAGE_DRAW, [ctx(c -> X.scale(c, factor, factor))], NO_ACTIONS)
Translate(x, y) = Attribute("Translate", STAGE_DRAW, [ctx(c -> X.translate(c, x, y))], NO_ACTIONS)

# --- (pango) text attributes

function to_font(fd::FontDescription)
    s = string(fd)
    "Font(\"$s\")"
end

function print_fonts()
    fm = get_font_map_default()
    ff = list_families(fm)
    sort!(ff, by=(f-> get_name(f)))
    for f in ff
        m = is_monospace(f) ? "[monospace]" : ""
        println("\n  $(get_name(f)) $m")
        fc = list_faces(f)
        for c in fc
            d = describe(c)
            @printf("  %20s %s\n", "$(get_name(c)):", to_font(d))
        end
    end
end

const STYLES = Dict("normal" => PANGO_STYLE_NORMAL,
                    "italic" => PANGO_STYLE_ITALIC,
                    "oblique" => PANGO_STYLE_OBLIQUE)
parse_style(style::AbstractString) = lookup("font style", STYLES, style)
parse_style(style::Integer) = style

const VARIANTS = Dict("normal" => PANGO_VARIANT_NORMAL,
                      "smallcaps" => PANGO_VARIANT_SMALL_CAPS)
parse_variant(variant::AbstractString) = lookup("font variant", VARIANTS, variant)
parse_variant(variant::Integer) = variant

const WEIGHTS = Dict("ultralight" => PANGO_WEIGHT_ULTRALIGHT,
                     "light" => PANGO_WEIGHT_LIGHT,
                     "normal" => PANGO_WEIGHT_NORMAL,
                     "bold" => PANGO_WEIGHT_BOLD,
                     "ultrabold" => PANGO_WEIGHT_ULTRABOLD,
                     "heavy" => PANGO_WEIGHT_HEAVY)
parse_weight(weight::AbstractString) = lookup("font weight", WEIGHTS, weight)
parse_weight(weight::Integer) = weight

const STRETCHES = Dict("ultracondensed" => PANGO_STRETCH_ULTRA_CONDENSED,
                       "extracondensed" => PANGO_STRETCH_EXTRA_CONDENSED,
                       "condensed" => PANGO_STRETCH_CONDENSED,
                       "semicondensed" => PANGO_STRETCH_SEMI_CONDENSED,
                       "normal" => PANGO_STRETCH_NORMAL,
                       "semiexpanded" => PANGO_STRETCH_SEMI_EXPANDED,
                       "expanded" => PANGO_STRETCH_EXPANDED,
                       "extraexpanded" => PANGO_STRETCH_EXTRA_EXPANDED,
                       "ultraexpanded" => PANGO_STRETCH_ULTRA_EXPANDED)
parse_stretch(stretch::AbstractString) = lookup("font stretch", STRETCHES, stretch)
parse_stretch(stretch::Integer) = stretch

const GRAVITY = Dict("auto" => PANGO_GRAVITY_AUTO,
                     "north" => PANGO_GRAVITY_NORTH,
                     "east" => PANGO_GRAVITY_EAST,
                     "south" => PANGO_GRAVITY_SOUTH,
                     "west" => PANGO_GRAVITY_WEST)
parse_gravity(gravity::AbstractString) = lookup("font gravity", GRAVITY, gravity)
parse_gravity(gravity::Integer) = gravity

set_via(setter, value) = (c, a) -> setter(c.font[end].fd, value)

function Font(desc; size=nothing, style=nothing, variant=nothing, weight=nothing, stretch=nothing, gravity=nothing)
    Attribute("Font", STAGE_DRAW,
              vcat(desc != nothing ? [(c, a) -> c.font[end].fd = desc] : [],
                   size != nothing ? [(c, a) -> c.font[end].size = size] : [],
                   style != nothing ? [set_via(set_style, parse_style(style))] : [],
                   variant != nothing ? [set_via(set_variant, parse_variant(variant))] : [],
                   weight != nothing ? [set_via(set_weight, parse_weight(weight))] : [],
                   stretch != nothing ? [set_via(set_stretch, parse_stretch(stretch))] : [],
                   gravity != nothing ? [set_via(set_gravity, parse_gravity(gravity))] : []),
              NO_ACTIONS)
end

Font(; size=nothing, style=nothing, variant=nothing, weight=nothing, stretch=nothing, gravity=nothing) = Font(nothing; size=size, style=style, weight=weight, stretch=stretch, gravity=gravity)

# --- paths

macro lift(new, old)
    :($(esc(new))(args...) = $(esc(old))(current_context(), args...))
end

@lift(move, X.move_to)
@lift(line, X.line_to)

# if (x, y) is the bottom left point of the box, then this returns an
# (x,y) that should be the new bottom-left point, such that the box
# appears aligned correctly relative to the given point(!)
# to make things more exciting, height should be negated when using the
# 'y goes down' convention (ie from pango), and then everything is
# relative to top left.
function alignment(x, y, align, width, height)
    xalign = x - ((align - 1) % 3) * width / 2
    if height > 0
        yalign = y - height + round((align - 2) / 3) * height / 2
    else
        yalign = y + round((align - 2) / 3) * height / 2
    end
    xalign, yalign
end

# TODO - ellipse + rotate
# TODO - more useful arcs joining lines etc
function circle(radius; from=0, to=360)
    c = current_context()
    x, y = get_current_point(c)
    X.new_sub_path(c)
    X.arc(c, x, y, radius, deg2rad(from), deg2rad(to))
    X.move_to(c, x, y)  # don't change current point
end

function get_current_scale(c)
    zero = X.device_to_user(c, 0, 0)
    one = X.device_to_user(c, 1, 1)
    dx, dy = one[1] - zero[1], one[2] - zero[2]
    sqrt(dx*dx + dy*dy) / sqrt(2)
end

function text(s; align=1, rotate=0)

    t = thread_context
    f = t.font[end]
    c = get(t.context)
    l = Layout(c.layout)    # our pango routines use wrapper

    # todo - maybe provide a flag to skip this transform and let the
    # user suffer whatever self-inflicted pain they want?
    sc = get_current_scale(c)
    X.save(c)
    X.reset_transform(c)
    X.rotate(c, -deg2rad(rotate))
    x, y = get_current_point(c)

    try

        update_layout(c, l)
        if f.size > 0
            set_absolute_size(f.fd, f.size / sc)
        end
        set_description(l, f.fd)

        set_text(l, s)

        ink, log = get_pixel_extents(l)
        # -height because y axis is inverted
        X.move_to(c, alignment(x-log.x, y-log.y, align, log.width, -log.height)...)

        show_path(c, l)

    finally
        X.move_to(c, x, y)
        X.restore(c)
    end
end

function rectangle(width, height; align=1, rotate=0)
    c = current_context()
    x, y = get_current_point(c)
    X.save(c)
    try
        a, b = alignment(x, y, align, width, height)
        X.rotate(c, deg2rad(rotate))
        X.move_to(c, a, b)
        X.line_to(c, a+width, b)
        X.line_to(c, a+width, b+height)
        X.line_to(c, a, b+height)
        X.line_to(c, a, b)
    finally
        X.restore(c)
        X.move_to(c, x, y)
    end
end

function square(side; align=1, rotate=0)
    rectangle(side, side; align=align, rotate=rotate)
end

# TODO - more shapes, including general polyogns and triangls.

# --- defaults during startup (working through the stages)

DEFAULT_ATTRIBUTES = Dict(STAGE_OUTPUT => TK(300, 200),
                          STAGE_PAPER => Paper(WHITE),
                          STAGE_AXES => Axes(),)

end


text:
#
# Helper module for Remote Terminal application. It
# contains shorthands (macroses) for more simple
# remote server configuring and requests.
#
# @author DeadbraiN
# TODO: add somehow REPL help tips for these functions
module Terminal
  import RpcApi
  import Config
  import Client
  import Helper
  import Connection
  import Creature

  # TODO: check it this module contains all functions from ManagerRpc
  export init
  export termGetRegion
  export termCreateOrganisms
  export termCreateOrganism
  export termSetConfig
  export termGetConfig
  export termSetQuiet
  export termMutate
  export termGetIps
  export termGetOrganism
  export termGetAmount
  export termGetOrganisms
  export termSetEnergy
  export termSetRandomEnergy
  export termBackup
  export termGetStatistics
  export termGetBest
  export termMarkOrganism
  #
  # Contains connection object for shorthands
  #
  type TerminalData
  	con::Client.ClientConnection
    cfg::Config.ConfigData

  	TerminalData() = new()
    TerminalData(con::Client.ClientConnection, cfg::Config.ConfigData) = new(con, cfg)
  end
  #
  # This function should be called first before any macros
  # will be called. It sets connection object as a global.
  # @param con Connection object
  #
  function init(con::Client.ClientConnection)
    _termData.con = con
    _termData.cfg = Config.create()
  end
  #
  # Shorthand for RpcApi.RPC_GET_REGION. See ManagerRpc.getRegion
  # for details
  #
  function termGetRegion(x::Int = 1, y::Int = 1, x1::Int = 0, y1::Int = 0)
    Client.request(_termData.con, RpcApi.RPC_GET_REGION, x, y, x1, y1)
  end
  #
  # Shorthand for RpcApi.RPC_CREATE_ORGANISMS. See
  # ManagerRpc.createOrganisms for details
  #
  function termCreateOrganisms()
    Client.request(_termData.con, RpcApi.RPC_CREATE_ORGANISMS)
  end
  #
  # Shorthand for RpcApi.RPC_CREATE_ORGANISM. See
  # ManagerRpc.createOrganism for details
  #
  function termCreateOrganism(pos::Helper.Point = Helper.Point(0, 0))
    Client.request(_termData.con, RpcApi.RPC_CREATE_ORGANISM, pos)
  end
  #
  # Shorthand for RpcApi.RPC_SET_CONFIG. See
  # ManagerRpc.setConfig for details
  #
  function termSetConfig(name::Symbol, value::Any)
    Client.request(_termData.con, RpcApi.RPC_SET_CONFIG, name, value)
  end
  #
  # Shorthand for RpcApi.RPC_GET_CONFIG. See
  # ManagerRpc.getConfig for details
  #
  function termGetConfig(name::Symbol)
    Client.request(_termData.con, RpcApi.RPC_GET_CONFIG, name)
  end
  #
  # Shorthand for RpcApi.RPC_SET_QUIET. See
  # ManagerRpc.setQuiet for details
  #
  function termSetQuiet(mode::Bool = true)
    Client.request(_termData.con, RpcApi.RPC_SET_QUIET, mode)
  end
  #
  # Shorthand for RpcApi.RPC_MUTATE. See
  # ManagerRpc.mutate for details
  #
  function termMutate(orgId::UInt, amount::Int = 1)
    Client.request(_termData.con, RpcApi.RPC_MUTATE, orgId, amount)
  end
  #
  # Shorthand for RpcApi.RPC_GET_IPS. See
  # ManagerRpc.getIps for details
  #
  function termGetIps()
    Client.request(_termData.con, RpcApi.RPC_GET_IPS)
  end
  #
  # Shorthand for RpcApi.RPC_GET_ORGANISM. See
  # ManagerRpc.getOrganism for details
  #
  function termGetOrganism(id::UInt)
    Client.request(_termData.con, RpcApi.RPC_GET_ORGANISM, id)
  end
  #
  # Shorthand for RpcApi.RPC_GET_AMOUNT. See
  # ManagerRpc.getAmount for details
  #
  function termGetAmount()
    Client.request(_termData.con, RpcApi.RPC_GET_AMOUNT)
  end
  #
  # Shorthand for RpcApi.RPC_GET_ORGANISMS. See
  # ManagerRpc.getOrganisms for details
  #
  function termGetOrganisms()
    Client.request(_termData.con, RpcApi.RPC_GET_ORGANISMS)
  end
  #
  # Shorthand for RpcApi.RPC_SET_ENERGY. See
  # ManagerRpc.setEnergy for details
  #
  function termSetEnergy(x::Int, y::Int, energy::UInt16)
    Client.request(_termData.con, RpcApi.RPC_SET_ENERGY, x, y, energy)
  end
  #
  # Shorthand for RpcApi.RPC_SET_ENERGY_RND. See
  # ManagerRpc.setEnergy for details
  #
  function termSetRandomEnergy(amount::Int = _termData.cfg.worldStartEnergyDots, energy::UInt16 = _termData.cfg.worldStartEnergyInDot)
    Client.request(_termData.con, RpcApi.RPC_SET_ENERGY_RND, amount, energy)
  end
  #
  # Shorthand for RpcApi.RPC_BACKUP. See
  # ManagerRpc.setEnergy for details
  #
  function termBackup()
    Client.request(_termData.con, RpcApi.RPC_BACKUP)
  end
  #
  # Shorthand for RpcApi.RPC_GET_STATISTICS. See
  # ManagerRpc.getStatistics for details
  #
  function termGetStatistics()
    Client.request(_termData.con, RpcApi.RPC_GET_STATISTICS)
  end
  #
  # Returns best organisms (with maximum energy). See
  # ManagerRpc.getBest for details
  #
  function termGetBest(amount::Int)
    Client.request(_termData.con, RpcApi.RPC_GET_BEST, amount)
  end
  #
  # Marks specified organism by color. See
  # ManagerRpc.markOrganism for details
  #
  function termMarkOrganism(orgId::UInt32, colorIndex::Int)
    Client.request(_termData.con, RpcApi.RPC_MARK_ORGANISM, orgId, colorIndex)
  end
  #
  # This object is created only for passing connection
  # data to macroses.
  #
  global _termData = TerminalData()
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
mutable struct Future <: AbstractFuture
    value::Any
    value_id::ValueId
    mutated::Bool
    stale::Bool

    function Future(value::Any, value_id::ValueId, mutated::Bool, stale::Bool)
        new_future = new(value, value_id, mutated, stale)

        # Create finalizer and register
        finalizer(new_future) do fut
            try
                record_request(DestroyRequest(fut.value_id))
            catch e
                # `record_request` will fail if there isn't any job to add the
                # request to. So we just continue silently.
                # @warn "Failed to destroy value $(fut.value_id) because job has stopped: $e"
            end
        end

        new_future
    end
end

isview(f::F) where F <: AbstractFuture = false

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
@testset "SurvivalSignature" begin

    A = zeros(6, 6)
    A[1, [2, 3]] .= 1.0
    A[2, [4, 5]] .= 1.0
    A[3, [4, 6]] .= 1.0
    A[4, [5, 6]] .= 1.0

    types = Dict(1 => [1, 2, 5], 2 => [3, 4, 6])

    φ = s_t_connectivity([1:6;], [1], [5, 6])

    @testset "exactentry" begin
        @test SurvivalSignature.exactentry(CartesianIndex(1, 1), A, types, φ) == 0.0
        @test SurvivalSignature.exactentry(CartesianIndex(1, 2), A, types, φ) == 0.0
        @test SurvivalSignature.exactentry(CartesianIndex(1, 3), A, types, φ) == 0.0
        @test SurvivalSignature.exactentry(CartesianIndex(1, 4), A, types, φ) == 0.0
        @test SurvivalSignature.exactentry(CartesianIndex(2, 1), A, types, φ) == 0.0
        @test SurvivalSignature.exactentry(CartesianIndex(2, 2), A, types, φ) == 0.0
        @test SurvivalSignature.exactentry(CartesianIndex(2, 3), A, types, φ) == 1.0 / 9.0
        @test SurvivalSignature.exactentry(CartesianIndex(2, 4), A, types, φ) == 3.0 / 9.0
        @test SurvivalSignature.exactentry(CartesianIndex(3, 1), A, types, φ) == 0.0
        @test SurvivalSignature.exactentry(CartesianIndex(3, 2), A, types, φ) == 0.0
        @test SurvivalSignature.exactentry(CartesianIndex(3, 3), A, types, φ) == 4.0 / 9.0
        @test SurvivalSignature.exactentry(CartesianIndex(3, 4), A, types, φ) == 6.0 / 9.0
        @test SurvivalSignature.exactentry(CartesianIndex(4, 1), A, types, φ) == 1.0
        @test SurvivalSignature.exactentry(CartesianIndex(4, 2), A, types, φ) == 1.0
        @test SurvivalSignature.exactentry(CartesianIndex(4, 3), A, types, φ) == 1.0
        @test SurvivalSignature.exactentry(CartesianIndex(4, 4), A, types, φ) == 1.0
    end
end


text:
"""
xxx
# Examples
```@example

```
"""
struct LTISource
  position::Vector{Float64}
  transmission ::Function
end

# Method
function (𝚽::LTISource)(𝛏₀::Vector{Float64}, t₀::Float64)
   𝐩ₛ, p = 𝚽.position, 𝚽.transmission
   delay = distBetween(𝐩ₛ,𝛏₀)/lightSpeed
   return A(delay) * p(t₀-delay)
end

# DISPLAY
Base.show(io::IO, x::LTISource) = print(io, "LTI omnidirectional source")


text:
using LinearAlgebra
mutable struct Symmetric_KKT_solver <: abstract_KKT_system_solver
    # abstract_KKT_system_solver
    ls_solver::abstract_linear_system_solver
    factor_it::Class_iterate
    delta_x_vec::Array{Float64,1}
    delta_s_vec::Array{Float64,1}
    rhs::System_rhs
    dir::Class_point
    kkt_err_norm::Class_kkt_error
    rhs_norm::Float64
    pars::Class_parameters
    schur_diag::Array{Float64,1}
    true_x_diag::Array{Float64,1}

    ready::Symbol #

    # Symmetric_KKT_solver only
    #true_diag::Array{Float64,1}
    Q::SparseMatrixCSC{Float64,Int64}
    #M::SparseMatrixCSC{Float64,Int64}
    #K::SparseMatrixCSC{Float64,Int64}
    #current_it::Class_iterate
    #reduct_factors::Class_reduction_factors

    function Symmetric_KKT_solver()
        ## reworks system to form of
        ## [[H A']]
        ## [[A D]].
        ##
      return new()
    end
end

function form_system!(kkt_solver::Symmetric_KKT_solver, iter::Class_iterate, timer::class_advanced_timer)
    start_advanced_timer(timer, "symmetric/form_system");
    J = get_jac(iter) #.cache.J
    J_T = get_jac_T(iter) #.cache.J_T
    #B = spdiagm(-iter.point.s ./ iter.point.y);
    B = sparse(Diagonal(-iter.point.s ./ iter.point.y))
    M = [[get_lag_hess(iter) J_T]; [J B]];

    pause_advanced_timer(timer, "symmetric/form_system");

    start_advanced_timer(timer, "symmetric/allocate");
    kkt_solver.Q = M;
    kkt_solver.factor_it = iter;
    kkt_solver.schur_diag = compute_schur_diag(iter)
    kkt_solver.true_x_diag = diag(M)[1:dim(iter)]
    kkt_solver.ready = :system_formed
    pause_advanced_timer(timer, "symmetric/allocate");

end

function factor_implementation!(kkt_solver::Symmetric_KKT_solver, timer::class_advanced_timer)
    return ls_factor!(kkt_solver.ls_solver, kkt_solver.Q, dim(kkt_solver.factor_it), ncon(kkt_solver.factor_it), timer)
end

function compute_direction_implementation!(kkt_solver::Symmetric_KKT_solver, timer::class_advanced_timer)
    start_advanced_timer(timer, "symmetric")
    y_org = get_y(kkt_solver.factor_it)
    s_org = get_s(kkt_solver.factor_it)
    rhs = kkt_solver.rhs;

    symmetric_rhs = [rhs.dual_r; rhs.primal_r + rhs.comp_r ./ y_org];
    dir_x_and_y = ls_solve(kkt_solver.ls_solver, symmetric_rhs, timer);

    shift_vector = [ones(dim(kkt_solver.factor_it));
    zeros(ncon(kkt_solver.factor_it))] * get_delta(kkt_solver.factor_it);

    dir = kkt_solver.dir;
    dir.x = dir_x_and_y[1:length(rhs.dual_r)];
    dir.y = -dir_x_and_y[(1+length(rhs.dual_r)):end];
    #dir.s = ( rhs.comp_r - dir.y .* s_org ) ./ y_org
    dir.s = eval_jac_prod(kkt_solver.factor_it, dir.x) - rhs.primal_r

    check_for_nan(dir)

    start_advanced_timer(timer, "symmetric/kkt_err");
    update_kkt_error!(kkt_solver, Inf, timer)
    pause_advanced_timer(timer, "symmetric/kkt_err");
    pause_advanced_timer(timer, "symmetric")
end

function update_delta_vecs!(kkt_solver::Symmetric_KKT_solver, delta_x_vec::Array{Float64,1}, delta_s_vec::Array{Float64,1}, timer::class_advanced_timer)
    start_advanced_timer(timer, "symmetric")
    start_advanced_timer(timer, "symmetric/delta_vecs")
    kkt_solver.delta_x_vec = delta_x_vec
    kkt_solver.delta_s_vec = delta_s_vec

    if sum(abs.(delta_s_vec)) > 0.0
        error("not implemented")
    else
        for i = 1:length(delta_x_vec)
            kkt_solver.Q[i,i] = kkt_solver.true_x_diag[i] + delta_x_vec[i]
        end
    end

    kkt_solver.ready = :delta_updated
    pause_advanced_timer(timer, "symmetric/delta_vecs")
    pause_advanced_timer(timer, "symmetric")
end


text:
# Autogenerated wrapper script for ONNX_jll for aarch64-linux-musl-cxx11
export libonnx, libonnxifi, libonnxifi_dummy

using protoc_jll
JLLWrappers.@generate_wrapper_header("ONNX")
JLLWrappers.@declare_file_product(libonnx)
JLLWrappers.@declare_library_product(libonnxifi, "libonnxifi.so")
JLLWrappers.@declare_library_product(libonnxifi_dummy, "libonnxifi_dummy.so")
function __init__()
    JLLWrappers.@generate_init_header(protoc_jll)
    JLLWrappers.@init_file_product(
        libonnx,
        "lib/libonnx.a",
    )

    JLLWrappers.@init_library_product(
        libonnxifi,
        "lib/libonnxifi.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libonnxifi_dummy,
        "lib/libonnxifi_dummy.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Return start of commn substring in s1, start of common substring in s2, and length of substring
# Indexes refer to character number, not index (differ for Unicode strings)
function longest_common_substring(s1::AbstractString, s2::AbstractString)
    if length(s1) > length(s2)
        start2, start1, len = longest_common_substring(s2, s1)
    else
        start1, start2, len = 0, 0, 0
        p = zeros(Int, length(s2))
        i1 = 0
        for ch1 in s1
            i1 += 1
            i2 = 0
            oldp = 0
            for ch2 in s2
                i2 += 1
                newp = 0
                if ch1 == ch2
                    newp = oldp > 0 ? oldp : i2
                    currentlength = (i2 - newp + 1)
                    if currentlength > len
                        start1, start2, len = i1 - currentlength + 1, newp, currentlength
                    end
                end
                p[i2], oldp = newp, p[i2]
            end
        end
    end
    return start1, start2, len
end

function matching_blocks!(x::Set{Tuple{Int, Int, Int}}, s1::AbstractString, s2::AbstractString, start1::Integer, start2::Integer)
    a = longest_common_substring(s1, s2)
    if a[3] > 0
        push!(x, (a[1] + start1 - 1, a[2] + start2 - 1, a[3]))
        s1before = SubString(s1, firstindex(s1), nextind(s1, 0, a[1]) - 1)
        s2before = SubString(s2, firstindex(s2), nextind(s2, 0, a[2]) - 1)
        matching_blocks!(x, s1before, s2before, start1, start2)
        if ((a[1] + a[3]) <= lastindex(s1)) & ((a[2] + a[3]) <= lastindex(s2))
            s1after = SubString(s1, nextind(s1, 0, a[1] + a[3]), lastindex(s1))
            s2after = SubString(s2, nextind(s2, 0, a[2] + a[3]), lastindex(s2))
            matching_blocks!(x, s1after, s2after, start1 + a[1] + a[3] - 1, start2 + a[2] + a[3] - 1)
        end
    end
end

function matching_blocks(s1::AbstractString, s2::AbstractString)
    x = Set{Tuple{Int, Int, Int}}()
    matching_blocks!(x, s1, s2, 1, 1)
    return x
end

struct RatcliffObershelp <: PreMetric end

function evaluate(dist::RatcliffObershelp, s1::AbstractString, s2::AbstractString)
    n_matched = 0
    for x in matching_blocks(s1, s2)
        n_matched += x[3]
    end
    1.0 - 2 * n_matched / (length(s1) + length(s2))
end


text:
# Bhattacharyya distances. Much like for KLDivergence we assume the vectors to
# be compared are probability distributions, frequencies or counts rather than
# vectors of samples. Pre-calc accordingly if you have samples.

struct BhattacharyyaDist <: SemiMetric end

struct HellingerDist <: Metric end

# Bhattacharyya coefficient

function bhattacharyya_coeff(a, b)
    n = length(a)
    if n != length(b)
        throw(DimensionMismatch("first argument has length $n which does not match the length of the second, $(length(b))."))
    end
    sqab, asum, bsum = _bhattacharyya_coeff(a, b)
    # We must normalize since we cannot assume that the vectors are normalized to probability vectors.
    return sqab / sqrt(asum * bsum)
end

@inline function _bhattacharyya_coeff(a, b)
    Ta = _eltype(a)
    Tb = _eltype(b)
    T = typeof(sqrt(zero(promote_type(Ta, Tb))))
    sqab = zero(T)
    asum = zero(Ta)
    bsum = zero(Tb)

    for (ai, bi) in zip(a, b)
        sqab += sqrt(ai * bi)
        asum += ai
        bsum += bi
    end
    return sqab, asum, bsum
end
@inline function _bhattacharyya_coeff(a::AbstractVector{Ta}, b::AbstractVector{Tb}) where {Ta<:Number,Tb<:Number}
    T = typeof(sqrt(oneunit(Ta)*oneunit(Tb)))
    sqab = zero(T)
    asum = zero(Ta)
    bsum = zero(Tb)

    @simd for i in eachindex(a, b)
        @inbounds ai = a[i]
        @inbounds bi = b[i]
        sqab += sqrt(ai * bi)
        asum += ai
        bsum += bi
    end
    return sqab, asum, bsum
end

# Faster pair- and column-wise versions TBD...

# Bhattacharyya distance
(::BhattacharyyaDist)(a, b) = -log(bhattacharyya_coeff(a, b))
const bhattacharyya = BhattacharyyaDist()

# Hellinger distance
(::HellingerDist)(a, b) = sqrt(1 - bhattacharyya_coeff(a, b))
const hellinger = HellingerDist()


text:
"""
    FelsensteinFunction(tree_postorder::Vector{N}, pi_::T, rates::Vector{Float64}, data::Array{Float64,3}, n_c::Int64, blv::Vector{Float64}) where {T<:Real, N<:GeneralNode}

This function calculates the log-likelihood of an evolutiuonary model using the
Felsensteins pruning algorithm.
The function is written such that it is differentiable by Zygote 0.5.3.

Returns log-likelihood as a Real number.

* `tree_postorder` : Vector of Nodes.

* `pi` : Real number used in calculation.

* `rates` : Vector of Floats, used in calculation.

* `data` : Array of Floats, currently a placeholder.

* `n_c` : Int64, nchar value derived from tree and dataset.

* `blv` : Vector of Floats, branchlength vector derived from tree.
"""
function FelsensteinFunction(tree_postorder::Vector{N}, pi_::Array{Float64}, rates::Float64,
                             data::Array{Float64,3}, n_c::Int64, blv::Array{Float64,1})::Float64 where {N<:GeneralNode}

    mu::Float64 =  1.0 / (2.0 * prod(pi_))
    mutationArray::Vector{Array{Float64, 2}} = calc_trans.(blv, pi_[1], mu, rates)
    root_node::N = last(tree_postorder)
    rns::Array{Float64} = zeros(Float64, 1, n_c)

    @views for node in tree_postorder
        if node.nchild > 0
            res = node_loop(node, mutationArray)
            if !node.root
                scaler::Array{Float64} = maximum(res, dims=1)
                rns = rns .+ log.(scaler)
                res = res ./ scaler
            end #if
            node.data = res
        end #if
    end # for

    return likelihood_root(root_node.data, pi_, rns)
end # function

@inline function likelihood_root(root::Array{Float64,2}, pi_::Array{Float64}, rns::Array{Float64,2})::Float64
    sum(log.(sum(root .* pi_, dims=1))+rns)
end

function node_loop(node::N, mutationArray::Vector{Array{Float64, 2}})::Array{Float64,2} where {N<:GeneralNode, S<:Real}
    # creating a new array is necessary because Zygote can not differentiate otherwise
    out = ones(size(node.data))
    @inbounds @views for child in node.children
            out = out .* BLAS.gemm('N','N',mutationArray[child.num], child.data)
    end
    out
end

function calc_trans(time::Float64, pi_::S, mu::Float64, r::Float64)::Array{Float64,2} where {S<:Real}

    v = exp(-r*time*mu)
    v1 = pi_ - (pi_ - 1)*v
    v2 = pi_ - pi_* v
    v3 = (pi_ - 1)*(v - 1)
    v4 = pi_*(v - 1) + 1
    mm::Array{Float64}=[v1 v3;
                        v2 v4]
    return mm
end

function calc_trans(time::Float64, pi_::Array{Float64}, mu::Float64, r::Float64)::Array{Float64,2}
    res::Array{Float64} = repeat(reshape(pi_, (1,2)), 2)
    di = Diagonal(ones(2))
    v = exp(-r*time*mu)
    res = res .*(1.0-v)
    res .+ (di .* v)
end

"""
    FelsensteinFunction(tree_postorder::Vector{N}, pi_::Array{Float64}, rates::Array{Float64},
                         U::Array{Float64,2}, D::Array{Float64}, Uinv::Array{Float64,2},
                         data::Array{Float64,4}, c_grad::Bool = true) where {N<:GeneralNode}

This function calculates the log-likelihood of an evolutiuonary model using the
Felsensteins pruning algorithm. If `c_grad` equals `true` (default) the analytic gradient
regarding the branch lengths of the tree gets computed as well.
"""
function FelsensteinFunction(tree_postorder::Vector{N}, pi_::Array{Float64}, rates::Array{Float64},
                     U::Array{M,2}, D::Array{M}, Uinv::Array{M,2}, mu::Float64,
                     data::Array{Float64,4}, c_grad::Bool = true) where {N<:GeneralNode, M<:Number}
    Nbases, Nsites, Nrates, Nnodes = size(data)
    mutationArray::Array{Float64,4} = Array{Float64,4}(undef, Nbases, Nbases, Nrates, Nnodes-1)
    mutationArray .= 0.0
    grv::Vector{Float64} = Vector{Float64}(undef, Nnodes-1)
    Down::Array{Float64,4} = similar(data)
    pre_order_partial::Array{Float64,4} = similar(data)
    ll = fels_ll(tree_postorder, data, D, U, Uinv, rates, mu, Nrates, Nsites, Down, pi_, mutationArray)
    if c_grad
        grv = fels_grad(tree_postorder, data, D, U, Uinv, rates, mu, Nrates,
                        Nsites, Nnodes, Down, pi_, mutationArray, pre_order_partial)
    end

    return ll, grv, data, Down, pre_order_partial

end

function fels_grad(tree_postorder::Vector{N}, data::Array{Float64,4},
         D::Array{M,1}, U::Array{M,2}, Uinv::Array{M,2},
         rates::Array{Float64,1}, mu::Float64, Nrates::Int64, Nsites::Int64, Nnodes::Int64, Down::Array{Float64,4},
         pi_::Array{Float64}, mutationArray::Array{Float64,4}, pre_order_partial::Array{Float64,4})::Vector{Float64} where {N <: GeneralNode, M <: Number}

    root_node::N = last(tree_postorder)
    pre_order_partial[:, :, :, root_node.num] .= pi_
    scaler::Array{Float64, 2} = Array{Float64,2}(undef, 1, Nsites)
    gradi::Array{Float64, 2} = Array{Float64,2}(undef, 1, Nsites)
    grv::Vector{Float64} = Vector{Float64}(undef, Nnodes-1)
    ptg::Array{Float64,2} = similar(mutationArray[:, :, 1, 1])

    @inbounds @views for node in reverse(tree_postorder)[2:end]
           mother::N = node.mother
           @inbounds @views for r in 1:Nrates
               pre_order_partial[:, :, r, node.num] .= pre_order_partial[:, :, r, mother.num]
               @inbounds @views for child in mother.children
                   if child.num != node.num
                       pre_order_partial[:, :, r, node.num] .*= Down[:, :, r, child.num]
                   end
               end
           end
           tg::Float64 = 0.0
           @inbounds @views for r in 1:Nrates
               BLAS.gemm!('N', 'N', 1.0, BLAS.symm('R', 'L', (D .* (rates[r] * mu)) .* diagm(exp.(D .* (rates[r]*mu*node.inc_length))), U), Uinv, 0.0, ptg)
               gradi = sum(pre_order_partial[:, :, r, node.num] .* BLAS.gemm('N', 'N', ptg , data[:, :, r, node.num]), dims=1)
               pre_order_partial[:, :, r, node.num] .= BLAS.gemm('T', 'N', mutationArray[:, :, r, node.num] , pre_order_partial[:, :, r, node.num])
               tg += sum(gradi ./ sum(pre_order_partial[:, :, r, node.num] .* data[:, :, r, node.num], dims=1))
               if node.nchild > 0
                   scaler = maximum(pre_order_partial[:, :, r, node.num], dims = 1)
                   pre_order_partial[:, :, r, node.num] ./= scaler
               end
           end
           grv[node.num] = tg
    end # for
    grv
end

function fels_ll(tree_postorder::Vector{N}, data::Array{Float64,4},
         D::Array{Float64,1}, U::Array{Float64,2}, Uinv::Array{Float64,2},
         rates::Array{Float64,1}, mu::Float64, Nrates::Int64, Nsites::Int64, Down::Array{Float64,4},
         pi_::Array{Float64}, mutationArray::Array{Float64,4})::Float64 where {N <: GeneralNode}

    scaler::Array{Float64, 2} = Array{Float64,2}(undef, 1, Nsites)
    ll::Float64 = 0.0
    root_node::N = last(tree_postorder)
    @inbounds @views for node in tree_postorder
        if node.nchild > 0
            data[:, :, :, node.num] .= 1.0
            for child in node.children
                @inbounds @views for r in 1:Nrates
                    BLAS.gemm!('N', 'N', 1.0, BLAS.symm('R', 'L', diagm(exp.(D .* (rates[r]*mu*child.inc_length))), U), Uinv, 0.0, mutationArray[:, :, r, child.num])
                    BLAS.gemm!('N','N', 1.0, mutationArray[:, :, r, child.num], data[:, :, r, child.num], 0.0, Down[:, :, r, child.num])
                    data[:, :, r, node.num] .*= Down[:, :, r, child.num]

                end

            end

            if !node.root
                @inbounds @views for r in 1:Nrates
                    scaler .= maximum(data[:, :, r, node.num], dims=1)
                    data[:, :, r, node.num] ./= scaler
                    ll += sum(log.(scaler))
                end
            end
        end #if
    end # for

    @inbounds @views @simd for r in 1:Nrates
        ll += sum(log.(sum(data[:, :, r, root_node.num] .* pi_, dims=1)))
    end
    ll
end

function fels_ll(tree_postorder::Vector{N}, data::Array{Float64,4},
         D::Array{ComplexF64,1}, U::Array{ComplexF64,2}, Uinv::Array{ComplexF64,2},
         rates::Array{Float64,1}, mu::Float64, Nrates::Int64, Nsites::Int64, Down::Array{Float64,4},
         pi_::Array{Float64}, mutationArray::Array{Float64,4})::Float64 where {N <: GeneralNode}

    scaler::Array{Float64, 2} = Array{Float64,2}(undef, 1, Nsites)
    ll::Float64 = 0.0
    root_node::N = last(tree_postorder)
    @inbounds @views for node in tree_postorder
        if node.nchild > 0
            data[:, :, :, node.num] .= 1.0
            for child in node.children
                @inbounds @views for r in 1:Nrates
                    mutationArray[:, :, r, child.num] .= abs.(BLAS.gemm('N', 'N', one(ComplexF64), BLAS.symm('R', 'L', diagm(exp.(D .* (rates[r]*mu*child.inc_length))), U), Uinv))
                    BLAS.gemm!('N','N', 1.0, mutationArray[:, :, r, child.num], data[:, :, r, child.num], 0.0, Down[:, :, r, child.num])
                    data[:, :, r, node.num] .*= Down[:, :, r, child.num]

                end

            end

            if !node.root
                @inbounds @views for r in 1:Nrates
                    scaler .= maximum(data[:, :, r, node.num], dims=1)
                    data[:, :, r, node.num] ./= scaler
                    ll += sum(log.(scaler))
                end
            end
        end #if
    end # for

    @inbounds @views @simd for r in 1:Nrates
        ll += sum(log.(sum(data[:, :, r, root_node.num] .* pi_, dims=1)))
    end
    ll
end


text:
@testset "Infinity Arithmetic" begin
    posinf = ECC.Infinity(1)
    neginf = ECC.Infinity(-1)
    @testset "Equal" begin
        @test posinf != 3
        @test posinf != -7
        @test neginf != 0
    end
    @testset "Addition" begin
        @test +posinf == posinf
        @test +neginf == neginf
        @test posinf + 3 == posinf
        @test posinf + 3 + 12 + 1 == posinf
        @test neginf + 2 == neginf
        @test neginf + 112 - 7 == neginf
    end
    @testset "Substraction" begin
        @test -posinf == neginf
        @test -neginf == posinf
        @test posinf - 3 == posinf
        @test posinf - 3 - 12 - 1 == posinf
        @test neginf - 2 == neginf
        @test neginf - 112 - 7 == neginf
    end
    @testset "Multiplication" begin
        @test 0 * posinf == 0
        @test 3 * posinf == posinf
        @test -3 * posinf == neginf
        @test 0 * neginf == 0
        @test 7 * neginf == neginf
        @test -2 * neginf == posinf
    end
    @testset "Exponentiation" begin
        @test posinf^0  == posinf
        @test neginf^0  == neginf
        @test posinf^2  == posinf
        @test neginf^5  == neginf
        @test posinf^-7 == posinf
        @test neginf^-1 == neginf
    end
end


text:
# Performance tests for InterestRates module

println("Running perftests...")

using InterestRates

dt_curve = Date(2015,08,08)

vert_x = [11, 15, 19, 23, 25, 200, 500]
vert_y = [0.10, 0.12, 0.20, 0.19, 0.21, 0.21, 0.22]
NS_params = [0.1, 0.2, 0.3, 0.5]
SVEN_params = [0.1, 0.2, 0.3, 0.4, 0.5, 0.8]

days_vec = [Dates.Day(i) for i=1:1500]

mat_vec = dt_curve + days_vec

curve_linear = InterestRates.IRCurve("dummy-linear", InterestRates.Actual360(),
	InterestRates.ExponentialCompounding(), InterestRates.Linear(), dt_curve,
	vert_x, vert_y)

curve_flatforward = InterestRates.IRCurve("dummy-flatforward", InterestRates.Actual360(),
	InterestRates.ExponentialCompounding(), InterestRates.FlatForward(), dt_curve,
	vert_x, vert_y)

curve_ns = InterestRates.IRCurve("dummy-ns", InterestRates.Actual360(),
	InterestRates.ExponentialCompounding(), InterestRates.NelsonSiegel(), dt_curve,
	NS_params)

curve_sven = InterestRates.IRCurve("dummy-sven", InterestRates.Actual360(),
	InterestRates.ExponentialCompounding(), InterestRates.Svensson(), dt_curve,
	SVEN_params)

curve_spline_rates = InterestRates.IRCurve("dummy-spline_rates", InterestRates.Actual360(),
	InterestRates.ExponentialCompounding(), InterestRates.CubicSplineOnRates(), dt_curve,
	vert_x, vert_y)

curve_spline_df = InterestRates.IRCurve("dummy-spline_df", InterestRates.Actual360(),
	InterestRates.ExponentialCompounding(), InterestRates.CubicSplineOnDiscountFactors(), dt_curve,
	vert_x, vert_y)

curve_step = InterestRates.IRCurve("dummy-step", InterestRates.Actual360(),
	InterestRates.ExponentialCompounding(), InterestRates.StepFunction(), dt_curve,
	vert_x, vert_y)

c_array = [curve_linear, curve_flatforward, curve_ns, curve_sven, curve_spline_rates, curve_spline_df, curve_step]

# Warm up

for c in c_array
	zero_rate(c, mat_vec)
	ERF(c, mat_vec)
	discountfactor(c, mat_vec)
end

sp = InterestRates.splinefit(vert_x, vert_y)
InterestRates.splineint(sp, convert(Vector{Int}, 1:30))

# Perftests

for c in c_array
	println("$(InterestRates.curve_get_method(c))")
	@time for i=1:1000 zero_rate(c, mat_vec) end
	@time for i=1:1000 ERF(c, mat_vec) end
	@time for i=1:1000 discountfactor(c, mat_vec) end
end

println("Perftests end")


text:
# To analyze the result file
function getlas(ARGS)
    fname = ARGS[1] # filename
    trains = Any[]; devs = Any[]; isdone=false;
    for line in eachline(fname)
        if contains(line, "acc")
            toparse = split(line)
            if length(toparse) < 9
                d0 = parse(toparse[end]); t0 = parse(toparse[end-2])
                if length(ARGS) > 1
                    println("Initial tracc $t0 devacc $d0")
                end
            else
                tracc = parse(toparse[end-3]); push!(trains, tracc)
                devacc = parse(toparse[end]); push!(devs, devacc)
            end
        end
        contains(line, "done") && (isdone=true)
    end
    if isempty(trains)
        println("No enough epoch"); return;
    end
    trmax, devmax = findmax(trains), findmax(devs)
    if isdone
        println("tr$(trmax[2]) $(trmax[1]) | dev$(devmax[2]) $(devmax[1]) DONE.")
    else
        println("tr$(trmax[2]) $(trmax[1]) | dev$(devmax[2]) $(devmax[1]) :epoch$(length(trains)+1)")
    end
end
!isinteractive() && getlas(ARGS)


text:
using LinearAlgebra
using Test

# 1
# a
A = [
   3 -1 -1  0 -1
  -1  3  0 -1 -1
  -1  0  3 -1 -1
   0 -1 -1  3 -1
  -1 -1 -1 -1  4
]

AᵀA = A' * A

# b
@test tr(AᵀA) == 68

@test rank(AᵀA) == 4

values, = eigen(AᵀA)

@test round(sum(values)) == 68

@test round(reduce(*, values)) == 0

# c
A₁ = [
   1 -1 0 0 0
  -1  1 0 0 0
   0  0 0 0 0
   0  0 0 0 0
   0  0 0 0 0
]

A₁ᵀA₁ = A₁' * A₁

@test A₁ᵀA₁ == [
  2  -2  0  0  0
 -2   2  0  0  0
  0   0  0  0  0
  0   0  0  0  0
  0   0  0  0  0
]

#d
b = ones(8, 1)

@test_throws SingularException inv(A)


text:
using MLStyle

@as_record Moment
@as_record CircuitGate
@as_record ControlledGate
@as_record HadamardGate
@as_record XGate
@as_record YGate
@as_record ZGate
@as_record SwapGate
@as_record RxGate
@as_record RyGate
@as_record RzGate
@as_record SGate
@as_record TGate
@as_record SdagGate
@as_record TdagGate

# function trans_bitarray(ba::BitArray{1})
#     name_register = "creg creg" * string(classical_register_count) * "[" * string(length(ba)) * "];\n"
#     global classical_register_count += 1
#     return name_register
# end

function trans_g(g::AbstractGate)
    @match g begin
        HadamardGate() => "h"
        XGate() => "x"
        YGate() => "y"
        ZGate() => "z"
        SGate() => "s"
        SdagGate() => "sdg"
        TGate() => "t"
        TdagGate() => "tdg"
        SwapGate() => "swap"
        RxGate(θ=θ) => "rx(" * string(θ[1]) * ")"
        RyGate(θ=θ) => "ry(" * string(θ[1]) * ")"
        RzGate(θ=θ) => "rz(" * string(θ[1]) * ")"
        ControlledGate(U=U
            ) => let U=trans_g(U)
                    return "c" * U
                end
    end
end

integer_pattern = r"[0-9]+$"
function trans_cg(cg::CircuitGate{M,G}) where {M,G}
    @match cg begin
        CircuitGate(iwire=iwire,
            gate=gate
            ) => let id = trans_g(gate)
                    return string(id) * " " * join( "qregister[" .* string.(reverse(iwire.-1)) .* "]" ,",") * ";"
                end
    end
end

function trans_moment(m::Moment)
    @match m begin
        Moment(
            gates=gates
        ) => let gates=gates
                expr = join(trans_cg.(gates), "\n")
                return expr
            end
    end
end

"""
    cgc2qasm(c::Circuit{N}) where {N}

converts Circuit{N} object to OpenQASM 2.0 representation
"""
function cgc2qasm(c::Circuit{N}) where {N}
    global classical_register_count = 1
    header = """OPENQASM 2.0;\n
    include "qelib1.inc";\n
    """
    header *= "qreg qregister[" * string(N) * "];\n"
    header = header .* join(trans_moment.(c.moments), "\n")
    return header
end


text:
export CG_HS

function CG_HS(nlp :: AbstractNLPModel;
               x   :: Vector{T}=copy(nlp.meta.x0),
               stp :: NLPStopping = NLPStopping(nlp,
                                                NLPAtX(nlp.meta.x0)),
               scaling   :: Bool = true,
               LS_algo   :: Function = bracket,
               LS_logger :: AbstractLogger = Logging.NullLogger(),
               kwargs...) where T

    return CG_generic(nlp;
                      x = x,
                      stp=stp,
                      scaling = scaling,
                      LS_algo = LS_algo,
                      LS_logger = LS_logger,
                      CG_formula  = formula_HS,
                      kwargs...)
end

max_stars_repo_name: vepiteski/LSDescent

max_issues_repo_name: vepiteski/LSDescent

max_forks_repo_name: vepiteski/LSDescent


text:
# Note: this file contains ideas/palns for future interface

# cases:
# sphere : Sphere or Vector(Sphere)
# wavelenght: single or vector

# =====================================================================================================

# Structure
s = Sphere(radius, epsilon)

# epsilon_media: real, positive
sol = solve(s, x, epsilon_media, MieAuto())
# sol = solve(s, wavelenght, epsilon_media, MieAuto()) ??
sol = solve(s, x, Rayleigh())
sol = solve(s, x, MieLarge())
sol = solve(s, x, MieConductive())

qsca = scattering(sol)
qext = extinction(sol)

# =====================================================================================================
# Structure
s = Sphere(radius, epsilon)

# epsilon_media: real, positive
model = MieAuto(s, wavelenght, epsilon_media) # ??

# model = Rayleigh(s, x)
# model = MieLarge(s, x)
# model = MieConductive(s, x)

qsca = scattering(model)
qext = extinction(model)

# =====================================================================================================
# Structure
s = Sphere(radius, epsilon)

# epsilon_media: real, positive
m, x = mx(s, wavelenght, epsilon_media)

model = Mie(m, x)

qsca = scattering(model)
qext = scattering(model)

max_stars_repo_name: fekad/jlmie

max_issues_repo_name: fekad/jlmie

max_forks_repo_name: fekad/jlmie


text:
@testset "footnotes" begin
    set_curpath("index.md")
    st = """
        A[^1] B[^blah] C
        """
    @test isapproxstr(st |> seval, """

A[1]
              B[2]
              C
""")

    st = """
        A[^1] B[^blah]
        C
        [^1]: first footnote
        [^blah]: second footnote
        """
    @test isapproxstr(st |> seval, """

                A
                [1]
                B[2]
                C

                    [1]
                    first footnote

                        [2]
                        second footnote

""")
end

@testset "Fn in code" begin
    s = raw"""
        ```markdown
       this has[^1]

       [^1]: def
       ```
       blah
       """
    @test isapproxstr(s |> fd2html_td, """

this has[^1]

        [^1]: def

blah
""")
end


text:
mutable struct SpanningTreeLightGraphsPrimSolver <: SpanningTreeSolver
  graph::SimpleGraph
  weights_matrix::Matrix{Float64}

  function SpanningTreeLightGraphsPrimSolver()
    return new(Graph(0), zeros(0, 0))
  end
end

has_lp_formulation(::SpanningTreeLightGraphsPrimSolver) = false
approximation_ratio(::SpanningTreeLightGraphsPrimSolver) = 1.0
approximation_term(::SpanningTreeLightGraphsPrimSolver) = 0.0
approximation_ratio_budgeted(::SpanningTreeLightGraphsPrimSolver) = 1.0
approximation_term_budgeted(::SpanningTreeLightGraphsPrimSolver) = 0.0
supports_solve_budgeted_linear(::SpanningTreeLightGraphsPrimSolver) = false
supports_solve_all_budgeted_linear(::SpanningTreeLightGraphsPrimSolver) = false

function build!(solver::SpanningTreeLightGraphsPrimSolver, graph::SimpleGraph)
  n = nv(graph)
  solver.graph = graph
  solver.weights_matrix = zeros(n, n)
end

function solve_linear(solver::SpanningTreeLightGraphsPrimSolver, reward::Dict{Tuple{Int, Int}, Float64})
  # Make up a rewards matrix by copying the input dictionary into the right data structure for the shortest path computations.
  fill!(solver.weights_matrix, 0.0)
  for (k, v) in reward
    solver.weights_matrix[k[1], k[2]] = v
  end
  solver.weights_matrix = maximum(solver.weights_matrix) .- solver.weights_matrix

  # Compute the maximum spanning tree.
  mst = prim_mst(solver.graph, solver.weights_matrix)

  # Transform edges into tuples.
  return _mst_solution_normalise(reward, [(src(e), dst(e)) for e in mst])
end


text:
"""
During optimization, a `SolverLog` is used to store convergence information.

A `SolverLog` can log the following:
- `x::Array{Float64,2}`: the history of the decision variables for every iteration.
- `abs_feas_error::Array{Float64,1}`: the absolute feasibility error evolution.
- `abs_opt_error::Array{Float64,1}`: the absolute optimality error evolution.
- `obj_value::Array{Float64,1}`: the evolution of the objective function value.
- `fc_evals::Array{Float64,1}`: the evolution of the number of function evaluations.
- `ga_evals::Array{Float64,1}`: the evolution of the number of gradient evaluations.
- `nStatus::Int`: the solution status return code.

The actual values stored in a `SolverLog` depend on the implementation of the callback used by each solver for every new point.
Using Ipopt, the callback is set with [`Ipopt.setIntermediateCallback`](https://github.com/jump-dev/Ipopt.jl/blob/d9e9176620a9b527a08991a3d41062fa948867f7/test/C_wrapper.jl#L123-L131).
Using Knitro, the callback is set with [`KNITRO.KN_set_newpt_callback`](https://github.com/jump-dev/KNITRO.jl/blob/43de7d38f01d7aa7c418bb24d0c1ab89b61e527a/examples/nlp2.jl#L201-L204).
"""
mutable struct SolverLog
    x::Array{Float64,2}
    abs_feas_error::Array{Float64,1}
    abs_opt_error::Array{Float64,1}
    obj_value::Array{Float64,1}
    fc_evals::Array{Float64,1}
    ga_evals::Array{Float64,1}
    nStatus::Int  # solution status return code

    @doc """
        SolverLog(N)

    Create an empty NLP solver log.

    # Arguments
    - `N::Int`: the total number of decision variables of the NLP.

    See also: [`update!`](@ref)
    """
    function SolverLog(N::Int)
        new(Array{Float64,2}(undef, N, 0),
            Array{Float64,1}(undef, 0),
            Array{Float64,1}(undef, 0),
            Array{Float64,1}(undef, 0),
            Array{Float64,1}(undef, 0),
            Array{Float64,1}(undef, 0),
            -999)
    end
end

"""
    update!(log; <keyword arguments>)

Update a [`SolverLog`](@ref) object.

# Keyword arguments
- `x=zeros(size(log.x, 1))`: an intermediate iterate point.
- `abs_feas_error=0`: the absolute feasibility error at the current point.
- `abs_opt_error=0`: the absolute optimality error at the current point.
- `obj_value=0`: the value of the objective function at the current point.
- `fc_evals=0`: the number of function evaluations requested by the NLP solver so far.
- `ga_evals=0`: the number of gradient evaluations requested by the NLP solver so far.
"""
function update!(log::SolverLog;
                 x=zeros(size(log.x, 1)),
                 abs_feas_error=0,
                 abs_opt_error=0,
                 obj_value=0,
                 fc_evals=0,
                 ga_evals=0)
    log.x = hcat(log.x, x)
    log.abs_feas_error = vcat(log.abs_feas_error, abs_feas_error)
    log.abs_opt_error = vcat(log.abs_opt_error, abs_opt_error)
    log.obj_value = vcat(log.obj_value, obj_value)
    log.fc_evals = vcat(log.fc_evals, fc_evals)
    log.ga_evals = vcat(log.ga_evals, ga_evals)
end

"""
    length(log) -> Integer

Return the length of a [`SolverLog`](@ref).

The length of a [`SolverLog`](@ref) is the same as the number of iterations taken to solve the problem.
"""
function length(log::SolverLog)
    num_iters = size(log.x, 2)
    @assert all(==(num_iters), length.([log.abs_feas_error, log.abs_opt_error, log.obj_value, log.fc_evals, log.ga_evals]))
    return num_iters
end

"""
    save(log, file)

Save a [`SolverLog`](@ref) to an `.npz` file.

See also: [`load!`](@ref)
"""
function save(log::SolverLog, file::String)
    npzwrite(file, Dict(
        "x" => log.x,
        "abs_feas_error" => log.abs_feas_error,
        "abs_opt_error" => log.abs_opt_error,
        "obj_value" => log.obj_value,
        "fc_evals" => log.fc_evals,
        "ga_evals" => log.ga_evals,
        "nStatus" => log.nStatus,
    ))
end

"""
    load!(log, file)

Load an `.npz` file and store it in a [`SolverLog`](@ref) object.

See also: [`save`](@ref)
"""
function load!(log::SolverLog, file::String)
    data = npzread(file)

    log.x = data["x"]
    log.abs_feas_error = data["abs_feas_error"]
    log.abs_opt_error = data["abs_opt_error"]
    log.obj_value = data["obj_value"]
    log.fc_evals = data["fc_evals"]
    log.ga_evals = data["ga_evals"]
    log.nStatus = data["nStatus"]

    return log
end

export
    SolverLog,
    update!,
    length,
    save,
    load!