text:
using Plots
using Plots.PlotMeasures
using Random
using LaTeXStrings
using StatsPlots
include("shuffle_config.jl")
include("initialize.jl")
include("sro_calc.jl")

"""
plot_configs(my_sros, filename="my_config_plot.png")

"""
function plot_configs(my_sros, filename="my_config_plot.png")
    #ticks = ["aaaaaaaabbbbbbbb", "bbbbbbbbaaaaaaaa"]
    x_axis = [i for i = 1:1:length(my_sros)]
    my_plt = plot(title = "SRO values",
            xlabel = "Configuration",
            ylabel = "Objective function",
            xticks = (1:length(my_sros)-1:length(my_sros)),#, ticks),
            tickfontsize = 20,
            legend = false,
            palette = :rainbow,
            bottom_margin = 10mm,
            left_margin = 10mm,
            right_margin = 35mm,
            size = (2000, 1300),
            guidefontsize = 25,
            titlefontsize = 35,
            )

    plot!(my_plt, x_axis, my_sros)
    display(my_plt)
end

"""
xyzplot(Backpack, best_sequence)

Just a little script for displaying the minimum structure. If you do not like my
chosen color palette, feel free to change it ;-)
"""
function xyzplot(X::Backpack, seq, cam_angle::Float64)
    plt = scatter(title = "Supercell", xlabel = "X in Å", ylabel = "Y in Å", gridlinewidth = 1, zlabel = "Z in Å", camera = (cam_angle, 15))
    my_colors = palette(:Accent_8) #Set3_12
    k = 0
    for (idx, val) in X.element_dict
        xyz = zeros(length(seq), 3)
        j = 0
        k += 1
        for i = 1:length(seq)
            if seq[i] == val
                j += 1
                xyz[i, :] =  X.supercell.lattice*X.supercell.fcoords[i, :]
            end
        end
        scatter!(plt, xyz[:, 1], xyz[:, 2], xyz[:, 3], markersize = 11, markercolor = my_colors[k],  label=idx)
    end
    return(plt)
end

"""
xyz_animation(X::Backpack, seq, gifname="supercell.gif", fps=15)

Animates a ittle gif of the supercell for better viewability
"""
function xyz_animation(X::Backpack, seq, gifname="supercell.gif", fps=15)
    n = 100
    t1 = [i*90/n for i=1:n]
    t = vcat(t1, reverse(t1))
    anim = @animate for i=1:length(t)
        xyzplot(X, seq, t[i])
    end
    return gif(anim, gifname, fps = 15)
end

max_stars_repo_name: ulipoto/WCSRO

max_issues_repo_name: ulipoto/WCSRO

max_forks_repo_name: ulipoto/WCSRO


text:
@testset "DisplayRow" begin

    row = DS.DisplayRow(40, background='.')

    @test repr(row) == "DisplayRow(" *
        "size=40, " *
        "background char=Char(46)" *
    ")"

    @test lastindex(row) == length(row)
    @test textwidth(Base.join(row.content)) == length(row)
    DS.render(row)
    @test T.read_stream(T.out_stream)  ==
        "........................................"

    row[10] = '字'
    DS.render(row)
    @test T.read_stream(T.out_stream)  ==
        ".........字............................."

    @test DS.get_element_index(row, 10) == (i=10, pre=0, post=1)

    row[3:7] = "兩bit"
    DS.render(row)
    @test T.read_stream(T.out_stream)  ==
        "..兩bit..字............................."

    @test DS.get_element_index(row, 4) == (i=3, pre=1, post=0)
    @test textwidth(Base.join(row.content)) == length(row)

    row[16:20] = "很棒"
    DS.render(row, pos=(5, 6))
    @test T.read_stream(T.out_stream)  ==
        "\e[5;6H" *
        "..兩bit..字....很棒....................."

    @test DS.get_element_index(row, 17) == (i=14, pre=1, post=0)
    @test textwidth(Base.join(row.content)) == length(row)

    deleteat!(row, 5) # 'b'
    @test DS.get_element_index(row, 30) == (i=26, pre=0, post=0)
    @test textwidth(Base.join(row.content)) == length(row)
    DS.render(row)
    @test T.read_stream(T.out_stream)  ==
        "..兩.it..字....很棒....................."

    deleteat!(row, 4) # '兩'
    @test DS.get_element_index(row, 30) == (i=27, pre=0, post=0)
    @test textwidth(Base.join(row.content)) == length(row)
    DS.render(row)
    @test T.read_stream(T.out_stream)  ==
        ".....it..字....很棒....................."

    # pos
    DS.render(row, pos=(5, 6))
    @test T.read_stream(T.out_stream)  ==
        "\e[5;6H" *
        ".....it..字....很棒....................."

end

ext: jl


text:
using CxxWrap
using Test

@testset "$(basename(@__FILE__)[1:end-3])" begin

let s = StdString("test")
  println("This prints a test string: ", s)
  @test s == "test"
  stref = CxxRef(s)
  @test stref == s
  @test stref == stref
end

let s = "šČô_φ_привет_일보"
  @show StdWString(s)
  @test StdWString(s) == s
  @test String(StdWString(s)) == s
end

let s = "😄😈😼"
  @show StdWString(s)
  @test StdWString(s) == s
  @test String(StdWString(s)) == s
end

let s = "café"
  @show StdString(s)
  @test StdString(s) == s
  @test String(StdString(s)) == s
end

let s = StdString("foo")
  @test String(s) == "foo"
  sref = CxxRef(s)
  @test sref == "foo"
  @test String(sref) == "foo"
  @test unsafe_string(CxxWrap.StdLib.c_str(s)) == "foo"
  @test unsafe_string(CxxWrap.StdLib.c_str(s),2) == "fo"
end

stvec = StdVector(Int32[1,2,3])
@test all(stvec .== [1,2,3])
push!(stvec,1)
@test all(stvec .== [1,2,3,1])
resize!(stvec,2)
@test all(stvec .== [1,2])
append!(stvec,Int32[2,1])
@test all(stvec .== [1,2,2,1])
empty!(stvec)
@test isempty(stvec)

@test all(StdVector([1,2,3]) .== [1,2,3])
@test all(StdVector([1.0,2.0,3.0]) .== [1,2,3])
@test all(StdVector([true, false, true]) .== [true, false, true])
bvec = StdVector([true, false, true])
append!(bvec, [true])
@test all(bvec .== [true, false, true, true])

cxxstrings = StdString["one", "two", "three"]
svec = StdVector(CxxRef.(cxxstrings))
@test all(svec .== ["one", "two", "three"])
push!(svec, StdString("four"))
@test all(svec .== ["one", "two", "three", "four"])
cxxappstrings = StdString["five", "six"]
append!(svec, CxxRef.(cxxappstrings))
@test all(svec .== ["one", "two", "three", "four", "five", "six"])
empty!(svec)
@test isempty(svec)

stvec = StdVector(Int32[1,2,3])
for vref in (CxxRef(stvec), CxxPtr(stvec))
  s = 0
  for x in vref
    s += x
  end
  @test s == sum(stvec)
  @test all(vref .== [1,2,3])
  @test vref[1] == 1
  vref[3] = 10
  @test vref[3] == 10
  vref[3] = 3
end

let
  valarr1 = StdValArray{Float64}()
  @test length(valarr1) == 0
  valarr2 = StdValArray([1.0, 2.0, 3.0])
  @test valarr2 == [1.0, 2.0, 3.0]
  valarr2[2] = sum(valarr2)
  @show valarr2
  @test valarr2[2] == 6
end

end


text:
"Default inference algorithm"
defalg(args...) = FailUnsat

"Default callbacks"
defcb(args...) = donothing

"Single sample form from `x`"
function Base.rand(x::RandVar;
                   alg::SamplingAlgorithm = defalg(x),
                   ΩT = defΩ(alg),
                   kwargs...)
  first(rand(Random.GLOBAL_RNG, x, 1, alg; ΩT = ΩT, kwargs...))
end

function Base.rand(rng::AbstractRNG,
                   x::RandVar;
                   alg::SamplingAlgorithm = defalg(x),
                   ΩT::Type{OT} = defΩ(alg),
                   kwargs...) where {OT <: Ω}
  first(rand(rng, x, 1, alg; ΩT = ΩT, kwargs...))
end

function Base.rand(rng::AbstractRNG,
                   x::RandVar,
                   n::Integer,
                   alg::SamplingAlgorithm;
                   ΩT::Type{OT} = defΩ(alg),
                   kwargs...) where {OT <: Ω}
  logdensity = Omega.mem(logerr(indomainₛ(x)))
  ωsamples = rand(rng, ΩT, logdensity, n, alg; kwargs...)
  map(ω -> applynotrackerr(Omega.mem(x), ω), ωsamples)
end

function Base.rand(x::RandVar,
                   n::Integer;
                   alg::SamplingAlgorithm = defalg(x),
                   ΩT::Type{OT} = defΩ(alg),
                   kwargs...) where {OT <: Ω}
  rand(Random.GLOBAL_RNG, x, n, alg; ΩT = ΩT, kwargs...)
end

function Base.rand(rng::AbstractRNG,
                   x::RandVar,
                   n::Integer;
                   alg::SamplingAlgorithm = defalg(x),
                   ΩT::Type{OT} = defΩ(alg),
                   kwargs...) where {OT <: Ω}
  rand(rng, x, n, alg; ΩT = ΩT, kwargs...)
end

function Base.rand(x::RandVar,
                   n::Integer,
                   alg::SamplingAlgorithm;
                   ΩT::Type{OT} = defΩ(alg),
                   kwargs...) where {OT <: Ω}
  rand(Random.GLOBAL_RNG, x, n, alg; ΩT = ΩT, kwargs...)
end

# Condition (deprecate this?) x | y

Base.rand(x::RandVar, y::RandVar, n; kwargs...) = rand(Random.GLOBAL_RNG, cond(x, y), n; kwargs...)
Base.rand(x::RandVar, y::RandVar; kwargs...) = rand(Random.GLOBAL_RNG, cond(x, y); kwargs...)
Base.rand(rng::AbstractRNG, x::RandVar, y::RandVar, n; kwargs...) = rand(rng, cond(x, y), n; kwargs...)
Base.rand(rng::AbstractRNG, x::RandVar, y::RandVar; kwargs...) = rand(rng, cond(x, y); kwargs...)

# Return Tuple

Base.rand(x::UTuple{RandVar}, n::Integer; kwargs...) = rand(randtuple(x), n; kwargs...)
Base.rand(x::UTuple{RandVar}; kwargs...) = rand(randtuple(x); kwargs...)

Base.rand(x::UTuple{RandVar}, y::RandVar, n::Integer; kwargs...) = rand(randtuple(x), y, n; kwargs...)
Base.rand(x::UTuple{RandVar}, y::RandVar; kwargs...) = rand(randtuple(x), y; kwargs...)

# Return Ω object

"Returns Ω object, .e.g `x = normal(0, 1); rand(Ω, x > 0)`"
function ld(rng::AbstractRNG,
            ΩT::Type{OT},
            x::RandVar,
            n::Integer,
            alg::IsSoft{ALG};
            kwargs...) where {ALG, OT <: Ω}
  logdensity = Omega.mem(logerr(indomainₛ(x)))
  ωsamples = rand(rng, ΩT, logdensity, n, ALG(); kwargs...)
end

function ld(rng::AbstractRNG,
            ΩT::Type{OT},
            x::RandVar,
            n::Integer,
            alg::IsHard{ALG};
            kwargs...) where {ALG, OT <: Ω}
  pred = Omega.indomain(x)
  ωsamples = rand(rng, ΩT, pred, n, ALG())
end

function ld(rng::AbstractRNG,
            ΩT::Type{OT},
            x::RandVar,
            n::Integer,
            alg::T;
            kwargs...) where {T, OT <: Ω}
  ld(rng, ΩT, x, n, softhard(T); kwargs...)
end

"Returns Ω object, e.g. `x = normal(0, 1); rand(Ω, cond(x, x >ₛ 2.0), 10)`"
Base.rand(rng::AbstractRNG, ::Type{Ω}, x::RandVar, n::Integer; alg::SamplingAlgorithm = defalg(x), kwargs...) =
  ld(rng, x, defΩ(alg), n, alg; kwargs...)

Base.rand(::Type{Ω}, x::RandVar, n::Integer; alg::SamplingAlgorithm = defalg(x), kwargs...) =
  ld(Random.GLOBAL_RNG, defΩ(alg), x, n, alg; kwargs...)


text:
import JuLIP
using JuLIP: JVec, Atoms
using JuLIP.Potentials: @pot, MPairPotential, SZList, ZList, z2i, i2z
using LinearAlgebra: dot

import JuLIP.Potentials: evaluate, evaluate_d

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

export PolyPairPot

struct PolyPairPot{T,TJ,NZ} <: MPairPotential
   coeffs::Vector{T}
   J::TJ
   zlist::SZList{NZ}
   bidx0::SMatrix{NZ,NZ,Int16}
end

@pot PolyPairPot

PolyPairPot(pB::PolyPairBasis, coeffs::Vector) =
            PolyPairPot(coeffs, pB.J, pB.zlist, pB.bidx0)

JuLIP.MLIPs.combine(pB::PolyPairBasis, coeffs::AbstractVector) =
            PolyPairPot(pB, collect(coeffs))

JuLIP.cutoff(V::PolyPairPot) = cutoff(V.J)

==(V1::PolyPairPot, V2::PolyPairPot) =
            ( (V1.J == V2.J) && (V1.coeffs == V2.coeffs) &&  V1.zlist == V2.zlist )

Dict(V::PolyPairPot) = Dict(
      "__id__" => "PolyPairPots_PolyPairPot",
      "coeffs" => V.coeffs,
      "J" => Dict(V.J),
      "zlist" => V.zlist.list
      )

function PolyPairPot(D::Dict)
   J = TransformedJacobi(D["J"])
   zlist = ZList(Int16.(D["zlist"]), static = true)
   return PolyPairPot( Vector{Float64}(D["coeffs"]),
                       J, zlist, get_bidx0(J, zlist) )
end

convert(::Val{:PolyPairPots_PolyPairPot}, D::Dict) = PolyPairPot(D)

alloc_temp(V::PolyPairPot{T}, N::Integer) where {T} =
      ( J = alloc_B(V.J),
        R = zeros(JVec{T}, N),
        Z = zeros(Int16, N) )

alloc_temp_d(V::PolyPairPot{T}, N::Integer) where {T} =
      ( J = alloc_B(V.J),
       dJ = alloc_dB(V.J),
       dV = zeros(JVec{T}, N),
        R = zeros(JVec{T}, N),
        Z = zeros(Int16, N) )

function _dot_zij(V, B, z, z0)
   i0 = _Bidx0(V, z, z0)
   return sum( V.coeffs[i0 + n] * B[n]  for n = 1:length(V.J) )
end

evaluate!(tmp, V::PolyPairPot, r::Number, z, z0) =
      _dot_zij(V, evaluate!(tmp.J, nothing, V.J, r), z, z0)

evaluate_d!(tmp, V::PolyPairPot, r::Number, z, z0) =
      _dot_zij(V, evaluate_d!(tmp.J, tmp.dJ, nothing, V.J, r), z, z0)

function evaluate!(tmp, V::PolyPairPot, r::Number)
   @assert length(V.zlist) == 1
   z = V.zlist.list[1]
   return evaluate!(tmp, V::PolyPairPot, r::Number, z, z)
end

function evaluate_d!(tmp, V::PolyPairPot, r::Number)
   @assert length(V.zlist) == 1
   z = V.zlist.list[1]
   return evaluate_d!(tmp, V::PolyPairPot, r::Number, z, z)
end

evaluate(V::PolyPairPot, r::Number) = evaluate!(alloc_temp(V, 1), V, r)
evaluate_d(V::PolyPairPot, r::Number) = evaluate_d!(alloc_temp_d(V, 1), V, r)


text:
using Test, TuringMachine, DataStructures

# test load_input
input_file = "test_input_1.txt"
input_1 = collect("111000011")
input = load_input(input_file)
@test typeof(load_input(input_file)) == Array{Char,1}
@test input == input_1

# test load_program
program_file = "test_program_1.txt"
program, init_state, halt_state = load_program(program_file)
@test typeof(program) == Dict{Any,Any}
@test init_state == "q0"
@test halt_state == "qAccept"
@test in(("q2", "_") => ("qReject", "_", "-"), program)
@test in(("q0", "_") => ("qAccept", "_", "-"), program)
@test in(("q2", "1") => ("q2", "1", ">"), program)

# test set_up
state, tape_left, tape_right = set_up(init_state, input)
@test state == "q0"
@test isempty(tape_left)
@test length(tape_right) == length(input_1)
@test front(tape_right) == parse(Int, input_1[1])
@test back(tape_right) == parse(Int, input_1[end])

# test simulate
state, tape_left, tape_right = simulate(state, program, tape_left, tape_right)
@test state == "q1"
@test length(tape_left) == 1
@test front(tape_left) == parse(Int, input_1[1])
@test length(tape_right) == length(input_1)-1
@test front(tape_right) == parse(Int, input_1[2])
@test back(tape_right) == parse(Int, input_1[end])
while true
        global state, tape_left, tape_right = simulate(state, program, tape_left, tape_right)
        state ==("qAccept") && break
        state ==("qReject") && break
end
@test state == "qReject"
@test isempty(tape_right)
@test length(tape_left) == length(input_1)
@test front(tape_left) == parse(Int, input_1[1])
@test back(tape_left) == parse(Int, input_1[end])


text:
# SPDX-License-Identifier: MIT

"""
    struct VectorFSM{T} <: AbstractMutableFSM{T}
        states::Vector{State{T}}
        arcs::Dict{State, Vector{Arc{T}}}
    end

Mutable FSM that store the states and arcs into vectors.
"""
struct VectorFSM{T} <: AbstractMutableFSM{T}
    states::Vector{State{T}}
    arcs::Dict{State, Vector{Arc{T}}}
end

VectorFSM{T}() where T = VectorFSM{T}([], Dict())

states(fsm::VectorFSM) = fsm.states
arcs(fsm::VectorFSM{T}, state) where T = get(fsm.arcs, state, Arc{T}[])

function addstate!(fsm::VectorFSM{T}, label; initweight = zero(T),
                   finalweight = zero(T)) where T
    s = State(length(fsm.states)+1, label, initweight, finalweight)
    push!(fsm.states, s)
    s
end

function addarc!(fsm::VectorFSM{T}, src, dest, weight::T = one(T)) where T
    arclist = get(fsm.arcs, src, Arc{T}[])
    arc = Arc{T}(dest, weight)
    push!(arclist, arc)
    fsm.arcs[src] = arclist
    arc
end

max_forks_repo_name: BUTSpeechFIT/ProbabilisticFSMs


text:
module PkgDeps

using PlotRecipes

# const _pkgs = Pkg.available()
# const _idxmap = Dict(p=>i for (i,p) in enumerate(_pkgs))
# const _alist = [Int[] for i=1:length(_pkgs)]

# for pkg in _pkgs
# 	i = _idxmap[pkg]
# 	for dep in Pkg.dependents(pkg)
# 		if !haskey(_idxmap, dep)
# 			push!(_pkgs, dep)
# 			push!(_alist, [])
# 			_idxmap[dep] = length(_pkgs)
# 		end
# 		j = _idxmap[dep]
# 		push!(_alist[j], i)
# 	end
# end

@userplot DepsGraph
@recipe function f(g::DepsGraph)
	source, destiny, names = g.args
    arrow --> arrow()
    markersize --> 50
    markeralpha --> 0.2
    linealpha --> 0.4
    linewidth --> 2
    names --> names
    func --> :tree
    root --> :left
    PlotRecipes.GraphPlot((source, destiny))
end

# const args = (source, destiny, pkgs)

const all_pkgs = Pkg.available()
const deplists = Dict(pkg => Pkg.dependents(pkg) for pkg in all_pkgs)

const childlists = Dict(pkg => Set{String}() for pkg in all_pkgs)
for pkg in all_pkgs
	for dep in deplists[pkg]
		if haskey(childlists, dep)
			push!(childlists[dep], pkg)
		else
			warn("Package $dep wasn't in Pkg.available()")
			deplists[dep] = []
			childlists[dep] = Set([pkg])
		end
	end
end

function add_deps(pkg, deps = Set([pkg]))
	for dep in deplists[pkg]
		if !(dep in deps)
			push!(deps, dep)
			add_deps(dep, deps)
		end
	end
	deps
end

function add_children(pkg, children = Set([pkg]))
	for child in childlists[pkg]
		if !(child in children)
			push!(children, child)
			add_children(child, children)
		end
	end
	children
end

function plotdeps(pkg)
	pkgs = unique(union(add_deps(pkg), add_children(pkg)))
	idxmap = Dict(p=>i for (i,p) in enumerate(pkgs))

	source, destiny = Int[], Int[]
	for pkg in pkgs
		i = idxmap[pkg]
		for dep in deplists[pkg]
			# if !haskey(_idxmap, dep)
			# 	push!(pkgs, dep)
			# 	push!(_alist, [])
			# 	_idxmap[dep] = length(pkgs)
			# end
			j = idxmap[dep]
			push!(source, j)
			push!(destiny, i)
			# push!(_alist[j], i)
		end
	end
	depsgraph(source, destiny, pkgs, root=:bottom)
end

	# # pkgs = Set([pkg])
	# idx = _idxmap[pkg]
	# source, destiny = Int[], Int[]
	# for j in _alist[i]
	# 	push!(pkgs, _pkgs[j])
	# 	push!(source, j)
	# 	push!(destiny, i)
	# end

# to use:
# depsgraph(PkgDeps.args...)

end # module


text:
import Voronoi

#using Random
include("Draw.jl")

function randf(start, finish, n)
   v = rand(n)
   return map(x -> start + x*(finish-start), v)
end

function demoLloyd()
   WIDTH = 100.0
   HEIGHT = 100.0
   #n = 10

   #seed = 100
   #println("Random seed: ", seed)
   #Random.seed!(seed)

   command = nothing
   for n in 2:100
      points = convert(Array{Tuple{Real, Real}}, collect(zip(randf(1, WIDTH-1, n), randf(1, HEIGHT-1, n))))
      #n = 55
      #points = convert(Array{Tuple{Real, Real}}, [(90.0, 90.0),# (80.0, 50), (70.0, 55.0), (30.0, 85.0), (30, 45), (35.0, 44.0),
      #                                            #(1, 90),
      #                                            #(5, 90), (10, 90), (15, 90),
      #                                            (10, 5), (10,10), (10, 20), (80, 20), (30, 20),
      #                                            (1, 1.1), (2, 1.2), (3, 1.3),
      #                                           ])
      #points = convert(Array{Tuple{Real, Real}}, [ (100*0.2219056800797507, 100*0.7000658172436546), (100*0.49169960385155936, 100*0.6636482987025588), (100*0.7618642131758435, 100*0.7432736698473899), (100*0.8846773855975993, 100*0.6317667493965284) ])

      max = 20*n
      if max > 100
         max = 100
      end
      for i in 1:max
         # test for repeated points
         repeat = false
         for i in 1:length(points)
            for j in 1:length(points)
               if j == i
                  continue
               end

               if points[i][1] == points[j][1] && points[i][2] == points[j][2]
                  repeat = true
               end
            end
         end
         if repeat
            println("REPEATED POINTS!")
            continue
         end

         V = Voronoi.Fortune.compute(points, WIDTH, HEIGHT)

         Voronoi.Intersect.intersect(V, Voronoi.Intersect.Rectangle(WIDTH, HEIGHT))

         gradient = Voronoi.Optimization.∇f(V, points)
         println(Voronoi.Optimization.f(V), "\t\t\t", Voronoi.Optimization.norm2(gradient))

         Draw.init(WIDTH, HEIGHT)
         Draw.voronoiDiagram(V)

         points, areas = Voronoi.Diagram.centroidsAndAreas(V)

         if command != "a"
            println("Press Return for a step or enter \"a\" to animate until the end.")
            command = readline(stdin)
         else
            sleep(0.0000001) # make sure stuff is drawn
         end
         n += 1
      end
   end
end

demoLloyd()

max_stars_repo_name: arthursantana/coveringWithCircles

max_issues_repo_name: arthursantana/coveringWithCircles

max_forks_repo_name: arthursantana/coveringWithCircles


text:
# Tests the routine kickfast jacobian.  This routine
# computes the impulse gradient after Dehnen & Hernandez (2017).

# Next, try computing three-body Keplerian Jacobian:

@testset "kickfast" begin
n = 3
t0 = 7257.93115525
h  = 0.05
tmax = 600.0
dlnq = big(1e-15)

elements = readdlm("elements.txt",',')
elements[2,1] = 1.0
elements[3,1] = 1.0

m =zeros(n)
x0=zeros(3,n)
v0=zeros(3,n)
# Define which pairs will have impulse rather than -drift+Kepler:
pair = ones(Bool,n,n)  # This does impulses
# We want Keplerian between star & planets, and impulses between
# planets.  Impulse is indicated with 'true', -drift+Kepler with 'false':
for i=2:n
  pair[1,i] = false  # This does a Kepler + drift
  # We don't need to define this, but let's anyways:
  pair[i,1] = false
end

jac_step = zeros(7*n,7*n)

for k=1:n
  m[k] = elements[k,1]
end
m0 = copy(m)

x0,v0 = init_nbody(elements,t0,n)

# Tilt the orbits a bit:
x0[2,1] = 5e-1*sqrt(x0[1,1]^2+x0[3,1]^2)
x0[2,2] = -5e-1*sqrt(x0[1,2]^2+x0[3,2]^2)
x0[2,3] = -5e-1*sqrt(x0[1,2]^2+x0[3,2]^2)
v0[2,1] = 5e-1*sqrt(v0[1,1]^2+v0[3,1]^2)
v0[2,2] = -5e-1*sqrt(v0[1,2]^2+v0[3,2]^2)
v0[2,3] = -5e-1*sqrt(v0[1,2]^2+v0[3,2]^2)

# Take a step:
dh17!(x0,v0,h,m,n,pair)

# Now, copy these to compute Jacobian (so that I don't step
# x0 & v0 forward in time):
x = copy(x0); v = copy(v0); m = copy(m0)
xerror = zeros(3,n); verror = zeros(3,n)
# Compute jacobian exactly:
dqdt_kick = zeros(7*n)
kickfast!(x,v,xerror,verror,h,m,n,jac_step,dqdt_kick,pair)
# Add in identity matrix:
for i=1:7*n
  jac_step[i,i] += 1
end

# Now compute numerical derivatives, using BigFloat to avoid
# round-off errors:
jac_step_num = zeros(BigFloat,7*n,7*n)
# Save these so that I can compute derivatives numerically:
xsave = big.(x0)
vsave = big.(v0)
msave = big.(m0)
hbig = big(h)
xerror = zeros(3,n); verror = zeros(3,n)
# Carry out step using BigFloat for extra precision:
kickfast!(xsave,vsave,hbig,msave,n,pair)
xbig = big.(x0)
vbig = big.(v0)
mbig = big.(m0)
# Compute numerical derivatives wrt time:
dqdt_num = zeros(BigFloat,7*n)
# Vary time:
kickfast!(xbig,vbig,hbig,mbig,n,pair)
# Initial positions, velocities & masses:
xbig .= big.(x0)
vbig .= big.(v0)
mbig .= big.(m0)
hbig = big(h)
dq = dlnq * hbig
hbig += dq
kickfast!(xbig,vbig,hbig,mbig,n,pair)
# Now x & v are final positions & velocities after time step
for i=1:n, k=1:3
  dqdt_num[(i-1)*7+  k] = (xbig[k,i]-xsave[k,i])/dq
  dqdt_num[(i-1)*7+3+k] = (vbig[k,i]-vsave[k,i])/dq
end
hbig = big(h)
# Vary the initial parameters of planet j:
for j=1:n
  # Vary the initial phase-space elements:
  for jj=1:3
  # Initial positions, velocities & masses:
    xbig .= big.(x0)
    vbig .= big.(v0)
    mbig .= big.(m0)
    dq = dlnq * xbig[jj,j]
    if xbig[jj,j] != 0.0
      xbig[jj,j] +=  dq
    else
      dq = dlnq
      xbig[jj,j] = dq
    end
    kickfast!(xbig,vbig,hbig,mbig,n,pair)
  # Now x & v are final positions & velocities after time step
    for i=1:n
      for k=1:3
        jac_step_num[(i-1)*7+  k,(j-1)*7+jj] = (xbig[k,i]-xsave[k,i])/dq
        jac_step_num[(i-1)*7+3+k,(j-1)*7+jj] = (vbig[k,i]-vsave[k,i])/dq
      end
    end
    xbig .= big.(x0)
    vbig .= big.(v0)
    mbig .= big.(m0)
    dq = dlnq * vbig[jj,j]
    if vbig[jj,j] != 0.0
      vbig[jj,j] +=  dq
    else
      dq = dlnq
      vbig[jj,j] = dq
    end
    kickfast!(xbig,vbig,hbig,mbig,n,pair)
    for i=1:n
      for k=1:3
        jac_step_num[(i-1)*7+  k,(j-1)*7+3+jj] = (xbig[k,i]-xsave[k,i])/dq
        jac_step_num[(i-1)*7+3+k,(j-1)*7+3+jj] = (vbig[k,i]-vsave[k,i])/dq
      end
    end
  end
# Now vary mass of planet:
  xbig .= big.(x0)
  vbig .= big.(v0)
  mbig .= big.(m0)
  dq = mbig[j]*dlnq
  mbig[j] += dq
  kickfast!(xbig,vbig,hbig,mbig,n,pair)
  for i=1:n
    for k=1:3
      jac_step_num[(i-1)*7+  k,j*7] = (xbig[k,i]-xsave[k,i])/dq
      jac_step_num[(i-1)*7+3+k,j*7] = (vbig[k,i]-vsave[k,i])/dq
    end
    # Mass unchanged -> identity
    jac_step_num[7*i,7*i] = big(1.0)
  end
end
jac_step_num = convert(Array{Float64,2},jac_step_num)

jacmax = 0.0
for i=1:7, j=1:3, k=1:7, l=1:3
  if jac_step[(j-1)*7+i,(l-1)*7+k] != 0
    # Compute the fractional error and absolute error, and take the minimum of the two:
    diff = minimum([abs(jac_step_num[(j-1)*7+i,(l-1)*7+k]/jac_step[(j-1)*7+i,(l-1)*7+k]-1.0);jac_step_num[(j-1)*7+i,(l-1)*7+k]-jac_step[(j-1)*7+i,(l-1)*7+k]])
    if diff > jacmax
      jacmax = diff
    end
  end
end
#println("jac_step: ",jac_step," jac_step-jac_step_num: ",jac_step-jac_step_num)

println("Maximum jac_step kickfast error: ",jacmax)
dqdt_num = convert(Array{Float64,1},dqdt_num)
println("Maximum dqdt_kick kickdfast error: ",maximum(abs.(dqdt_kick-dqdt_num)))

@test isapprox(jac_step,jac_step_num;norm=maxabs)
@test isapprox(dqdt_kick,dqdt_num;norm=maxabs)
end


text:
export
    TestMode, TrainMode,
    FluxOptApp, OptimOptApp, SciMLOptApp

abstract type Flows end
abstract type NormalizingFlows <: Flows end
abstract type ContinuousNormalizingFlows <: NormalizingFlows end
abstract type InfinitesimalContinuousNormalizingFlows <: ContinuousNormalizingFlows end

abstract type AbstractICNF{T} <: InfinitesimalContinuousNormalizingFlows where {T <: AbstractFloat} end
abstract type AbstractCondICNF{T} <: InfinitesimalContinuousNormalizingFlows where {T <: AbstractFloat} end

abstract type Mode end
struct TestMode <: Mode end
struct TrainMode <: Mode end

abstract type OptApp end
struct FluxOptApp <: OptApp end
struct OptimOptApp <: OptApp end
struct SciMLOptApp <: OptApp end

# MLJ interface

abstract type MLJICNF <: MLJModelInterface.Unsupervised end

# Distributions interface

abstract type ICNFDistribution <: ContinuousMultivariateDistribution end


text:
########################
###### Data Types ######
########################

# (1 + 2)::Int
# 3

function f(y::Int)
    x::Int8 = 100
    x + y
end

x = f(20) # 120

typeof(x) # Int8

############################
###### Abstract Types ######
############################

# abstract type «name» end
# abstract type «name» <: «supertype» end

abstract type Number end
abstract type Real     <: Number end
abstract type AbstractFloat <: Real end
abstract type Integer  <: Real end
abstract type Signed   <: Integer end
abstract type Unsigned <: Integer end

# The Number type is a direct child type of Any, and Real is its child.

#############################
###### Primitive Types ######
#############################

# primitive type «name» «bits» end
# primitive type «name» <: «supertype» «bits» end

# primitive type Float16 <: AbstractFloat 16 end
# primitive type Float32 <: AbstractFloat 32 end
# primitive type Float64 <: AbstractFloat 64 end

# primitive type Bool <: Integer 8 end
# primitive type Char <: AbstractChar 32 end

# primitive type Int8    <: Signed   8 end
# primitive type UInt8   <: Unsigned 8 end
# primitive type Int16   <: Signed   16 end
# primitive type UInt16  <: Unsigned 16 end
# primitive type Int32   <: Signed   32 end
# primitive type UInt32  <: Unsigned 32 end
# primitive type Int64   <: Signed   64 end
# primitive type UInt64  <: Unsigned 64 end
# primitive type Int128  <: Signed   128 end
# primitive type UInt128 <: Unsigned 128 end

####################
###### Struct ######
####################

struct Foo
    bar  # Any
    baz::Int
    qux::Float64
end

foo = Foo("Hello", 23, 1.5)

typeof(foo) # Foo

fieldnames(Foo) # (:bar, :baz, :qux)

foo.bar # "Hello"

# Composite objects declared with struct are immutable; they cannot be modified after construction.

mutable struct Bar
    baz
    qux::Float64
end

bar = Bar("Hello", 1.5)

bar.qux = 2.0

bar.baz = 1//2

#########################
###### Type Unions ######
#########################

IntOrString = Union{Int,AbstractString}

1 :: IntOrString

"Hello" :: IntOrString

# 1.0 :: IntOrString - ERROR

##############################
###### Parametric Types ######
##############################

# An important and powerful feature of Julia's type system is that it is parametric: types can
# take parameters, so that type declarations actually introduce a whole family of new types – one
# for each possible combination of parameter values.

struct Point{T}
    x::T
    y::T
end

# Point{Float64}
# Point{AbstractString}

# Point{Float64} <: Point
# true

# Point{AbstractString} <: Point
# true

# Point{Float64} <: Point{Real}
# false

# !!!! Even though Float64 <: Real we DO NOT have Point{Float64} <: Point{Real}. !!!!

# Define a method that accepts all arguments of type Point{T} where T is a subtype of Real
function norm(p::Point{<:Real})
    sqrt(p.x^2 + p.y^2)
end

p = Point{Float64}(1.0, 2.0)

typeof(p) # Point{Float64}(1.0, 2.0)

p1 = Point(1,2)

typeof(p1) # Point{Int64}(1, 2)

# p3 = Point(1, 2.5) !!!! Error !!!!

#######################################
###### Parametric Abstract Types ######
#######################################

abstract type Pointy{T} end

struct Point{T} <: Pointy{T}
    x::T
    y::T
end

# Point{Float64} <: Pointy{Float64}
# true

# Point{Float64} <: Pointy{Real}
# false

# Point{Float64} <: Pointy{<:Real}
# true

struct DiagPoint{T} <: Pointy{T}
    x::T
end

# Now both Point{Float64} and DiagPoint{Float64} are implementations of the Pointy{Float64} abstraction,
# and similarly for every other possible choice of type T. This allows programming to a common interface
# shared by all Pointy objects, implemented for both Point and DiagPoint

struct Point{T<:Real} <: Pointy{T}
    x::T
    y::T
end

# Real-world example of how all this parametric type machinery can be useful

struct Rational{T<:Integer} <: Real
    num::T
    den::T
end

# It only makes sense to take ratios of integer values, so the parameter type T is restricted to
# being a subtype of Integer, and a ratio of integers represents a value on the real number line,
# so any Rational is an instance of the Real abstraction.

#########################
###### Tuple Types ######
#########################

struct Tuple2{A,B}
    a::A
    b::B
end

typeof((1,"foo",2.5)) # Tuple{Int64, String, Float64}

# Vararg Tuple types

mytupletype = Tuple{AbstractString,Vararg{Int}}
# Tuple{AbstractString, Vararg{Int64, N} where N}

isa(("1",), mytupletype)
# true

# Named Tuple Types

typeof((a=1,b="hello"))
# NamedTuple{(:a, :b), Tuple{Int64, String}}

# The @NamedTuple macro provides a more convenient struct-like syntax for declaring NamedTuple types
# via key::Type declarations, where an omitted ::Type corresponds to ::Any.

@NamedTuple{a::Int, b::String}
# NamedTuple{(:a, :b), Tuple{Int64, String}}

@NamedTuple begin
    a::Int
    b::String
end

# NamedTuple{(:a, :b), Tuple{Int64, String}}

############################
###### UnionAll Types ######
############################

const T1 = Array{Array{T,1} where T, 1}
# Vector{Vector{T} where T} (alias for Array{Array{T, 1} where T, 1})

const T2 = Array{Array{T, 1}, 1} where T
# Array{Vector{T}, 1} where T

# Type T1 defines a 1-dimensional array of 1-dimensional arrays; each of the inner arrays consists
# of objects of the same type, but this type may vary from one inner array to the next. On the other hand,
# type T2 defines a 1-dimensional array of 1-dimensional arrays all of whose inner arrays must have the
# same type. Note that T2 is an abstract type, e.g., Array{Array{Int,1},1} <: T2, whereas T1 is a
# concrete type. As a consequence, T1 can be constructed with a zero-argument constructor a=T1() but
# T2 cannot.

Vector{T} = Array{T, 1}

# This is equivalent to const Vector = Array{T,1} where T.
# Writing Vector{Float64} is equivalent to writing Array{Float64,1}

####################################
###### Type{T} type selectors ######
####################################

# For each type T, Type{T} is an abstract parametric type whose only instance is the object T.

isa(Float64, Type{Float64})
# true

isa(Real, Type{Float64})
# false

isa(Real, Type{Real})
# true

# In other words, isa(A, Type{B}) is true if and only if A and B are the same object and that object is a type.

struct WrapType{T}
    value::T
end

WrapType(Float64) # default constructor, note DataType
# WrapType{DataType}(Float64)

WrapType(::Type{T}) where T = WrapType{Type{T}}(T)
# WrapType

WrapType(Float64) # sharpened constructor, note more precise Type{Float64}
# WrapType{Type{Float64}}(Float64)

####################################
###### Custom pretty-printing ######
####################################

struct Polar{T<:Real} <: Number
    r::T
    Θ::T
end

Polar(r::Real,Θ::Real) = Polar(promote(r,Θ)...)
# Polar

max_stars_repo_name: mauricioriva/learn-julia

max_issues_repo_name: mauricioriva/learn-julia

max_forks_repo_name: mauricioriva/learn-julia


text:
## save window size, pan positions, etc
## create folder if necessary

mutable struct Project
    name::String
    path::String
    files::Array{String, 1}
    scroll_position::Array{Float64, 1}
    ntbook_idx::Integer
    main_window::MainWindow

    Project(main_window::MainWindow, name::String) =
        new(name, "", String[], Float64[], 1, main_window)
end

#let's not serialize main_window
JSON.lower(w::Project) = Dict(
    "name" => w.name,
    "path" => w.path,
    "files" => w.files,
    "scroll_position" => w.scroll_position,
     "ntbook_idx" => w.ntbook_idx
 )

function update!(w::Project)
    editor = _editor(w.main_window)
    w.path = pwd()
    w.files = String[]
    w.scroll_position = Float64[]
    w.ntbook_idx = index(editor)

    for i=1:length(editor)
        t = get_tab(editor, i)
        if typeof(t) == EditorTab && t.filename != ""#in case we want to have something else in the editor

            adj = get_gtk_property(t, :vadjustment, GtkAdjustment)
            push!(w.scroll_position, get_gtk_property(adj, :value, AbstractFloat))

            push!(w.files, t.filename)
        end
    end
end

#upgrade smoothly from the old project system
function upgrade_project()
    !isdir( joinpath(HOMEDIR, "config", "projects") ) && mkdir( joinpath(HOMEDIR, "config", "projects") )
    if !isfile(joinpath(HOMEDIR, "config", "projects", "default.json"))
        isfile(joinpath(HOMEDIR, "config", "project")) &&
        cp(joinpath(HOMEDIR, "config", "project"), joinpath(HOMEDIR, "config", "projects", "default.json"))
    end
end

function save(w::Project)
    update!(w)
    !isdir( joinpath(HOMEDIR, "config", "projects") ) && mkdir( joinpath(HOMEDIR, "config", "projects") )
    open( joinpath(HOMEDIR, "config", "projects", "$(w.name).json") , "w") do io
        JSON.print(io, w, 4)
    end
end

function load(w::Project)
    !isdir( joinpath(HOMEDIR, "config", "projects") ) && mkdir( joinpath(HOMEDIR, "config", "projects") )
    pth = joinpath(HOMEDIR, "config", "projects", "$(w.name).json")
    if !isfile(pth)
        w.path = pwd()
        return
    end
#	println( joinpath(HOMEDIR, "config", "project"))
    j = JSON.parsefile(pth)

    if haskey(j, "name")#allow smooth upgrade
        w.name = j["name"]
    else
        w.name = "default"
    end

    w.path = j["path"]
    w.files = j["files"]
    w.scroll_position = j["scroll_position"]
    w.ntbook_idx = j["ntbook_idx"]
end


text:
module MarginalLogDensities

using Optim
using ForwardDiff
using LinearAlgebra
using SuiteSparse
using SparseArrays
using HCubature
using Distributions
using SparsityDetection
# using Symbolics
using SparseDiffTools
using NLSolversBase

export MarginalLogDensity,
    HessianConfig,
    AbstractMarginalizer,
    LaplaceApprox,
    Cubature,
    dimension,
    imarginal,
    ijoint,
    nmarginal,
    njoint,
    merge_parameters,
    split_parameters

abstract type AbstractMarginalizer end

struct LaplaceApprox <: AbstractMarginalizer end

struct Cubature{T} <: AbstractMarginalizer
    upper::AbstractVector{T}
    lower::AbstractVector{T}
end

"""
    `MarginalLogDensity(logdensity, n, imarginal, [method=LaplaceApprox()])`

Construct a callable object which wraps the function `logdensity` and
integrates over a subset of its arguments.
* `logdensity` : function taking a vector of `n` parameters and returning a positive
log-probability (e.g. a log-pdf, log-likelihood, or log-posterior).
* `imarginal` : Vector of indices indicating which arguments to `logdensity` to marginalize
* `method` : How to perform the marginalization.  Defaults to `LaplaceApprox()`; `Cubature()`
is also available.

If `length(imarginal) == m`, then the constructed `MarginalLogDensity` object it `mld`
can be called as `mld(θ)`, where `θ` is a vector with length `n-m`.  It can also be called
as `mld(u, θ)`, where `u` is a length-`m` vector of the marginalized variables.  In this
case, the return value is the same as the full conditional `logdensity` with `u` and `θ`
"""
struct MarginalLogDensity{TI<:Integer, TM<:AbstractMarginalizer,
        TV<:AbstractVector{TI}, TF, THP}
    logdensity::TF
    n::TI
    imarginal::TV
    ijoint::TV
    method::TM
    hessconfig::THP
end

function MarginalLogDensity(logdensity, n, im, method=LaplaceApprox(), forwarddiff_sparsity=false)
    ij = setdiff(1:n, im)
    hessconfig = HessianConfig(logdensity, im, ij, forwarddiff_sparsity)
    return MarginalLogDensity(logdensity, n, im, ij, method, hessconfig)
end

# Was using this for testing/troubleshooting, will probably delete later
# """
#     `num_hessian_sparsity(f, x, [δ=1.0])`
#
# Calculate the sparsity pattern of the Hessian matrix of function `f`. This is a brute-force
# approach, but more robust than the one in SparsityDetection
# """
# function num_hessian_sparsity(f, x, δ=1.0)
#     N = length(x)
#     g(x) = ForwardDiff.gradient(f, x)
#     y = g(x)
#     ii = Int[]
#     jj = Int[]
#     vv = Float64[]
#     for j in 1:N
#         x[j] += δ
#         yj = g(x)
#         di = findall(.! (yj .≈ y))
#         for i in di
#             push!(jj, j)
#             push!(ii, i)
#             push!(vv, 1.0)
#         end
#         x[j] -= δ
#     end
#     return sparse(ii, jj, vv)
# end

struct HessianConfig{THS, THC, TI<:Integer, TD, TG}
    Hsparsity::THS
    Hcolors::THC
    ncolors::TI
    D::TD
    Hcomp_buffer::TD
    G::TG
    δG::TG
end

function HessianConfig(logdensity, imarginal, ijoint, forwarddiff_sparsity=false)
    x = ones(length(imarginal) + length(ijoint))
    if forwarddiff_sparsity
        println("Detecting Hessian sparsity via ForwardDiff...")
        H = ForwardDiff.hessian(logdensity, x)
        Hsparsity = sparse(H)[imarginal, imarginal] .!= 0
    else
        println("Detecting Hessian sparsity via SparsityDetection...")
        Hsparsity = hessian_sparsity(logdensity, x)[imarginal, imarginal]
    end
    Hcolors = matrix_colors(Hsparsity)
    D = hcat([float.(i .== Hcolors) for i in 1:maximum(Hcolors)]...)
    Hcomp_buffer = similar(D)
    G = zeros(length(imarginal))
    δG = zeros(length(imarginal))
    return HessianConfig(Hsparsity, Hcolors, size(Hcolors, 2), D, Hcomp_buffer, G, δG)
end

function sparse_hessian!(H, f, g!, θ, hessconfig::HessianConfig, δ=sqrt(eps(Float64)))
    nc = hessconfig.ncolors
    for j in one(nc):nc
        g!(hessconfig.G, θ)
        g!(hessconfig.δG, θ + δ * @view hessconfig.D[:, j])
        hessconfig.Hcomp_buffer[:, j] .= (hessconfig.δG .- hessconfig.G) ./ δ
    end
    ii, jj, vv = findnz(hessconfig.Hsparsity)
    for (i, j) in zip(ii, jj)
        H[i, j] = hessconfig.Hcomp_buffer[i, hessconfig.Hcolors[j]]
    end
end

function sparse_hessian(f, g!, θ,  hessconfig::HessianConfig, δ=sqrt(eps(Float64)))
    i, j, v = findnz(hessconfig.Hsparsity)
    H = sparse(i, j, zeros(eltype(θ), length(v)))
    sparse_hessian!(H, f, g!, θ, hessconfig, δ)
    return H
end

function MarginalLogDensity(logdensity::Function, n::TI,
        imarginal::AbstractVector{TI}; method=LaplaceApprox(), forwarddiff_sparsity=false) where {TI<:Integer}
    ijoint = setdiff(1:n, imarginal)
    hessconfig = HessianConfig(logdensity, imarginal, ijoint, forwarddiff_sparsity)
    mld  = MarginalLogDensity(logdensity, n, imarginal, ijoint, method, hessconfig)
    return mld
end

dimension(mld::MarginalLogDensity) = mld.n
imarginal(mld::MarginalLogDensity) = mld.imarginal
ijoint(mld::MarginalLogDensity) = mld.ijoint
nmarginal(mld::MarginalLogDensity) = length(mld.imarginal)
njoint(mld::MarginalLogDensity) = length(mld.ijoint)

function merge_parameters(θmarg::AbstractVector{T1}, θjoint::AbstractVector{T2}, imarg, ijoint) where {T1,T2}
    N = length(θmarg) + length(θjoint)
    θ = Vector{promote_type(T1, T2)}(undef, N)
    θ[imarg] .= θmarg
    θ[ijoint] .= θjoint
    return θ
end

split_parameters(θ, imarg, ijoint) = (θ[imarg], θ[ijoint])

function (mld::MarginalLogDensity)(θmarg::AbstractVector{T1}, θjoint::AbstractVector{T2}) where {T1, T2}
    θ = merge_parameters(θmarg, θjoint, imarginal(mld), ijoint(mld))
    return mld.logdensity(θ)
end

function (mld::MarginalLogDensity)(θjoint::AbstractVector{T}, verbose=false) where T
    integral = _marginalize(mld, θjoint, mld.method, verbose)
    return integral
end

############################################################################################
# This section implements a `logabsdet` method for sparse LU decompositions, so we can
# calculate the Laplace approximation for sparse Hessians.  This is an import from the
# future; the method will be in the next Julia release but it's copy-pasted in here for now
# so that this package will work.

# compute the sign/parity of a permutation
function _signperm(p)
    n = length(p)
    result = 0
    todo = trues(n)
    while any(todo)
        k = findfirst(todo)
        todo[k] = false
        result += 1 # increment element count
        j = p[k]
        while j != k
            result += 1 # increment element count
            todo[j] = false
            j = p[j]
        end
        result += 1 # increment cycle count
    end
    return ifelse(isodd(result), -1, 1)
end

function LinearAlgebra.logabsdet(F::SuiteSparse.UMFPACK.UmfpackLU{T, TI}) where {T<:Union{Float64,ComplexF64}, TI<:Integer} # return log(abs(det)) and sign(det)
    n = LinearAlgebra.checksquare(F)
    issuccess(F) || return log(zero(real(T))), zero(T)
    U = F.U
    Rs = F.Rs
    p = F.p
    q = F.q
    s = _signperm(p)*_signperm(q)*one(real(T))
    P = one(T)
    abs_det = zero(real(T))
    @inbounds for i in 1:n
        dg_ii = U[i, i] / Rs[i]
        P *= sign(dg_ii)
        abs_det += log(abs(dg_ii))
    end
    return abs_det, s * P
end
############################################################################################

function _marginalize(mld::MarginalLogDensity, θjoint::AbstractVector{T},
        method::LaplaceApprox, verbose) where T
    N = nmarginal(mld)
    θmarginal0 = ones(T, N)
    f = (θmarginal) -> -mld(θmarginal, θjoint)
    gconfig = ForwardDiff.GradientConfig(f, θmarginal0)
    g! = (G, x) -> ForwardDiff.gradient!(G, f, x, gconfig)
    h! = (H, x) -> sparse_hessian!(H, f, g!, x, mld.hessconfig)
    H0 = sparse_hessian(f, g!, θmarginal0, mld.hessconfig)
    td = TwiceDifferentiable(f, g!, h!, θmarginal0, zero(T), zeros(T, N), H0)

    verbose && println("Optimizing...")
    opt = optimize(td, θmarginal0)
    verbose && println("Calculating Hessian at mode...")
    h!(H0, opt.minimizer)
    verbose && println("Laplace approximating...")
    integral = -opt.minimum + 0.5 * (log((2π)^N) - logdet(H0))
    return integral
end

function _marginalize(mld::MarginalLogDensity, θjoint, method::Cubature, verbose)
    f = θmarginal -> exp(mld(θmarginal, θjoint))
    int, err = hcubature(f, method.lower, method.upper)
    return log(int)
end

end # module


text:
module FunHPC

using Comm, Par, Refs

end


text:
using Test, FrameFun.Platforms, BasisFunctions

@testset "Generic Platforms" begin
    P = platform(Fourier(10))
    @test dictionary(P,param(Fourier(10))) == Fourier(10)
    @test P[10] == Fourier(10)
    @test operator(dualdictionary(P,10,measure(P))) ≈ operator(dual(Fourier(10)))
    @test DictionaryStyle(P)==BasisStyle()
    @test SamplingStyle(P) == InterpolationStyle()
    @test SolverStyle(P,SamplingStyle(P)) == InverseStyle()
    @test param_first(P) == 10
    @test param_double(P,10) == 20
    @test param_increment(P,10) == 11
    @test param_inbetween(P,10,20) == 15
    @test correctparamformat(P,10)
    @test !correctparamformat(P,(10,))

    @test elements(platform(Fourier(10)^2)) == (P,P)
    @test element(platform(Fourier(10)^2),1) == P

    P = platform(Fourier(10)^2)
    @test correctparamformat(P,(10,10))
    @test !correctparamformat(P,(10,))
    @test !correctparamformat(P,10)
    @test param_first(P) == (10,10)
    @test param_double(P,(10,20)) == (20,40)
    @test param_increment(P,(10,11)) == (11,12)
    @test param_inbetween(P,(10,10),(20,10)) == (15,10)
    @test ProductPlatform(platform(Fourier(10)),2) == P
    @test element(P,1) == platform(Fourier(10))
    length(elements(P))==2
    @test dictionary(P,10) == Fourier(10)^2
    @test dictionary(P,param(Fourier(10)^2)) == Fourier(10)^2
    @test P[10] == Fourier(10)^2
    @test dictionary(P,(11,12)) == Fourier(11)⊗Fourier(12)
    @test dualdictionary(P,(10,10), measure(Fourier(10)^2)) isa TensorProductDict
    isbasis(Fourier(10)^2)
    @test DictionaryStyle(P)==BasisStyle()
    @test SamplingStyle(P) == ProductSamplingStyle(InterpolationStyle(),InterpolationStyle())
    @test SolverStyle(P,SamplingStyle(P)) == ProductSolverStyle(InverseStyle(),InverseStyle())
end

@testset "Platform Styles" begin
    @test SamplingStyle(platform(Fourier(10))) == InterpolationStyle()
    @test SolverStyle(platform(Fourier(10)), InterpolationStyle()) == InverseStyle()
    @test SolverStyle(platform(Fourier(10)), OversamplingStyle()) == DirectStyle()
    @test SolverStyle(platform(Fourier(10)), GridStyle()) == DirectStyle()
    @test SolverStyle(platform(Fourier(10)), GenericSamplingStyle()) == DirectStyle()

    @test SolverStyle(FrameStyle(), platform(Fourier(10)), InterpolationStyle()) == AZStyle()
    @test SolverStyle(FrameStyle(), platform(Fourier(10)), OversamplingStyle()) == AZStyle()
    @test SolverStyle(FrameStyle(), platform(Fourier(10)), GridStyle()) == AZStyle()
    @test SolverStyle(FrameStyle(), platform(Fourier(10)), GenericSamplingStyle()) == AZStyle()

    @test DictionaryStyle(platform(Fourier(10))) == BasisStyle()
    @test DictionaryStyle(platform(MultiDict((Fourier(10),Fourier(10))))) == UnknownDictionaryStyle()
    @test elements(SolverStyle(platform(Fourier(10)^2), SamplingStyle(platform(Fourier(10)^2)))) == (InverseStyle(), InverseStyle())
    @test ProblemStyle(platform(Fourier(10)))  == DictionaryOperatorStyle()

end


text:
# This Example app uses Salsa to build an incremental Spreadsheet, which is correct by
# construction.
#
# Incremental means that whenever you change the contents of any one cell, the system
# performs _only the minimum necessary computations_ to update the displayed output. That
# is, if you update cell A1, and it is used in the formula for B1, only those two cells'
# values will be recomputed and re-displayed, but any other computations needed to compute
# the display values for other cells will be remembered from when they were last updated,
# and will not need to be recomputed. This is achieved via Salsa's memoization of the
# results of `@derived` functions, and its automatic dependency detection to determine
# when to invalidate those cached values.
#
# We say it's correct by construction because we did not have to manually write _any code_
# to achieve this incrementality. Instead, below we write simple, straight-line code that
# computes the values to display for the _entire spreadsheet_, relying on Salsa to
# automatically determine which intermediate values to keep and which to recompute whenever
# a cell is modified.
#
# The code below stores the text values _entered by users_ as Salsa Inputs, and provides
# functions to _derive_ the output display strings for the entire spreadsheet, given only
# those inputs.
# In order to guarantee correctness, we only have to follow these simple rules:
#   - Do not mutate values stored in Salsa. To make this easier, this example _only uses
#     immutable types_.
#   - Derived functions must be stateless, pure functions (depend on no outside state, and
#     have no side-effects), which the below functions achieve by only depending on the
#     `@input`s on the Spreadsheet.
# In order to achieve good incremental performance, we simply have to consider which
# functions we mark as `@derived` functions, and how we split up our computations to provide
# good caching and invaldiation points. You'll see some comments on this inline, below.

using Salsa

# TODO: Change this to just use a Tuple{Symbol,Int}, since we're not using a sparse
# representation, not a matrix, anyway.
# A Cell is identified by its position. We use this identifier as the key in all the
# derived functions, below, so it's good to use the simplest possible identifier.
const CellId = Tuple{Int,Int}

# --- Salsa Inputs ------------------------------------------------------------------------
# Below we declare the Salsa inputs that we will use for our spreadsheet app. These inputs
# store values that come from users, and are the only part of the Salsa state that we can
# manipulate directly. When we make changes to the cells' source text, those changes will
# ripple through the Salsa derived functions that ultimately compute the displayed values.

# NOTE: The "Manifest Pattern"
# Notice that we have a single input that stores all the cells that the user has
# entered, and a _separate_ input that maps from those cell IDs to their text values.
# This is called the "manifest pattern," because you have one input that acts as the
# manifest (valid_cells), containing all the keys to the other map input(s).
# This is important for two reasons:
#   1. This is a syntactic requirement, because Salsa does not track changes to the
#      keys in a map, it only tracks changes to _values_.
#   2. There might be many separate things that you want to record for a given Cell:
#      Currently, we store only its text, but in the future we may also want to store
#      its color, it's formatting (bold, italics, ...), etc. For the best incremental
#      performance, you will want to keep these inputs _separate_ ("normalized"), so you
#      keep one manifest that contains the keys to all the other maps.
# NOTE: Storing Immutable Values
# For correctness, you should only use _immutable values_ in Salsa. However, for now
# for performance here we use a mutable Set. It would be better to switch to an immutable
# Set once there is better support for immutable collections. For example, we could use a
# proper immutable collection datastructure, such as those provided by
# https://github.com/JuliaCollections/FunctionalCollections.jl.
@declare_input valid_cells(rt) :: Set{CellId}
@declare_input stored_cell_text(rt, cid::CellId) :: String

# Create a struct to represent our Salsa state, which encapsulates the Salsa state,
# and initializes it with valid initial values. A new instance of a Spreadsheet
# will have a clean salsa state, with nothing in the inputs and no derived values cached.
struct Spreadsheet
    rt::Salsa.DefaultRuntime
    # The default Constructor initializes the manifest to be empty (a common part of the
    # "Manifest Pattern").
    function Spreadsheet()
        rt = Runtime()
        set_valid_cells!(rt, Set{CellId}())
        new(rt)
    end
end

function set_cell_text!(ss::Spreadsheet, args...)
    set_cell_text!(ss.rt, args...)
end

# The manifest and the attributes (valid_cells and cell_text) must be maintained together,
# so we should always use this Setter function to update the cell text.
function set_cell_text!(salsa::Runtime, id::CellId, text::String)
    if isempty(text)
        current_valid_cells = valid_cells(salsa)
        if id in current_valid_cells
            delete_stored_cell_text!(salsa, id)
            # Invalidate the cell
            set_valid_cells!(salsa, filter(x->x!=id, current_valid_cells))
        end
    else
        new_valid_cells = copy(valid_cells(salsa))
        push!(new_valid_cells, id)
        set_valid_cells!(salsa, new_valid_cells)
        set_stored_cell_text!(salsa, id, text)
    end
end

# This is the first derived function we've seen. It acts like a filter, which generates
# empty text ("") for strings that the user hasn't set directly, which allows us to use a
# sparse representation to store the user-entered cell text values.
# Marking this function as `@derived` causes Salsa to memoize its return values based on
# the input arguments, and also to track any _other_ derived functions or Salsa inputs that
# this function accesses, so that if the values from those functions _change_, this
# function's cached value is invalidated and it will be recomputed the next time it's
# called.
# NOTE: On Tuning the `@derived` annotation
# From a performance perspective, it might not actually make sense to make this a derived
# function, since it's body is very cheap, but we did so for explanation purposes. You will
# need to make trade-off decisions between cache size and cpu performance, motivated by
# running performance measurements. This tuning is common to working with Salsa.
@derived function cell_text(rt::Runtime, id::CellId)::String
    if cell_is_set(rt, id)
        stored_cell_text(rt,id)
    else
        ""
    end
end
# NOTE: "Firewall Pattern"
# Notice that we have a derived function here for what seems like a very cheap operation.
# This is called the "Firewall Pattern" (or the "Selector Pattern"). This derived function
# acts as a recompuation barrier by limiting the exposure to a potentially quite commonly
# modified input: the `valid_cells` manifest.
# Consider: _every time_ a new cell is added or removed, this variable will change. Without
# this barrier, that means that any function that checks if its `id` is in the manifest
# would have to rerun. In this case, if this function weren't marked `@derived`, it would
# mean that whenever you add or remove a CellId, _all the values stored in the `cell_text()`
# function would be invalidated, and would have to rerun. This barrier function prevents
# that.
@derived cell_is_set(rt::Runtime, id::CellId)::Bool = id in valid_cells(rt)

# --- Cell Type Computations -------------------------------------------------------------
# Here we have some derived functions which determine what type of cell we have, based
# on the cell's text. They look up the cell's contents based on its ID, then do some
# computation to determine the cell's type.
#
# Notice that the functions all take the *CellId* as arguments, rather than the cell's text.
# This is important for cache performance: when the cell's contents change, the _key_ for
# the Salsa cache (which method was called, and what were its arguments) remains the same,
# so that the new returned value _overwrites_ the previous one, rather than accumulating
# beside it in the cache.
#
# Notice also that we include return-type annotations on these functions, because Salsa
# uses them to construct strongly-typed Dictionaries for storing the cached values.

@derived cell_is_empty(rt::Runtime, id::CellId)::Bool = isempty(cell_text(rt, id))

@derived function cell_is_formula(rt::Runtime, id::CellId)::Bool
    text = cell_text(rt, id)
    !isempty(text) && first(text) === '='
end

# --- Cell Value Computations -------------------------------------------------------------
# Finally, we have derived functions to determine the value and display text for a given
# cell.

# AbstractUserFacingException are caught and displayed, rather than allowed to propogate and
# terminate execution. These are supposed to be in response to an error in the
# _user-supplied_ logic, rather than an error in our Julia program.
abstract type AbstractUserFacingException <: Base.Exception end

# An Exception that we throw when we encounter a bug in the user's program.
struct UserError <: AbstractUserFacingException
    err
end
# Define pretty-print method to print the UserError in the UI.
Base.print(io::IO, err::UserError) = print(io, err.err)

# This is the main derived function that runs the whole Spreadsheet App. It parses the
# contents of the user's cell into a valid Julia type, and if the cell contains a Formula,
# it evaluates the formula to compute a value.
# Notice that this function can be _recursive_ thanks to evaluating formulas:
#  - The call to `replace_varnames` will find valid `CellId`s in the user's expression,
#    and recursively call `call_value(rt, var_id)` to compute _those_ values. The cached
#    values from any intermediate computations we've written above will be automatically
#    reused.
#  - Salsa will dynamically track the dependencies between the calls to `cell_value()`,
#    which automatically provides the dependencies between cells, based on the user-provided
#    values. Whenever formulas change, they'll be recomputed and the dependencies will be
#    automatically updated.
@derived function cell_value(rt::Runtime, id::CellId)::Any
    if cell_is_empty(rt, id)
        ""
    elseif cell_is_formula(rt, id)
        try
            expr = Meta.parse(cell_text(rt, id)[2:end])
            expr = replace_varnames(rt, expr)
            Meta.eval(expr)
        catch e
            if e isa AbstractUserFacingException
                e  # Return the exception
            elseif !(e isa InterruptException)
                UserError(e)
            else
                rethrow()
            end
        end
    else
        try
            Meta.parse(cell_text(rt, id))
        catch
            # If it doesn't parse, treat it as a string! :)
            cell_text(rt, id)
        end
    end
end

function replace_varnames(rt::Runtime, expr)
    # Recursively walk through the expression replacing references to cell ids and ranges
    # with their values. The recursive walk!() function is defined below and called on
    # `expr` at the very end.

    walk!(node) = node
    # Convert variables that look like :A1 to the value stored at cell-id "A1"
    function walk!(var::Symbol)
        var_id = cell_id_from_name(var)
        if var_id === nothing
            # If it's not a cell id, we assume it's a builtin symbol, like :sum or :+
            var
        else
            # NOTE: Here is the recursive call to cell_value(rt, var_id). Without any
            # explicit work on our part, this creates a dependency from the current CellId
            # being computed to the
            return QuoteNode(cell_value(rt, var_id))
        end
    end
    # Convert Range expressions containing variable names :(A1:B1) to a 2D matrix of the
    # values in the cells stored in that range.
    function walk!(e::Expr)
        if e.head === :call && e.args[1] === Symbol(":")
            length(e.args) === 3 || throw(UserError("BAD RANGE: `$e`"))
            r_start = cell_id_from_name(e.args[2])
            r_end = cell_id_from_name(e.args[3])
            # Return a 2D Matrix of the values
            tasks = Task[Threads.@spawn cell_value(rt, (r,c))
             for r in r_start[1]:r_end[1], c in r_start[2]:r_end[2]
            ]
            map(fetch, tasks)
        else
            # Recursive walk on all elements of the Expr
            Expr(e.head, walk!.(e.args)...)
        end
    end

    # Walk the Expr and do the replacements.
    walk!(expr)
end

# Finally, we have the function that converts a cell's _actual value_ into a String.
# This is also a `@derived` function because printing the value to a string might be
# expensive, so it might be worth caching. Again, this tradeoff should be made after
# performance measurements.
#
# The final missing piece of our spreadsheet app is in the UI.jl file, where the display
# functions loop over the entire range of the spreadsheet currently in-view from the UI, and
# call this function on each cell.
@derived function cell_display_str(rt::Runtime, id::CellId)::String
    try
        v = cell_value(rt, id)
        if v isa Exception
            "#ERR#"
        else
            string(v)
        end
    catch e
        "#CRSH"
    end
end

function column_idx_to_name(i)
    out = ""
    while i > 0
        l = (i-1)%26+1
        out = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"[l] * out
        i ÷= 26
    end
    out
end
function column_name_to_idx(name)
    name = uppercase(name)
    i = 0
    l = 0
    while l < length(name)
        i += 26^l*findfirst(name[end-l], "ABCDEFGHIJKLMNOPQRSTUVWXYZ")[1]
        l += 1
    end
    i
end

function cell_id_from_name(var::Symbol)
    s = String(var)
    s = uppercase(s)
    if (m = match(r"([A-Z]+)([0-9]+)", s); m !== nothing)
        col,row = m.captures
        row,col = parse(Int,row),column_name_to_idx(col)
        return (row,col)
    else
        nothing
    end
end


text:
"""
    Box([origin::AbstractVector{T}=zeros(),] dimensions::AbstractVector{T})

A [`Box`](@ref) for 2D and 3D with axis aligned sides, defined by dimensions and origin (at the center).
"""
struct Box{T,Dim} <: Shape{T,Dim}
    origin::SVector{Dim,T}
    dimensions::SVector{Dim,T}
end

# Alternative constructor, where bottomleft and topright corners are specified
function Box(origin::AbstractVector{T}, dimensions::AbstractVector{T}) where T
    Dim = length(origin)
    Box{T,Dim}(origin, dimensions)
end

function Box(origin::NTuple{Dim,T}, dimensions::NTuple{Dim,T}) where {T,Dim}
    Box{T,Dim}(origin, dimensions)
end

function Box(corners::Vector{S}) where S<:AbstractVector
    centre = mean(corners)
    cs = hcat([abs.(c - centre) for c in corners]...)
    dimensions = 2 .* abs.(corners[1] - centre)

    if dimensions != 2 .* abs.(corners[2] - centre)
        @error "expected $corners to be a vector of corners of a $(length(corners[1])) dimensional box, with sides aligned with the coordinate axis."
    end

    Box(centre,dimensions)
end

Box(dimensions::AbstractVector{T}) where T = Box(zeros(T,length(dimensions)),dimensions)

Rectangle(bottomleft::Union{AbstractVector{T},NTuple{2,T}}, topright::Union{AbstractVector{T},NTuple{2,T}}) where T = Box([bottomleft,topright])

name(shape::Box) = "Box"
name(shape::Box{T,2}) where T = "Rectangle"

outer_radius(r::Box) = sqrt(sum(r.dimensions .^ 2))
volume(box::Box) = prod(box.dimensions)

import Base.issubset
function issubset(inner::Box, outer::Box)
    corners = box_corners(inner)
    return all(c ∈ outer for c in corners)
end

import Base.in
function in(x::AbstractVector, b::Box{T})::Bool where T
    all(abs.(x .- b.origin) .<= b.dimensions ./ T(2))
end

import Base.(==)
function ==(r1::Box, r2::Box)
    r1.origin == r2.origin &&
    r1.dimensions  == r2.dimensions
end

import Base.isequal
function isequal(r1::Box, r2::Box)
    isequal(r1.origin, r2.origin) &&
    isequal(r1.dimensions , r2.dimensions )
end

function iscongruent(r1::Box, r2::Box)
    r1.dimensions  == r2.dimensions
end

function congruent(r::Box, x)
    Box(x, r.dimensions)
end

# Box bounds itself
bounding_box(b::Box) = b

"Returns a vector of SVector with the coordinates of the corners of the box"
function box_corners(b::Box{T,Dim}) where {T,Dim}
    d1 = origin(b) - b.dimensions ./ T(2)
    d2 = origin(b) + b.dimensions ./ T(2)
    ds =[d1,d2]

    corners = map(0:(2^Dim-1)) do i
        switches = 1 .+ [parse(Int,c)  for c in bitstring(i)[end-Dim+1:end]]
        [ds[switches[j]][j] for j in eachindex(switches)]
    end

    return corners
end

"Return SVector with the coordinates of the bottom left of a rectangle"
bottomleft(rect::Box{T,2}) where T = origin(rect) .- SVector(rect.dimensions) ./ 2

"Return SVector with the coordinates of the top right of a rectangle"
topright(rect::Box{T,2}) where T = origin(rect) .+ SVector(rect.dimensions) ./ 2

function boundary_functions(rect::Box{T,2}) where T
    w, h = rect.dimensions
    function x(t)
        check_boundary_coord_range(t)
        if     t <= 1//4 x = bottomleft(rect)[1] + 4t * w
        elseif t <= 2//4 x = topright(rect)[1]
        elseif t <= 3//4 x = topright(rect)[1] - 4*(t-2//4) * w
        else             x = bottomleft(rect)[1]
        end
    end

    function y(t)
        check_boundary_coord_range(t)
        if     t <= 1//4 x = bottomleft(rect)[2]
        elseif t <= 2//4 x = bottomleft(rect)[2] + 4*(t-1//4) * h
        elseif t <= 3//4 x = topright(rect)[2]
        else             x = topright(rect)[2] - 4*(t-3//4) * h
        end
    end

    return x, y
end


text:
"""
$(TYPEDEF)

A system of ordinary differential equations.

# Fields
$(FIELDS)

# Example

```
using ModelingToolkit

@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

de = ODESystem(eqs,[x,y,z],[σ,ρ,β])
```
"""
struct ODESystem <: AbstractODESystem
    """The ODEs defining the system."""
    eqs::Vector{Equation}
    """Independent variable."""
    iv::Variable
    """Dependent (state) variables."""
    states::Vector{Variable}
    """Parameter variables."""
    ps::Vector{Variable}
    """
    Time-derivative matrix. Note: this field will not be defined until
    [`calculate_tgrad`](@ref) is called on the system.
    """
    tgrad::RefValue{Vector{Expression}}
    """
    Jacobian matrix. Note: this field will not be defined until
    [`calculate_jacobian`](@ref) is called on the system.
    """
    jac::RefValue{Matrix{Expression}}
    """
    `Wfact` matrix. Note: this field will not be defined until
    [`generate_factorized_W`](@ref) is called on the system.
    """
    Wfact::RefValue{Matrix{Expression}}
    """
    `Wfact_t` matrix. Note: this field will not be defined until
    [`generate_factorized_W`](@ref) is called on the system.
    """
    Wfact_t::RefValue{Matrix{Expression}}
    """
    Name: the name of the system
    """
    name::Symbol
    """
    systems: The internal systems
    """
    systems::Vector{ODESystem}
end

function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
                   systems = ODESystem[],
                   name=gensym(:ODESystem))
    iv′ = convert(Variable,iv)
    dvs′ = convert.(Variable,dvs)
    ps′ = convert.(Variable,ps)
    tgrad = RefValue(Vector{Expression}(undef, 0))
    jac = RefValue(Matrix{Expression}(undef, 0, 0))
    Wfact   = RefValue(Matrix{Expression}(undef, 0, 0))
    Wfact_t = RefValue(Matrix{Expression}(undef, 0, 0))
    ODESystem(deqs, iv′, dvs′, ps′, tgrad, jac, Wfact, Wfact_t, name, systems)
end

var_from_nested_derivative(x) = var_from_nested_derivative(x,0)
var_from_nested_derivative(x::Constant) = (missing, missing)
var_from_nested_derivative(x,i) = x.op isa Differential ? var_from_nested_derivative(x.args[1],i+1) : (x.op,i)
iv_from_nested_derivative(x) = x.op isa Differential ? iv_from_nested_derivative(x.args[1]) : x.args[1].op
iv_from_nested_derivative(x::Constant) = missing

function ODESystem(eqs; kwargs...)
    ivs = unique(skipmissing(iv_from_nested_derivative(eq.lhs) for eq ∈ eqs))
    length(ivs) == 1 || throw(ArgumentError("one independent variable currently supported"))
    iv = first(ivs)

    dvs = unique(skipmissing(var_from_nested_derivative(eq.lhs)[1] for eq ∈ eqs))
    ps = filter(vars(eq.rhs for eq ∈ eqs)) do x
        isparameter(x) & !isequal(x, iv)
    end |> collect
    ODESystem(eqs, iv, dvs, ps; kwargs...)
end

Base.:(==)(sys1::ODESystem, sys2::ODESystem) =
    _eq_unordered(sys1.eqs, sys2.eqs) && isequal(sys1.iv, sys2.iv) &&
    _eq_unordered(sys1.states, sys2.states) && _eq_unordered(sys1.ps, sys2.ps)
# NOTE: equality does not check cached Jacobian

function rename(sys::ODESystem,name)
    ODESystem(sys.eqs, sys.iv, sys.states, sys.ps, sys.tgrad, sys.jac, sys.Wfact, sys.Wfact_t, name, sys.systems)
end


text:
"""
    CALCEPH

  This module is a wrapper of CALCEPH, IMCCE planetary ephemeris access
  library. It supports INPOPxx, JPL DExxx and SPICE ephemeris.

  https://www.imcce.fr/inpop/calceph

"""
module CALCEPH

# To make Pkg aware that this dependency
# will be injected by BinaryProvider.
using Libdl

deps = abspath(joinpath(splitdir(@__FILE__)[1], "..", "deps", "deps.jl"))
if isfile(deps)
    include(deps)
else
    error("libcalceph was not found. Please run 'Pkg.build(\"CALCEPH\").")
end

struct CALCEPHException <: Exception
    msg::String
end

include("ephem.jl")
export Ephem, prefetch, CALCEPHException

include("compute.jl")
export compute

include("timespan.jl")
export timespan

include("bodyId.jl")
export naifId

include("units.jl")
export unitAU, unitKM, unitDay, unitSec, unitRad, useNaifId, outputEulerAngles, outputNutationAngles

include("orient.jl")
export orient

include("rotAngMom.jl")
export rotAngMom

include("constants.jl")
export constants

include("introspection.jl")
export timeScale, positionRecords, orientationRecords

include("fivePointStencil.jl")

include("errorHandling.jl")

end # module


text:
mutable struct Bee{T}
    sol::T
    fit::Float64
    t::Int
end

Bees = Array{Bee, 1}

Bee(solution) = Bee(solution, fit(solution.f), 0)

"""
    get_position(bee)

Get the position vector of a bee when optimize using ABC algorithm.
"""
get_position(bee::Bee) = bee.sol.x

"""
    fval(solution)

Get the fitness of a bee when optimize using ABC algorithm.
"""
fval(bee::Bee) = bee.sol.f

Base.show(io::IO, ::MIME"text/html", pop::Vector{Bee{xf_indiv}}) = show(io, "text/html", [p.sol for p in pop])

minimum(st::State{Bee{xf_indiv}}) = st.best_sol.sol.f
minimum(st::State{Bee{xfgh_indiv}}) = st.best_sol.sol.f
minimizer(st::State{Bee{xf_indiv}}) = st.best_sol.sol.x
minimizer(st::State{Bee{xfgh_indiv}}) = st.best_sol.sol.x

@inline function updateFit!(bee::Bee)
    bee.fit = fit(bee.sol.f)
end

@inline function updateFit!(bees::Bees)
    for bee in bees
        updateFit!(bee)
    end
end

function fit(fx)
    if fx >= 0.0
        return 1.0/(1.0+fx)
    end

    1.0 + abs(fx)
end

function updateBee!(bee, bee2, problem)
    D = length(bee.sol.x)
    ϕ = -1.0 + 2.0rand()

    v = ϕ*(bee.sol.x - bee2.sol.x)

    x_new = bee.sol.x + v
    replace_with_random_in_bounds!(x_new, problem.bounds)

    new_sol = create_solution(x_new, problem)

    if is_better(new_sol, bee.sol) #fx < bee.sol.f
        bee.sol = new_sol
        bee.fit = fit(new_sol.f)
        bee.t = 0
    else
        bee.t += 1
    end
end

function getK(i, N)
    j = rand(1:N)
    while j == i
        j = rand(1:N)
    end
    j
end

function employedPhase!(bees, problem, Ne)
    N = length(bees)
    for i in randperm(N)[1:Ne]
        updateBee!(bees[i], bees[getK(i, N)], problem)
    end
end

function roulettSelect(bees, sum_f)
    r = rand()
    rs = 0.0
    for i in 1:length(bees)
        rs += bees[i].fit / sum_f
        if r <= rs
            return i
        end
    end

    return length(bees)
end

function outlookerPhase!(bees, problem, No::Int)
    N = length(bees)
    sum_f = sum(map(x->x.fit, bees))

    for i=1:No
        j = roulettSelect(bees, sum_f)
        updateBee!(bees[j], bees[getK(j, N)], problem)
    end
end

function scoutPhase!(bees, problem, genPos::Function, limit::Int)
    bees_scout = filter(x->x.t >= limit, bees)
    for i in 1:length(bees_scout)
        bees_scout[i].sol = create_solution(genPos(), problem)
        bees_scout[i].t = 0
    end

    return length(bees_scout)
end

function getBestBee(bees)
    best = bees[1]
    for bee in bees
        if bee.sol.f < best.sol.f
            best = bee
        end
    end

    return best
end

function chooseBest(bees, best)
    bee_cand = getBestBee(bees)
    if bee_cand.sol.f < best.sol.f
        return deepcopy(bee_cand)
    end

    return best
end

function initialbees(N, problem)
     P = generate_population(N, problem)

     return [ Bee(sol) for sol in P ]
end


text:
import Base
using Printf

include("transition.jl")

DIAMOND = "◇"

"""
This type represents a (*nondeterministic*) *finite automaton*, defined as a tuple ``A=(N,T,transitions,q0,F)``
    Automaton(N::Iterable, T::Iterable, transitions::Iterable, q0::AbstractString, F::Iterable)
Builds an [`Automaton`](@ref) obtained from the given components.
"""
struct Automaton
    "The states of the automaton."
    N::Set
    "The transition labels."
    T::Set
    "The transition of the automata."
    transitions::Array
    "The starting state of the automation"
    q0::AbstractString
    "The set of *final* states."
    F::Set
    function Automaton(N::Iterable, T::Iterable, transitions::Iterable, q0::AbstractString, F::Iterable)
        N = isa(N, Set) ? N : Set(N)
        T = isa(T, Set) ? T : Set(T)
        F = isa(F, Set) ? F : Set(F)
        transitions = if all(x->isa(x, Transition), transitions) collect(transitions) else throw(ArgumentError("There are non-transitions in the transitions set")) end
        if ~isempty(N ∩ T) throw(ArgumentError("The set of states and input symbols are not disjoint, but have" * string(collect(N ∩ T)) * "in common.")) end
        if q0 ∉ N throw(ArgumentError("The specified q0 (" * q0 * ") is not a state.")) end
        if F ⊈ N throw(ArgumentError("The accepting states " * string(collect(F - N)) * " in F are not states.")) end
        bad_trans = [t for t ∈ transitions if t.from ∉ N || t.to ∉ N || t.label ∉ (T ∪ Set(["ε"]))]
        if ~isempty(bad_trans) throw(ArgumentError("The following transitions contain states or symbols that are neither states nor input symbols: " * string(bad_trans))) end
        return new(N, T, transitions, q0, F)
    end
end

"""
    Automaton(transitions::AbstractString, F::Union{Nothing,Set} = nothing, q0::Union{Nothing, AbstractString} = nothing)::Automaton
Builds an [`Automaton`](@ref) obtained from the given transitions.
"""
function Automaton(transitions::AbstractString, F::Union{Nothing,Set} = nothing, q0::Union{Nothing, AbstractString} = nothing)::Automaton
    transitions = parsetransitions(transitions)
    if q0 === nothing q0 = transitions[1].from end
    if F === nothing F = Set() end

    N = Set(map(x::Transition -> x.from, transitions) ∪ map(x::Transition -> x.to, transitions))
    T = Set(map(x::Transition -> x.label, transitions)) - "ε"
    return Automaton(N, T, transitions, q0, F)
end

"""
    Automaton(G::Grammar)::Automaton
Builds the [`Automaton`](@ref) corresponding to the given *regular* [`Grammar`](@ref).
"""
function Automaton(G::Grammar)::Automaton
    transitions = []

    for P ∈ G.P
        if length(P.right) > 2 throw(ArgumentError("Production " * string(P) * " has more than two symbols on the left-hand side")) end
        if length(P.right) == 2
            A = P.left
            a = P.right[1]
            B = P.right[2]
            if ~(a ∈ G.T && B ∈ G.N) throw(ArgumentError("Production " * string(P) * " right-hand side is not of the aB form")) end
            push!(transitions, Transition(A, a, B))
        elseif P.right[1] ∈ G.N
            push!(transitions, Transition(P.left, "ε", P.right))
        else
            push!(transitions, Transition(P.left, P.right[1], DIAMOND))
        end
    end

    return Automaton(G.N ∪ Set([DIAMOND]), G.T, transitions, G.S, Set([DIAMOND]))
end

### Functions ###

"""
    δ(a::Automaton, X, x)
The transition function.

This function returns the set of states reachable from the given state and input symbol.
"""
δ(a::Automaton, X, x) = Set([t.to for t ∈ a.transitions if t.from == X && t.label == x])

### Operators ###

Base.show(io::IO, a::Automaton) = Base.show(io, @sprintf "Automaton(N=%s, T=%s, transitions=%s, F=%s, q0=%s)" a.N a.T a.F a.transitions a.q0) #TODO


text:
using Pseudospectral
using Test

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


text:
using CGP
using Logging

inputs = Dict("cart_pole"=>4, "mountain_car"=>2, "pendulum"=>3)
outputs = Dict("cart_pole"=>2, "mountain_car"=>1, "pendulum"=>1)

function cart_pole(c::Chromosome)
    gravity = 9.81; m_cart = 1.0; m_pole = 0.1; tau = 0.02; force_mag = 10.0
    theta_threshold = 12 * 2 * pi / 360; x_threshold = 2.4; pole_length = 0.5
    m_total = m_cart + m_pole
    limits = [x_threshold, force_mag, theta_threshold, 2*pi]
    reward = 0
    state = 0.1*rand(4)-0.05
    x, x_dot, theta, theta_dot = state

    while ((-x_threshold <= x <= x_threshold) &&
           (-theta_threshold <= theta <= theta_threshold) &&
           reward < 100)
        # get action
        inps = [-1.0; max.(min.(state ./ limits, 1.0), -1.0)[:]]
        action = indmax(process(c, inps))
        force = force_mag * (2 * (action - 1) - 1)
        reward += 1

        # update
        temp = (force + m_pole * pole_length * theta_dot * theta_dot * sin(theta))/m_total
        theta_acc = ((gravity * sin(theta) - cos(theta) * temp) /
                     (pole_length * (4.0/3 - m_pole * cos(theta) * cos(theta) / m_total)))
        x_acc = temp - m_pole * pole_length * theta_acc * cos(theta) / m_total
        x += tau * x_dot
        x_dot += tau * x_acc
        theta += tau * theta_dot
        theta_dot += tau * theta_acc
        state = [x, x_dot, theta, theta_dot]
        # println(reward, " ", action, " ", state, " ", inps)
    end
    reward / 100
end

function mountain_car(c::Chromosome)
    min_x = -1.2; max_x = 0.6; max_v = 0.07
    goal_position = 0.45; power = 0.0015
    steps = 0
    x = -0.6 + 0.2 * rand(); v = 0.0
    best_x = x
    while ((x <= goal_position) && (steps < 100))
        inps = [0.0, 0.0, 2*((x-min_x)/(max_x-min_x))-1.0, v / max_v, 0.0]
        force = process(c, inps)[1]
        v = min(max(v + force * power - 0.0025 * cos(3 * x), -max_v), max_v)
        x = min(max(x + v, min_x), max_x)
        if (x == min_x) && (v < 0)
            v = 0
        end
        steps += 1
        best_x = max(best_x, x)
        # println(steps, " ", inps, " ", force, " ", x, " ", v, " ")
    end
    min(max(best_x / goal_position, 0.0), 1.0)
end

function pendulum(c::Chromosome)
    max_speed = 8; max_torque = 2.0; dt = 0.5;
    g = 10.; m = 1.; l = 1.

    th = 2*pi*rand() - pi
    thdot = 2*rand() - 1.0

    cost = 0.0; steps = 0;
    while steps < 100
        inps = [1.0, cos(th), sin(th), thdot / max_speed, 0.0]
        u = max_torque * process(c, inps)[2]

        angle = (((th+pi) % (2*pi)) - pi)
        cost -= angle^2 + (.1*thdot)^2 + 0.001*(u^2)

        thdot = min(max(thdot + (-3*g/(2*l) * sin(th +pi) + 3 ./(m*l^2)*u) * dt,
                        -max_speed), max_speed)
        th = th + thdot * dt
        steps += 1
        # println(steps, " ", inps, " ", u, " ", cost)
    end
    max(min((cost + 4000) / 3850, 1.0), 0.0)
    # min((-100-cost)/-100, 1.0)
end

function combined(c::Chromosome)
    score = cart_pole(c); CGP.reset!(c)
    score += mountain_car(c); CGP.reset!(c)
    score += pendulum(c)
    score
end

function repeat(c::Chromosome, func::Function)
    score = 0
    for i = 1:5
        score += func(c)
    end
    score/5
end

seed = 0
log = "log"
fitness = "pendulum"
if length(ARGS) > 0; seed = parse(Int64, ARGS[1]); end
if length(ARGS) > 1; log = ARGS[2]; end
if length(ARGS) > 2; fitness = ARGS[3]; end

# CGP.Config.init("cfg/base.yaml")
# CGP.Config.init("cfg/classic.yaml")
CGP.Config.init("cfg/test.yaml")

Logging.configure(filename=log, level=INFO)
nin = inputs[fitness]; nout = outputs[fitness]
fit = x->repeat(x, eval(parse(fitness)))

include("param_sweep.jl")
param_sweep()


text:
using Test
using ManifoldDiffEq
using Manifolds
using OrdinaryDiffEq: OrdinaryDiffEq, alg_order
using LinearAlgebra
using DiffEqBase

function test_solver_frozen(manifold_to_alg; expected_order = nothing, adaptive = false)
    expected_order !== nothing && @testset "alg_order" begin
        alg = manifold_to_alg(Sphere(2))
        @test alg_order(alg) == expected_order
    end

    @testset "Sphere" begin
        M = Sphere(2)
        A = ManifoldDiffEq.FrozenManifoldDiffEqOperator{Float64}() do u, p, t
            return cross(u, [1.0, 0.0, 0.0])
        end
        u0 = [0.0, 1.0, 0.0]
        alg = manifold_to_alg(M)
        prob = ManifoldDiffEq.ManifoldODEProblem(A, u0, (0, 2.0), M)
        sol1 = if adaptive
            solve(prob, alg)
        else
            solve(prob, alg, dt = 1 / 8)
        end

        @test sol1(0.0) ≈ u0
        @test is_point(M, sol1(1.0))
    end

    @testset "Product manifold" begin
        M = ProductManifold(Sphere(2), Euclidean(3))
        A = ManifoldDiffEq.FrozenManifoldDiffEqOperator{Float64}() do u, p, t
            return ProductRepr(cross(u.parts[1], [1.0, 0.0, 0.0]), u.parts[2])
        end
        u0 = ProductRepr([0.0, 1.0, 0.0], [1.0, 0.0, 0.0])
        alg = manifold_to_alg(M)
        prob = ManifoldDiffEq.ManifoldODEProblem(A, u0, (0, 2.0), M)
        sol1 = if adaptive
            solve(prob, alg)
        else
            solve(prob, alg, dt = 1 / 8)
        end

        @test isapprox(M, sol1(0.0), u0)
        @test is_point(M, sol1(1.0))
    end

end

function test_solver_lie(manifold_to_alg; expected_order = nothing)
    expected_order !== nothing && @testset "alg_order" begin
        action = RotationAction(Euclidean(3), SpecialOrthogonal(3))
        alg = manifold_to_alg(Sphere(2), action)
        @test alg_order(alg) == expected_order
    end

    @testset "Sphere" begin
        M = Sphere(2)
        action = RotationAction(Euclidean(3), SpecialOrthogonal(3))

        A = ManifoldDiffEq.LieManifoldDiffEqOperator{Float64}() do u, p, t
            return hat(SpecialOrthogonal(3), Matrix(I(3)), cross(u, [1.0, 0.0, 0.0]))
        end
        u0 = [0.0, 1.0, 0.0]
        alg = manifold_to_alg(M, action)
        prob = ManifoldDiffEq.ManifoldODEProblem(A, u0, (0, 2.0), M)
        sol1 = solve(prob, alg, dt = 1 / 8)

        @test sol1(0.0) ≈ u0
        @test is_point(M, sol1(1.0))
    end

end

# constructing tableaus for comparison with DiffEq on the Euclidean space
function constructME(T::Type = Float64)
    A = fill(0, 1, 1)
    c = [0]
    α = [1]
    A = map(T, A)
    α = map(T, α)
    c = map(T, c)
    return (DiffEqBase.ExplicitRKTableau(A, c, α, 1))
end

function constructCG2(T::Type = Float64)
    A = [
        0 0
        1//2 0
    ]
    c = [0, 1 // 2]
    α = [0, 1]
    A = map(T, A)
    α = map(T, α)
    c = map(T, c)
    return (DiffEqBase.ExplicitRKTableau(A, c, α, 2))
end

function constructCG3(T::Type = Float64)
    A = [
        0 0 0
        3//4 0 0
        119//216 17/108 0
    ]
    c = [0, 3 // 4, 17 // 24]
    α = [13 // 51, -2 // 3, 24 // 17]
    A = map(T, A)
    α = map(T, α)
    c = map(T, c)
    return (DiffEqBase.ExplicitRKTableau(A, c, α, 3))
end

function constructCG4a(T::Type = Float64)
    c2 = 0.8177227988124852
    c3 = 0.3859740639032449
    c4 = 0.3242290522866937
    c5 = 0.8768903263420429
    a21 = 0.8177227988124852
    a31 = 0.3199876375476427
    a41 = 0.9214417194464946
    a51 = 0.3552358559023322
    a32 = 0.0659864263556022
    a42 = 0.4997857776773573
    a52 = 0.2390958372307326
    a43 = -1.0969984448371582
    a53 = 1.3918565724203246
    a54 = -1.1092979392113465
    b1 = 0.1370831520630755
    b2 = -0.0183698531564020
    b3 = 0.7397813985370780
    b4 = -0.1907142565505889
    b5 = 0.3322195591068374

    A = [
        0 0 0 0 0
        a21 0 0 0 0
        a31 a32 0 0 0
        a41 a42 a43 0 0
        a51 a52 a53 a54 0
    ]
    c = [0, c2, c3, c4, c5]
    α = [b1, b2, b3, b4, b5]
    A = map(T, A)
    α = map(T, α)
    c = map(T, c)
    return (DiffEqBase.ExplicitRKTableau(A, c, α, 4))
end

function constructRKMK4(T::Type = Float64)
    A = [
        0 0 0 0
        1//2 0 0 0
        0 1//2 0 0
        0 0 1 0
    ]
    c = [0, 1 // 2, 1 // 2, 1]
    α = [1 // 6, 1 // 3, 1 // 3, 1 // 6]
    A = map(T, A)
    α = map(T, α)
    c = map(T, c)
    return (DiffEqBase.ExplicitRKTableau(A, c, α, 3))
end

function compare_with_diffeq_frozen(manifold_to_alg, tableau)
    M = Euclidean(2)
    # damped harmonic oscillator d^2u/dt^2 + c du/dt + ku = 0
    k = -0.25
    c = 1
    f(u, p, t) = [-c * u[1] - k * u[2], u[1]]
    A = ManifoldDiffEq.FrozenManifoldDiffEqOperator{Float64}(f)
    u0 = [-1.0, 1.0]
    alg = manifold_to_alg(M)
    tspan = (0, 2.0)
    dt = 1 / 8
    prob_frozen = ManifoldDiffEq.ManifoldODEProblem(A, u0, tspan, M)
    sol_frozen = solve(prob_frozen, alg, dt = dt)

    alg_diffeq = OrdinaryDiffEq.ExplicitRK(tableau)
    prob_diffeq = ODEProblem(A, u0, tspan)
    sol_diffeq = solve(prob_diffeq, alg_diffeq; dt = dt, adaptive = false)

    @test isapprox(sol_frozen.u, sol_diffeq.u)
end

function compare_with_diffeq_lie(manifold_to_alg, tableau)
    M = Euclidean(2)
    # damped harmonic oscillator d^2u/dt^2 + c du/dt + ku = 0
    k = -0.25
    c = 1
    f(u, p, t) = [-c * u[1] - k * u[2], u[1]]
    action = TranslationAction(Euclidean(2), TranslationGroup(2))

    A = ManifoldDiffEq.LieManifoldDiffEqOperator{Float64}(f)
    u0 = [-1.0, 1.0]
    alg = manifold_to_alg(M, action)
    tspan = (0, 2.0)
    dt = 1 / 8
    prob_lie = ManifoldDiffEq.ManifoldODEProblem(A, u0, tspan, M)
    sol_lie = solve(prob_lie, alg, dt = dt)

    alg_diffeq = OrdinaryDiffEq.ExplicitRK(tableau)
    prob_diffeq = ODEProblem(A, u0, tspan)
    sol_diffeq = solve(prob_diffeq, alg_diffeq; dt = dt, adaptive = false)

    @test isapprox(sol_lie.u, sol_diffeq.u)
end

@testset "ManifoldDiffEq" begin
    manifold_to_alg_me = M -> ManifoldDiffEq.ManifoldEuler(M, ExponentialRetraction())
    test_solver_frozen(manifold_to_alg_me; expected_order = 1)
    compare_with_diffeq_frozen(manifold_to_alg_me, constructME())

    manifold_to_alg_cg2 = M -> ManifoldDiffEq.CG2(M, ExponentialRetraction())
    test_solver_frozen(manifold_to_alg_cg2; expected_order = 2)
    compare_with_diffeq_frozen(manifold_to_alg_cg2, constructCG2())

    manifold_to_alg_cg23 = M -> ManifoldDiffEq.CG2_3(M, ExponentialRetraction())
    test_solver_frozen(manifold_to_alg_cg23; expected_order = 2, adaptive = true)

    manifold_to_alg_cg3 = M -> ManifoldDiffEq.CG3(M, ExponentialRetraction())
    test_solver_frozen(manifold_to_alg_cg3; expected_order = 3)
    compare_with_diffeq_frozen(manifold_to_alg_cg3, constructCG3())

    manifold_to_alg_cg4a = M -> ManifoldDiffEq.CG4a(M, ExponentialRetraction())
    test_solver_frozen(manifold_to_alg_cg4a; expected_order = 4)
    compare_with_diffeq_frozen(manifold_to_alg_cg4a, constructCG4a())

    manifold_to_alg_lie_euler =
        (M, action) -> ManifoldDiffEq.ManifoldLieEuler(M, ExponentialRetraction(), action)
    test_solver_lie(manifold_to_alg_lie_euler; expected_order = 1)
    compare_with_diffeq_lie(manifold_to_alg_lie_euler, constructME())

    manifold_to_alg_rkmk4 =
        (M, action) -> ManifoldDiffEq.RKMK4(M, ExponentialRetraction(), action)
    test_solver_lie(manifold_to_alg_rkmk4; expected_order = 4)
    compare_with_diffeq_lie(manifold_to_alg_rkmk4, constructRKMK4())
end


text: run(`make`)

max_stars_repo_name: vOptSolver/vOptSpecific

max_issues_repo_name: vOptSolver/vOptSpecific

max_forks_repo_name: vOptSolver/vOptSpecific