text:
import JLD

function write_cache(doc::WeaveDoc, cache_path)
    cache_dir = "$(doc.cwd)/$cache_path"
    isdir(cache_dir) || mkpath(cache_dir)
    JLD.save("$cache_dir/$(doc.basename).jld", Dict("doc" => doc))
    return nothing
end

function read_cache(doc::WeaveDoc, cache_path)
    name = "$(doc.cwd)/$cache_path/$(doc.basename).jld"
    isfile(name) || return nothing
    return JLD.load(name)["doc"]
end

function restore_chunk(chunk::CodeChunk, cached)
    chunks = filter(x -> x.number == chunk.number &&
            string(typeof(x)) == "Weave.CodeChunk", cached.chunks)

    #Chunk types, don't match after loading. Fix by constructing chunks
    #from loaded content
    new_chunks = Any[]
    for c in chunks
      newc = CodeChunk(c.content, c.number, c.start_line, c.optionstring, c.options)
      newc.result_no = c.result_no
      newc.figures = c.figures
      newc.result = c.result
      newc.output = c.output
      newc.rich_output = c.rich_output
      push!(new_chunks, newc)
    end
    return new_chunks
end

#Could be used to restore inline code in future
function restore_chunk(chunk::DocChunk, cached)
      return chunk
end


text:
using LinearMaps
using Magnus
using PyPlot

include("plot_convergence.jl")

b = -1.
ω = 2π
f = t -> sin(ω*t)
c = 10.
y0 = [1.0]

tmax = π
N = 200

t = linspace(0,tmax,N+1)
y = exp.(-(b*t - (cos.(ω*t)-1)/ω * c))*y0[1]

function plot_approx(t_approx,y_approx,subs,label)
    st = length(t_approx) > 100 ? "-" : ".-"
    subplot(subs[1])
    plot(t_approx, y_approx, st, label=label)
    subplot(subs[2])
    semilogy(t_approx, abs(y_approx-y), st, label=label)
end

Exp = ScalarExponentiator(y0)
midpoint = MidpointPropagator(t -> LinearMap(b + f(t)*c), -1, Exp)
cfet = CFET4BfCPropagator(LinearMap(b), f, LinearMap(c), -1, Exp)

figure("scalar",figsize=(10,14))
clf()
subplot(411)
for (prop,name) in [(midpoint,"Midpoint"),
                    (cfet, "CFET (25)")]
    result = integrate(y0, tmax, N, prop,
                       save_intermediate = true,
                       verbose = true)
    plot_approx(t, vec(result[:V]),
                [411,412], name)
end
legend(framealpha=0.75)
margins(0,1)
xlabel(L"t")
ylabel("Error")
subplot(411)
plot(t, y, "k--", label="Exact")
ylabel(L"y(t)")
margins(0,0.1)
gca()[:set_xticklabels]([])

Nl = 0,5
plot_convergence(y[end], y0, tmax, Nl..., midpoint, [413,414], "Midpoint")
plot_convergence(y[end], y0, tmax, Nl..., cfet, [413,414], "CFET (25)")
subplot(413)
Nl = logspace(Nl...,10)
plot(Nl,1./Nl, ":", label="Linear")
plot(Nl,1./Nl.^2, ":", label="Quadratic")
plot(Nl,1./Nl.^3, ":", label="Cubic")
plot(Nl,1./Nl.^4, ":", label="Quartic")
plot(Nl,1./Nl.^5, ":", label="Quintic")
plot(Nl,1./Nl.^6, ":", label="Sextic")

xscale("log")
yscale("log")
xlabel("Number of time steps")
ylabel("Error")
legend(framealpha=0.75,loc=4,ncol=2)

subplot(414)
xscale("log")
yscale("log")
xlabel("Error")
ylabel("Execution time [s]")
legend()
tight_layout()

ext: jl


text:
mutable struct NewickXMLElement <: MyXMLElement
    el::XMLOrNothing
    newick::String
    fix_tree::Bool

    NewickXMLElement(newick::String) = new(nothing, newick, true)
end

function make_xml(nl::NewickXMLElement)

    el = new_element(bn.NEWICK)
    set_attribute(el, bn.ID, bn.DEFAULT_TREE_NAME)
    if nl.fix_tree
        set_attribute(el, bn.USING_HEIGHTS, bn.TRUE)
        set_attribute(el, bn.USING_DATES, bn.FALSE)
    end
    set_content(el, "$(nl.newick);")
    nl.el = el
    return el
end


text:
"""
    `generate_vtk_cuboid(filename::String, x1::Real, x2::Real, x3::Real;
                                    O::Array{T1, 1}=zeros(Float64, 3),
                                    Oaxis::Array{T1, 1}=Array{Float64}(I, 3, 3),
                                    optargs...)`

    Returns a cuboid of sides `x1`, `x2`, `x3` with origin `O` and orientation
`Oaxis`. Vertex [0,0,0] is placed at the origin. `Oaxis[i, :]` contains the
unit vector of the i-th axis.
"""
function generate_vtk_cuboid(filename::String, x1::Real, x2::Real, x3::Real;
                                O::Array{T1, 1}=zeros(Float64, 3),
                                Oaxis::Array{T2, 2}=Array{Float64}(I, 3, 3),
                                optargs...) where {T1<:Real, T2<:Real}

    check_coord_sys(Oaxis)
    invOaxis = collect(Oaxis)'

    # Defining points
    p0 = countertransform([0,0,0], invOaxis, O)
    p1 = countertransform([x1,0,0], invOaxis, O)
    p2 = countertransform([x1,x2,0], invOaxis, O)
    p3 = countertransform([0,x2,0], invOaxis, O)
    p4 = countertransform([0,0,x3], invOaxis, O)
    p5 = countertransform([x1,0,x3], invOaxis, O)
    p6 = countertransform([x1,x2,x3], invOaxis, O)
    p7 = countertransform([0,x2,x3], invOaxis, O)
    points = [p0, p1, p2, p3, p4, p5, p6 ,p7]

    # Defining cells
    c0 = [0, 3, 2, 1]
    c1 = [4, 5, 6, 7]
    c2 = [0, 4, 7, 3]
    c3 = [1, 2, 6, 5]
    c4 = [0, 1, 5, 4]
    c5 = [2, 3, 7, 6]

    cells = [c0, c1, c2, c3, c4, c5]

    # Generate vtk file
    generateVTK(filename, points; cells=cells, optargs...)
end

"""
    Returns a cylinder of radius `r` and height `h` with origin `O` at the
center of the lower face, and orientation `Oaxis` with the centerline along
the third axis. `Oaxis[i, :]` contains the unit vector of the i-th axis.
"""
function generate_vtk_cyl(filename::String, r::Real, h::Real;
                                np::Int64=36,
                                O::Array{T1, 1}=zeros(Float64, 3),
                                Oaxis::Array{T2, 2}=Array{Float64}(I, 3, 3),
                                optargs...) where {T1<:Real, T2<:Real}

    check_coord_sys(Oaxis)
    invOaxis = collect(Oaxis')
    cen = countertransform([0, 0, h], invOaxis, zeros(3)) # Centerline to upper face

    # Generate lower face
    points = [countertransform([r*cos(a),r*sin(a),0], invOaxis, O)
                                    for a in range(0, 2*pi, length=np+1)[1:end-1]]
    # Generate upper face and add centerline points
    points = vcat([O], points, [O+cen], [p+cen for p in points])

    # Defining lower face cells
    cells = [[0, i, i%np + 1] for i in 1:np]

    # Defining side face cells
    cells = vcat(cells, [[i, np+1+i, np+1+(i%np+1), i%np+1] for i in 1:np])

    # Defining uper face cells
    cells = vcat(cells, [[np+1, np+1 + i%np + 1, np+1 + i] for i in 1:np])

    # Generate vtk file
    generateVTK(filename, points; cells=cells, optargs...)
end


text:
# an example on mixed integer programming
#
#   maximize x + 2 y + 5 z
#
#   s.t.  x + y + z <= 10
#         x + 2 y + z <= 15
#
#         x is continuous: 0 <= x <= 5
#         y is integer: 0 <= y <= 10
#         z is binary
#
# z = 1, y= 7, x = 0
using Xpress, Compat.Test
using Compat.SparseArrays

@testset "MIP 4" begin

    model = Xpress.Model( "mip_01", :maximize)

    add_cvar!(model, 1., 0., 5.)  # x
    add_ivar!(model, 2., 0, 10)   # y
    add_bvar!(model, 5.)          # z

    add_constr!(model, ones(3), '<', 10.)
    add_constr!(model, [1., 2., 1.], '<', 15.)

    @test Xpress.get_obj(model) == [1, 2, 5]
    @test Xpress.get_constrmatrix(model) == sparse([1. 1. 1.; 1. 2. 1.])
    @test Xpress.get_rhs(model) == [10, 15]
    @test Xpress.get_lb(model) == [0, 0, 0]
    @test Xpress.get_ub(model) == [5, 10, 1]
    @test Xpress.get_sense(model) == [:<=,:<=]

    optimize(model)

    ans = get_optiminfo(model)
    @test ans.status_mip == :mip_optimal

    @test isapprox(get_solution(model), [0, 7, 1]; atol = 1e-3)

    @test isapprox(get_objval(model), 19; atol = 1e-3)
end


text:
using SPRING
using Test


text:
#
# XPA.jl --
#
# Implement XPA communication via the dynamic library.
#
#------------------------------------------------------------------------------
#
# This file is part of XPA.jl released under the MIT "expat" license.
# Copyright (C) 2016-2020, Éric Thiébaut (https://github.com/JuliaAstro/XPA.jl).
#

module XPA

using XPA_jll

using Base: ENV, @propagate_inbounds

include("cdefs.jl")
include("types.jl")
include("misc.jl")
include("client.jl")
include("server.jl")

end # module


text:
## Promotions and conversions

# Promote to Basket
function Base.promote_rule(
        ::Type{Basket},
        ::Type{B}) where B<:AbstractMonetary
    Basket
end

# Convert to basket
Base.convert(::Type{Basket}, m::AbstractMonetary) = Basket((m,))

# Convert to Monetary
function Base.convert(::Type{T}, b::Basket) where T<:Monetary
    len = length(b)
    if len == 0
        zero(T)
    elseif len == 1
        convert(T, first(b))
    else
        throw(InexactError(:convert, T, b))
    end
end

## Functions that promote
b::AbstractMonetary == c::AbstractMonetary = iszero(-(promote(b, c)...))

b::Basket            + c::AbstractMonetary = +(promote(b, c)...)
b::AbstractMonetary  + c::Basket           = c + b
b::AbstractMonetary  - c::AbstractMonetary = b + (-c)


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(MachO.r_extern, (MachO.relocation_info,))
    precompile(MachO.getindex, (MachO.Sections, Int64,))
    precompile(MachO.unsafe_string, (Base.AbstractIOBuffer{Array{UInt8, 1}}, UInt32,))
    precompile(MachO.Type, (Type{MachO.LoadCmds}, MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, Void, Int64,))
    precompile(MachO.getindex, (MachO.Relocations, Int64,))
    precompile(MachO.unpack, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, Type{MachO.load_command},))
    precompile(MachO.r_type, (MachO.relocation_info,))
    precompile(MachO.getindex, (MachO.Sections, UInt32,))
    precompile(MachO.getindex, (MachO.FatMachOHandle, Int64,))
    precompile(MachO.getindex, (MachO.Symbols, UInt32,))
    precompile(MachO.next, (MachO.LoadCmds, Tuple{UInt64, Int64, Int64},))
    precompile(MachO.unpack, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, Type{UInt32},))
    precompile(MachO.readloadcmd, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}},))
    precompile(MachO.unpack, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, Type{MachO.dyld_info_command},))
    precompile(MachO.debugsections, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}},))
    precompile(MachO.Type, (Type{MachO.FatMachOHandle}, Base.AbstractIOBuffer{Array{UInt8, 1}}, Int64,))
    precompile(MachO.unsafe_string, (MachO.small_fixed_string,))
    precompile(MachO.debugsections, (MachO.LoadCmd{MachO.segment_command_64},))
    precompile(MachO.readheader, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}},))
    precompile(MachO.Type, (Type{MachO.LoadCmd}, MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, UInt64, UInt32, MachO.source_version_command,))
    precompile(MachO.unpack, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, Type{MachO.routines_command_64},))
    precompile(MachO.Type, (Type{MachO.SectionRef{MachO.section_64}}, MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, MachO.section_64,))
    precompile(MachO.Type, (Type{MachO.LoadCmd}, MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, UInt64, UInt32, MachO.version_min_macosx_command,))
    precompile(MachO.unpack, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, Type{MachO.section_64},))
    precompile(MachO.unpack, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, Type{MachO.uuid_command},))
    precompile(MachO.Type, (Type{MachO.LoadCmd}, MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, UInt64, UInt32, MachO.uuid_command,))
    precompile(MachO.unpack, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, Type{MachO.linkedit_data_command},))
    precompile(MachO.Type, (Type{MachO.LoadCmd}, MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, UInt64, UInt32, MachO.linkedit_data_command,))
    precompile(MachO.Type, (Type{MachO.LoadCmd}, MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, UInt64, UInt32, MachO.thread_command,))
    precompile(MachO.unpack, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, Type{MachO.relocation_info},))
    precompile(MachO.Type, (Type{MachO.LoadCmd}, MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, UInt64, UInt32, MachO.segment_command_64,))
    precompile(MachO.readuntil, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, Char,))
    precompile(MachO.__init__, ())
    precompile(MachO.Type, (Type{MachO.sub_framework_command}, UInt32, String,))
    precompile(MachO.unpack, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, Type{MachO.version_min_macosx_command},))
    precompile(MachO.Type, (Type{MachO.LoadCmd}, MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, UInt64, UInt32, MachO.sub_framework_command,))
    precompile(MachO.unpack, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, Type{MachO.segment_command},))
    precompile(MachO.unpack, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, Type{MachO.source_version_command},))
    precompile(MachO.Type, (Type{MachO.LoadCmd}, MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, UInt64, UInt32, MachO.dyld_info_command,))
    precompile(MachO.Type, (Type{MachO.LoadCmd}, MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, UInt64, UInt32, MachO.dylinker_command,))
    precompile(MachO.unpack, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, Type{MachO.dysymtab_command},))
    precompile(MachO.Type, (Type{MachO.LoadCmd}, MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, UInt64, UInt32, MachO.symtab_command,))
    precompile(MachO.Type, (Type{MachO.LoadCmd}, MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, UInt64, UInt32, MachO.sub_client_command,))
    precompile(MachO.r_length, (MachO.relocation_info,))
    precompile(MachO.seek, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, UInt64,))
    precompile(MachO.readmeta, (Base.AbstractIOBuffer{Array{UInt8, 1}}, Type{MachO.MachOHandle},))
    precompile(MachO.Type, (Type{MachO.LoadCmd}, MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, UInt64, UInt32, MachO.dysymtab_command,))
    precompile(MachO.unpack, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, Type{MachO.symtab_command},))
    precompile(MachO.Type, (Type{MachO.LoadCmd}, MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, UInt64, UInt32, MachO.routines_command_64,))
    precompile(MachO.unpack, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, Type{MachO.segment_command_64},))
    precompile(MachO.position, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}},))
    precompile(MachO.Type, (Type{MachO.LoadCmd}, MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, UInt64, UInt32, MachO.dylib_command,))
    precompile(MachO.compute_X86_64, (MachO.RelocationRef, UInt64,))
    precompile(MachO.Type, (Type{MachO.LoadCmds}, MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}},))
    precompile(MachO.unpack, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, Type{MachO.mach_header_64},))
    precompile(MachO.Type, (Type{MachO.LoadCmd}, MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, UInt64, UInt32, MachO.rpath_command,))
    precompile(MachO.sectionsize, (MachO.section_64,))
    precompile(MachO.sizeof, (Type{MachO.section_64},))
    precompile(MachO.Type, (Type{MachO.Sections}, MachO.LoadCmd{MachO.segment_command_64},))
    precompile(MachO.Type, (Type{MachO.dylib_command}, UInt32, UInt32, UInt32, UInt32, String,))
    precompile(MachO.print, (Base.AbstractIOBuffer{Array{UInt8, 1}}, MachO.small_fixed_string,))
    precompile(MachO.seek, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, UInt32,))
    precompile(MachO.sectionname, (MachO.SectionRef{MachO.section_64},))
    precompile(MachO.Type, (Type{MachO.sub_client_command}, String,))
    precompile(MachO.Type, (Type{MachO.Symbols}, MachO.LoadCmd{MachO.symtab_command},))
    precompile(MachO.load_strtab, (MachO.SectionRef{MachO.section_64},))
    precompile(MachO.sectionoffset, (MachO.section_64,))
    precompile(MachO.Type, (Type{MachO.dylinker_command}, String,))
    precompile(MachO.seek, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, Int64,))
    precompile(MachO.sectionname, (MachO.section_64,))
    precompile(MachO.Type, (Type{MachO.Sections}, MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}},))
    precompile(MachO.deref, (MachO.LoadCmd{MachO.uuid_command},))
    precompile(MachO.Type, (Type{ObjFileBase.Sections}, MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}},))
    precompile(MachO.deref, (MachO.LoadCmd{MachO.thread_command},))
    precompile(MachO.Type, (Type{MachO.rpath_command}, UInt32, String,))
    precompile(MachO.deref, (MachO.LoadCmd{MachO.segment_command_64},))
    precompile(MachO.deref, (MachO.SectionRef{MachO.section_64},))
    precompile(MachO.handle, (MachO.SectionRef{MachO.section_64},))
    precompile(MachO.Type, (Type{ObjFileBase.StrTab}, MachO.SectionRef{MachO.section_64},))
    precompile(MachO.unpack, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, Type{MachO.dylib_command}, UInt32,))
    precompile(MachO.compute_value, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, MachO.Symbols, MachO.Sections, Base.Dict{Symbol, Any}, MachO.RelocationRef,))
    precompile(MachO.unpack_lcstr, (MachO.MachOHandle{Base.AbstractIOBuffer{Array{UInt8, 1}}}, UInt32, Int64, UInt32,))
    precompile(MachO.compute_relocation, (MachO.mach_header_64, MachO.RelocationRef, UInt64,))
    precompile(MachO.Type, (Type{ObjFileBase.Sections}, MachO.LoadCmd{MachO.segment_command_64},))
end


text:
using Test
using NQCDynamics
using FiniteDiff
using Random
using OrdinaryDiffEq
using DiffEqDevTools
using NQCDynamics: DynamicsMethods, DynamicsUtils
using NQCDynamics.DynamicsMethods.MappingVariableMethods
Random.seed!(1)

@test MappingVariableMethods.NRPMD{Float64}(10, 0.0) isa MappingVariableMethods.NRPMD
atoms = Atoms(1)
sim = RingPolymerSimulation{NRPMD}(atoms, DoubleWell(), 10; temperature=1e-3)

v = randn(size(sim))
r = randn(size(sim))
u = DynamicsVariables(sim, v, r, PureState(2))

@testset "Population correlation" begin
    K = NQCModels.nstates(sim)
    out = zeros(K, K)
    n = 1e4
    correlation = TimeCorrelationFunctions.PopulationCorrelationFunction(sim, Diabatic())
    normalisation = TimeCorrelationFunctions.evaluate_normalisation(sim, correlation)

    for i=1:n
        u = DynamicsVariables(sim, v, r, PureState(1))
        K = Estimators.initial_diabatic_population(sim, u)
        Kinv = Estimators.diabatic_population(sim, u)
        out += normalisation * K * Kinv'
    end
   @test out ./ n ≈ [1 0; 0 0] atol=0.1
end

@testset "motion! obeys Hamilton's equations" begin
    function test_motion!(sim, u)
        f(x) = DynamicsUtils.classical_hamiltonian(sim, x)

        grad = FiniteDiff.finite_difference_gradient(f, u)

        du = zero(u)
        DynamicsMethods.motion!(du, u, sim, 0.0)

        @test DynamicsUtils.get_positions(du) ≈ DynamicsUtils.get_velocities(grad) ./ sim.atoms.masses' rtol=1e-3
        @test DynamicsUtils.get_velocities(du) ≈ -DynamicsUtils.get_positions(grad) ./ sim.atoms.masses' rtol=1e-3
        @test MappingVariableMethods.get_mapping_positions(du) ≈ MappingVariableMethods.get_mapping_momenta(grad) rtol=1e-3
        @test MappingVariableMethods.get_mapping_momenta(du) ≈ -MappingVariableMethods.get_mapping_positions(grad) rtol=1e-3
    end

    test_motion!(sim, u)
end

algs = (DynamicsMethods.IntegrationAlgorithms.RingPolymerMInt(), Tsit5())
@testset "Energy conservation $alg" for alg in algs
    sol = run_trajectory(u, (0, 10.0), sim; output=(:hamiltonian, :population), dt=1e-2, algorithm=alg, abstol=1e-8, reltol=1e-8)
    @test sol.hamiltonian[1] ≈ sol.hamiltonian[end] rtol=1e-2
end

@testset "Algorithm comparison" begin
    sol = run_trajectory(u, (0, 10.0), sim; dt=1e-2, algorithm=DynamicsMethods.IntegrationAlgorithms.RingPolymerMInt())
    sol1 = run_trajectory(u, (0, 10.0), sim; algorithm=Tsit5(), reltol=1e-10, abstol=1e-10, saveat=sol.t)
    @test sol.u ≈ sol1.u rtol=1e-2
end

@testset "MInt algorithm convergence" begin
    tspan=(0, 10.0)
    prob = DynamicsMethods.create_problem(u, tspan, sim)
    dts = 1 .// 2 .^(8:-1:4)

    alg = DynamicsMethods.IntegrationAlgorithms.RingPolymerMInt()
    test_alg = Vern9()
    setup = Dict(:alg => test_alg, :adaptive=>true, :abstol=>1e-8, :reltol=>1e-8)
    res = analyticless_test_convergence(dts, prob, alg, setup)
    @test res.𝒪est[:final] ≈ 2 atol=0.1
end


text:
using Bridge, StaticArrays, Distributions, PyPlot
using Test
import Base.Math.gamma # just to use the name
#import Bridge: b, σ, a, transitionprob
using Bridge: runmean, sqmahal, Gaussian

const percentile = 3.0
const SV = SVector{2,Float64}
const SM = SMatrix{2,2,Float64,4}

kernel(x, a=0.0025) = exp(Bridge.logpdfnormal(x, a*I))

refine(tt, n) =  first(tt):(Base.step(tt)/n):last(tt)

TEST = false
CLASSIC = false

using Bridge.outer

# Define a diffusion process
if ! @isdefined(Target)
struct Target  <: ContinuousTimeProcess{SV}
    c::Float64
    κ::Float64
end
end

if ! @isdefined(Linear)
struct Linear  <: ContinuousTimeProcess{SV}
    T::Float64
    v::SV
    b11::Float64
    b21::Float64
    b12::Float64
    b22::Float64
end
end

g(t, x) = sin(x)
gamma(t, x) = 1.2 - sech(x)/2

# define drift and sigma of Target

Bridge.b(t, x, P::Target) = SV(P.κ*x[2] - P.c*x[1],  -P.c*x[2] + g(t, x[2]))::SV
Bridge.σ(t, x, P::Target) = SM(0.5, 0.0, 0.0, gamma(t, x[2]))
Bridge.a(t, x, P::Target) = SM(0.25, 0, 0, outer(gamma(t, x[2])))
Bridge.constdiff(::Target) = false

# define drift and sigma of Linear approximation

Bridge.b(t, x, P::Linear) = SV(P.b11*x[1] + P.b12*x[2], P.b21*x[1] + P.b22*x[2] + g(P.T, P.v[2]))
Bridge.B(t, P::Linear) = SM(P.b11, P.b21, P.b12, P.b22)

Bridge.β(t, P::Linear) = SV(0, g(P.T, P.v[2]))

Bridge.σ(t, x, P::Linear) = SM(0.5, 0, 0, gamma(P.T, P.v[2]))
Bridge.a(t, x, P::Linear) = SM(0.25, 0, 0, outer(gamma(P.T, P.v[2])))
Bridge.a(t, P::Linear) = SM(0.25, 0, 0, outer(gamma(P.T, P.v[2])))
Bridge.constdiff(::Linear) = false

c = 0.0
κ = 3.0

# times

t = 0.7
T = 1.5
S = (t + T)/2

# grid

n = 401
dt = (T-t)/(n-1)
tt = t:dt:T
S = tt[n÷2 + 1]
tt1 = t:dt:S
tt2 = S:dt:T

m = 1_800_000

Ti = n
Si = n÷2

# observations

Σ = 0.8 # observation noise
Q = Normal(0, sqrt(Σ)) # noise distribution
L = @SMatrix [1.0 0.0]

xt = @SVector [0.1, 0.0]
vS = @SVector [0.1]
xS = @SVector [0.1, -0.2]
xT = @SVector [0.3, -0.4]

# processes

P = Target(c, κ)
P̃ = Pt = Linear(T, xT, -c-0.1, -0.1, κ-0.1, -c/2)

# parameters

B = Bridge.B(0, Pt)
β = Bridge.β(0, Pt)
a = Bridge.a(0, Pt)
σ = sqrt(Bridge.a(0, Pt))

W = sample(tt, Wiener{SV}())
W1 = sample(tt1, Wiener{SV}())
W2 = sample(tt2, Wiener{SV}())

# Target and log probability density (forward simulation)

YT = SV[]
YS = Float64[]
p = Float64[]
X = SamplePath(tt, zeros(SV, length(tt)))
Xs = SamplePath(tt, zeros(SV, length(tt)))
best = Inf
for i in 1:m
    W = sample!(W, Wiener{SV}())
    Bridge.solve!(Euler(), X, xt, W, P)
    push!(YT, X.yy[end])
    push!(YS, X.yy[end÷2][2]) # depends on L

    eta = rand(Q)
    nrm = norm(xT - X.yy[Ti]) + norm(vS - eta - L*X.yy[Si])
    l = kernel(xT - X.yy[Ti]) * kernel(vS - L*X.yy[Si] - eta)
    push!(p, l)
    if nrm < best
        best = nrm
        Xs.yy .= X.yy
    end
end

lphat = log(mean(p))

# Proposal log probability density (forward simulation)

YtT = SV[]
p̃ = pt = Float64[]
Xt = SamplePath(tt, zeros(SV, length(tt)))
l = 0.0
for i in 1:m
    W = sample!(W, Wiener{SV}())
    Bridge.solve!(Euler(), Xt, xt, W, Pt)
    push!(YtT, Xt.yy[end])
    eta = rand(Q)
    l =  kernel(xT - Xt.yy[Ti]) * kernel(vS - L*Xt.yy[Si] - eta) # likelihood
    push!(pt, l)
end
lpthat = log(mean(pt))

@show lpthat

# Proposal

Z = Float64[]
Xo1 = SamplePath(tt1, zeros(SV, length(tt)))
Xo2 = SamplePath(tt2, zeros(SV, length(tt)))

GP2 = GuidedBridge(tt2, P, Pt, xT)
H♢, V = Bridge.gpupdate(GP2, L, Σ, vS)
GP1 = GuidedBridge(tt1, P, Pt, V, H♢)

@time for i in 1:m

    sample!(W1, Wiener{SV}())
    sample!(W2, Wiener{SV}())
    y = Bridge.bridge!(Xo1, xt, W1, GP1)
    Bridge.bridge!(Xo2, y, W2, GP2)

    ll = llikelihood(LeftRule(), Xo1, GP1) + llikelihood(LeftRule(), Xo2, GP2)
    push!(Z, exp(ll))
end

@show log(mean(exp.(lpthat))), lphat

@show lpthat
#@show lpt
@show lptilde(GP1, xt) - Bridge.traceB(tt2, Pt)

PhiS = Bridge.fundamental_matrix(tt1, Pt)
PhiT = Bridge.fundamental_matrix(tt, Pt)
PhiTS = Bridge.fundamental_matrix(tt2, Pt)

KS = Bridge.gpK(tt1, zero(SM), Pt)
KT = Bridge.gpK(tt, zero(SM), Pt)
KTS = Bridge.gpK(tt2, zero(SM), Pt)

muS = Bridge.gpmu(tt1, xt, Pt)
muT = Bridge.gpmu(tt, xt, Pt)
muTS = Bridge.gpmu(tt2, xS, Pt)

hS = Bridge.gpmu(tt1, 0*xt, Pt)
hT = Bridge.gpmu(tt, 0*xt, Pt)
hTS = Bridge.gpmu(tt2, 0*xS, Pt)

@test norm(PhiTS\(xT - muTS) - (GP2.V[1] - xS)) < 1e-9

#@test norm(muT - mean(YtT)) < 1/sqrt(m)

H♢ = GP1.H♢[1]
H♢S_ = GP1.H♢[end]
H♢S = GP2.H♢[1]
V = GP1.V[1]
VS = GP2.V[1]
VS_ = GP1.V[end]

u = [vS - L*muS; xT - muT]

Upsilon = [
L*KS*L' + Σ     L*KS*PhiTS'
PhiTS*KS*L'     KT
]

Upsilon2 = [
    L*(PhiS*H♢*PhiS' - GP1.H♢[end])*L' + Σ     L*(PhiS*H♢*PhiT' - GP1.H♢[end]*PhiTS')
    (PhiT*H♢*PhiS' - PhiTS*GP1.H♢[end])*L'     KT
    ]

Upsilon2 = [
    L*KS*L' + Σ     L*KS
    KS*L'     KS + GP2.H♢[1]
    ]
@test norm(cat((1,2),eye(1),PhiTS)*Upsilon2*cat((1,2),eye(1),PhiTS)' - Upsilon)<1e-8
U = inv(Upsilon)
#0 = logdet(Upsilon2) + 2*logdet(PhiTS) + logdet(U)

A = [L*PhiS; PhiT ]

z = u + A*xt

lpt = Bridge.logpdfnormal(u, Upsilon)

@test norm(logdet(U) + Bridge.logdetU(GP1, GP2, L, Σ)) < 1e-8

@test norm(PhiTS*xS + (hTS - muTS)) < 1e-7

@test norm(PhiTS*(VS-xS) - ( xT-PhiTS*xS - hTS ) ) < 1e-7

@test_broken norm(sqmahal(Gaussian( A*xt, Upsilon), z) - sqmahal(Gaussian(GP1.V[1], GP1.H♢[1]) , xt)) < 1e-8
@show lpt
@show lpthat

@show logdet(Upsilon)
@show logdet(KT) + logdet(Σ) + 2*logdet(PhiTS)
@test -logdet(U) ≈ logdet( L*KS*L' + Σ - L*KS*PhiTS'*inv(KT)*PhiTS*KS*L' ) +
    logdet(KT)
@test -logdet(U) ≈ logdet( KT - PhiTS*KS*L'*inv(L*KS*L' + Σ)*L*KS*PhiTS' ) +
    logdet(L*KS*L' + Σ )

@test 1/det(U) ≈  det(L*KS*L' + Σ - 1)*det(KT) +  det(KT -  PhiTS*KS*L'*L*KS*PhiTS')

@test norm(logdet(KTS) - (logdet(GP2.H♢[1]) + 2logdet(PhiTS))) < 1e-8

@test norm(A'*U*A - inv(GP1.H♢[1])) < 1e-7

@test norm(PhiTS\(xT - muTS) - (GP2.V[1] - xS)) < 1e-9
@test norm(inv(GP1.H♢[end]) - inv(GP2.H♢[1]) -  L'*inv(Σ)*L) < 1e-10
@test norm(KS - PhiS*H♢*PhiS' + GP1.H♢[end]) < 1e-8
@test norm(KT - PhiT*H♢*PhiT' - PhiTS*(GP2.H♢[1]- GP1.H♢[end])*PhiTS' ) < 1e-8

@test norm(xt'*inv(H♢)*xt - xt'*A'*U*A*xt) < 1e-7
@test norm(xt'*inv(H♢)*V - xt'*A'*U*z) < 1e-7

@test norm(PhiT*(V-xt) + PhiTS*(VS-VS_)  - ( xT-PhiT*xt - hT ) ) < 1e-7

mh1 = [vS - L*muS;xT - muT]'*U*[vS - L*muS;xT - muT]
mh2 = (V - xt)'*inv(H♢)*(V - xt)

r1, r2 = H♢\(V - xt), A'*U*[vS - L*muS;xT - muT]
@test norm(r1 - r2)  < 1e-7

# Plot best "bridge"

clf()
subplot(121)
plot(Xs.tt, first.(Xs.yy), label="X₁˟")
plot(Xs.tt, last.(Xs.yy), label="X₂˟")
plot.(t, xt, "ro")
plot(S, vS, "ro")
plot.(T, xT, "ro")
legend()

subplot(122)
step = 10
plot(mean(pt)*runmean(Z)[1:step:end], label="Ψ p̃")
plot(runmean(p)[1:step:end], label="p")
plot(runmean(pt)[1:step:end], label="p̃")
plot(fill(exp(lpt),length(1:step:m)), label="p̃ theor.")
legend()
axis([1, div(m,step), 0, 3*exp(lpthat)])

error("done")

figure()
subplot(411)
plot(Xs.tt, Xs.yy, label="X*")
legend()

subplot(413)
plot(Xo.tt, Xo.yy, label="Xo")
legend()

subplot(414)
step = 10
plot(runmean(exp.(Z))[1:step:end], label="Xo")
plot(runmean(kernel.(Yv))[1:step:end], label="X")
plot(runmean(kernel.(Ytv))[1:step:end], label="Xt")
legend()
axis([1, div(m,step), 0, 2*exp(lpt)])

ex = Dict(
"u" => u,
"v" => v,
"Xo" => Xo,
"Xs" => Xs,
"Xts" => Xts,
"Yt" => Yt,
"Y" => Y,
"Z" => Z,
"lpt" => lpt
)


text:
module URIParser

export URI
export defrag, userinfo, path_params, query_params, isvalid
export escape, escape_form, escape_with, unescape, unescape_form

using Unicode
using Base.Unicode: isascii

import Base: isequal, isvalid, show, print, (==)

include("parser.jl")
include("esc.jl")
include("utils.jl")
include("precompile.jl")
_precompile_()

end # module


text:
using Test,Revise
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\utils\\CUDAGpuUtils.jl")

includet("C:/GitHub/GitHub/NuclearMedEval/src/distanceMetrics/SimplerHousdorff.jl")
using ..SimplerHousdorff, ..CUDAGpuUtils
using CUDA

# @testset "getBlockNumb" begin
# arrGold = CUDA.zeros(1023,200,800);
#     @test  ..SimplerHousdorff.getBlockNumb(size(arrGold))==(800, 1023  )
#     CUDA.reclaim()# just to destroy from gpu our dummy data
#     arrGold = CUDA.zeros(1028,200,2000);
#     @test  ..SimplerHousdorff.getBlockNumb(size(arrGold))==(4000,1024)
#     CUDA.reclaim()# just to destroy from gpu our dummy data
# end #

@testset "allocateMomory" begin
    arrGold = CUDA.zeros(1024,200,800);
    @test  size(..SimplerHousdorff.allocateMomory(size(arrGold))) == (32, cld(200,32), cld(800,32),6)
    CUDA.reclaim()# just to destroy from gpu our dummy data

end #


text:
export PermCycs
export randpermcycs

immutable PermCycs{T<:Real} <: AbstractPerm{T}
    data::Array{Array{T,1},1}
end

## Construct PermCycs objects ##

# construct permutation from list of disjoint cycles
function permcycs(T::Type, cycs...)
    length(cycs) == 0 && return PermCycs()
    pc = PermCycs(PermPlain.canoncycles([collect(T,c) for c in cycs]))
    isperm(pc) || error("Trying to construct PermCycs from illegal or non-disjoint cycles.")
    return pc
end

PermCycs(a::Tuple) = PermCycs(PermPlain.tupcollect(a))
PermCycs{T}(::Type{T}, a::Tuple) = PermCycs(PermPlain.tupcollect(T,a))
PermCycs{T<:Real}(::Type{T}) = PermCycs(Array(Array{T,1},0))

function permcycs(cycs...)
    length(cycs) == 0 && return PermCycs()
    permcycs(typeof(cycs[1][1]),cycs...)
end

PermCycs() = PermCycs(Array(Array{Int,1},0))
#PermCycs{T}(::Type{T}) = PermCycs(Array(Array{T,1},0))
permcycs() = PermCycs()
idpermcyc() = PermCycs()
idpermcyc{T}(::Type{T}) = PermCycs(T)
one(::Type{PermCycs}) = PermCycs()
one{T}(::Type{PermCycs{T}}) = PermCycs(T)
idperm{T}(c::PermCycs{T}) = PermCycs(T)
plength(c::PermCycs) = greatestmoved(c)

## Copying, indexing, ... ##

copy(c::PermCycs) = PermCycs(copy(c.data))
length(c::PermCycs) = length(c.data)
getindex(c::PermCycs, k::Real) = c.data[k]
setindex!(c::PermCycs, ci::Real, k::Real) = c.data[k] = ci

## Compare, test, and/or return properties ##

#ispermcycs(c::PermCycs) = ispermcycs(c.data) # probably not useful/wise
isperm(c::PermCycs) = PermPlain.isperm(c.data)
isid(c::PermCycs) = length(c) == 0
# commute()
# distance()
==(c1::PermCycs, c2::PermCycs) = c1.data == c2.data # c1 and c2 must be in canonical order
sign(c::PermCycs) = PermPlain.permsgn(c.data)
order(c::PermCycs) = PermPlain.permorder(c.data)
cyclelengths(c::PermCycs) = PermPlain.cyclelengths(c.data)
cycletype(c::PermCycs) = PermPlain.cycletype(c.data)

## Apply permutation, and permutation operations ##
getindex{T}(v::Array{T,1},c::PermCycs{Bool}) = error("Silence compiler. You don't want Bool, anyway")
getindex(v::Array, c::PermCycs) = v[list(c).data]
getindex(v::String, c::PermCycs) = v[list(c).data] # How to define this for everything?

# Whether to return as-is or sort ?
# hmm, sort flag is overkill, since wrapping call in sort! is easy
function flatten(c::PermCycs; sort = false ) # a bit inefficient
    length(c) == 0 && return Array(eltype(c),0)
    outa = Array(eltype(c),0)
    for cyc in c.data
        append!(outa,cyc)
    end
    sort ? sort!(outa) : nothing
    return outa
end

# sort, to agree with PermList, and for fixed()
support(c::PermCycs) = flatten(c,sort=true)
supportsize(c::PermCycs) = sum([length(c1) for c1 in c.data])

# must be an idiomatic way to do this
# Not sure this is useful. fixed wrt what set ?
function fixed(c::PermCycs)
    outa = Array(eltype(c),0)
    mvd = support(c)
    ci = 1
    mi = mvd[ci]
    for i in 1:mvd[end]
        if i != mi
            push!(outa,i)
        elseif ci < length(mvd)
            ci += 1
            mi = mvd[ci]
        end
    end
    return outa
end

greatestmoved(cs::PermCycs) = (isempty(cs) ? zero(eltype(cs)) : maximum([ maximum(c) for c in cs.data ]))
leastmoved(cs::PermCycs) = minimum([ minimum(c) for c in cs.data ])

## Output ##

function _print(io::IO, c::PermCycs)
    print(io,"(")
    for cyc in c.data PermPlain.cycleprint(io, cyc) end
    print(io,")")
end

show(io::IO, c::PermCycs) = _print(io,c)
# This is needed to avoid trying to print PermList with showarray and failing in 1000 ways
writemime(io::IO, ::MIME"text/plain", p::PermCycs) = _print(io,p)


text:
using DataFrames, CategoricalArrays, GLM, Distributions, PrettyTables, QuasiGLM, Test

#------------- Run package tests --------------

@testset "QuasiGLM.jl" begin

    #-------------
    # Quasipoisson
    #-------------

    # Define some data

    dobson = DataFrame(Counts = [18,17,15,20,10,20,25,13,12], Outcome = categorical([1,2,3,1,2,3,1,2,3]), Treatment = categorical([1,1,1,2,2,2,3,3,3]))

    # Fit Poisson model

    gm = fit(GeneralizedLinearModel, @formula(Counts ~ Outcome + Treatment), dobson, Poisson())

    # Correct standard errors using quasi correction

    testOutputs = AdjustQuasiGLM(gm, dobson; level=0.95)
    @test testOutputs isa DataFrames.DataFrame

    #--------------
    # Quasibinomial
    #--------------

    # Set up data and divide percentage by 100 to get proportion

    blotchData = DataFrame(blotch = [0.05,0.00,1.25,2.50,5.50,1.00,5.00,5.00,17.50,0.00,0.05,1.25,0.50,1.00,5.00,0.10,10.00,25.00,0.00,0.05,2.50,0.01,6.00,5.00,5.00,5.00,42.50,0.10,0.30,16.60,3.00,1.10,5.00,5.00,5.00,50.00,0.25,0.75,2.50,2.50,2.50,5.00,50.00,25.00,37.50,0.05,0.30,2.50,0.01,8.00,5.00,10.00,75.00,95.00,0.50,3.00,0.00,25.00,16.50,10.00,50.00,50.00,62.50,1.30,7.50,20.00,55.00,29.50,5.00,25.00,75.00,95.00,1.50,1.00,37.50,5.00,20.00,50.00,50.00,75.00,95.00,1.50,12.70,26.25,40.00,43.50,75.00,75.00,75.00,95.00], variety = categorical(repeat([1,2,3,4,5,6,7,8,9], inner=1, outer=10)), site = categorical(repeat([1,2,3,4,5,6,7,8,9,10], inner=9, outer=1)))

    blotchData.blotch = blotchData.blotch ./ 100

    # Fit binomial model

    gm2 = fit(GeneralizedLinearModel, @formula(blotch ~ variety + site), blotchData, Binomial())

    # Correct standard errors using quasi correction

    testOutputs2 = AdjustQuasiGLM(gm2, blotchData; level=0.95)
    @test testOutputs2 isa DataFrames.DataFrame
end


text:
export quadratic_estimate

"""

    quadratic_estimate(ds::DataSet, which; wiener_filtered=true)
    quadratic_estimate((ds1::DataSet,ds2::DataSet), which; wiener_filtered=true)

Compute quadratic estimate of ϕ given data.

The `ds` or `(ds1,ds2)` tuple contain the DataSet object(s) which houses the
data and covariances used in the estimate. Note that only the Fourier-diagonal
approximations for the beam, mask, and noise,, i.e. `ds.B̂`, `ds.M̂`, and
`ds.Cn̂`, are accounted for. To account full operators (if they are not actually
Fourier-diagonal), you should compute the impact using Monte Carlo.

If a tuple is passed in, the result will come from correlating the data from
`ds1` with that from `ds2`, which can be useful for debugging / isolating
various noise terms.

An optional keyword argument `AL` can be passed in in case the QE normalization
was already computed, in which case it won't be recomputed during the
calculation.

Returns a NamedTuple `(ϕqe, AL, Nϕ)` where `ϕqe` is the (possibly Wiener
filtered, depending on `wiener_filtered` option) quadratic estimate, `AL` is the
normalization (which is already applied to ϕqe, it does not need to be applied
again), and `Nϕ` is the analytic N0 noise bias (Nϕ==AL if using unlensed
weights, currently only Nϕ==AL is always returned, no matter the weights)
"""
function quadratic_estimate((ds1,ds2)::Tuple{DataSet{F1},DataSet{F2}}, which; wiener_filtered=true, AL=nothing, weights=:unlensed) where {F1,F2}
    @assert weights in (:lensed, :unlensed) "weights should be :lensed or :unlensed"
    @assert (which in [:TT, :EE, :EB]) "which='$which' not implemented"
    @assert (ds1.Cf===ds2.Cf && ds1.Cf̃===ds2.Cf̃ && ds1.Cn̂===ds2.Cn̂ && ds1.Cϕ===ds2.Cϕ && ds1.B̂===ds2.B̂) "operators in `ds1` and `ds2` should be the same"
    @assert spin(F1)==spin(F2)
    check_hat_operators(ds1)
    @unpack Cf, Cf̃, Cn̂, Cϕ, B̂, M̂ = ds1()
    (quadratic_estimate_func, pol) = @match which begin
        :TT => (quadratic_estimate_TT, :I)
        :EE => (quadratic_estimate_EE, :P)
        :EB => (quadratic_estimate_EB, :P)
        _   => error("`which` argument to `quadratic_estimate` should be one of (:TT, :EE, :EB)")
    end
    quadratic_estimate_func((ds1.d[pol], ds2.d[pol]), Cf[pol], Cf̃[pol], Cn̂[pol], Cϕ, (M̂*B̂)[pol], wiener_filtered, weights, AL)
end
quadratic_estimate(ds::DataSet, which; kwargs...) = quadratic_estimate((ds,ds), which; kwargs...)
quadratic_estimate(ds::DataSet{<:FlatS0}; kwargs...) = quadratic_estimate(ds, :TT; kwargs...)
quadratic_estimate(ds::DataSet{<:FlatS2}; kwargs...) = quadratic_estimate(ds, :EB; kwargs...) # somewhat arbitraritly make default P estimate be EB

quadratic_estimate_TT(d::Field, args...) = quadratic_estimate_TT((d,d), args...)
quadratic_estimate_EB(d::Field, args...) = quadratic_estimate_EB((d,d), args...)
quadratic_estimate_EE(d::Field, args...) = quadratic_estimate_EE((d,d), args...)

@doc doc"""

All of the terms in the quadratic estimate and normalization expressions look like

    C * l[i] * l̂[j] * l̂[k] * ...

where C is some field or diagonal covariance. For example, there's a term in the EB
estimator that looks like:

    (CE * (CẼ+Cn) \ d[:E])) * l[i] * l̂[j] * l̂[k]

(where note that `l̂[j]` and `l̂[k]` are unit vectors, but `l[i]` is not).  The
function `get_term_memoizer` returns a function `term` which could be called in
the following way to compute this term:

    term((CE * (CẼ+Cn) \ d[:E])), [i], j, k)

(note that the fact that `l[i]` is not a unit vector is specified by putting the
`[i]` index in brackets).

Additionally, all of these terms are symmetric in their indices, i.e. in
`(i,j,k)` in this case. The `term` function is smart about this, and is memoized
so that each unique set of indices is only computed once. This leads to a pretty
drastic speedup for terms with many indices like those that arize in the EE and
EB normalizations, and lets us write code which is both clear and fast without
having to think too hard about these symmetries.

"""
function get_term_memoizer(f)
    Ðf = Ð(f)
    Lfactors() = 1
    Lfactors(inds...) = broadcast!(*, similar(Ðf), getproperty.(getindex.(Ref(∇),inds),:diag)...)
    term(C::Diagonal, inds...) = term(C.diag, inds...)
    term(C::Field, inds...) = term(count((x->x isa Int),inds)/2f0, C, first.(inds)...)
    @sym_memo term(n, C::Field, @sym(inds...)) = Map(nan2zero.(C .* Lfactors(inds...) ./ ∇².diag.^n))
    term
end
ϵ(x...) = levicivita([x...])
inds(n) = collect(product(repeated(1:2,n)...))[:]

function quadratic_estimate_TT((d1,d2)::NTuple{2,FlatS0}, Cf, Cf̃, Cn, Cϕ, TF, wiener_filtered, weights, AL=nothing)

    term = get_term_memoizer(d1)
    ΣTtot = TF^2 * Cf̃ + Cn
    CT = (weights==:unlensed) ? Cf : Cf̃
    @unpack Ωpix = fieldinfo(d1)

    # unnormalized estimate
    ϕqe_unnormalized = @subst -sum(∇[i] * Fourier(term($(ΣTtot\(TF*d1))) * term($(CT*(ΣTtot\(TF*d2))), [i])) for i=1:2)

    # normalization
    if AL == nothing
        AL = @subst begin
            A(i,j) = (
                term($(@. TF^2 * CT^2 / ΣTtot), [i], [j]) * term($(@. TF^2      / ΣTtot)     )
              + term($(@. TF^2 * CT   / ΣTtot), [i]     ) * term($(@. TF^2 * CT / ΣTtot), [j])
            )
            Ωpix * pinv(Diagonal(sum(∇[i].diag .* ∇[j].diag .* Fourier(A(i,j)) for (i,j) in inds(2))))
        end
    end
    Nϕ = AL # true only for unlensed weights

    ϕqe = (wiener_filtered ? (Cϕ*pinv(Cϕ+Nϕ)) : 1) * (AL*ϕqe_unnormalized)
    @namedtuple ϕqe AL Nϕ

end

function quadratic_estimate_EE((d1,d2)::NTuple{2,FlatS2}, Cf, Cf̃, Cn, Cϕ, TF, wiener_filtered, weights, AL=nothing)

    term = get_term_memoizer(d1[:E])
    TF² = TF[:E]^2
    ΣEtot = TF² * Cf̃[:E] + Cn[:E]
    CE = ((weights==:unlensed) ? Cf : Cf̃)[:E]
    @unpack Ωpix = fieldinfo(d1)

    # unnormalized estimate
    ϕqe_unnormalized = @subst begin
        I(i) = -(
            2sum(term($(CE * (ΣEtot \ (TF*d1)[:E])), [i], j, k) * term($((ΣEtot \ (TF*d2)[:E])), j, k) for (j,k) in inds(2))
               - term($(CE * (ΣEtot \ (TF*d1)[:E])), [i]      ) * term($((ΣEtot \ (TF*d2)[:E]))      )
        )
        sum(∇[i] * Fourier(I(i)) for i=1:2)
    end

    # normalization
    if AL == nothing
        AL = @subst begin
            A1(i,j) = -4 * sum( ϵ(m,p,3) * ϵ(n,q,3) * (
                  term($(@. TF² * CE^2 / ΣEtot), [i], [j], k, l, m, n) * term($(@. TF²      / ΣEtot),      k, l, p, q)
                + term($(@. TF² * CE   / ΣEtot), [i],      k, l, m, n) * term($(@. TF² * CE / ΣEtot), [j], k, l, p, q))
                for (k,l,m,n,p,q) in inds(6)
            )
            A2(i,j) = (
                  term($(@. TF² * CE^2 / ΣEtot), [i], [j]) * term($(@. TF²      / ΣEtot)     )
                + term($(@. TF² * CE   / ΣEtot), [i]     ) * term($(@. TF² * CE / ΣEtot), [j])
            )
            Ωpix * pinv(Diagonal(sum(∇[i].diag .* ∇[j].diag .* Fourier(A1(i,j) + A2(i,j)) for i=1:2,j=1:2)))
        end
    end
    Nϕ = AL # true only for unlensed weights

    ϕqe = (wiener_filtered ? (Cϕ*pinv(Cϕ+Nϕ)) : 1) * (AL*ϕqe_unnormalized)
    @namedtuple ϕqe AL Nϕ

end

function quadratic_estimate_EB((d1,d2)::NTuple{2,FlatS2}, Cf, Cf̃, Cn, Cϕ, TF, wiener_filtered, weights, AL=nothing; zeroB=false)

    term = get_term_memoizer(d1[:E])
    CE, CB = getindex.(Ref((weights==:unlensed) ? Cf : Cf̃),(:E,:B))
    TF²E, TF²B = TF[:E]^2, TF[:B]^2
    ΣEtot = TF²E * Cf̃[:E] + Cn[:E]
    ΣBtot = TF²B * Cf̃[:B] + Cn[:B]
    @unpack Ωpix = fieldinfo(d1)

    # unnormalized estimate
    ϕqe_unnormalized = @subst begin
        I(i) = 2 * sum(  ϵ(k,l,3) * (
                           term($(CE * (ΣEtot \ (TF*d1)[:E])), [i], j, k) * term($(     (ΣBtot \ (TF*d2)[:B])),      j, l)
            - (zeroB ? 0 : term($(      ΣEtot \ (TF*d1)[:E]),       j, k) * term($(CB * (ΣBtot \ (TF*d2)[:B])), [i], j, l)))
            for (j,k,l) in inds(3)
        )
        sum(∇[i] * Fourier(I(i)) for i=1:2)
    end

    # normalization
    if AL == nothing
        AL = @subst begin
            @sym_memo A(@sym(i,j)) = 4 * sum( ϵ(m,p,3) * ϵ(n,q,3) * (
                                 term($(@. TF²E * CE^2 / ΣEtot), [i], [j], k, l, m, n) * term($(@. TF²B        / ΣBtot),           k, l, p, q)
                + (zeroB ? 0 : -2term($(@. TF²E * CE   / ΣEtot), [i],      k, l, m, n) * term($(@. TF²B * CB   / ΣBtot),      [j], k, l, p, q))
                + (zeroB ? 0 :   term($(@. TF²E        / ΣEtot),           k, l, m, n) * term($(@. TF²B * CB^2 / ΣBtot), [i], [j], k, l, p, q)))
                for (k,l,m,n,p,q) in inds(6)
            )
            Ωpix * pinv(Diagonal(sum(∇[i].diag .* ∇[j].diag .* Fourier(A(i,j)) for i=1:2,j=1:2)))
        end
    end
    Nϕ = AL # true only for unlensed weights

    ϕqe = (wiener_filtered ? (Cϕ*pinv(Cϕ+Nϕ)) : 1) * (AL * ϕqe_unnormalized)
    @namedtuple ϕqe AL Nϕ

end


text:
using DynamicPolynomials
using MultivariateMoments
using SemialgebraicSets
using CSDP

# -----------------------------------------------------------------------------
# Jiawang Nie, Zi Yang, & Guangming Zhou. (2018).
# The Saddle Point Problem of Polynomials.
# -----------------------------------------------------------------------------

@testset "Bad Saddle Example 6.1 (I)" begin
        @polyvar x[1:3] y[1:3]
        p = x[1]x[2] + x[2]x[3] + x[3]y[1] + x[1]y[3] + y[1]y[2] + y[2]y[3]
        Sx = @set sum(x) == 1 && x[1] >= 0 && x[2] >= 0 && x[3] >= 0
        Sy = @set sum(y) == 1 && y[1] >= 0 && y[2] >= 0 && y[3] >= 0
        opt = CSDP.Optimizer

        expected_measure_x = AtomicMeasure(
                variables(Sx),
                [WeightedDiracMeasure([0, 1, 0], 1)],
        )
        expected_measure_y = AtomicMeasure(
                variables(Sy),
                [
                        WeightedDiracMeasure([0.02, 0.5, 0.48], 0.3),
                        WeightedDiracMeasure([0.25, 0.5, 0.25], 0.3),
                        WeightedDiracMeasure([0.48, 0.5, 0.02], 0.3),
                ],
        )
        expected_value = 0.25

        pvalue, pmeasure_y, pmeasure_x =
                solve_game(+p, Sy, Sx, opt, iteration = 0x2)
        dvalue, dmeasure_x, dmeasure_y =
                solve_game(-p, Sx, Sy, opt, iteration = 0x1)

        eps = 1e-3
        @test isapprox(expected_value, pvalue, atol = eps)
        @test isapprox(dvalue, -pvalue, atol = eps)
        @test isapprox(expected_measure_y, pmeasure_y, atol = 1e-1) # WARNING!
        @test isapprox(expected_measure_x, dmeasure_x, atol = eps)
end

@testset "Bad Saddle Example 6.2 (I)  " begin
        @polyvar x[1:2] y[1:2]
        p =
                (x[1] + x[2] + y[1] + y[2] + 1)^2 -
                4 * (x[1]x[2] + x[2]y[1] + y[1]y[2] + y[2] + x[1])
        Sx = @set x[1] >= 0 && x[2] >= 0 && 1 - x[1] >= 0 && 1 - x[2] >= 0
        Sy = @set y[1] >= 0 && y[2] >= 0 && 1 - y[1] >= 0 && 1 - y[2] >= 0
        opt = CSDP.Optimizer

        expected_measure_y =
                AtomicMeasure(variables(Sy), [WeightedDiracMeasure([1, 0], 1)])
        expected_value = 4

        pvalue, pmeasure_y, pmeasure_x =
                solve_game(+p, Sy, Sx, opt, iteration = 0x1)
        dvalue, dmeasure_x, dmeasure_y =
                solve_game(-p, Sx, Sy, opt, iteration = 0x1)

        eps = 1e-3
        @test isapprox(expected_value, pvalue, atol = eps)
        @test isapprox(dvalue, -pvalue, atol = eps)
        @test isapprox(expected_measure_y, pmeasure_y, atol = eps)
end

@testset "Unknown Saddle Example 6.2 (II)  " begin
        @polyvar x[1:3] y[1:3]
        p =
                sum([(x[i] + y[i]) for i in 1:3]) - prod([(x[i] - y[i]) for i in 1:3])
        Sx = @set x[1] >= 0 &&
             x[2] >= 0 &&
             x[3] >= 0 &&
             1 - x[1] >= 0 &&
             1 - x[2] >= 0 &&
             1 - x[3] >= 0
        Sy = @set y[1] >= 0 &&
             y[2] >= 0 &&
             y[3] >= 0 &&
             1 - y[1] >= 0 &&
             1 - y[2] >= 0 &&
             1 - y[3] >= 0
        opt = CSDP.Optimizer

        expected_measure_x = AtomicMeasure(
                variables(Sx),
                [
                        WeightedDiracMeasure([0, 1, 0], 1 / 3),
                        WeightedDiracMeasure([1, 0, 0], 1 / 3),
                        WeightedDiracMeasure([0, 0, 1], 1 / 3),
                ],
        )
        expected_measure_y = AtomicMeasure(
                variables(Sy),
                [WeightedDiracMeasure([1, 1, 1], 1)],
        )
        expected_value = 4

        pvalue, pmeasure_y, pmeasure_x =
                solve_game(+p, Sy, Sx, opt, iteration = 0x1)
        dvalue, dmeasure_x, dmeasure_y =
                solve_game(-p, Sx, Sy, opt, iteration = 0x1)

        eps = 1e-3
        @test isapprox(expected_value, pvalue, atol = eps)
        @test isapprox(dvalue, -pvalue, atol = eps)
        @test isapprox(expected_measure_y, pmeasure_y, atol = eps)
        @test isapprox(expected_measure_x, dmeasure_x, atol = eps)
end

#=
#@testset "Saddle Example 6.4 (II)  " begin
#	no saddle point - useless test
#	@polyvar x[1:3] y[1:3]
#	p = x[1]^2y[1]^2 + x[2]^2y[2]^2 + x[3]^2y[3]^2 + x[1]^2y[2]y[3] + x[2]^2y[1]y[3] + x[3]^2y[1]y[2] + y[1]^2x[2]x[3] + y[2]^2x[1]x[3] + y[3]^2x[1]x[2]
#	Sx = @set x'x == 1
#	Sy = @set y'y == 1
#       opt = CSDP.Optimizer
#
#	expected_value_y = 0.7666
#
#	actual_value_x, actual_measure_x = solve_game(-p, Sy, Sx, opt, iteration=1)
#	actual_value_y, actual_measure_y = solve_game(p, Sx, Sy, opt, iteration=1)
#
#	eps = 1e-3
#	@test isapprox(actual_value_x, -actual_value_y, atol=eps)
#	@test isapprox(expected_value_y, actual_value_y, atol=eps)
#end

=#


text:
"""
   update_check_regularity!(state, agent_idx, config)

Update an agent's check regularity depending on its perceived public opinion

# Arguments
- `state`: a tuple of the current graph and agent_list
- `agent_idx`: agent index
- `config`: Config object as provided by Config()

See also: [`Config`](@ref), [`update__perceiv_publ_opinion!`](@ref)
"""
function update_check_regularity!(
    state::Tuple{AbstractGraph, AbstractArray},
    agent_idx::Integer,
    config::Config
)

    agent_list = state[2]
    this_agent = agent_list[agent_idx]

    if (
        abs(this_agent.opinion - this_agent.perceiv_publ_opinion)
        > config.opinion_threshs.check_unease
    )
        this_agent.check_regularity = (
            config.agent_props.check_decrease * this_agent.check_regularity
        )
    else
        this_agent.check_regularity = 1.0
    end

    return state
end

# suppress output of include()
;


text:
_update(::Val{false}, args...) = update(args...)
_update(::Val{true}, args...) = update_training(args...)
_getindex(is_training::Nothing, s) = false
_getindex(is_training, s) = is_training[s]

"""
    stopping_time(criterion, losses; verbosity=0)
    stopping_time(criterion, losses, is_training; verbosity=0)

Determine the stopping time for the iterator `losses`, given
`stopping::StoppingCriterion`. Include the `Bool` vector `is_training`
of matching length, when there is a distinction between
"out-of-sample" losses and "training" losses.

If losses completes before a stop, then `0` is returned.

```
julia> stopping_time(NotANumber(), [10.0, 3.0, NaN, 4.0])
3

julia> stopping_time(NotANumber(), [10.0, 3.0, 5.0, 4.0])
0
```
"""
function stopping_time(criterion, losses, training; verbosity=0)

    t_stop = 0 # meaning no stop
    t = 0 # counts regular `update` calls but ignores `update_training` calls
    s = 0 # counter for iteration over `losses`

    is_training = training === nothing ?
        nothing : collect(training)
    global state

    for loss in losses
        s += 1
        _is_training = _getindex(is_training, s)
        state = if s == 1
            _update(Val(_is_training), criterion, loss)
        else
            _update(Val(_is_training), criterion, loss, state)
        end
        _is_training || (t += 1)
        verbosity < 1 || begin
            @info "loss updates: $t"
            @info "state: $state"
        end
        if !_is_training && done(criterion, state)
            t_stop = t
            break
        end
    end

    return t_stop

end

stopping_time(criterion, losses; kwargs...) =
    stopping_time(criterion, losses, nothing; kwargs...)


text:
"""
METROPOLIS-HASTING

The following codes explore the use of MH for Bayesian Inference.
"""
immutable MH
  logposterior::Function
  proposal    ::Function
  init_est    ::Array{Float64}
  d           ::Int64

  function default_proposal(theta::Array{Float64}, d::Int64)
    sigmas = ones(d);
    random = Float64[]

    for i in 1:length(theta)
      push!(random, rand(Normal(theta[i], sigmas[i])))
    end

    return random
  end

  MH(logposterior = Union{}; proposal::Function = default_proposal, init_est::Array{Float64} = [0.; 0.], d::Int64 = 2) = new(logposterior, proposal, init_est, d)
end

function mcmc(parameters::MH; r = 1000)
  chain = zeros(r, parameters.d)
  chain[1, :] = parameters.init_est

  for i in 1:(r - 1)
    propose = parameters.proposal(chain[i, :], parameters.d)
    probab = exp(parameters.logposterior(propose) - parameters.logposterior(chain[i, :]))

    if rand(Uniform()) < probab
      chain[i + 1, :] = propose
    else
      chain[i + 1, :] = chain[i, :]
    end
  end

  return chain

end


text:
using SchurPolynomials
using Test

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


text:
# This is an image comparison via time domain
# subtraction in a sliding wavelet sized window.
#
# author: GRAM 5/1/16

function compare(s_p, s_n, s_o, sei, i, j, k, wvln, layer, verbose)

	# split the interface sample
	ind_1 = Int(layer[k + 1] + ceil(wvln / 2))
	ind_2 = Int(layer[k + 1] - floor(wvln / 2))

	if verbose == 2
		print("\nwindow goes $ind_2 to $ind_1\n")
	end

	# find error
	dif_p = abs(sei[ind_2 : ind_1, i, j] - s_p[ind_2 : ind_1])
	dif_n = abs(sei[ind_2 : ind_1, i, j] - s_n[ind_2 : ind_1])
	dif_o = abs(sei[ind_2 : ind_1, i, j] - s_o[ind_2 : ind_1])

	# sum the difference vectors
	return [sum(dif_p),sum(dif_n),sum(dif_o)]

end


text:
# version of library to download
const version = v"3.1.2"
const glfw = "glfw-$version"

# TODO: if the latest version is already installed, don't bother with any of this

# download and compile the library from source
@linux_only begin
	mkpath("downloads")
	const tarball = "downloads/$glfw.tar.gz"
	if !isfile(tarball)
		info("Downloading GLFW $version source tarball")
		download("https://github.com/glfw/glfw/archive/$version.tar.gz", tarball)
	end
	mkpath("src")
	run(`tar xzf $tarball -C src`)
	map(mkpath, ("builds/$glfw", "usr$WORD_SIZE"))
	try
		info("Building GLFW $version library from source")
		cd("builds/$glfw") do
			options = map(x -> "-D$(x[1])=$(x[2])", [
				("BUILD_SHARED_LIBS",    "ON"),
				("CMAKE_INSTALL_PREFIX", "../../usr$WORD_SIZE"),
				("GLFW_BUILD_DOCS",      "OFF"),
				("GLFW_BUILD_EXAMPLES",  "OFF"),
				("GLFW_BUILD_TESTS",     "OFF")
			])
			run(`cmake $options ../../src/$glfw`)
			run(`make install`)
		end
	catch e
		warn(e)
		warn("""
On Linux, you might not have xorg-dev and libglu1-mesa-dev installed, which you can get via the commandline with for example
sudo apt-get install  xorg-dev libglu1-mesa-dev
If that doesn't help, try to install GLFW manually
(see http://www.glfw.org/download.html for more information)""")
	end
end

# download a pre-compiled binary (built by Bintray for Homebrew)
@osx_only begin
	const osx_version = convert(VersionNumber, readall(`sw_vers -productVersion`))
	if osx_version >= v"10.11"
		codename = "el_capitan"
	elseif osx_version >= v"10.10"
		codename = "yosemite"
	else
		codename = "mavericks"
	end
	const tarball = "glfw3-$version.$codename.bottle.tar.gz"
	if (!isfile("downloads/$tarball"))
		info("Downloading pre-compiled GLFW $version binary")
		mkpath("downloads")
		download("https://homebrew.bintray.com/bottles-versions/$tarball", "downloads/$tarball")
	end
	mkpath("builds")
	run(`tar xzf downloads/$tarball -C builds`)
	mkpath("usr64")
	cp("builds/glfw3/$version/lib", "usr64/lib", remove_destination=true)
end

# download a pre-compiled binary (built by GLFW)
@windows_only begin
	mkpath("downloads")
	for (sz, suffix) in ((32, ""), (64, "-w64"))
		const build = "$glfw.bin.WIN$sz"
		const archive = "downloads/$build.zip"
		if !isfile(archive)
			info("Downloading pre-compiled GLFW $version $sz-bit binary")
			download("https://github.com/glfw/glfw/releases/download/$version/$build.zip", archive)
		end
		run(`7z x -obuilds -y $archive`)
		dllpath = "builds/$build/lib-mingw$suffix"
		mkpath("usr$sz")
		cp(dllpath, "usr$sz/lib", remove_destination=true)
	end
end


text:
@testset "@map macro(ex)" begin
    n = 3
    A = rand(2, 2);
    b = rand(2)
    identitymap = @map x -> I(n)*x
    @test identitymap == MathematicalSystems.IdentityMap(n)
    linearmap = @map x -> A*x
    @test linearmap == MathematicalSystems.LinearMap(A)
    affinemap = @map x -> A*x + b
    @test affinemap == MathematicalSystems.AffineMap(A, b)
end


text:
module LazyTaylorSeries

using StaticArrays

export Taylor1, Taylor, tt
export variable, variables, constant, degree, evaluate!
export reset!

import Base:
    +, -, *,
    exp,
    getindex

include("powers.jl")
include("taylor1.jl")
include("taylorN.jl")

variable(T) = Taylor1(T, (t, i) -> (i == 1) ? one(T) : zero(T))

end

max_stars_repo_name: dpsanders/LazyTaylorSeries

max_issues_repo_name: dpsanders/LazyTaylorSeries

max_forks_repo_name: dpsanders/LazyTaylorSeries


text:
__precompile__()

module ECC

import Base: +, -, *, ^, /, ==, inv, sqrt, show, div
export S256Point, Signature, PrivateKey,
       point2sec, sec2point, verify, pksign, sig2der, der2sig,
       ∞, int2bytes, bytes2int

include("helper.jl")
include("primefield.jl")
include("infinity.jl")
include("point.jl")
include("signature.jl")
include("privatekey.jl")
include("scep256k1.jl")

end # module


text:
using ONNXRunTime
using Documenter

DocMeta.setdocmeta!(ONNXRunTime, :DocTestSetup, :(using ONNXRunTime); recursive=true)

cp(joinpath(@__DIR__, "..", "README.md"), joinpath(@__DIR__, "src", "index.md"), force=true)

makedocs(;
    modules=[ONNXRunTime],
    authors="Jan Weidner <jw3126@gmail.com> and contributors",
    repo="https://github.com/jw3126/ONNXRunTime.jl/blob/{commit}{path}#{line}",
    sitename="ONNXRunTime.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jw3126.github.io/ONNXRunTime.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "API" => "api.md",
        "Explanation" => "explanation.md",
    ],
)

deploydocs(;
    repo="github.com/jw3126/ONNXRunTime.jl",
    devbranch = "main",
)


text:
export SkipBlock

"""
$(TYPEDEF)

TODO (vtjeng)

## Fields:
$(FIELDS)
"""
struct SkipBlock <: Layer
    layers::Array{<:Layer}
end

function Base.show(io::IO, p::SkipBlock)
    println(io, "SkipBlock")
    for (index, value) in enumerate(p.layers)
        println(io, "    ($index) $value")
    end
end

# TODO (vtjeng): better typing on xs
function apply(p::SkipBlock, xs::Array)
    num_layers = length(p.layers)
    inputs = xs[end-num_layers+1:end]
    outputs = map((f, x) -> f(x), p.layers, inputs)
    return reduce((x1, x2) -> x1 .+ x2, outputs)
end

(p::SkipBlock)(xs::Array) = apply(p, xs)


text:
# Hello world gtk-server demo in Julia
#

using GtkServer

gtk_init()

gtk_send("gtk_init NULL NULL")
win = gtk_send("gtk_window_new 0")
gtk_send("gtk_window_set_title $win 'Julia GTK-Server demo'")
gtk_send("gtk_window_set_default_size $win 400 200")
tbl = gtk_send("gtk_table_new 10 10 1")
gtk_send("gtk_container_add $win $tbl")
info = gtk_send("gtk_button_new_with_label 'Click here!'")
gtk_send("gtk_table_attach_defaults $tbl $info 1 4 1 9")
but = gtk_send("gtk_button_new_with_label '...or click to quit.'")
gtk_send("gtk_table_attach_defaults $tbl $but 6 9 1 9")
gtk_send("gtk_widget_show_all $win")
dia = gtk_send("u_dialog Information '\"Hello, World!\"' 200 130")

# Main loop
event=0
while true
    event = gtk_send("gtk_server_callback WAIT")
    # The user clicked on the Hello World button.
    if event == dia
        gtk_send("gtk_widget_hide $dia")
    # The user clicked to see the dialog.
    elseif event == info
        gtk_send("gtk_widget_show_all $dia")
    # The user click on the quit button.
    elseif event == but || event == win
        gtk_exit()
        break
    end
end


text:
# The goal of this file is to provide tools that implements the MFK algorithm from these two references:
# @article{furman2017,
#   title = {On {{Log}}-{{Normal Convolutions}}: {{An Analytical}}-{{Numerical Method With Applications}} to {{Economic Capital Determination}}},
#   shorttitle = {On {{Log}}-{{Normal Convolutions}}},
#   author = {Furman, Edward and Hackmann, Daniel and Kuznetsov, Alexey},
#   year = {2017},
#   issn = {1556-5068},
#   abstract = {We put forward an efficient algorithm for approximating the sums of independent and lognormally distributed random variables. Namely, by merging tools from probability theory and numerical analysis, we are able to compute the cumulative distribution functions of the just-mentioned sums with any desired precision. Importantly, our algorithm is fast and can tackle equally well sums with just a few or thousands of summands. We illustrate the effectiveness of the new method in the contexts of the individual and collective risk models, aggregate economic capital determination, and economic capital allocation.},
#   journal = {SSRN Electronic Journal},
#   language = {en}
# }

# @article{miles2019,
#   ids = {mile,milesa},
#   title = {Risk Aggregation: {{A}} General Approach via the Class of {{Generalized Gamma Convolutions}}},
#   author = {Miles, Justin and Furman, Edward and Kuznetsov, Alexey},
#   year = {2019},
#   journal = {Variance},
#   keywords = {gamma}
# }

"""

    compute_g(dist,n, integrator)

This function will compute the n first theretical (-1)-exponentially shifted moments of a distribution,
using takashi-mori tanh-sinh exponentials formulas from the DoubleExponentialFormulas package.

This might take a long time. As an integrator, you should pass the result of `QuadDE(BigFloat)` after setting enough precision.

"""
function compute_g(dist,n, integrator)
    g = Array{BigFloat}(undef, 1,2n+1)
    residuals_g = deepcopy(g)
    for i in 0:(2n)
        g[i+1],residuals_g[i+1] = integrator(x -> (-x)^(i) * Distributions.pdf(dist,x) * exp(-x), 0, +Inf)
        print("g_{",i,"} = ",Float64(g[i+1]),", rez = ",Float64(residuals_g[i+1]),"\n")
    end
    return g
end

"""

    E_from_g(g)

Simply computes the empirical laguerre coefficients from the exponentially shifted moments,
if you computed those moments from a theoretical distribution.

"""
function E_from_g(g)
    # this function should compute laguerre coefficients from g.
    a = deepcopy(g)
    for i in 1:length(a)
        a[i] = sqrt(2) * sum(binomial(big(i-1),big(k))*big(2)^(k)/factorial(big(k)) * g[k+1] for k in 0:(i-1))
    end
    return a
end

function build_s!(s,g,facts)
    s[1] = -g[2]/g[1]
    for k in 1:(length(s)-1)
        s[k+1] = g[k+2] / facts[k+1]
        for i in 0:(k-1)
            s[k+1] += s[i+1] * g[k-i+1] / facts[k-i+1]
        end
        s[k+1] = - s[k+1]/g[1]
    end
end

function s_from_k(k)
    if isodd(length(k))
        n = Int((length(k)-1)/2)
    else
        n = Int(length(k)/2) - 1
    end
    # n = Int(round(length(k)/2))
    s = k[2:2n+1] .* (-1) .^(2:(2n+1))
    return s
end

function MFK_end(s::Vector{T}) where T
    n = Int(length(s)//2)
    S = zeros(T,(n,n))
    for i in 0:(n-1)
        for j in 0:(n-1)
            S[i+1,j+1] = s[i+j+1]
        end
    end

    sol_b = LinearAlgebra.Symmetric(S) \ (-s[(n+1):end])
    b = deepcopy(sol_b)
    append!(b,1)
    b = reverse(b)
    b_deriv = reverse(sol_b) .* (1:n)

    a = zeros(T,n)
    a[1] = s[1]
    for k in 1:(n-1)
        a[k+1] = s[k+1]
        for i in 0:(k-1)
            a[k+1] = a[k+1] + b[i+2] * s[k-i]
        end
    end

    z = real.(PolynomialRoots.roots(b, polish=true))

    beta = -z .-1
    alpha = deepcopy(beta)

    for i in 1:length(alpha)
        rez_num = 0
        rez_denom = 0
        for k in 1:length(a)
            rez_num += a[k] * z[i]^(k-1)
        end
        for k in 1:length(b_deriv)
            rez_denom += b_deriv[k] * z[i]^(k-1)
        end
        alpha[i] = rez_num/rez_denom
    end
    theta = 1 ./beta
    return alpha,theta
end

"""

    MFK_Projection(g_integrals,n_gammas)

From a set of g_integrals computed by the compute_g function, this function runs the algorithm from Miles, Furman & Kuznetsov
to produce a univariate gamma convolution.

This algorithm works only if the distribution you computed the g_integrals from is, indeed, a generalized gamma convolution.
otherwise, it might fail and return garbage.

"""
function MFK_Projection(g_integrals::Vector{T},n) where T
    facts = factorial.(T.(0:2n))
    s = Array{T}(undef, 2n)
    build_s!(s,g_integrals,facts)
    alpha,beta = MFK_end(s,n)
    return ThorinDistributions.UnivariateGammaConvolution(Float64.(alpha), Float64.(1 ./ beta))
end


text:
"""
Pluto2Jupyter module converts Pluto Notebooks to Jupyter Notebooks and viceversa

jl2nb("myfile.jl")
to convert pluto Notebook to  "myfile.ipynb" Jupyter Notebook

nb2jl("myfile.ipynb")
to convert Jupyter Notebook to "myfile.jl" Pluto Notebook
"""
module Pluto2Jupyter
export nb2jl, jl2nb
using JSON, UUIDs, Pluto

"""
    nb2jl(path::AbstractString)

Convert Jupyter Notebook to Pluto Notebook
"""
function nb2jl(path::AbstractString)

    nb = open(JSON.parse, path, "r")

    # test if notebook is acceptable
    nb["nbformat"] == 4 || error("unrecognized version", nb["nbformat"])
    lang = lowercase(nb["metadata"]["language_info"]["name"])
    lang == "julia" || error("notebook is for unregognized language $lang")

    # Shell or help mode
    sh_mode = r"^\s*[;?]"

    nbcells = Array{Pluto.Cell}(undef,0)
    for cell in nb["cells"]
        if cell["cell_type"] == "markdown"
            code = string("md\"\"\"",join(cell["source"]),"\n\"\"\"")
        elseif cell["cell_type"] == "code"
            code = join(cell["source"])
            occursin(sh_mode, code) && continue  # skip ? and ;
        end
        push!(nbcells,Pluto.Cell(uuid1(),code))
    end

    savepath = string(path[1:end-5],"jl")
    plutonb = Pluto.Notebook(nbcells,savepath)
    Pluto.save_notebook(plutonb,savepath);
end

"""
    jl2nb(path::AbstractString)

Convert Pluto Notebook to Jupyter Notebook
"""
function jl2nb(path::AbstractString)
    plutonb = Pluto.load_notebook_nobackup(path)
    nbcells = []
    counter = 1
    for cell in plutonb.cells
        tempcell = Dict("cell_type"=> "code",
        "metadata"=> Dict(),
        "source"=>[])

        if occursin(r"^(md\"\"\")",cell.code)
            tempcell["cell_type"]= "markdown"
            push!(tempcell["source"],cell.code[6:end-3])
        else
            push!(tempcell["source"],cell.code)
            tempcell["outputs"]=[]
            tempcell["execution_count"]= counter
            counter +=1
        end
        push!(nbcells,tempcell)
    end

    # Todo: get version and spec from the system
    language_info = Dict("file_extension" => ".jl",
    "mimetype" => "application/julia",
    "name" => "julia",
    "version" => "1.5.3")
    kernelspec = Dict(
        "display_name"=> "Julia 1.5.3",
        "language"=> "julia",
        "name"=> "julia-1.5"
       )
    metadata = Dict("language_info"=>language_info,
        "kernelspec"=>kernelspec)

    nb = Dict("cells"=>nbcells,
            "metadata" => metadata,
            "nbformat" => 4,
            "nbformat_minor" => 2)

    savepath = string(path[1:end-2],"ipynb")
    open(savepath,"w") do f
        write(f,JSON.json(nb))
    end
end

end


text:
abstract type AbstractIndirectOptimizationProblem end

# Not implementing this for now
# May want to allow users to only suppply bvps without STM and use AD for STM computation
struct UserProvidedIndirectOptimizationProblem <: AbstractIndirectOptimizationProblem
end

struct IndirectOptimizationProblem{IBVPF, BVPF, BVPSTMF} <: AbstractIndirectOptimizationProblem
    # BVP for initialization phase (likely to include additional callbacks for early termination)
    initBVPFunc::IBVPF

    # BVP without STM
    BVPFunc::BVPF

    # BVP with STM
    BVPWithSTMFunc::BVPSTMF

    # Initial conditions
    iConds::Vector{Float64}

    # Final conditions
    fConds::Vector{Float64}

    # Flag to indicate if homotopy continuation will be employed when solving BVP.
    # If true, BVP functions must accept continuation parameter
    homotopy::Bool
end

# TSPAN in units of days!!! Scaled for respective problems!
function IndirectOptimizationProblem(scenario::String, iConds::AbstractVector, fConds::AbstractVector, tspan::Tuple;
    homotopy = true, homotopyFlag::HomotopyFlag = MEMF(), initializationFlag::InitializationFlag = Initialization())

    # Check for scenario type
    if occursin("CR3BP", scenario)
        # Initialize parameters
        ps = initCR3BPIndirectParams(scenario)
        psSTM = initCR3BPIndirectWithSTMParams(scenario)

        # Scale tspan
        daysToTU = 86400/ps.crp.TU
        tspanScalled = (tspan[1]*daysToTU, tspan[2]*daysToTU)

        # Initialize functions
        initBVPFunc(y0) = integrate(y0, tspanScalled, ps, CR3BP(), initializationFlag, homotopyFlag;
                                    copyParams = true, termCallbacks = true)
        if homotopy == true
            BVPFunc(y0, ϵ)        = integrateWithHomotopy(y0, tspanScalled, ϵ, ps, CR3BP(), Solving(), homotopyFlag;
                                                          copyParams = true)
            BVPWithSTMFunc(z0, ϵ) = integrateWithHomotopy(z0, tspanScalled, ϵ, psSTM, CR3BP(), SolvingWithSTM(),
                                                          homotopyFlag; copyParams = true)
        else
            BVPFunc(y0)         = integrate(y0, tspanScalled, ps, CR3BP(), Solving(), homotopyFlag; copyParams = true)
            BVPWithSTMFunc(z0)  = integrate(z0, tspanScalled, psSTM, CR3BP(), SolvingWithSTM(), homotopyFlag; copyParams = true)
        end
    else
        throw(ArgumentError("Only CR3BP scenarios are implemented now."))
    end

    IndirectOptimizationProblem{typeof(initBVPFunc), typeof(BVPFunc), typeof(BVPWithSTMFunc)}(
        initBVPFunc, BVPFunc, BVPWithSTMFunc, iConds, fConds, homotopy
    )
end


text:
# This function calculates the free kurtosis of a matrix
# The function is uniform for both Hermitian or rectangular X
"""
    κ₄(X)

This function return the free kurtosis of a matrix X

# Arguments
-   `X`: a matrix, Hermitian or rectangular

# Outputs
-   `fk`: a scalar, the free kurtosis of X
"""
function κ₄(X)
    # get dimension
    N, M = size(X)

    # compute the free kurtosis
    fk = tr((X*X')^2)/N - (1 + N/M) * (tr(X*X')/N)^2
    return fk
end

"""
    neg_abs_sum_free_kurt(Z)

This function calculates the sum of negative absolute value of free kurtosis of components of Z
"""
function neg_abs_sum_free_kurt(Z)
    #--------------------------------------------------------------------------
    # Syntax:       F = neg_abs_sum_free_kurt(Z)
    #
    # Input:        Z: a array of Hermitian or rectangular matrices
    #
    # Outputs:      F: a scalr, it is the summation of negative absolute values
    #               of the free kurtosis of Z[i]
    #
    # Description:  This function return the summation of negative absolute value
    #               of free kurtosis of input array of matrices
    #
    # Authors:      Hao Wu
    #               lingluanwh@gmail.com
    #
    # Date:         Jan 20, 2019
    #--------------------------------------------------------------------------

    return -sum(abs.(κ₄.(Z)))
end

"""
    free_ent(X; mat = "her")

This function return the free entropy of a matrix X

# Arguments
-   `X`: a matrix, Hermitian or rectangular
-   `mat`: the type of X, valid options are "her" and "rec"

# Outputs
-   `chi`: a scalr, it is the free entropy of X
"""
function free_ent(X; mat = "her")
    # get dimension
    N, M = size(X)

    # compute the entropy
    # self-adjoint case
    if mat == "her"
        Λ = eigvals(Hermitian(X))
        chi = mean(log.(abs.([Λ[i] - Λ[j] for i = 1:N for j = (i + 1):N])))
    end

    # rectangular case
    if mat == "rec"
        a, b = N/(N + M), M/(N + M)
        Λ = svdvals(X*X')
        chi = a^2 * mean(log.(abs.([Λ[i] - Λ[j] for i = 1:N for j = (i + 1):N]))) + a*(b - a)*mean(log.(abs.(Λ)))
    end

    return chi
end

"""
    sum_free_ent(Z; mat = "her")

This function calculates the sum of free entropy of components of Z
"""
function sum_free_ent(Z; mat = "her")
    return sum(free_ent.(Z, mat = mat))
end


text:
# # Soil types

## load packages
using PlantHydraulics
using PlotPlants
FT = Float32;
#------------------------------------------------------------------------------

## prescribed soil types
_sh1 = BrooksCorey{FT}();
_sh2 = VanGenuchten{FT}();
_STS = ["Sand", "Loamy Sand", "Sandy Loam", "Loam", "Sandy Clay Loam",
        "Silt Loam", "Silt", "Clay Loam", "Silty Clay Loam", "Sandy Clay",
        "Silty Clay", "Clay"];

function compare_soil_vc(sn::Int)
    _bc = create_soil_VC(_sh1, _STS[sn]);
    _vg = create_soil_VC(_sh2, _STS[sn]);
    @show _bc.b, _bc.ϕs;

    _fig,_axes = create_canvas(_STS[sn]; ncol=2);
    _ax1,_ax2  = _axes;

    _Θs = range(_bc.Θr+FT(1e-3); stop=_bc.Θs, length=50);
    _P1 = -1 .* soil_p_25_swc.([_bc], _Θs);
    _P2 = -1 .* soil_p_25_swc.([_vg], _Θs);
    _K1 = soil_k_ratio_swc.([_bc], _Θs);
    _K2 = soil_k_ratio_swc.([_vg], _Θs);

    _ax1.plot(_Θs, _P1, "b-", label="Brooks Corey" );
    _ax1.plot(_Θs, _P2, "r-", label="van Genuchten");
    _ax2.plot(_Θs, _K1, "b-");
    _ax2.plot(_Θs, _K2, "r-");
    _ax1.set_yscale("log");
    _ax1.legend(loc="upper right");
    set_xlabels!(_axes, ["Soil Water Content (-)" for i in 1:2]);
    set_ylabels!(_axes, ["Metric potential (-MPa)", "Relative K"]);
    _fig.set_tight_layout(true);

    return _fig
end
#------------------------------------------------------------------------------

# ## Pre-defined soil types
# ### Sand
compare_soil_vc(1)
#------------------------------------------------------------------------------

# ### Loamy Sand
compare_soil_vc(2)
#------------------------------------------------------------------------------

# ### Sandy Loam
compare_soil_vc(3)
#------------------------------------------------------------------------------

# ### Loam
compare_soil_vc(4)
#------------------------------------------------------------------------------

# ### Sandy Clay Loam
compare_soil_vc(5)
#------------------------------------------------------------------------------

# ### Silt Loam
compare_soil_vc(6)
#------------------------------------------------------------------------------

# ### Silt
compare_soil_vc(7)
#------------------------------------------------------------------------------

# ### Clay Loam
compare_soil_vc(8)
#------------------------------------------------------------------------------

# ### Silty Clay Loam
compare_soil_vc(9)
#------------------------------------------------------------------------------

# ### Sandy Clay
compare_soil_vc(10)
#------------------------------------------------------------------------------

# ### Silty Clay
compare_soil_vc(11)
#------------------------------------------------------------------------------

# ### Clay
compare_soil_vc(12)
#------------------------------------------------------------------------------


text:
module TestUtils

export random_complex_matrix, random_real_matrix, random_hermitian_matrix
export random_complex_sparse_matrix, random_real_sparse_matrix
export random_hermitian_sparse_matrix, random_state_vector

using Random
using Distributions
using LinearAlgebra
using SparseArrays

"""Construct a random complex matrix of size N×N with spectral radius ρ.

```julia
random_complex_matrix(N, ρ)
```
"""
function random_complex_matrix(N, ρ)
    σ = 1/√N
    d = Normal(0.0, σ)
    H = ρ * (rand(d, (N, N)) + rand(d, (N, N)) * 1im) / √2
end

"""Construct a random real-valued matrix of size N×N with spectral radius ρ.

```julia
random_real_matrix(N, ρ)
```
"""
function random_real_matrix(N, ρ)
    σ = 1/√N
    d = Normal(0.0, σ)
    H = ρ * rand(d, (N, N))
end

"""Construct a random Hermitian matrix of size N×N with spectral radius ρ.

```julia
random_hermitian_matrix(N, ρ)
```
"""
function random_hermitian_matrix(N, ρ)
    σ = 1/√N
    d = Normal(0.0, σ)
    X = rand(d, (N, N))
    H = ρ * (X + X') / (2*√2)
end

"""Construct a random sparse complex matrix.

```julia
random_complex_sparse_matrix(N, ρ, sparsity)
```

returns a matrix of size N×N with spectral radius ρ and the given sparsity
(number between zero and one that is the approximate fraction of non-zero
elements).
"""
function random_complex_sparse_matrix(N, ρ, sparsity)
    σ = 1/√(sparsity * N)
    d = Normal(0.0, σ)
    Hre = sprand(N, N, sparsity, (dims...) -> rand(d, dims...))
    Him = sprand(N, N, sparsity, (dims...) -> rand(d, dims...))
    H = ρ * (Hre + Him * 1im) / √2
end

"""Construct a random sparse real-valued matrix.

```julia
random_real_sparse_matrix(N, ρ, sparsity)
```

returns a matrix of size N×N with spectral radius ρ and the given sparsity
(number between zero and one that is the approximate fraction of non-zero
elements).
"""
function random_real_sparse_matrix(N, ρ, sparsity)
    σ = 1/√(sparsity * N)
    d = Normal(0.0, σ)
    H = ρ * sprand(N, N, sparsity, (dims...) -> rand(d, dims...))
end

"""Construct a random sparse Hermitian matrix.

```julia
random_hermitian_sparse_matrix(N, ρ, sparsity)
```

returns a matrix of size N×N with spectral radius ρ and the given sparsity
(number between zero and one that is the approximate fraction of non-zero
elements).
"""
function random_hermitian_sparse_matrix(N, ρ, sparsity)
    σ = 1/√(sparsity * N)
    d = Normal(0.0, σ)
    H = sprand(N, N, sparsity, (dims...) -> rand(d, dims...))
    return 0.5ρ * (H + H') / √2
end

"""Return a random, normalized Hilbert space state vector of dimension `N`.

```julia
random_state_vector(N)
```
"""
function random_state_vector(N)
    Ψ = rand(N) .* exp.((2π * im) .* rand(N))
    Ψ ./= norm(Ψ)
    return Ψ
end

end


text:
using PackageCompiler
create_sysimage([:Pluto, :PlutoUI,
                 :Plots, :LaTeXStrings,
                 :DataFrames, :Unitful,
                 :LinearAlgebra, :Markdown,
                 :InteractiveUtils];
                #precompile_execution_file = "warmup.jl",
                precompile_statements_file = "combined_trace.jl",
                #replace_default = true,
                sysimage_path="/home/jovyan/sysimage.so",
                cpu_target = PackageCompiler.default_app_cpu_target())


text:
@testset "High-level operations on distributed data" begin

    W = addprocs(2)
    @everywhere using GigaSOM

    Random.seed!(1)
    dd = rand(11111, 5)
    buckets = rand([1, 2, 3], 11111)

    di1 = distribute_array(:test1, dd, W[1:1])
    di2 = distribute_array(:test2, dd, W)
    buckets1 = distribute_array(:buckets1, buckets, W[1:1])
    buckets2 = distribute_array(:buckets2, buckets, W)

    @testset "Distribution works as expected" begin
        @test distributed_collect(di1) == dd
        @test distributed_collect(di2) == dd
    end

    subsel = [1, 3, 4]
    dd = dd[:, subsel]
    dselect(di1, subsel)
    dselect(di2, subsel)

    @testset "dapply_cols, dapply_rows" begin
        diV = distribute_array(:testV, dd, W[1:1])
        diH = distribute_array(:testH, Matrix(dd'), W[1:1])
        dapply_cols(diV, (x, _) -> x ./ sum(x), Vector(1:3))
        dapply_rows(diH, x -> x ./ sum(x))
        @test isapprox(distributed_collect(diV), Matrix(distributed_collect(diH)'))
        undistribute(diV)
        undistribute(diH)
    end

    @testset "dselect" begin
        @test distributed_collect(di1) == dd
        @test distributed_collect(di2) == dd
    end

    di3 = dcopy(di2, :test3)

    @testset "dcopy" begin
        @test di2.workers == di3.workers
        @test dd == distributed_collect(di3)
    end

    undistribute(di3)

    (means1, sds1) = dstat(di1, [1, 3])
    (means2, sds2) = dstat(di2, [1, 3])

    @testset "dstat" begin
        @test length(means1) == 2
        @test length(sds1) == 2
        @test isapprox(means1, means2, atol = 1e-8)
        @test isapprox(sds1, sds2, atol = 1e-8)
        @test all([isapprox(0.5, m, atol = 0.1) for m in means1])
        @test all([isapprox(0.29, s, atol = 0.1) for s in sds1])
    end

    medians1 = dmedian(di1, [1, 3])
    medians2 = dmedian(di2, [1, 3])

    @testset "dmedian" begin
        @test length(medians1) == 2
        @test isapprox(medians1, medians2, atol = 1e-4)
        @test all([isapprox(0.5, m, atol = 0.1) for m in medians1])
    end

    (bmeans1, bsds1) = dstat_buckets(di1, 3, buckets1, [1, 3])
    (bmeans2, bsds2) = dstat_buckets(di2, 3, buckets2, [1, 3])

    @testset "bucketed dstat" begin
        @test size(bmeans1) == (3, 2)
        @test size(bsds1) == (3, 2)
        @test isapprox(bmeans1, bmeans2, atol = 1e-8)
        @test isapprox(bsds1, bsds2, atol = 1e-8)
        @test all([isapprox(0.5, m, atol = 0.1) for m in bmeans1])
        @test all([isapprox(0.29, s, atol = 0.1) for s in bsds1])
    end

    medians1 = dmedian_buckets(di1, 3, buckets1, [1, 3])
    medians2 = dmedian_buckets(di2, 3, buckets2, [1, 3])

    @testset "dmedian" begin
        @test size(medians1) == (3, 2)
        @test isapprox(medians1, medians2, atol = 1e-4)
        @test all([isapprox(0.5, m, atol = 0.1) for m in medians1])
    end

    dscale(di1, [2, 3])
    dscale(di2, [2, 3])
    (means1, sds1) = dstat(di1, [1, 2, 3])
    (means2, sds2) = dstat(di2, [1, 2, 3])

    @testset "dscale" begin
        @test isapprox(means1, means2, atol = 1e-8)
        @test isapprox(sds1, sds2, atol = 1e-8)
        @test isapprox(0.5, means1[1], atol = 0.1)
        @test all([isapprox(0.0, m, atol = 1e-3) for m in means1[2:3]])
        @test all([isapprox(1.0, s, atol = 1e-3) for s in sds1[2:3]])
    end

    dc = distributed_collect(di1)
    dc[:, 1:2] = asinh.(dc[:, 1:2] ./ 1.23)
    dtransform_asinh(di1, [1, 2], 1.23)
    dtransform_asinh(di2, [1, 2], 1.23)

    @testset "dtransform_asinh" begin
        @test isapprox(distributed_collect(di1), dc)
        @test isapprox(distributed_collect(di2), dc)
    end

    undistribute(di1)
    undistribute(di2)
    rmprocs(W)

    d = ones(Float64, 2, 2)
    di = distribute_array(:test, d, [1])
    dscale(di, [1, 2])
    @testset "dscale does not produce NaNs on sdev==0" begin
        @test distributed_collect(di) == zeros(Float64, 2, 2)
    end
    undistribute(di)

end


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

function tag_list(repo::GitRepo)
    with(StrArrayStruct()) do sa
        sa_ref = Ref(sa)
        @check ccall((:git_tag_list, :libgit2), Cint,
                      (Ptr{StrArrayStruct}, Ptr{Void}), sa_ref, repo.ptr)
        convert(Vector{AbstractString}, sa_ref[])
    end
end

function tag_delete(repo::GitRepo, tag::AbstractString)
    @check ccall((:git_tag_delete, :libgit2), Cint,
                  (Ptr{Void}, Cstring, ), repo.ptr, tag)
end

function tag_create(repo::GitRepo, tag::AbstractString, commit::Union{AbstractString,Oid};
                    msg::AbstractString = "",
                    force::Bool = false,
                    sig::Signature = Signature(repo))
    oid_ptr = Ref(Oid())
    with(get(GitCommit, repo, commit)) do commit_obj
        commit_obj === nothing && return oid_ptr[] # return empty oid
        with(convert(GitSignature, sig)) do git_sig
            @check ccall((:git_tag_create, :libgit2), Cint,
                 (Ptr{Oid}, Ptr{Void}, Cstring, Ptr{Void}, Ptr{SignatureStruct}, Cstring, Cint),
                  oid_ptr, repo.ptr, tag, commit_obj.ptr, git_sig.ptr, msg, Cint(force))
        end
    end
    return oid_ptr[]
end

function name(tag::GitTag)
    str_ptr = ccall((:git_tag_name, :libgit2), Cstring, (Ptr{Void}, ), tag.ptr)
    str_ptr == C_NULL && throw(Error.GitError(Error.ERROR))
    return unsafe_string(str_ptr)
end

function target(tag::GitTag)
    oid_ptr = ccall((:git_tag_target_id, :libgit2), Ptr{Oid}, (Ptr{Void}, ), tag.ptr)
    oid_ptr == C_NULL && throw(Error.GitError(Error.ERROR))
    return Oid(oid_ptr)
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
"""
    OneSkeleton

A realization of the one-skeleton of a filtered simplicial complex at time t as an
`AbstractGraph`. Some edges of the graph may be removed by placing them into the `removed`
set. This will ignore said edges when looking for neighbours or shortest paths.
Weights are usually simplex birth times, but can be customized.
"""
struct OneSkeleton{T,F<:AbstractFiltration,S<:AbstractCell{1},A} <: AbstractGraph{Int}
    filtration::F
    threshold::T
    weights::A
    removed::Set{S}
end
function OneSkeleton(
    filtration::F,
    thresh::T=threshold(filtration),
    removed=(),
    weights=distance_matrix(filtration),
) where {F,T}
    S = simplex_type(filtration, 1)
    removed = Set{S}(removed)
    return OneSkeleton{T,F,S,typeof(weights)}(filtration, thresh, weights, removed)
end

_birth_or_value(σ::AbstractCell) = birth(σ)
_birth_or_value(σ) = σ
_in(σ::S, g::OneSkeleton{S}) where {S} = !isnothing(σ) && σ ≤ g.threshold && σ ∉ g.removed
_in(σ, g::OneSkeleton) = !isnothing(σ) && birth(σ) ≤ g.threshold && σ ∉ g.removed

# Convert linear indices to filtration indices. Applicable to Cubical because it has
# vertices of type CartesianIndex.
_linear(g, i) = LinearIndices(vertices(g.filtration))[i]
_inv_linear(g, i) = vertices(g.filtration)[i]

Graphs.edgetype(::OneSkeleton) = Edge{Int}

function Graphs.has_edge(g::OneSkeleton, u::Integer, v::Integer)
    if u ∉ vertices(g) || v ∉ vertices(g)
        return false
    else
        σ = simplex(g.filtration, Val(1), _inv_linear.(Ref(g), (v, u)))
        return _in(σ, g)
    end
end

function Graphs.edges(g::OneSkeleton)
    result = edgetype(g)[]
    for sx in Ripserer.edges(g.filtration)
        if _in(sx, g)
            u, v = _linear.(Ref(g), sx)
            push!(result, Edge(u, v))
        end
    end
    return result
end
function Graphs.outneighbors(g::OneSkeleton, u::Integer)
    root = simplex(g.filtration, Val(0), (_inv_linear(g, u),))
    neighbors = Int[]
    for sx in Ripserer.coboundary(g.filtration, root)
        if _in(sx, g)
            v, w = _linear.(Ref(g), sx)
            if v == u
                push!(neighbors, w)
            else
                push!(neighbors, v)
            end
        end
    end
    return neighbors
end
Graphs.inneighbors(g::OneSkeleton, u::Integer) = outneighbors(g, u)

Base.eltype(::OneSkeleton) = Int
Graphs.has_vertex(g::OneSkeleton, u) = 1 ≤ u ≤ nv(g)
Graphs.ne(g::OneSkeleton) = length(edges(g))
Graphs.nv(g::OneSkeleton) = nv(g.filtration)
Graphs.vertices(g::OneSkeleton) = Base.OneTo(nv(g))
Graphs.is_directed(::OneSkeleton) = false
Graphs.is_directed(::Type{<:OneSkeleton}) = false

# need default dist for filtrations...
Graphs.weights(g::OneSkeleton) = g.weights

_heuristic(filtration, src, distances) = dst -> distances[dst, src]
function _path_length(dists, cyc)
    result = zero(eltype(dists))
    for edge in cyc
        i, j = src(edge), dst(edge)
        result += dists[i, j]
    end
    return result
end

"""
    _find_cycle(g, dists)

Find the shortest cycle in `g` that has exactly one edge from `g.removed`.
"""
function _find_cycle(g, dists)
    # Idea:
    # best_weight is the current best length of shortest cycle.
    # best_path is the current best candidate for shortest cycle.
    # best_simplex completes best_path to a cycle.
    # We go through all edges in g.removed and find the shortest path through its
    # endpoints. The shortest among those is the shortest cycle.
    best_weight = missing
    best_path = edgetype(g)[]
    best_sx = first(g.removed)
    for sx in g.removed
        u, v = _linear.(Ref(g), sx)
        path = a_star(g, u, v, dists, _heuristic(g.filtration, v, dists))
        weight = _path_length(dists, path) + dists[u, v]
        if !isempty(path) && isless(weight, best_weight)
            best_weight = weight
            best_path = path
            best_sx = sx
        end
    end
    if ismissing(best_weight)
        error("no cycle found!")
    else
        result = [best_sx]
        # Convert the type of best_path from `Edge`s to the simplex type of the filtration.
        for e in best_path
            u, v = _inv_linear.(Ref(g), (src(e), dst(e)))
            push!(result, simplex(g.filtration, Val(1), (u, v)))
        end
        return result
    end
end

"""
    reconstruct_cycle(filtration, interval[, t]; distances=distance_matrix(filtration))

Reconstruct the shortest representative cycle for the first homology group of given
`interval`. The optional argument `t` sets the time at which the cycle is to be computed. It
defaults to interval birth time, which gives a cycle similar to a representative cycle
computed from homology. In general, higher times will yield nicer cycles. `t` can be a
simplex or a number.

The optional `distances` keyword argument can be used to change the distance matrix used for
determining edge lengths.

This method uses the representative _co_cycle to compute the cycle. As such, the interval
must be computed with the default cohomology algorithm and must include a representative. To
get such an interval, run `ripserer` with the keyword argument `reps=true` or `reps=1`.

Note that this method only works in the first dimension, as it is based on finding shortest
paths in a graph.

!!! warning
    This feature is still experimental.
"""
function reconstruct_cycle(
    filtration::AbstractFiltration{<:Any,T},
    interval,
    r=birth_simplex(interval);
    distances=distance_matrix(filtration),
) where {T}
    if !hasproperty(interval, :representative)
        throw(
            ArgumentError("interval has no representative! Run `ripserer` with `reps=true`")
        )
    elseif !(eltype(interval.representative) <: AbstractChainElement{<:AbstractCell{1}})
        throw(
            ArgumentError("cycles can only be reconstructed for 1-dimensional intervals.")
        )
    elseif !(birth(interval) ≤ _birth_or_value(r) < death(interval))
        return simplex_type(filtration, 1)[]
    else
        reps = filter!(simplex.(representative(interval))) do sx
            birth(sx) ≤ _birth_or_value(r)
        end
        g = OneSkeleton(filtration, r, reps)
        return _find_cycle(g, distances)
    end
end


text:
wl = filter(!isempty, split("""You will rejoice to hear that no disaster has accompanied the
    commencement of an enterprise which you have regarded with such evil
    forebodings.""", r"\W+"))

println("Original list:\n - ", join(wl, "\n - "))
sort!(wl; by=x -> (-length(x), lowercase(x)))
println("\nSorted list:\n - ", join(wl, "\n - "))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using BenchmarkTools
using unfold,StatsModels,MixedModels,DataFrames
using Random: randn,seed!

const SUITE = BenchmarkGroup()
seed!(3)

# Preparatory code
include("../test/test_utilities.jl"); # to load the simulated data

data, evts = loadtestdata("testcase6","test/")
categorical!(evts,:subject);
categorical!(evts,:stimulus);
evts.subjectB = evts.subject;
evts1 = evts[evts.condA.==0,:]
evts2 = evts[evts.condA.==1,:]
data = data.+ 0.5*randn(size(data)) # we have to add minimal noise, else mixed models crashes.

f1_lmm  = @formula 0~1+condB+(1|subject)
f2_lmm  = @formula 0~1+condB+(1|subjectB)

f1  = @formula 0~1+condB
f2  = @formula 0~1+condB

data_r = reshape(data,(1,:)) #make ch x time
data_epochs,times = unfold.epoch(data=data_r,tbl=evts,τ=(-0.4,0.8),sfreq=100); #epoch
evts,data_epochs = unfold.dropMissingEpochs(evts,data_epochs) # remove missing

ba1 = firbasis(τ=(0,1),sfreq = 10,name="evts1")
ba2 = firbasis(τ=(0,1),sfreq = 10,name="evts2")

X1_lmm  = designmatrix(UnfoldLinearMixedModel,f1_lmm,evts1,ba1)
X2_lmm  = designmatrix(UnfoldLinearMixedModel,f2_lmm,evts2,ba2)

X1   = designmatrix(UnfoldLinearModel,f1,evts1,ba1)
X2   = designmatrix(UnfoldLinearModel,f2,evts2,ba2)

SUITE = BenchmarkGroup()
SUITE["nodc"] = BenchmarkGroup(
    ["nodc"])
SUITE["dc"] = BenchmarkGroup(["dc"])
# designmatrix generation
SUITE["dc"]["X_gen_lin"] = @benchmarkable designmatrix(UnfoldLinearModel,$f1,$evts1,$ba1)
SUITE["dc"]["X_gen_mix"] = @benchmarkable designmatrix(UnfoldLinearMixedModel,$f1_lmm,$evts1,$ba1)

# designmat concatination
SUITE["dc"]["X_cat_lin"] = @benchmarkable $X1+$X2
SUITE["dc"]["X_cat_mix"] = @benchmarkable $X1_lmm+$X2_lmm

# Model Fit
SUITE["nodc"]["fit_lin"] = @benchmarkable unfold.fit(UnfoldLinearModel,$f1,$evts1,$data_epochs,$times)
#SUITE["nodc"]["fit_mix"] = @benchmarkable unfold.fit(UnfoldLinearMixedModel,$f1_lmm,$evts1,$data_epochs,$times)
SUITE["dc"]["fit_lin"] = @benchmarkable unfold.unfoldfit(UnfoldLinearModel,$X1+$X2,$data_r);
SUITE["dc"]["fit_mix"] = @benchmarkable unfold.unfoldfit(UnfoldLinearMixedModel,$X1_lmm+$X2_lmm,$data_r);


text:
include("cg_sem.jl")
using LinearAlgebra: eigen

function boundary_forcing_test(
    N,
    K,
    z0,
    z1;
    δw0 = t -> 0,
    δw1 = t -> 0,
    cfl = 1,
    tspan = (0, 1),
)
    elem = element_operators(N)
    mesh = create_mesh(K, elem, z0, z1)

    problem = Isothermal{Float64}()

    # Matrices for upwinding the boundary conditions
    A0 = A_(problem, z0)
    (λ, V) = eigen(A0)
    A0⁺ = V * Diagonal(max.(0, λ)) / V
    A0⁻ = V * Diagonal(min.(0, λ)) / V

    A1 = -A_(problem, z1)
    (λ, V) = eigen(A1)
    A1⁺ = V * Diagonal(max.(0, λ)) / V
    A1⁻ = V * Diagonal(min.(0, λ)) / V

    bc = (
        (A = A0, A⁺ = A0⁺, A⁻ = A0⁻, δw = δw0),
        (A = A1, A⁺ = A1⁺, A⁻ = A1⁻, δw = δw1),
    )

    # Calculate the wave speed
    c̄ = maximum(c̄_(problem, mesh.zdg))

    # estimate time step
    dt0 = cfl * mesh.Δz / (N^2 * c̄)

    # Initial conditions
    q = (
        δρ = fill!(similar(mesh.zdg), 0),
        δw = fill!(similar(mesh.zdg), 0),
        δE = fill!(similar(mesh.zdg), 0),
    )

    # advance in time
    rhs!(∂q, q, t) = tendency!(∂q, q, t, problem, mesh, elem, bc)

    for step in 1:length(tspan)
        if step > 1
            ts = (tspan[step - 1], tspan[step])
            steps = ceil(Int, (ts[2] - ts[1]) / dt0)
            dt = (ts[2] - ts[1]) / steps
            timestep!(q, rhs!, dt, ts)
        end

        # Calculate the L2 error
        L2 = zeros(3)
        for i in 1:3
            L2[i] = sqrt(sum(mesh.J .* elem.ω .* q[i] .^ 2))
        end
        t = tspan[step]
        @show t, L2
    end
end