text:
using SFML
using AnimatedPlots

t = 4*pi/4
a = 0*pi/8
td = 0.1 * rand()
ad = 0.2 * rand()
g = 9.8
delta = 1/60
torque = 0

# These functions are used for plotting
function ft(x)
	td/2
end

function fa(x)
	ad/2
end

# Create the SFML render window
settings = ContextSettings()
settings.antialiasing_level = 3
window = RenderWindow(VideoMode(800, 600), "Double Pendulum", settings, window_defaultstyle)

# Create the animated graphs
t_graph = AnimatedGraph(ft)
t_graph.accuracy = 3
a_graph = AnimatedGraph(fa, color=SFML.blue)
a_graph.accuracy = 3

# Create the plot window using the render window we created earlier
plotwindow = create_window(window)
# Add the graphs
add_graph(plotwindow, t_graph)
add_graph(plotwindow, a_graph)

# Follow one of them (they will go at the same speed)
follow(t_graph)

set_vsync_enabled(window, true)

# Get the views (these are used for drawing
view = get_default_view(window)
plotview = plotwindow.view

# You can change the viewports if you want
# set_viewport(view, FloatRect(0, 0, 1.0, 1.0))
# set_viewport(plotwindow.view, FloatRect(0.0, 0.75, 1.0, 0.25))

set_framerate_limit(window, 60)
event = Event()
circles = CircleShape[]

for i = 1:3
	circle = CircleShape()
	set_radius(circle, 10)
	set_fillcolor(circle, SFML.red)
	set_origin(circle, Vector2f(10, 10))
	push!(circles, circle)
end

rectangles = RectangleShape[]
for i = 1:2
	rect = RectangleShape()
	set_size(rect, Vector2f(10, 100))
	set_fillcolor(rect, SFML.blue)
	set_origin(rect, Vector2f(5, 0))
	push!(rectangles, rect)
end
set_position(rectangles[1], Vector2f(400, 300))

# Used for creating gifs
# clock = Clock()
# gif_made = false

while isopen(window)
	sleep(0)
	while pollevent(window, event)
		if get_type(event) == EventType.CLOSED
			close(window)
		end
		if get_type(event) == EventType.KEY_PRESSED
			if get_key(event).key_code == KeyCode.LEFT
				torque = -1
			elseif get_key(event).key_code == KeyCode.RIGHT
				torque = 1
			end
		else
			torque = 0
		end
	end
	cycles = 2000
	for i = 1:cycles
		tdd = (1 / (2 - cos(a)*cos(a))) * (td*td*sin(a)*(1+cos(a)) + (2*td+ad)*ad*sin(a) + g*(cos(a)*sin(t+a)-2*sin(t)) - torque*(1 + cos(a)))
		add = -tdd - tdd*cos(a) - td*td*sin(a) - g*sin(t+a) + torque
		td += delta * tdd / cycles
		ad += delta * add / cycles
		t += delta * td / cycles
		a += delta * ad / cycles
	end

	# Create a gif after 10 seconds (this will take awhile)
	# if as_seconds(get_elapsed_time(clock)) >= 10 && !gif_made
	# 	gif_made = true
	# 	println("Making gif")
	# 	make_gif(plotwindow, 600, 450, 2, "double_pendulum.gif", 0.04)
	# end

	x1 = 400 + 100 * sin(t)
	y1 = 300 + 100 * cos(t)
	x2 = 400 + 100 * (sin(t) + sin(a + t))
	y2 = 300 + 100 * (cos(t) + cos(a + t))
	set_position(circles[3], Vector2f(400, 300))
	set_position(circles[1], Vector2f(x1, y1))
	set_position(circles[2], Vector2f(x2, y2))

	set_rotation(rectangles[1], -180*t/pi )
	set_rotation(rectangles[2], -180*(t+a)/pi )
	set_position(rectangles[2], Vector2f(x1, y1))

	# Calculate the energy
	# T = td*td + 0.5*(td+ad)^2 + td*(td+ad)*cos(a)
	# V = -g*(2*cos(t)+cos(t+a))
	# E = T+V
	# println(E)

	clear(window, SFML.white)
	# Set the view for drawing the double pendulum
	set_view(window, view)
	for i = 1:length(rectangles)
		draw(window, rectangles[i])
	end
	for i = 1:length(circles)
		draw(window, circles[i])
	end
	# Set the view for drawing the plots
	set_view(window, plotview)
	# Update the plotwindow
	redraw(plotwindow)
	# Draw the plots
	draw(plotwindow)
	display(window)
end


text:
abstract type AeroAcousticType end

# TODO: Check that length of arr and fc match!
struct FreqArray{T, N, AA<:AbstractVector} <: AbstractArray{T, N}
    arr::Array{T, N}
    fc::AA
end

# Standard array functions and indexing
Base.size(A::FreqArray) = size(A.arr)
Base.axes(A::FreqArray) = axes(A.arr)
Base.getindex(A::FreqArray, i::Int) = A.arr[i]
Base.getindex(A::FreqArray{T, N}, I::Vararg{Int, N}) where {T, N} = A.arr[I...]
Base.setindex!(A::FreqArray,v,i) = A.arr[i] = v

"""
    Environment

The Environment struct is required for most methods and defines the geometrical
setup, constants, and stores the relevant data together. The microphone array is assumed to be at the center of the coordinate system. The fields are

# Arguments

*Required:*

- `micgeom::Array{T,2}`: (x,y,z) coordinates for microhone array.
- `z0::Real`: distance to source plane from microphone array.
- `CSM::FreqArray`: Cross-spectral matrix, size `M x M x Nf`, where `M` is the number of microphones and `Nf` is the number of frequency bins.

*Optional:*

- `flim::Tuple=extrema(CSM.fc)`: Frequency limits (fmin,fmax).
- `Nx::Integer=21`: Number of computational gridpoint in x direction.
- `Ny::Integer=21`: Number of computational gridpoint in y direction.
- `xlim::Tuple=(-1.,1.)`: Cartesian x-coordinate limits.
- `ylim::Tuple=(-1.,1.)`: Cartesian y-coordinate limits.
- `wv=ones(size(micgeom,2))`: Microphones on/off.
- `wc::Bool=false`: Microphones coherence weighting.
- `dr::Bool=false`: Diagonal removal of CSM
- `w=ones(M,Nf)`: Microphone weights. If `wc=true` it calculates coherence weighting.
- `shear::Bool = false`: Amiet phase correction
- `ampcorr::Bool = shear`: Amiet amplitude correction (only applies when shear = true)
- `c::Real=343.`: Speed of sound [m/s].
- `Ma::Real=0.0`: Mach number (sign determines flow direction)
- `h::Real=0.0`: Distance from array center to shear layer (Amiet correction) should be supplied when `Ma != 0`.

"""
@with_kw mutable struct Environment <: AeroAcousticType
    micgeom::Matrix{<:AbstractFloat}
    z0::Real
    CSM::FreqArray
    flim::NTuple{2,Real} = extrema(CSM.fc)
    Nx::Int = 21
    Ny::Int = 21
    xlim::NTuple{2,Real} = (-1.,1.)
    ylim::NTuple{2,Real} = (-1.,1.)
    shear::Bool = false
    ampcorr::Bool = shear
    c::Real = 343.0 # Speed of sound
    Ma::Real = 0.0 # Flow Mach speed (in positive x-direction) TODO: Generalize to NTuple{3,Real}
    h::Real = 0.0 # Distance from array center to shear layer
    ### Compute extra parameters
    dr::Bool = false
    Cinds = (CSM.fc.>=flim[1]) .& (CSM.fc.<=flim[2])
    fn = CSM.fc[Cinds]
    wv::Vector = ones(size(micgeom,2)) # weights vector
    micgeom_s = micgeom[:,Bool.(wv)]
    CSM_w = FreqArray(CSM[Bool.(wv),Bool.(wv),Cinds],fn)
    CSM_s = ifelse(dr,DR(CSM_w),CSM_w)
    M::Int = size(micgeom_s,2) #count((!iszero).(wv))
    Nf::Int = length(fn)
    wc::Bool = false
    w = ifelse(wc,coherence_weights(CSM_w),ones(M,Nf))
    N = Nx*Ny
    rx = range(xlim[1],stop = xlim[2], length = Nx)
    ry = range(ylim[1],stop = ylim[2], length = Ny)
    Rxy = hcat([[x, y, z] for x in rx, y in ry, z in z0]...)
    D0 = colwise(Euclidean(), Rxy, [0,0,0]) # Distance from center of array to grid points
    D = pairwise(Euclidean(), Rxy, micgeom_s; dims=2) # Distance from each mic to grid points
    steeringvec = nothing
end

ext: jl


text:
push!(LOAD_PATH, ".")

module IntegerIndividual

include("utility.jl")

using Individual
using Parallel

abstract type AbstractIntegerIndividual{GeneType <: Integer} <: Individual.AbstractIndividual{GeneType} end

@define IntegerIndividual begin
  IntegerIndividual.Individual.@Individual
end

struct IntegerIndividualType{GeneType} <: AbstractIntegerIndividual{GeneType}
  @IntegerIndividual

  function IntegerIndividualType{GeneType}(x...) where {GeneType}
    args = build(x...)
    new(args[1], args[2], args[3])
  end
end

function build(x...)
  return Individual.build(x...)
end

function Individual.generateRandom!(this::AbstractIntegerIndividual)
  for i in eachindex(this)
    gene_lb, gene_ub = this.bounds[i]
    this[i] = Parallel.threadRand(gene_lb : gene_ub)
  end
end

end

max_stars_repo_name: Markhyz/Julia-Genetic-Algorithm

max_issues_repo_name: Markhyz/Julia-Genetic-Algorithm

max_forks_repo_name: Markhyz/Julia-Genetic-Algorithm


text:
type CVSystem # entire solution
    heart::CVModule.Heart
    branches::CVModule.ArterialBranches
    svc::CVModule.VenaCava
    ivc::CVModule.VenaCava
    lungs::CVModule.Lungs
    liver::CVModule.Liver
    hemo::CVModule.Hemorrhage
    cns::CVModule.CNS
    solverparams::CVModule.SolverParams
    t::Vector{Float64}
    arterialvolume::Float64
    peripheralvolume::Float64
    vcvolume::Float64
    heartvolume::Float64
    lungvolume::Float64
    livervolume::Float64
    initialvolume::Float64
    finalvolume::Float64

    function CVSystem(filename="test.csv",restart="no")
        this = new()
        if restart == "no"
            this.heart = CVModule.Heart();
            this.branches = CVModule.ArterialBranches(filename);
        elseif restart == "yes"
            vars = MAT.matread(filename);
            sys = vars["system"];
            heart = sys["heart"];
            branches = sys["branches"];
            this.heart = CVModule.Heart(heart,restart)
            this.branches = CVModule.ArterialBranches(filename,branches,restart)
        else
            error("Keyword restart must either be yes or no. Aborting.")
        end
        this.svc = CVModule.VenaCava();
        this.ivc = CVModule.VenaCava();
        this.lungs = CVModule.Lungs();
        this.liver = CVModule.Liver();
        this.cns = CVModule.CNS();
        this.hemo = CVModule.Hemorrhage();
        this.solverparams = CVModule.SolverParams();
        this.t = Vector{Float64}[];
        return this
    end
end

# build solution struct
function buildall(filename="test.csv";numbeatstotal=1,restart="no",injury="no")
    if restart == "no"
        system = CVModule.CVSystem(filename);
        system.solverparams.numbeatstotal = numbeatstotal;
        CVModule.calcbranchprops!(system);
        CVModule.discretizebranches!(system);
        CVModule.assignterminals!(system);
        CVModule.discretizeperiphery!(system);
        CVModule.discretizeheart!(system);
        CVModule.discretizelungs!(system);
        CVModule.discretizeliver!(system);
        CVModule.discretizecns!(system);
        CVModule.applybranchics!(system);
        CVModule.applyperipheryics!(system);
        CVModule.applyheartics!(system);
        CVModule.applylungics!(system);
        CVModule.applyliverics!(system);
        CVModule.applycnsics!(system);
        CVModule.applycustomics!(system);
        CVModule.applyhemoics!(system);
    elseif restart == "yes"
        vars = MAT.matread(filename);
        sys = vars["system"];
        branches = sys["branches"];
        term = branches["term"];
        heart = sys["heart"];
        lungs = sys["lungs"];
        liver = sys["liver"];
        cns = sys["cns"];
        system = CVModule.CVSystem(filename,restart);
        system.solverparams.numbeatstotal = numbeatstotal;
        CVModule.calcbranchprops!(system,branches,restart);
        CVModule.discretizebranches!(system,sys,restart);
        CVModule.assignterminals!(system,term,restart);
        CVModule.discretizeperiphery!(system);
        CVModule.discretizeheart!(system);
        CVModule.discretizelungs!(system);
        CVModule.discretizeliver!(system);
        CVModule.discretizecns!(system);
        CVModule.applybranchics!(system,sys,restart);
        CVModule.applyperipheryics!(system,sys,restart);
        CVModule.applyheartics!(system,heart,restart);
        CVModule.applylungics!(system,lungs,restart);
        CVModule.applyliverics!(system,liver,restart);
        CVModule.applycnsics!(system,cns,restart);
        CVModule.applyhemoics!(system,sys);
    end
    CVModule.updatevolumes!(system,0);
    return system
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function test_spirals(fname)
    Drawing(2000, 2000, fname)
    origin()
    setopacity(0.5)
    setline(2)
    for i in 0:pi/6:2pi
        rotate(i)
        randomhue()
        v1 = spiral(10, log(Base.MathConstants.golden)/(pi/2), log=true, stepby=pi/5, period=4pi, vertices=true)
        poly(v1, :fill)
        for j in 1:length(v1) - 3
            poly([v1[j], v1[j+1], v1[j+2]], :fill)
        end
    end
    for i in 0:pi/6:2pi
        rotate(i)
        v2 = spiral(30, 1, log=false, stepby=pi/10, period=4pi, vertices=true)
        poly(v2, :stroke)
        rand(Bool) ? sethue("white") : sethue("black")
        for j in 1:length(v2) - 3
            poly([v2[j], v2[j+1], v2[j+2]], :fill)
        end
    end
    @test finish() == true
end

fname = "spirals.pdf"
test_spirals(fname)
println("...finished test: output in $(fname)")


text:
module Day8

export get_inputs, get_solution1, get_solution2

using ..Machines: read_instructions, run_tape, Machine

## Input getting
function get_inputs()
    test_input1 = test_input2 = read_instructions(IOBuffer("""
        nop +0
        acc +1
        jmp +4
        acc +3
        jmp -3
        acc -99
        acc +1
        jmp -4
        acc +6
    """))
    test_output1 = 5
    test_output2 = 8
    data = read_instructions(joinpath(@__DIR__, "input.txt"))
    return (; test_input1, test_input2, test_output1, test_output2, data)
end

## Solution functions
# Part 1
infinite_loop(machine) = machine.current_line in machine.lines_visited

get_solution1(tape) = run_tape(tape, Machine(track=true); halt_condition=infinite_loop)[1].accumulator

# Part 2
function run_until_valid(tape)
    new_tape = copy(tape)

    for i in eachindex(tape)
        f, args = new_tape[i]

        if f == :nop
            new_tape[i] = (:jmp, args)
        elseif f == :jmp
            new_tape[i] = (:nop, args)
        end

        machine, halted = run_tape(new_tape, Machine(track=true); halt_condition=infinite_loop)
        new_tape[i] = tape[i]

        if !halted
            return machine, i
        end
    end

    return Machine(), -1
end

get_solution2(tape) = run_until_valid(tape)[1].accumulator

end


text:
abstract type LateralDriverModel end
get_name(::LateralDriverModel) = "???"
reset_hidden_state!(model::LateralDriverModel) = model # do nothing by default
observe!(model::LateralDriverModel, scene::Scene, roadway::Roadway, egoid::Int) = model  # do nothing by default
Base.rand(model::LateralDriverModel) = error("rand not implemented for model $model")
Distributions.pdf(model::LateralDriverModel, a_lon::Float64) = error("pdf not implemented for model $model")
Distributions.logpdf(model::LateralDriverModel, a_lon::Float64) = error("logpdf not implemented for model $model")

# TODO: Why does this struct live in the `lateral_driver_models.jl` file
"""
	ProportionalLaneTracker

A controller that executes the lane change decision made by the `lane change models`

# Constructors
	ProportionalLaneTracker(;σ::Float64 = NaN,kp::Float64 = 3.0,kd::Float64 = 2.0)

# Fields
- `a::Float64 = NaN` predicted acceleration
- `σ::Float64 = NaN` optional stdev on top of the model, set to zero or NaN for deterministic behavior
- `kp::Float64 = 3.0` proportional constant for lane tracking
- `kd::Float64 = 2.0` derivative constant for lane tracking
"""
mutable struct ProportionalLaneTracker <: LateralDriverModel
    a::Float64 # predicted acceleration
    σ::Float64 # optional stdev on top of the model, set to zero or NaN for deterministic behavior
    kp::Float64 # proportional constant for lane tracking
    kd::Float64 # derivative constant for lane tracking

    function ProportionalLaneTracker(;
        σ::Float64 = NaN,
        kp::Float64 = 3.0,
        kd::Float64 = 2.0,
        )

        retval = new()
        retval.a = NaN
        retval.σ = σ
        retval.kp = kp
        retval.kd = kd
        retval
    end
end
get_name(::ProportionalLaneTracker) = "ProportionalLaneTracker"
function track_lateral!(model::ProportionalLaneTracker, laneoffset::Float64, lateral_speed::Float64)
    model.a = -laneoffset*model.kp - lateral_speed*model.kd
    model
end
function observe!(model::ProportionalLaneTracker, scene::Scene, roadway::Roadway, egoid::Int)

    ego_index = findfirst(egoid, scene)
    veh_ego = scene[ego_index]
    t = veh_ego.state.posF.t # lane offset
    dt = veh_ego.state.v * sin(veh_ego.state.posF.ϕ) # rate of change of lane offset
    model.a = -t*model.kp - dt*model.kd

    model
end
function Base.rand(model::ProportionalLaneTracker)
    if isnan(model.σ) || model.σ ≤ 0.0
        model.a
    else
        rand(Normal(model.a, model.σ))
    end
end
function Distributions.pdf(model::ProportionalLaneTracker, a_lat::Float64)
    if isnan(model.σ) || model.σ ≤ 0.0
        Inf
    else
        pdf(Normal(model.a, model.σ), a_lat)
    end
end
function Distributions.logpdf(model::ProportionalLaneTracker, a_lat::Float64)
    if isnan(model.σ) || model.σ ≤ 0.0
        Inf
    else
        logpdf(Normal(model.a, model.σ), a_lat)
    end
end


text:
## export to JSON

"""
    toDict(s::FittedShape)

Convert `s` to a `Dict{String,Any}`. It's "type" is defined by [`strt`](@ref).
All fields of the struct is saved to the dict.
"""
function toDict(s::FittedShape)
    d = Dict{String,Any}("type"=>strt(s))
    for f in fieldnames(typeof(s))
        sf = String(f)
        merge!(d, Dict(sf=>getfield(s, f)))
    end
    d
end

toDict(sc::ExtractedShape) = toDict(sc.shape)

"""
    toDict(a::Vector{T}) where {T<:Union{FittedShape,ExtractedShape}}

Convert a vector of shapes to a `Dict{String,Any}`.
The top key is a "primitive", whose value is the array of the shapes.
See the documentation for examples.
"""
function toDict(a::Vector{T}) where {T<:Union{FittedShape,ExtractedShape}}
    ad = toDict.(a)
    return Dict("primitives"=>ad)
end

"""
    printJSON(io::IO, s, indent)

Print a `FittedShape`, `ExtractedShape` or a vector of them to `io` as a JSON string.
With `indent` given, it prints a representation with newlines and indents.

# Arguments:
- `io::IO`: must be specified, use `stdout` for interactive purposes.
- `s`: a `FittedShape`, `ExtractedShape` or a vector of one of them.
- `indent::Int`: indentation level.
"""
function exportJSON(io::IO, s, indent)
    JSON.print(io, toDict(s), indent)
end

"""
    printJSON(io::IO, s)

Print a `FittedShape`, `ExtractedShape`
or a vector of them to `io` as a compact JSON string.

# Arguments:
- `io::IO`: must be specified, use `stdout` for interactive purposes.
- `s`: a `FittedShape`, `ExtractedShape` or a vector of one of them.
"""
function exportJSON(io::IO, s)
    JSON.print(io, toDict(s))
end

## import yaml config file

"""
    readconfig(fname; toextend=DEFAULT_PARAMETERS, shapedict=DEFAULT_SHAPE_DICT)

Read a config file to a `NamedTuple`.
A "base" ntuple is expected, that gets overwritten/extended
with the values in the config file.
A `Dict{String,FittedShape}` dictionary is also expected,
that translates the string primitive types to julia types.

# Arguments

- `fname`: name of the config file.
- `toextend=DEFAULT_PARAMETERS`: a named tuple,
    that will be overwritten/extended with the values from the config file.
- `shapedict=DEFAULT_SHAPE_DICT`: a dictionary that
    translates the string primitive names to julia types.
"""
function readconfig(fname; toextend=DEFAULT_PARAMETERS, shapedict=DEFAULT_SHAPE_DICT)
    fio = open(fname, "r")
    fdict = YAML.load(fio)
    close(fio)

    newnd = ((;Symbol(k)=>dict2nt(v)) for (k,v) in pairs(fdict))
    for v in newnd
        if haskey(v[1], :shape_types)
            nv = merge(v[1], (shape_types=[shapedict[k] for k in v[1].shape_types],))
            toextend = ransacparameters(toextend; iteration=nv)
        else
            toextend = ransacparameters(toextend; v...)
        end
    end
    return toextend
end

function dict2nt(dict)
    dd = [NamedTuple{Tuple(Symbol.(keys(d)))}(values(d)) for d in dict]
    if length(dd) < 2
        return dd[1]
    end
    newd = dd[1]
    for i in eachindex(dd)
        i == 1 && continue
        newd = merge(newd, dd[i])
    end
    return newd
end


text:
require("ZMQ")
using ZMQ

println("Testing with ZMQ version $(ZMQ.version)")

ctx=Context(1)

@assert typeof(ctx) == Context

ZMQ.close(ctx)

#try to create socket with expired context
try
	Socket(ctx, PUB)
	@assert false
catch ex
	@assert typeof(ex) == StateError
end

ctx2=Context(1)
s=Socket(ctx2, PUB)
@assert typeof(s) == Socket
ZMQ.close(s)

#trying to close already closed socket
try
	ZMQ.close(s)
catch ex
	@assert typeof(ex) == StateError
end

s1=Socket(ctx2, REP)
if ZMQ.version.major == 2
	ZMQ.set_hwm(s1, 1000)
else
	ZMQ.set_sndhwm(s1, 1000)
end
ZMQ.set_linger(s1, 1)
ZMQ.set_identity(s1, "abcd")

@assert ZMQ.get_identity(s1)::String == "abcd"
if ZMQ.version.major == 2
	@assert ZMQ.get_hwm(s1)::Integer == 1000
else
	@assert ZMQ.get_sndhwm(s1)::Integer == 1000
end
@assert ZMQ.get_linger(s1)::Integer == 1
@assert ZMQ.ismore(s1) == false

s2=Socket(ctx2, REQ)
@assert ZMQ.get_type(s1) == REP
@assert ZMQ.get_type(s2) == REQ

ZMQ.bind(s1, "tcp://*:5555")
ZMQ.connect(s2, "tcp://localhost:5555")

ZMQ.send(s2, Message("test request"))
@assert (bytestring(ZMQ.recv(s1)) == "test request")
ZMQ.send(s1, Message("test response"))
@assert (bytestring(ZMQ.recv(s2)) == "test response")

# Test task-blocking behavior
c = Base.Condition()
msg_sent = false
@async begin
	global msg_sent
	sleep(0.5)
	msg_sent = true
	ZMQ.send(s2, Message("test request"))
	@assert (bytestring(ZMQ.recv(s2)) == "test response")
	notify(c)
end

# This will hang forver if ZMQ blocks the entire process since
# we'll never switch to the other task
@assert (bytestring(ZMQ.recv(s1)) == "test request")
@assert msg_sent == true
ZMQ.send(s1, Message("test response"))
wait(c)

ZMQ.send(s2, Message("another test request"))
msg = ZMQ.recv(s1)
o=convert(IOStream, msg)
seek(o, 0)
@assert (takebuf_string(o)=="another test request")

ZMQ.close(s1)
ZMQ.close(s2)
ZMQ.close(ctx2)


text:
# Implementation of the AbstractObjective interface
# (see file `AbstractObjectiveInterface.jl`)
@with_kw mutable struct VectorObjectiveFunction{
        SC <: SurrogateConfig,
        FT <: Function
    } <: AbstractObjective
    n_in :: Int = 0;
    n_out :: Int = 0;
    n_evals :: Int64 = 0;   # true function evaluations (also counts fdm evaluations)

    # max_evals :: Int64 = typemax(Int64);

    model_config :: SC

    function_handle :: FT
end

# Required methods (see file `AbstractObjectiveInterface.jl`)
function _wrap_func( ::Type{<:VectorObjectiveFunction},
        fn::Function, model_cfg::SurrogateConfig,
        n_vars :: Int, n_out :: Int
    ) :: VectorObjectiveFunction

    return VectorObjectiveFunction(;
        n_in = n_vars,
        n_out = n_out,
        function_handle = fn,
        model_config = model_cfg
    );
end

num_vars( objf :: VectorObjectiveFunction ) = objf.n_in;

num_evals( objf :: VectorObjectiveFunction ) = objf.n_evals;
function num_evals!( objf :: VectorObjectiveFunction, N :: Int )
    objf.n_evals = N;
    nothing
end

num_outputs( objf :: VectorObjectiveFunction ) = objf.n_out;

# can_batch( objf :: VectorObjectiveFunction ) = isa( objf.function_handle, BatchObjectiveFunction );

model_cfg( objf :: VectorObjectiveFunction ) = objf.model_config;

combinable(::Type{<:VectorObjectiveFunction}, ::Type{<:VectorObjectiveFunction}) = true
eval_handle(objf::VectorObjectiveFunction) = objf.function_handle;


text:
using PowerSystems
using JuMP
using PowerSimulations
using Ipopt
const PS = PowerSimulations

# ED Testing
base_dir = dirname(dirname(pathof(PowerSystems)))
ps_dict = PowerSystems.parsestandardfiles(abspath(joinpath(base_dir, "data/matpower/case5_dc.m")));
sys5 = PowerSystem(ps_dict);

solver = Ipopt.Optimizer

@test try

    ED = PS.PowerOperationModel(PS.EconomicDispatch,
        [(device = ThermalGen, formulation =PS.ThermalDispatch)],
        nothing,
        nothing,
        [(device=Line, formulation=PS.PiLine)],
        PS.CopperPlatePowerModel,
        nothing,
        sys5,
        Model(),
        false,
        nothing)

    PS.buildmodel!(sys5,ED)

    ED_sim = PS.buildsimulation!(sys5, ED)

    PS.run_simulations(ED_sim, solver, ps_dict)

true finally end


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

JLLWrappers.@generate_main_file_header("FreeType2")
JLLWrappers.@generate_main_file("FreeType2", UUID("d7e528f0-a631-5988-bf34-fe36492bcfd7"))
end  # module FreeType2_jll


text:
struct SDEDiffusionTermWrapper{iip,G,H}
  g::G
  h::H
end
function SDEDiffusionTermWrapper(g::G,h::H) where {G,H}
  SDEDiffusionTermWrapper{isinplace(g,5),G,H}(g,h)
end
(g::SDEDiffusionTermWrapper{true})(du, u, p, t) = g.g(du, u, g.h, p, t)
(g::SDEDiffusionTermWrapper{false})(u, p, t) = g.g(u, g.h, p, t)

struct SDEFunctionWrapper{iip,F,G,H,TMM,Ta,Tt,TJ,JVP,VJP,JP,SP,TW,TWt,TPJ,GG,S,TCV} <: DiffEqBase.AbstractRODEFunction{iip}
  f::F
  g::G
  h::H
  mass_matrix::TMM
  analytic::Ta
  tgrad::Tt
  jac::TJ
  jvp::JVP
  vjp::VJP
  jac_prototype::JP
  sparsity::SP
  Wfact::TW
  Wfact_t::TWt
  paramjac::TPJ
  ggprime::GG
  syms::S
  colorvec::TCV
end

(f::SDEFunctionWrapper{true})(du, u, p, t) = f.f(du, u, f.h, p, t)
(f::SDEFunctionWrapper{false})(u, p, t) = f.f(u, f.h, p, t)

function wrap_functions_and_history(f::SDDEFunction, g, h)
  gwh = SDEDiffusionTermWrapper{isinplace(g,5),typeof(g),typeof(h)}(g,h)

  if f.jac === nothing
    jac = nothing
  else
    if isinplace(f)
      jac = let f_jac = f.jac, h = h
        (J, u, p, t) -> f_jac(J, u, h, p, t)
      end
    else
      jac = let f_jac = f.jac, h = h
        (u, p, t) -> f_jac(u, h, p, t)
      end
    end
  end

  SDEFunctionWrapper{isinplace(f),typeof(f.f),typeof(gwh),typeof(h),typeof(f.mass_matrix),
                     typeof(f.analytic),typeof(f.tgrad),typeof(jac),typeof(f.jvp),
                     typeof(f.vjp),typeof(f.jac_prototype),typeof(f.sparsity),
                     typeof(f.Wfact),typeof(f.Wfact_t),typeof(f.paramjac),
                     typeof(f.ggprime),typeof(f.syms),typeof(f.colorvec)}(
                       f.f, gwh, h, f.mass_matrix, f.analytic, f.tgrad, jac,
                       f.jvp, f.vjp,
                       f.jac_prototype, f.sparsity, f.Wfact, f.Wfact_t,
                       f.paramjac, f.ggprime, f.syms, f.colorvec), gwh
end


text:
"""
"""
function FETerm(
  res::Function, trian::Triangulation, quad::CellQuadrature)
  NonlinearFETermWithAutodiff(res,trian,quad)
end

struct NonlinearFETermWithAutodiff <: FETerm
  res::Function
  trian::Triangulation
  quad::CellQuadrature
end

function get_cell_residual(t::NonlinearFETermWithAutodiff,uh,v)
  @assert is_a_fe_function(uh)
  @assert is_a_fe_cell_basis(v)
  _v = restrict(v,t.trian)
  _uh = restrict(uh,t.trian)
  integrate(t.res(_uh,_v),t.trian,t.quad)
end

function get_cell_jacobian(t::NonlinearFETermWithAutodiff,uh,du,v)
  @assert is_a_fe_function(uh)
  @assert is_a_fe_cell_basis(v)
  @assert is_a_fe_cell_basis(du)
  _v = restrict(v,t.trian)
  function uh_to_cell_residual(uh)
    _uh = restrict(uh,t.trian)
    integrate(t.res(_uh,_v),t.trian,t.quad)
  end
  cell_j = autodiff_cell_jacobian_from_residual(uh_to_cell_residual,uh,get_cell_id(t))
  cell_j
end

function get_cell_id(t::NonlinearFETermWithAutodiff)
  get_cell_id(t.trian)
end

struct FEEnergy <: FETerm
  ener::Function
  trian::Triangulation
  quad::CellQuadrature
end

function get_cell_residual(t::FEEnergy,uh,v)
  @assert is_a_fe_function(uh)
  @assert is_a_fe_cell_basis(v)
  function uh_to_cell_energy(uh)
    _uh = restrict(uh,t.trian)
    integrate(t.ener(_uh),t.trian,t.quad)
  end
  cell_r = autodiff_cell_residual_from_energy(uh_to_cell_energy,uh,get_cell_id(t))
  cell_r
end

function get_cell_jacobian(t::FEEnergy,uh,du,v)
  @assert is_a_fe_function(uh)
  @assert is_a_fe_cell_basis(v)
  @assert is_a_fe_cell_basis(du)
  function uh_to_cell_energy(uh)
    _uh = restrict(uh,t.trian)
    integrate(t.ener(_uh),t.trian,t.quad)
  end
  cell_j = autodiff_cell_jacobian_from_energy(uh_to_cell_energy,uh,get_cell_id(t))
  cell_j
end

function get_cell_id(t::FEEnergy)
  get_cell_id(t.trian)
end


text:
using Test
using System

# NOTE: tests require certain headers and libraries to be installed
#       reference `System.jl/.github/workflows/CI.yml` for commands

@testset "System.jl" begin
	@test @eval begin
		@sys using sys
		true
	end

	@test @eval begin
		@sys using libc
		true
	end

	if Sys.islinux()
		@test @eval begin
			@sys using linux
			true
		end

		@test @eval begin
			@sys using alsa
			true
		end
	end
end


text:
@testset "StoppingCondition" begin

    @testset "Defaults" begin
        p = DSProblem(4)
        @test typeof(p.stoppingconditions[1]) == DS.IterationStoppingCondition
        @test typeof(p.stoppingconditions[2]) == DS.MeshPrecisionStoppingCondition
        @test p.stoppingconditions[1].limit == 1000

        p = DSProblem(4; iteration_limit=100)
        @test p.stoppingconditions[1].limit == 100
    end

    @testset "AddStoppingCondition" begin
        p = DSProblem(4)

        test_stopping_cond = DS.IterationStoppingCondition(1234)
        DS.AddStoppingCondition(p, test_stopping_cond)
        @test p.stoppingconditions[3].limit == 1234
    end

    @testset "_check_stoppingconditions" begin
        p = DSProblem(4)
        c = DS.AbstractStoppingCondition[]
        @test DS._check_stoppingconditions(p, c) == true

        p.status.iteration = 100
        push!(c, DS.IterationStoppingCondition(200))
        @test DS._check_stoppingconditions(p, c) == true

        p.status.iteration = 200
        @test DS._check_stoppingconditions(p, c) == false

        p.status.iteration = 300
        @test DS._check_stoppingconditions(p, c) == false

        push!(c, DS.IterationStoppingCondition(400))
        @test DS._check_stoppingconditions(p, c) == false
    end

    @testset "setstatus" begin
        p = DSProblem(4)

        @test p.status.optimization_status == DS.Unoptimized
        @test p.status.optimization_status_string == "Unoptimized"

        struct test_sc <: DS.AbstractStoppingCondition end
        DS.StoppingConditionStatus(::test_sc) = "test status message"
        DS.CheckStoppingCondition(p::DSProblem, s::test_sc) = false

        DS.setstatus(p, test_sc())
        @test p.status.optimization_status_string == "test status message"

        #First false sc reached in the array is the reported status
        DS.AddStoppingCondition(p, test_sc())
        @test DS._check_stoppingconditions(p) == false
        @test p.status.optimization_status_string == "test status message"
        p.status.iteration = 1001
        @test DS._check_stoppingconditions(p) == false
        @test p.status.optimization_status_string == "Iteration limit"
    end

    @testset "StoppingConditionStatus" begin
        struct another_test_sc <: DS.AbstractStoppingCondition end
        DS.CheckStoppingCondition(p::DSProblem, s::another_test_sc) = false

        @test DS.StoppingConditionStatus(another_test_sc()) == "Unknown stopping condition status"
        p = DSProblem(4)
        DS.AddStoppingCondition(p, another_test_sc())
        DS._check_stoppingconditions(p)
        @test p.status.optimization_status == DS.OtherStoppingCondition
        @test p.status.optimization_status_string == "Unknown stopping condition status"
    end

end


text:
# The DH17 integrator WITHOUT derivatives.
"""
Carries out the DH17 mapping with compensated summation.
"""
function dh17!(x::Array{T,2},v::Array{T,2},xerror::Array{T,2},verror::Array{T,2},h::T,m::Array{T,1},n::Int64,pair::Array{Bool,2}) where {T <: Real}
    alpha = convert(T,alpha0)
    h2 = 0.5*h
    # alpha = 0. is similar in precision to alpha=0.25
    kickfast!(x,v,xerror,verror,h/6,m,n,pair)
    if alpha != 0.0
        phisalpha!(x,v,xerror,verror,h,m,alpha,n,pair)
    end
    #mbig = big.(m); h2big = big(h2)
    #xbig = big.(x); vbig = big.(v)
    drift!(x,v,xerror,verror,h2,n)
    #drift!(xbig,vbig,h2big,n)
    #xcompare = convert(Array{T,2},xbig); vcompare = convert(Array{T,2},vbig)
    #x .= xcompare; v .= vcompare
    #hbig = big(h)
    @inbounds for i=1:n-1
        @inbounds for j=i+1:n
            if ~pair[i,j]
                #      xbig = big.(x); vbig = big.(v)
                driftij!(x,v,xerror,verror,i,j,-h2)
                #      driftij!(xbig,vbig,i,j,-h2big)
                #      xcompare = convert(Array{T,2},xbig); vcompare = convert(Array{T,2},vbig)
                #      x .= xcompare; v .= vcompare
                #      xbig = big.(x); vbig = big.(v)
                #      keplerij!(mbig,xbig,vbig,i,j,h2big)
                #      xcompare = convert(Array{T,2},xbig); vcompare = convert(Array{T,2},vbig)
                keplerij!(m,x,v,xerror,verror,i,j,h2)
                #      x .= xcompare; v .= vcompare
            end
        end
    end
    # kick and correction for pairs which are kicked:
    phic!(x,v,xerror,verror,h,m,n,pair)
    if alpha != 1.0
        #  xbig = big.(x); vbig = big.(v)
        phisalpha!(x,v,xerror,verror,h,m,2.0*(1.0-alpha),n,pair)
        #  phisalpha!(xbig,vbig,hbig,mbig,2*(1-big(alpha)),n,pair)
        #  xcompare = convert(Array{T,2},xbig); vcompare = convert(Array{T,2},vbig)
        #  x .= xcompare; v .= vcompare
    end
    @inbounds for i=n-1:-1:1
        @inbounds for j=n:-1:i+1
            if ~pair[i,j]
                #      xbig = big.(x); vbig = big.(v)
                #      keplerij!(mbig,xbig,vbig,i,j,h2big)
                #x = convert(Array{T,2},xbig); v = convert(Array{T,2},vbig)
                #      xcompare = convert(Array{T,2},xbig); vcompare = convert(Array{T,2},vbig)
                keplerij!(m,x,v,xerror,verror,i,j,h2)
                #      x .= xcompare; v .= vcompare
                #      xbig = big.(x); vbig = big.(v)
                driftij!(x,v,xerror,verror,i,j,-h2)
                #      driftij!(xbig,vbig,i,j,-h2big)
                #      xcompare = convert(Array{T,2},xbig); vcompare = convert(Array{T,2},vbig)
                #      x .= xcompare; v .= vcompare
            end
        end
    end
    #xbig = big.(x); vbig = big.(v)
    drift!(x,v,xerror,verror,h2,n)
    #drift!(xbig,vbig,h2big,n)
    #xcompare = convert(Array{T,2},xbig); vcompare = convert(Array{T,2},vbig)
    #x .= xcompare; v .= vcompare
    if alpha != 0.0
        phisalpha!(x,v,xerror,verror,h,m,alpha,n,pair)
    end
    kickfast!(x,v,xerror,verror,h/6,m,n,pair)
    return
end

"""

Drifts bodies i & j with compensated summation:
"""
function driftij!(x::Array{T,2},v::Array{T,2},xerror::Array{T,2},verror::Array{T,2},i::Int64,j::Int64,h::T) where {T <: Real}
    for k=1:NDIM
        x[k,i],xerror[k,i] = comp_sum(x[k,i],xerror[k,i],h*v[k,i])
        x[k,j],xerror[k,j] = comp_sum(x[k,j],xerror[k,j],h*v[k,j])
    end
    return
end

"""

Carries out a Kepler step for bodies i & j with compensated summation:
"""
function keplerij!(m::Array{T,1},x::Array{T,2},v::Array{T,2},xerror::Array{T,2},verror::Array{T,2},i::Int64,j::Int64,h::T) where {T <: Real}
    # The state vector has: 1 time; 2-4 position; 5-7 velocity; 8 r0; 9 dr0dt; 10 beta; 11 s; 12 ds
    # Initial state:
    state0 = zeros(T,12)
    # Final state (after a step):
    state = zeros(T,12)
    delx = zeros(T,NDIM)
    delv = zeros(T,NDIM)
    #println("Masses: ",i," ",j)
    for k=1:NDIM
        state0[1+k     ] = x[k,i] - x[k,j]
        state0[1+k+NDIM] = v[k,i] - v[k,j]
    end
    gm = GNEWT*(m[i]+m[j])
    if gm == 0
        for k=1:NDIM
            #x[k,i] += h*v[k,i]
            x[k,i],xerror[k,i] = comp_sum(x[k,i],xerror[k,i],h*v[k,i])
            #x[k,j] += h*v[k,j]
            x[k,j],xerror[k,j] = comp_sum(x[k,j],xerror[k,j],h*v[k,j])
        end
    else
        # predicted value of s
        kepler_step!(gm, h, state0, state)
        for k=1:NDIM
            delx[k] = state[1+k] - state0[1+k]
            delv[k] = state[1+NDIM+k] - state0[1+NDIM+k]
        end
        # Advance center of mass:
        # Compute COM coords:
        mijinv =1.0/(m[i] + m[j])
        vcm = zeros(T,NDIM)
        for k=1:NDIM
            vcm[k] = (m[i]*v[k,i] + m[j]*v[k,j])*mijinv
        end
        centerm!(m,mijinv,x,v,xerror,verror,vcm,delx,delv,i,j,h)
    end
    return
end

"""

Advances the center of mass of a binary (any pair of bodies) with compensated summation:
"""
function centerm!(m::Array{T,1},mijinv::T,x::Array{T,2},v::Array{T,2},xerror::Array{T,2},verror::Array{T,2},vcm::Array{T,1},delx::Array{T,1},delv::Array{T,1},i::Int64,j::Int64,h::T) where {T <: Real}
    for k=1:NDIM
        #x[k,i] +=  m[j]*mijinv*delx[k] + h*vcm[k]
        x[k,i],xerror[k,i] =  comp_sum(x[k,i],xerror[k,i],m[j]*mijinv*delx[k])
        x[k,i],xerror[k,i] =  comp_sum(x[k,i],xerror[k,i],h*vcm[k])
        #x[k,j] += -m[i]*mijinv*delx[k] + h*vcm[k]
        x[k,j],xerror[k,j] = comp_sum(x[k,j],xerror[k,j],-m[i]*mijinv*delx[k])
        x[k,j],xerror[k,j] = comp_sum(x[k,j],xerror[k,j],h*vcm[k])
        #v[k,i] +=  m[j]*mijinv*delv[k]
        v[k,i],verror[k,i] =  comp_sum(v[k,i],verror[k,i],m[j]*mijinv*delv[k])
        #v[k,j] += -m[i]*mijinv*delv[k]
        v[k,j],verror[k,j] = comp_sum(v[k,j],verror[k,j],-m[i]*mijinv*delv[k])
    end
    return
end

"""

Takes a single kepler step, calling Wisdom & Hernandez solver
"""
function kepler_step!(gm::T,h::T,state0::Array{T,1},state::Array{T,1}) where {T <: Real}
    # compute beta, r0,  get x/v from state vector & call correct subroutine
    x0 = zeros(eltype(state0),3)
    v0 = zeros(eltype(state0),3)
    for k=1:3
        x0[k]=state0[k+1]
        v0[k]=state0[k+4]
    end
    #  x0=state0[2:4]
    r0 = sqrt(x0[1]*x0[1]+x0[2]*x0[2]+x0[3]*x0[3])
    #  v0 = state0[5:7]
    beta0 = 2*gm/r0-(v0[1]*v0[1]+v0[2]*v0[2]+v0[3]*v0[3])
    s0=state0[11]
    iter = kep_ell_hyp!(x0,v0,r0,gm,h,beta0,s0,state)
    #  if beta0 > zero
    #    iter = kep_elliptic!(x0,v0,r0,gm,h,beta0,s0,state)
    #  else
    #    iter = kep_hyperbolic!(x0,v0,r0,gm,h,beta0,s0,state)
    #  end
    return
end

"""

Solves equation (35) from Wisdom & Hernandez for the elliptic case.
"""
function kep_ell_hyp!(x0::Array{T,1},v0::Array{T,1},r0::T,k::T,h::T,
                      beta0::T,s0::T,state::Array{T,1}) where {T <: Real}
    # Now, solve for s in elliptical Kepler case:
    f = zero(T); g=zero(T); dfdt=zero(T); dgdt=zero(T); cx=zero(T);sx=zero(T);g1bs=zero(T);g2bs=zero(T)
    s=zero(T); ds = zero(T); r = zero(T);rinv=zero(T); iter=0
    if beta0 > zero(T) || beta0 < zero(T)
        s,f,g,dfdt,dgdt,cx,sx,g1bs,g2bs,r,rinv,ds,iter = solve_kepler!(h,k,x0,v0,beta0,r0,
                                                                       s0,state)
    else
        println("Not elliptic or hyperbolic ",beta0," x0 ",x0)
        r= zero(T); fill!(state,zero(T)); rinv=zero(T); s=zero(T); ds=zero(T); iter = 0
    end
    state[8]= r
    state[9] = (state[2]*state[5]+state[3]*state[6]+state[4]*state[7])*rinv
    # recompute beta:
    # beta is element 10 of state:
    state[10] = 2.0*k*rinv-(state[5]*state[5]+state[6]*state[6]+state[7]*state[7])
    # s is element 11 of state:
    state[11] = s
    # ds is element 12 of state:
    state[12] = ds
    return iter
end

"""

Solves elliptic Kepler's equation for both elliptic and hyperbolic cases.
"""
function solve_kepler!(h::T,k::T,x0::Array{T,1},v0::Array{T,1},beta0::T,
                       r0::T,s0::T,state::Array{T,1}) where {T <: Real}
    # Initial guess (if s0 = 0):
    r0inv = inv(r0)
    beta0inv = inv(beta0)
    signb = sign(beta0)
    sqb = sqrt(signb*beta0)
    zeta = k-r0*beta0
    eta = dot(x0,v0)
    sguess = 0.0
    if s0 == zero(T)
        # Use cubic estimate:
        if zeta != zero(T)
            sguess = cubic1(3eta/zeta,6r0/zeta,-6h/zeta)
        else
            if eta != zero(T)
                reta = r0/eta
                disc = reta^2+2h/eta
                if disc > zero(T)
                    sguess =-reta+sqrt(disc)
                else
                    sguess = h*r0inv
                end
            else
                sguess = h*r0inv
            end
        end
        s = copy(sguess)
    else
        s = copy(s0)
    end
    s0 = copy(s)
    y = zero(T); yp = one(T)
    iter = 0
    ds = Inf
    KEPLER_TOL = sqrt(eps(h))
    ITMAX = 20
    while iter == 0 || (abs(ds) > KEPLER_TOL && iter < ITMAX)
        xx = sqb*s
        if beta0 > 0
            sx = sin(xx); cx = cos(xx)
        else
            cx = cosh(xx); sx = exp(xx)-cx
        end
        sx *= sqb
        # Third derivative:
        yppp = zeta*cx - signb*eta*sx
        # First derivative:
        yp = (-yppp+ k)*beta0inv
        # Second derivative:
        ypp = signb*zeta*beta0inv*sx + eta*cx
        y  = (-ypp + eta +k*s)*beta0inv - h  # eqn 35
        # Now, compute fourth-order estimate:
        ds = calc_ds_opt(y,yp,ypp,yppp)
        s += ds
        iter +=1
    end
    if iter == ITMAX
        println("Reached max iterations in solve_kepler: h ",h," s0: ",s0," s: ",s," ds: ",ds)
    end
    #println("sguess: ",sguess," s: ",s," s-sguess: ",s-sguess," ds: ",ds," iter: ",iter)
    # Since we updated s, need to recompute:
    xx = 0.5*sqb*s
    if beta0 > 0
        sx = sin(xx); cx = cos(xx)
    else
        cx = cosh(xx); sx = exp(xx)-cx
    end
    # Now, compute final values:
    g1bs = 2.0*sx*cx/sqb
    g2bs = 2.0*signb*sx^2*beta0inv
    f = one(T) - k*r0inv*g2bs # eqn (25)
    g = r0*g1bs + eta*g2bs # eqn (27)
    for j=1:3
        # Position is components 2-4 of state:
        state[1+j] = x0[j]*f+v0[j]*g
    end
    r = sqrt(state[2]*state[2]+state[3]*state[3]+state[4]*state[4])
    rinv = inv(r)
    dfdt = -k*g1bs*rinv*r0inv
    dgdt = (r0-r0*beta0*g2bs+eta*g1bs)*rinv
    for j=1:3
        # Velocity is components 5-7 of state:
        state[4+j] = x0[j]*dfdt+v0[j]*dgdt
    end
    return s,f,g,dfdt,dgdt,cx,sx,g1bs,g2bs,r,rinv,ds,iter
end

"""

Computes quartic Newton's update to equation y=0 using first through 3rd derivatives.
Uses techniques outlined in Murray & Dermott for Kepler solver.
"""
function calc_ds_opt(y::T,yp::T,ypp::T,yppp::T) where {T <: Real}
# Rearrange to reduce number of divisions:
num = y*yp
den1 = yp*yp-y*ypp*.5
den12 = den1*den1
den2 = yp*den12-num*.5*(ypp*den1-third*num*yppp)
return -y*den12/den2
end


text:
module GoldenSequences
#using Random
using FixedPointNumbers
using Roots
using Base.Iterators
using Base.MathConstants
export GoldenIntSequence, GoldenSequence, GoldenCartesianSequence

"""
	nextodd(x)

Next odd integer or next fixpoint number with odd representation.
"""
nextodd(i::Integer) = i >> 1 << 1 + one(i)
nextodd(z::T) where {T<:Normed} = reinterpret(T, nextodd(z.i))
prevodd(i::Integer) = i >> 1 << 1 - one(i)
prevodd(z::T) where {T<:Normed} = reinterpret(T, prevodd(z.i))

"""
    closecoprime(r, q)

Find number coprime to `q` close to `r`.
"""
function closecoprime(r, q)
	i = zero(r)
	while true
		gcd(r + i, q) == 1 && return r + i
		gcd(r - i, q) == 1 && return r - i
	    i += 1
	end # ends at r == 1 or r == q
end

struct GoldenSequence{T}
	x0::T
	coeff::T
end
Iterators.eltype(::Type{GoldenSequence{T}}) where {T} = T
Iterators.IteratorEltype(::Type{GoldenSequence{T}}) where {T} = Base.HasEltype()
Iterators.IteratorSize(::Type{GoldenSequence{T}}) where {T} = Iterators.IsInfinite()

struct GoldenIntSequence{T}
	x0::T
	coeff::T
end
Iterators.eltype(::Type{GoldenIntSequence{T}}) where {T} = T
Iterators.IteratorEltype(::Type{GoldenIntSequence{T}}) where {T} = Base.HasEltype()
Iterators.IteratorSize(::Type{GoldenIntSequence{T}}) where {T} = Iterators.IsInfinite()

f(x, d) = x^(d+1) - (x + one(x))

"""
GoldenCartesianSequence(dims...)

Golden Cartesian sequence on `CartesianIndices(dims)`.
Iterates `CartesianIndex`.
Spacefilling with  period `prod(dims...)`
if `dims` are coprime.

"""
struct GoldenCartesianSequence{T,S}
	coeff::T
	R::CartesianIndices{S}
end
Iterators.eltype(::Type{GoldenCartesianSequence{T,S}}) where {T,S} = T
Iterators.IteratorSize(::Type{GoldenCartesianSequence{T,S}}) where {T,S} = Iterators.IsInfinite()
Iterators.IteratorEltype(::Type{GoldenCartesianSequence{T,S}}) where {T,S} = Base.HasEltype()

"""
Compute or look up solution to x^(d+1) = x + 1
"""
function Phi(d)
	if d <= length(phis)
		phis[d]
	else
		find_zero(x->f(x, d), (1.0, 2.0))
	end
end
Phi(_, d) = Phi(d)

function GoldenSequence(x0)
	phi = Phi(eltype(T), n)
	t = ntuple(n) do i
	       T(phi^(-i))
	end
	GoldenSequence{typeof(t)}(x0, false*x0 .+ t)
end
function GoldenSequence(x0::NTuple{N,T}) where {N, T}
	phi = Phi(T, N)
	t = ntuple(N) do i
	       T(phi^(-i))
	end
	GoldenSequence{typeof(x0)}(x0, t)
end
GoldenSequence(n::Int) = GoldenSequence(ntuple(n) do i; 0.0 end)
GoldenSequence(T, n) = GoldenSequence(ntuple(n) do i; zero(T) end)

function GoldenSequence(x0::T) where {T <: Number}
	GoldenSequence{T}(x0, convert(T, 1/golden))
end

function GoldenSequence(::Type{T}, n) where {T<:Normed}
	phi = Phi(T, n)
	t = ntuple(n) do i
	       prevodd(T(phi^(-i)))
	end
	GoldenSequence{typeof(t)}(t)
end

function GoldenIntSequence(::Type{T}, n) where {T<:Unsigned}
	phi = Phi(T, n)
	t = ntuple(n) do i
	       (round(T, (typemax(T)+1.0)*phi^(-i)))
	end
	GoldenIntSequence{typeof(t)}(zero.(t), t)
end

GoldenCartesianSequence(R) = GoldenCartesianSequence(CartesianIndices(R))
function GoldenCartesianSequence(R::CartesianIndices)
	n = ndims(R)
	phi = Phi(n)
	t = ntuple(n) do i
	    	closecoprime(round(Int, size(R, i)*phi^(-i)), size(R, i))
	end
	GoldenCartesianSequence(t, R)
end

function Base.getindex(g::GoldenSequence{<:Tuple}, ind)
    ntuple(length(g.coeff)) do i
         (g.x0[i] + g.coeff[i]*ind) % true
    end
end
function Base.getindex(g::GoldenSequence, ind)
     (g.x0 .+ g.coeff.*ind) % true
end

function Base.iterate(g::GoldenSequence, val=g.x0)
	 state = (g.coeff .+ val) .% true
	 val, state
end
function Base.iterate(g::GoldenIntSequence, val=g.x0)
	 state = g.coeff .+ val # intentional overflow
	 val, state
end

function Base.iterate(g::GoldenCartesianSequence, val=g.R[1])
	state = mod.(g.coeff .+ Tuple(val), g.R.indices)
	val, typeof(val)(state)
end

const phis = [
1.618033988749895 # MathConstants.golden
1.324717957244746 # plastic number
1.2207440846057596
1.1673039782614187
1.1347241384015196
1.1127756842787055
1.09698155779856
1.085070245491451
1.0757660660868373
1.0682971889208415
1.0621691678642553
1.0570505752212285
1.0527109201475584
1.0489849347570346
1.045751024156342
1.0429177323017866
1.0404149477818474
1.0381880194364501
1.0361937171306836
1.034397396133807
1.032770966441043
1.0312914124792476
1.0299396966705234
1.0286999355282511
1.0275587723930215
1.0265048941462673
1.0255286547645814
1.0246217791365844
1.0237771278614651
1.0229885088662907
1.022250525317837
1.0215584519243472
1.020908133630809
1.0202959021164726
1.0197185065485594
1.0191730558310543
1.0186569701821966
1.0181679403286723
1.0177038929544147
1.0172629613133726
1.0168434601276772
1.0164438640594171
1.0160627891762175
1.0156989769358962
1.015351280299597
1.0150186516505193
1.0147001322501537
1.0143948430084788
1.0141019763809709
1.0138207892351239
1.013550596553797
1.0132907658630492
1.0130407122890395
1.0127998941626677
1.0125678091024284
1.0123439905158556
1.0121280044682868
1.0119194468747255
1.0117179409765626
1.0115231350700071
1.0113347004574071
1.01115232959636
1.010975734424683
1.0108046448420556
1.0106388073314982
1.0104779837058908
1.0103219499664957
1.0101704952619777
1.0100234209377514
1.00988053966664
1.0097416746528447
1.009606658902116
1.0094753345517855
1.009347552255012
1.009223170614191
1.0091020556590062
1.0089840803650738
1.0088691242095444
1.008757072760392
1.00864781729645
1.0085412544555385
1.0084372859082966
1.008335818055544
1.0082367617472274
1.0081400320211664
1.0080455478599986
1.0079532319648552
1.0078630105444437
1.0077748131183248
1.0076885723332831
1.0076042237917846
1.007521705891604
1.007440959675782
1.007361928692145
1.0072845588616808
1.0072087983551326
1.0071345974772095
1.0070619085578796
1.006990685850238
1.0069208854344955
]

end


text:
"""
  function rxdetect(X::Matrix{T}, alpha::Float64 = 0.99)

Takes data in the form of matrix where first index correspond to realisations and
second to feratures (marginals).
Using the RX (Reed-Xiaoli) Anomaly Detection returns the array of Bool that
correspond to outlier realisations. alpha is the sensitivity parameter of the RX detector

```jldoctest
julia> Random.seed!(42);

julia> x = vcat(rand(8,2), 20*rand(2,2))
10×2 Array{Float64,2}:
  0.533183    0.956916
  0.454029    0.584284
  0.0176868   0.937466
  0.172933    0.160006
  0.958926    0.422956
  0.973566    0.602298
  0.30387     0.363458
  0.176909    0.383491
 11.8582      5.25618
 14.9036     10.059

julia> rxdetect(x, 0.95)
10-element Array{Bool,1}:
 false
 false
 false
 false
 false
 false
 false
 false
  true
  true
```
"""
function rxdetect(X::Matrix{T}, alpha::Float64 = 0.99) where T <: AbstractFloat
  t = size(X,1)
  outliers = fill(false, t)
  mu = mean(X,dims=1)[1,:]
  Kinv = inv(cov(X))
  d = Chisq(size(X,2))
  for i in 1:t
    if (X[i,:] - mu)'*Kinv*(X[i,:] - mu) > quantile(d, alpha)
      outliers[i] = true
    end
  end
  outliers
end

"""
hosvdstep(X::Matrix{T}, ls::Vector{Bool}, β::Float64, r::Int, cc::SymmetricTensor{T,4}) where T <: AbstractFloat

Returns Vector{Bool} - outliers form an itteration step of th hosvd algorithm
and Vector{Float64} vector of univariate kurtosis for data projected on specific
directions
"""
function hosvdstep(X::Matrix{T}, ls::Vector{Bool}, β::Float64, r::Int, cc::SymmetricTensor{T,4}) where T <: AbstractFloat
  bestls = copy(ls)
  M = Array(cum2mat(cc))
  # sorting eigenvalues for highest to lowest
  p = sortperm(eigvals(M), rev = true)
  W = eigvecs(M)[:,p[1:r]]

  Z = X*W
  mm = [mad(Z[ls,i]; center=median(Z[ls,i]), normalize=true) for i in 1:r]
  me = [median(Z[ls,i]) for i in 1:r]
  for i in findall(ls)
    if maximum(abs.(Z[i,:].-me)./mm) .> β
     bestls[i] = false
   end
 end
 bestls, norm([kurtosis(Z[bestls,i]) for i in 1:r])
end

"""
  function hosvdc4detect(X::Matrix{T}, β::Float64 = 4.1, r::Int = 3; b::Int = 4)

Takes data in the form of matrix where first index correspond to realisations and
second to feratures (marginals).
Using the HOSVD of the 4'th cumulant's tensor of data returns the array of Bool that
correspond to outlier realisations. β is the sensitivity parameter while r a
number of specific directions, data are projected onto. Parameter b is a size of
blocks in a SymmetricTensors structure

```jldoctest
julia> Random.seed!(42);

julia> x = vcat(rand(8,2), 20*rand(2,2))
10×2 Array{Float64,2}:
  0.533183    0.956916
  0.454029    0.584284
  0.0176868   0.937466
  0.172933    0.160006
  0.958926    0.422956
  0.973566    0.602298
  0.30387     0.363458
  0.176909    0.383491
 11.8582      5.25618
 14.9036     10.059

julia> rxdetect(x, 0.95)
10-element Array{Bool,1}:
 false
 false
 false
 false
 false
 false
 false
 false
  true
  true
```
"""
function hosvdc4detect(X::Matrix{T}, β::Float64 = 4.1, r::Int = 3; b::Int = 4) where T <: AbstractFloat
  X = X.-mean(X,dims=1)
  s = cov(X)
  X = X*Real.(sqrt(inv(s)))
  ls = fill(true, size(X,1))
  lsold = copy(ls)
  aold = 1000000000.
  ma = momentarray(X, 4, b)
  t = size(X,1)
  while count(ls) > div(size(X,1)+r+1, 2)
    ma, t = updatemoments(ma, t, X, ls, lsold)
    lsold = copy(ls)
    c = copy(ma)
    moms2cums!(c)
    ls, a = hosvdstep(X, lsold, β, r, c[4])
    if (aold-a)/a < 0.0001
      return .!lsold
    end
    aold = a
  end
  .!lsold
end

"""
 updatemoments(ma::Vector{SymmetricTensor{Float64,N} where N}, t::Int, X::Matix{T}, ls::Vector{Bool}, lsold::Vector{Bool})

Returns Array{SymmetricTensor} - an array of updated moments (after outliers removal)
 and t - number of realisations of data after outliers removal
"""
function updatemoments(ma::Vector{SymmetricTensor{Float64,N} where N}, t::Int, X::Matrix{T}, ls::Vector{Bool}, lsold::Vector{Bool}) where T <: AbstractFloat
  lstemp = (lsold.+ls) .== 1
  if true in lstemp
    dt = count(lstemp)
    ma = t/(t-dt).*ma .- dt/(t-dt).*momentarray(X[lstemp,:], 4, ma[1].bls)
    t = t - dt
  end
  ma, t
end


text:
export MakeTestData

function MakeTestData(;
        nz = 15,
        nψ = 13,
        num_i = 50,
        num_zt = 150,
        obs_per_well = 10:20,
        royalty_rates = RoyaltyRatesSmall(),
        reward    = DefaultDrillReward(),
        θ_royalty = Theta(RoyaltyModel()),
        θ_produce = Theta(ProductionModel()),
        θ_drill   = Theta(DefaultDrillReward()),
        wp = FullLeasedProblem(),
        anticipate = true,
        seed=1234,
        ddm_opts = (minmaxleases=1:1, nper_initial=40:40, tstart=1:50)
)

        Random.seed!(seed)

        # grids
        ψs = PsiSpace(nψ)
        u,v = randn(num_i), randn(num_i)
        z = TestPriceProcess(;nsim=num_zt, ngrid=nz)
        _zchars = series(z)
        log_ogip = logOGIP(num_i)

        # Royalty data
        data_roy = TestDataRoyalty(u, v, royalty_rates, θ_royalty, log_ogip)

        # endogenous geology + royalty
        _ichars = SDLParameters.ichars(log_ogip, observed_choices(data_roy))

        ddm   = TestDynamicDrillModel(z, anticipate; reward=reward)
        data_drill = DataDynamicDrill(u, v, _zchars, _ichars, ddm, θ_drill; ddm_opts...)
        data_produce = TestDataProduce(u, data_drill, obs_per_well, θ_produce, log_ogip)

        # big model
        dataset_full  = DataSetofSets(data_drill, data_roy, data_produce, CoefLinks(data_drill))
        theta_full = merge_thetas((θ_drill, θ_royalty, θ_produce), dataset_full)

        rwrd_c = ConstrainedProblem(reward)
        ddm_c = TestDynamicDrillModel(z, anticipate; reward=rwrd_c)
        data_drill_c = DataDrill(data_drill, ddm_c)
        dataset_drill = DataSetofSets(data_drill_c, EmptyDataSet(), EmptyDataSet())
        theta_drill = Theta(rwrd_c)

        return ((dataset_drill, theta_drill, ), (dataset_full, theta_full))
end


text:
import ..update!, ..select_action, ..reset!

using StaticArrays
using Random
using LinearAlgebra:Diagonal

using Ramnet.Models.AltDiscriminators:FunctionalDiscriminator
using Ramnet.Optimizers
import Ramnet.Loss

struct AltContinousActionPolicyLoss <: Loss.AbstractLoss end

# Policy approximator with only two admissible actions: {-1, 1}.
struct AltContinousActionPolicy{OS,AS,T <: Real,O <: StaticArray{Tuple{OS},T,1},A <: StaticArray{Tuple{AS},T,1},E <: AbstractEncoder{T}}
    opt::FunctionalOptimizer{AltContinousActionPolicyLoss}
    f::FunctionalDiscriminator{6,T,E}
    rng::MersenneTwister
end

function AltContinousActionPolicy(::Type{O}, ::Type{A}, n::Int, encoder::E; η::Float64=0.1, epochs::Int=1, partitioner::Symbol=:uniform_random, seed::Union{Nothing,Int}=nothing) where {AS,OS,T <: Real,O <: StaticArray{Tuple{OS},T,1},A <: StaticArray{Tuple{AS},T,1},E <: AbstractEncoder{T}}
    AltContinousActionPolicy{OS,AS,T,O,A,E}(
        FunctionalOptimizer(AltContinousActionPolicyLoss(); learning_rate=η, epochs),
        FunctionalDiscriminator{2 * AS}(OS, n, encoder, partitioner; seed),
        isnothing(seed) ? MersenneTwister() : MersenneTwister(seed)
    )
    end

function gaussian_params(p::AltContinousActionPolicy{OS,AS,T,O,A,E}, observation::O) where {OS,AS,T,O,A,E}
    params = predict(p.f, observation)
    # @show params
    μ = params[1:AS]
    σ = exp.(params[AS + 1:end]) .+ eps()

    return (μ, σ)
end

function Loss.grad(::AltContinousActionPolicyLoss, p::AltContinousActionPolicy{OS,AS,T,O,A,E}, observation::O, action::A, G) where {OS,AS,T,O,A,E}
    μ, σ = gaussian_params(p, observation)

    # @show μ, σ
    inv_cov = inv(Diagonal(σ)^2)
    rel_action = action - μ

    gradient = [-1 * inv_cov * rel_action; 1 .- inv_cov * rel_action.^2] * G

    @assert !any(isnan.(gradient)) "Something wrong with the gradient"

    return gradient
end

# TODO: Validate action (must be ether -1 or +1)
function update!(p::AltContinousActionPolicy{OS,AS,T,O,A,E}, observation::O, action::A, G) where {OS,AS,T,O,A,E}
    gradient = Optimizers.grad(p.opt, p, observation, action, G)

    train!(p.f, observation, learning_rate(p.opt) * SVector{2 * AS}(gradient))

    return nothing
end

# TODO: The clamping of the action is specific to the hopper agent!
function select_action(p::AltContinousActionPolicy{OS,AS,T,O,A,E}, observation::O) where {OS,AS,T,O,A,E}
    μ, σ = gaussian_params(p, observation)

# @show μ
    # @show σ

    cov = Diagonal(σ)^2

    return A(clamp.(cov * randn(p.rng, AS) + μ, -1.0f0, 1.0f0))
end

function reset!(p::AltContinousActionPolicy)
    Ramnet.reset!(p.f)

    return nothing
end


text:
using MonteCarloMarkovKernels
using StaticArrays
import Statistics.mean
using Random
using LinearAlgebra
using Test

function makelogMVN(μ::SVector{d, Float64}, Σ::SMatrix{d, d, Float64}) where d
  invΣ = inv(Σ)
  lognc = - 0.5 * d * log(2 * π) - 0.5 * logdet(Σ)
  function lpi(x::SVector{d, Float64})
    v = x - μ
    return lognc - 0.5*dot(v, invΣ * v)
  end
  return lpi
end

function testd(d::Int64)
  zerod = zeros(d)
  Szerod = SVector{d, Float64}(zerod)
  μ = randn(d)
  A = randn(d, d)
  Σ = A * A'

  logtarget = makelogMVN(SVector{d, Float64}(μ), SMatrix{d, d, Float64}(Σ))

  # P_AM = makeAMKernel(logtarget, d, Random.GLOBAL_RNG, d^2)
  P_AM = makeAMKernel(logtarget, Matrix{Float64}(I, d, d), Random.GLOBAL_RNG, d^2)
  chain = simulateChainProgress(P_AM, Szerod, d*2^21)

  @test maximum(abs.(mean(chain) - μ)) < 0.1
  @test maximum(abs.(MonteCarloMarkovKernels.cov(chain) - Σ)) < 0.1

  @test 0.234 < P_AM(:acceptanceRate) < 0.45
  @test P_AM(:meanEstimate) ≈ mean(chain) atol=0.01
  @test maximum(abs.(MonteCarloMarkovKernels.cov(chain) - P_AM(:covEstimate))) < 0.01

end

Random.seed!(12345)

for d in 1:3
  testd(d)
end


text:
### A Pluto.jl notebook ###
# v0.18.4

using Markdown
using InteractiveUtils

# ╔═╡ d7967150-5a59-4199-83ad-4e79476c0d44
using PlutoUI: with_terminal

# ╔═╡ 1a45c71a-4204-461b-8f85-61d9a79f8cf8
using BenchmarkTools

# ╔═╡ bdacad08-a9b1-11ec-1013-45393b4e19c9
# hideall
title = "Benchmarking";

# ╔═╡ 031e949e-e784-4d9d-b767-0a3e2d29d3a2
"""
+++
title = "$title"
+++
""" |> Base.Text

# ╔═╡ 8a859941-0f2d-45fb-b56b-a3f58f97a36a
md"""
# $title

This notebook shows some `@benchmark` and `@code_warntype` output.
Both outputs are shown via the `with_terminal` from [PlutoUI.jl](https://github.com/JuliaPluto/PlutoUI.jl), see below.

We define some function `double` and a dictionary of numbers in order to show type inference problems via `@code_warntype`:
"""

# ╔═╡ 90cb7685-09fc-4f2a-88c9-64fada9e50d0
numbers = Dict(:one => 1f0, :two => 2.0)

# ╔═╡ 3f0e2049-8597-4dac-b499-4d7a8a35978e
function double(mapping, key::Symbol)
    return 2 * mapping[key]
end;

# ╔═╡ faef1ae5-eeac-4e23-ba75-23ae2f1ca629
md"""
Now, the code works.
"""

# ╔═╡ a3c923e5-2873-48a2-ad0f-0087585136a7
double(numbers, :one)

# ╔═╡ a16c4a67-8666-4dc5-b100-4c2c7c51f981
double(numbers, :two)

# ╔═╡ 187e21bc-db38-48c9-97cd-7894685e1b43
md"""
But the `@code_warntype` shows some big red warnings:
"""

# ╔═╡ a3b4bcc6-4aa1-48b2-b598-0613794beab2
with_terminal() do
    @code_warntype double(numbers, :one)
end

# ╔═╡ 8c147812-f0e7-4b14-9082-3ac4d5f07000
md"""
We can fix this by forcing all elements in the dictionary to have the same type.
Specifically, to we force all elements to be of type `Float64`:
"""

# ╔═╡ 0f78717a-ccbd-4af1-a150-adcce25cddb4
typednumbers = Dict{Symbol, Float64}(:one => 1f0, :two => 2.0)

# ╔═╡ 48d6f2fa-04e8-475c-92b6-43418f6a2c6c
md"""
This gets rid of all the type warnings:
"""

# ╔═╡ c678d759-9a68-439e-9cb3-3cfe94088794
with_terminal() do
    @code_warntype double(typednumbers, :one)
end

# ╔═╡ 3cf79c82-9d43-41d3-b968-e2d283bc32ea
md"""
And makes the method more quick:
"""

# ╔═╡ b52f7449-aed4-4eba-905f-3c078fba4f3f
md"""
## Appendix
"""

# ╔═╡ 064395ff-2cc6-4be1-a6f0-0c3ccefdf8ce
function with_benchmark_terminal(f)
    out = sprint(show, "text/plain", f())
    with_terminal() do
        print(out)
    end
end;

# ╔═╡ bc815bd5-1f5e-4641-81f4-80e0306398bf
with_benchmark_terminal() do
    @benchmark double(numbers, :one)
end

# ╔═╡ b8f4d737-4647-47f1-8c45-c3ed30994e1f
with_benchmark_terminal() do
    @benchmark double(typednumbers, :one)
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
BenchmarkTools = "~1.3.1"
PlutoUI = "~0.7.37"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.BenchmarkTools]]
deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
git-tree-sha1 = "4c10eee4af024676200bc7752e536f858c6b8f93"
uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
version = "1.3.1"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.3"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "bf0a1121af131d9974241ba53f601211e9303a9e"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.37"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.Profile]]
deps = ["Printf"]
uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═bdacad08-a9b1-11ec-1013-45393b4e19c9
# ╠═031e949e-e784-4d9d-b767-0a3e2d29d3a2
# ╠═8a859941-0f2d-45fb-b56b-a3f58f97a36a
# ╠═90cb7685-09fc-4f2a-88c9-64fada9e50d0
# ╠═3f0e2049-8597-4dac-b499-4d7a8a35978e
# ╠═faef1ae5-eeac-4e23-ba75-23ae2f1ca629
# ╠═a3c923e5-2873-48a2-ad0f-0087585136a7
# ╠═a16c4a67-8666-4dc5-b100-4c2c7c51f981
# ╠═187e21bc-db38-48c9-97cd-7894685e1b43
# ╠═d7967150-5a59-4199-83ad-4e79476c0d44
# ╠═a3b4bcc6-4aa1-48b2-b598-0613794beab2
# ╠═8c147812-f0e7-4b14-9082-3ac4d5f07000
# ╠═0f78717a-ccbd-4af1-a150-adcce25cddb4
# ╠═48d6f2fa-04e8-475c-92b6-43418f6a2c6c
# ╠═c678d759-9a68-439e-9cb3-3cfe94088794
# ╠═3cf79c82-9d43-41d3-b968-e2d283bc32ea
# ╠═1a45c71a-4204-461b-8f85-61d9a79f8cf8
# ╠═bc815bd5-1f5e-4641-81f4-80e0306398bf
# ╠═b8f4d737-4647-47f1-8c45-c3ed30994e1f
# ╠═b52f7449-aed4-4eba-905f-3c078fba4f3f
# ╠═064395ff-2cc6-4be1-a6f0-0c3ccefdf8ce
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: rikhuijzer/JuliaTutorialsTemplate

max_issues_repo_name: rikhuijzer/JuliaTutorialsTemplate

max_forks_repo_name: rikhuijzer/JuliaTutorialsTemplate


text:
module UpdateIndex

export @update!, updateindex!

macro update!(ex)
    @assert isa(ex.head, Symbol)
    opstr = string(ex.head)
    @assert length(opstr) == 2
    op = Symbol(opstr[1:1])
    @assert Symbol(opstr[2:2]) == :(=)
    arr = ex.args[1].args[1]
    inds = ex.args[1].args[2:end]
    s = ex.args[2]
    return :(updateindex!($(esc(arr)), $(esc(op)), $(esc(s)), $(map(esc, inds)...)))
end

updateindex!(A, f, s, i) = A[i] = f(A[i], s)
updateindex!(A, f, s, i, j) = A[i, j] = f(A[i, j], s)
updateindex!(A, f, s, i, j, k) = A[i, j, k] = f(A[i, j, k], s)
updateindex!(A, f, s, i, j, k, l) = A[i, j, k, l] = f(A[i, j, k, l], s)
updateindex!(A, f, s, I::Vararg) = A[I...] = f(A[I...], s)

# SparseMatrix updateindex!
import Base.Order.Forward

function updateindex!{T,Ti}(A::SparseMatrixCSC{T,Ti}, f, s, i0::Integer, i1::Integer)
    i0 = convert(Ti, i0)
    i1 = convert(Ti, i1)
    if !(1 <= i0 <= A.m && 1 <= i1 <= A.n); throw(BoundsError()); end
    r1 = Int(A.colptr[i1])
    r2 = Int(A.colptr[i1+1]-1)
    i = (r1 > r2) ? r1 : searchsortedfirst(A.rowval, i0, r1, r2, Forward)

    if (i <= r2) && (A.rowval[i] == i0)
        A.nzval[i] = f(A.nzval[i], s)
    else
        insert!(A.rowval, i, i0)
        insert!(A.nzval, i, f(zero(T), s))
        @simd for j = (i1+1):(A.n+1)
            @inbounds A.colptr[j] += 1
        end
    end
    return A
end

end # module


text:
let
	bn = BayesNet()
	@test length(bn) == 0

	push!(bn, StaticCPD(:a, Normal(0.0,1.0)))

	@test length(bn) == 1
	@test name(get(bn, :a)) == :a
	@test name(get(bn, 1)) == :a
	@test names(bn) == [:a]
	@test parents(bn, :a) == NodeName[]
	@test children(bn, :a) == NodeName[]
	@test neighbors(bn, :a) == NodeName[]
	@test descendants(bn, :a) == NodeName[]
	@test pdf(bn, Assignment(:a=>0.0)) == pdf(Normal(0.0,1.0), 0.0)
	@test logpdf(bn, Assignment(:a=>0.0)) == logpdf(Normal(0.0,1.0), 0.0)
	@test isapprox(pdf(bn, DataFrame(a=[0.0, 1.0])), pdf(Normal(0.0,1.0), 0.0) * pdf(Normal(0.0,1.0), 1.0))
	@test isapprox(logpdf(bn, DataFrame(a=[0.0, 1.0])), logpdf(Normal(0.0,1.0), 0.0) + logpdf(Normal(0.0,1.0), 1.0))

	# b = 2a + 1
	push!(bn, LinearGaussianCPD(:b, [:a], [2.0], 1.0, 1.0))
	@test length(bn) == 2
	@test name(get(bn, :b)) == :b
	@test name(get(bn, 2)) == :b
	@test has_edge(bn, :a, :b)
	@test names(bn) == [:a, :b]
	@test parents(bn, :a) == NodeName[]
	@test children(bn, :a) == [:b]
	@test neighbors(bn, :a) == [:b]
	@test descendants(bn, :a) == [:b]
	@test parents(bn, :b) == [:a]
	@test children(bn, :b) == NodeName[]
	@test neighbors(bn, :b) == [:a]
	@test descendants(bn, :b) == NodeName[]
	@test pdf(bn, Assignment(:a=>1.0, :b=>2.0)) == pdf(Normal(0.0,1.0), 1.0) * pdf(Normal(3.0, 1.0), 2.0)
	@test logpdf(bn, Assignment(:a=>1.0, :b=>2.0)) == logpdf(Normal(0.0,1.0), 1.0) + logpdf(Normal(3.0, 1.0), 2.0)

	# test removal
	delete!(bn, :b)
	@test length(bn) == 1
	@test name(get(bn, :a)) == :a
	@test name(get(bn, 1)) == :a
	@test !has_edge(bn, :a, :b)
	@test names(bn) == [:a]
	@test parents(bn, :a) == NodeName[]
	@test children(bn, :a) == NodeName[]
end

let
	bn = BayesNet()
	push!(bn, StaticCPD(:A, Normal(0.0,1.0)))
	push!(bn, StaticCPD(:B, [:A], Normal(0.0,1.0)))
	push!(bn, StaticCPD(:D, [:A], Normal(0.0,1.0)))
	push!(bn, StaticCPD(:C, [:A, :B, :D], Normal(0.0, 1.0)))
	push!(bn, StaticCPD(:F, Normal(0.0,1.0)))
	push!(bn, StaticCPD(:E, [:F], Normal(0.0,1.0)))
	@test is_independent(bn, [:B], [:D], [:A])
	@test !is_independent(bn, [:B], [:D], [:C])
	@test !is_independent(bn, [:B], [:D], [:E])
	@test is_independent(bn, [:B], [:C], [:A])
end

let
	bn = BayesNet()
	push!(bn, StaticCPD(:B, Normal(0.0,1.0)))
	push!(bn, StaticCPD(:S, Normal(0.0,1.0)))
	push!(bn, StaticCPD(:E, [:B, :S], Normal(0.0,1.0)))
	push!(bn, StaticCPD(:D, [:E], Normal(0.0,1.0)))
	push!(bn, StaticCPD(:C, [:E], Normal(0.0,1.0)))
	@test !is_independent(bn, [:B], [:S], [:E, :D])
	@test !is_independent(bn, [:B], [:S], [:C])
end


text:
using Catlab.CategoricalAlgebra, Catlab.Graphs, Catlab.Present, Catlab.Graphics, Catlab.Theories
using Catlab.CategoricalAlgebra.FinCats: FinCatGraphEq

using Distributions: Exponential, Geometric, cdf
using Random: randsubseq
using StatsBase: sample

using Plots

# include("utils.jl")

# Bernoulli sample of vector `m` (all same probability)
function sample_matches(m::AbstractVector{T}, r::AbstractFloat, Δt) where {T}
    p = cdf(Exponential(), r * Δt)
    randsubseq(m, p)
end

# update matches; return nothing if invalid
function postcompose_partial(kg::ACSetTransformation, kh::ACSetTransformation, m::ACSetTransformation)
    d = Dict()
    for (k,vs) in pairs(components(m))
      vs_ = Int[]
      for v in collect(vs)
        kv = findfirst(==(v), collect(kg[k]))
        if isnothing(kv)
          return nothing
        else
          push!(vs_, kh[k](kv))
        end
      end
      d[k] = vs_
    end
    return ACSetTransformation(dom(m), codom(kh); d...)
  end

  # stuff to store stuff
  struct Rule
      L::ACSetTransformation
      R::ACSetTransformation
  end

  mutable struct MatchedRule
      rule::Rule
      match::Union{Nothing, ACSetTransformation}
  end

# schema
@present ThSIR(FreeSchema) begin
    (S,I,R)::Ob
end

@acset_type SIR(ThSIR)

# infection rules
I = @acset SIR begin I=1 end
I2 = @acset SIR begin I=2 end
SI = @acset SIR begin S=1; I=1 end
L_infect =  ACSetTransformation(I, SI; I=[1]) # fn from dom I to codom SI, I is the injective function mapping stuff in dom to codom
R_infect = ACSetTransformation(I, I2; I=[1])

# recovery rules
R = @acset SIR begin R=1 end
L_recovery = ACSetTransformation(SIR(), I)
R_recovery = ACSetTransformation(SIR(), R)

# structs
rule_inf = Rule(L_infect, R_infect)
rule_rec = Rule(L_recovery, R_recovery)

N = 500
I0 = 10
S0 = N - I0

Δt = 0.1
tmax = 100
steps = Int(tmax/Δt)
γ = 1/10 # recovery rate
R0 = 2.5
β = R0 * γ # R0 for corresponding ODEs

state = @acset SIR begin S=S0; I=I0; R=N-S0-I0 end

# simulation loop
out = fill(-1, (steps, 3))

for t = 1:steps
    # possible events
    infections = homomorphisms(SI, state)
    recoveries = homomorphisms(I, state)
    # sample occurances
    infections = sample_matches(infections, β/N, Δt)
    recoveries = sample_matches(recoveries, γ, Δt)
    # apply events
    events = MatchedRule[[MatchedRule(rule_inf, m) for m in infections]..., [MatchedRule(rule_rec, m) for m in recoveries]...]
    while length(events) > 0
        ev = pop!(events)
        # apply event
        _, kg, _, kh = rewrite_match_maps(ev.rule.L, ev.rule.R, ev.match)
        state = codom(kh)
        # update the remaining matches post rewrite
        for j in 1:length(events)
            events[j].match = postcompose_partial(kg, kh, events[j].match)
        end
        events = filter((e) -> !isnothing(e.match), events)
    end
    # write output
    out[t, :] = [nparts(state, x) for x in [:S, :I, :R]]
end

plot(
    (1:steps) * Δt,
    out,
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number"
)

# # --------------------------------------------------------
# match = homomorphism(SI, state)
# match = homomorphisms(SI, state)

# match = homomorphism(I, state)
# match = homomorphisms(I, state)

# # rewrite_match(L_recovery, R_recovery, match[1])

# # sample_matches(match, γ, Δt)

# state

# rewrite(L_recovery, R_recovery, state)

# # not sure how to update morphisms after the first one was applied.
# nparts(state, :R)
# state = apply_matches(state,match,L_recovery,R_recovery)
# nparts(state, :R)


text:
"""
    mat_C1(q,N)

Generate the matrix corresponding to even, π-periodic solutions to the
angular Mathieu equation.

* `q`: real parameter of Mathieu’s equation
* `N`: truncated matrix dimension (matrix is N×N)
"""
mat_C1(q::Number,N::Integer) = mat_C1(blasfloat(q),N)
function mat_C1(q::T,N::Integer) where T<:BlasReal
    dd = [T(2k)^2 for k = 0:N-1]
    du = fill(q,N-1)
    du[1] *= √2
    return SymTridiagonal(dd,du)
end

"""
    mat_C2(q,N)

Generate the matrix corresponding to even, 2π-periodic (π-antiperiodic)
solutions to the angular Mathieu equation.

* `q`: real parameter of Mathieu’s equation
* `N`: truncated matrix dimension (matrix is N×N)
"""
mat_C2(q::Number,N::Integer) = mat_C2(blasfloat(q),N)
function mat_C2(q::T,N::Integer) where T<:BlasReal
    dd = [T(2k+1)^2 for k = 0:N-1]
    dd[1] += q
    du = fill(q,N-1)
    return SymTridiagonal(dd,du)
end

"""
    mat_C3(q,N)

Generate the matrix corresponding to odd, 2π-periodic (π-antiperiodic)
solutions to the angular Mathieu equation.

* `q`: real parameter of Mathieu’s equation
* `N`: truncated matrix dimension (matrix is N×N)
"""
mat_C3(q::Number,N::Integer) = mat_C3(blasfloat(q),N)
function mat_C3(q::T,N::Integer) where T<:BlasReal
    dd = [T(2k+1)^2 for k = 0:N-1]
    dd[1] -= q
    du = fill(q,N-1)
    return SymTridiagonal(dd,du)
end

"""
    mat_C4(q,N)

Generate the matrix corresponding to odd, π-periodic solutions to the
angular Mathieu equation.

* `q`: real parameter of Mathieu’s equation
* `N`: truncated matrix dimension (matrix is N×N)
"""
mat_C4(q::Number,N::Integer) = mat_C4(blasfloat(q),N)
function mat_C4(q::T,N::Integer) where T<:BlasReal
    dd = [T(2k+2)^2 for k = 0:N-1]
    du = fill(q,N-1)
    return SymTridiagonal(dd,du)
end


text:
function michaelis_menten(u, p, t)
    [0.6 - 1.5u[1]/(0.3+u[1])]
end

u0 = [0.5]
tspan = (0.0, 4.0)

problem_1 = ODEProblem(michaelis_menten, u0, tspan)
solution_1 = solve(problem_1, Tsit5(), saveat = 0.1)
problem_2 = ODEProblem(michaelis_menten, 2*u0, tspan)
solution_2 = solve(problem_2, Tsit5(), saveat = 0.1)
X = [solution_1[:,:] solution_2[:,:]]
ts = [solution_1.t; solution_2.t]

DX = similar(X)
for (i, xi) in enumerate(eachcol(X))
    DX[:, i] = michaelis_menten(xi, [], ts[i])
end

@parameters t
@variables u[1:2]
h = [monomial_basis(u[1:1], 4)...]
basis = Basis([h; h .* u[2]], u)

@testset "Ideal data" begin

    prob = ContinuousDataDrivenProblem(X, ts, DX)

    opts = [ImplicitOptimizer(5e-1);ImplicitOptimizer(0.4:0.1:0.7)]
    for opt in opts
        res = solve(prob, basis, opt, normalize = false, denoise = false, maxiter = 10000)
        m = metrics(res)
        @test m.Error < 5e-1
        @test m.AICC < 23.0
        @test m.Sparsity == 4
    end

    for opt in [ADM(5e-1); ADM(4e-1:1e-3:5e-1)]
        res = solve(prob, basis, opt, normalize = false, denoise = false, maxiter = 1000)
        m = metrics(res)
        @test m.Error < 8e-1
        @test m.AICC < 23.0
        @test m.Sparsity == 4
    end
end

Random.seed!(2345)
X = X .+ 1e-3*randn(size(X))

@testset "Noisy data" begin

    prob = ContinuousDataDrivenProblem(X, ts, GaussianKernel())

    for opt in [ImplicitOptimizer(4e-1);ImplicitOptimizer([0.3; 0.4; 0.5])]
        res = solve(prob, basis, opt, normalize = true, denoise = true)
        m = metrics(res)
        @test m.Error < 3e-1
        @test m.AICC < 35.0
        @test m.Sparsity == 4
    end

    for opt in [ADM(0.01:0.01:4e-1)]
        res = solve(prob, basis, opt, normalize = false, denoise = false)
        m = metrics(res)
        @test m.Error < 5e-1
        @test m.AICC < 12.0
        @test m.Sparsity == 4
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct UpdateEncryptionConfig <: SwaggerModel
    keyVaultMetaInfo::Any # spec type: Union{ Nothing, UpdateKeyVaultMetaInfo } # spec name: keyVaultMetaInfo

    function UpdateEncryptionConfig(;keyVaultMetaInfo=nothing)
        o = new()
        validate_property(UpdateEncryptionConfig, Symbol("keyVaultMetaInfo"), keyVaultMetaInfo)
        setfield!(o, Symbol("keyVaultMetaInfo"), keyVaultMetaInfo)
        o
    end
end # type UpdateEncryptionConfig

const _property_map_UpdateEncryptionConfig = Dict{Symbol,Symbol}(Symbol("keyVaultMetaInfo")=>Symbol("keyVaultMetaInfo"))
const _property_types_UpdateEncryptionConfig = Dict{Symbol,String}(Symbol("keyVaultMetaInfo")=>"UpdateKeyVaultMetaInfo")
Base.propertynames(::Type{ UpdateEncryptionConfig }) = collect(keys(_property_map_UpdateEncryptionConfig))
Swagger.property_type(::Type{ UpdateEncryptionConfig }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_UpdateEncryptionConfig[name]))}
Swagger.field_name(::Type{ UpdateEncryptionConfig }, property_name::Symbol) =  _property_map_UpdateEncryptionConfig[property_name]

function check_required(o::UpdateEncryptionConfig)
    true
end

function validate_property(::Type{ UpdateEncryptionConfig }, name::Symbol, val)
end


text:
# This is part of Calendars.jl. See the copyright note there.
# ============================ Tests ========================

using Calendars
using Dates: yearmonthday, now
using Test

# Planetary week
sunday    = 0
monday    = 1
tuesday   = 2
wednesday = 3
thursday  = 4
friday    = 5
saturday  = 6
weeklen   = 7

function TestDateGregorian()
    SomeDateGregorian = [
        (CE,    1, 1,  1)::CDate,
        (CE, 1756, 1, 27)::CDate,
        (CE, 2022, 1,  1)::CDate,
        (CE, 2022, 9, 26)::CDate
    ]

    @testset "DateGregorian" begin
    for date in SomeDateGregorian
        dn = DayNumberFromDate(date)
        gdate = DateFromDayNumber("Gregorian", dn)
        @test gdate == (CE, date[2], date[3], date[4])
        println(CDateStr(date), " -> ", CDateStr(EN, dn), " -> ", CDateStr(gdate))
    end
    end
end

function TestDateJulian()
    SomeDateJulian = [
        (JD,    1,  1,  3)::CDate,
        (JD, 1756,  1, 16)::CDate,
        (JD, 2021, 12, 19)::CDate,
        (JD, 2022,  9, 13)::CDate
    ]

    @testset "DateJulian" begin
    for date in SomeDateJulian
        dn = DayNumberFromDate(date)
        jdate = DateFromDayNumber("Julian", dn)
        @test jdate == (JD, date[2], date[3], date[4])
        println(CDateStr(date), " -> ", CDateStr(EN, dn), " -> ", CDateStr(jdate))
    end
    end
end

function TestDateHebrew()
    SomeDateHebrew  = [
        (AM, 5516, 11, 25)::CDate,
        (AM, 5782,  7,  1)::CDate,
        (AM, 5782, 10, 27)::CDate,
        (AM, 5782,  6, 29)::CDate,
        (AM, 5783,  7,  1)::CDate
    ]

    @testset "DateHebrew" begin
    for date in SomeDateHebrew
        dn = DayNumberFromDate(date)
        hdate = DateFromDayNumber("Hebrew", dn)
        @test hdate == (AM, date[2], date[3], date[4])
        println(CDateStr(date), " -> ", CDateStr(EN, dn), " -> ", CDateStr(hdate))
    end
    end
end

function TestDateIslamic()
    SomeDateIslamic = [
        (AH,    1, 1,  1)::CDate,
        (AH, 1756, 1, 27)::CDate,
        (AH, 2022, 1,  1)::CDate,
        (AH, 2022, 9, 26)::CDate
    ]

    @testset "DateIslamic" begin
    for date in SomeDateIslamic
        dn = DayNumberFromDate(date)
        idate = DateFromDayNumber("Islamic", dn)
        @test idate == (AH, date[2], date[3], date[4])
        println(CDateStr(date), " -> ", CDateStr(EN, dn), " -> ", CDateStr(idate))
    end
    end
end

function TestDateIso()
    SomeDateIso = [
        (ID,    1,  1, 1)::CDate,
        (ID, 1756,  5, 2)::CDate,
        (ID, 2021, 52, 6)::CDate,
        (ID, 2022, 39, 1)::CDate
    ]

    @testset "DateIso" begin
    for date in SomeDateIso
        dn = DayNumberFromDate(date)
        idate = DateFromDayNumber("IsoDate", dn)
        @test idate == (ID, date[2], date[3], date[4])
        println(CDateStr(date), " -> ", CDateStr(EN, dn), " -> ", CDateStr(idate))
    end
    end
end

function TestConversions()

    @testset "Conversions" begin
    @test (AM, 9543, 7,24) == ConvertDate("Gregorian", 5782, 10, 28, "Hebrew")
    @test (AH, 5319, 8,25) == ConvertDate("Gregorian", 5782, 10, 28, "Islamic")
    @test (AM, 5516,11,25) == ConvertDate("Gregorian", 1756,  1, 27, "Hebrew")
    @test (AH,    1, 1, 1) == ConvertDate("Gregorian", 622,   7, 19, "Islamic")

    println("=")
    @test (AH, 1443, 5,27) == ConvertDate("Hebrew", 5782, 10, 28, "Islamic")
    @test (CE, 5782,10,28) == ConvertDate("Hebrew", 9543,  7, 24, "Gregorian")
    @test (AH, 4501, 8, 3) == ConvertDate("Hebrew", 8749, 11, 03, "Islamic")

    println("==")
    @test (AM, 7025,11,23) == ConvertDate("Islamic", 2724, 08, 23, "Hebrew")
    @test (CE, 7025,11,29) == ConvertDate("Islamic", 6600, 11, 21, "Gregorian")
    @test (AM, 9992,12,27) == ConvertDate("Islamic", 5782, 10, 28, "Hebrew")

    @test (AM, 9543, 7,24) == ConvertDate(CE, 5782, 10, 28, AM)
    @test (AH, 5319, 8,25) == ConvertDate(CE, 5782, 10, 28, AH)
    @test (AM, 5516,11,25) == ConvertDate(CE, 1756,  1, 27, AM)
    @test (AH,    1, 1, 1) == ConvertDate(CE, 622,   7, 19, AH)

    println("===")
    @test (AH, 1443, 5,27) == ConvertDate(AM, 5782, 10, 28, AH)
    @test (CE, 5782,10,28) == ConvertDate(AM, 9543,  7, 24, CE)
    @test (AH, 4501, 8, 3) == ConvertDate(AM, 8749, 11, 03, AH)

    println("====")
    @test (AM, 7025,11,23) == ConvertDate(AH, 2724, 08, 23, AM)
    @test (CE, 7025,11,29) == ConvertDate(AH, 6600, 11, 21, CE)
    @test (AM, 9992,12,27) == ConvertDate(AH, 5782, 10, 28, AM)

    end

    println("\nThe following 3 conversions test warnings, so worry only if you do not see them.\n")
    println("This Hebrew date is before the Julian epoch.")
    ConvertDate("Hebrew", 2022,  1,  1, "Gregorian") |> println
    ConvertDate("Hebrew", 5782, 10, 28, "Georgian")  |> println
    ConvertDate("Homebrew", 2022, 1, 1, "Gregorian") |> println
end

function TestIso()

    ISO = [
        ((CE, 1977,01,01 ), (ID, 1976,53,6 )),
        ((CE, 1977,01,02 ), (ID, 1976,53,7 )),
        ((CE, 1977,12,31 ), (ID, 1977,52,6 )),
        ((CE, 1978,01,01 ), (ID, 1977,52,7 )),
        ((CE, 1978,01,02 ), (ID, 1978,01,1 )),
        ((CE, 1978,12,31 ), (ID, 1978,52,7 )),
        ((CE, 1979,01,01 ), (ID, 1979,01,1 )),
        ((CE, 1979,12,30 ), (ID, 1979,52,7 )),
        ((CE, 1979,12,31 ), (ID, 1980,01,1 )),
        ((CE, 1980,01,01 ), (ID, 1980,01,2 )),
        ((CE, 1980,12,28 ), (ID, 1980,52,7 )),
        ((CE, 1980,12,29 ), (ID, 1981,01,1 )),
        ((CE, 1980,12,30 ), (ID, 1981,01,2 )),
        ((CE, 1980,12,31 ), (ID, 1981,01,3 )),
        ((CE, 1981,01,01 ), (ID, 1981,01,4 )),
        ((CE, 1981,12,31 ), (ID, 1981,53,4 )),
        ((CE, 1982,01,01 ), (ID, 1981,53,5 )),
        ((CE, 1982,01,02 ), (ID, 1981,53,6 )),
        ((CE, 1982,01,03 ), (ID, 1981,53,7 ))
    ]

    @testset "ISODates" begin
    for iso in ISO
        a = DayNumberFromDate(iso[1])
        b = DayNumberFromDate(iso[2])
        @test a == b
        println(a, " = ", b)
    end
    end
end

function TestDateTables()

TestDates = [
((EN,0,0,405735),(CE,1111,11,11),(JD,1111,11,4), (EC,1111,11,4),(ID,1111,45,6),(AM,4872,9,2),  (AH, 505,4,29), (JN,0,0,2127158))::DateTable,
((EN,0,0,811258),(CE,2222,2,22), (JD,2222,2,7),  (EC,2222,2,22),(ID,2222,8,5), (AM,5982,12,10),(AH,1649,9,10), (JN,0,0,2532681))::DateTable,
((EN,0,0,641029),(CE,1756,1,27), (JD,1756,1,16), (EC,1756,1,27),(ID,1756,5,2), (AM,5516,11,25),(AH,1169,4,24), (JN,0,0,2362452))::DateTable,
((EN,0,0,738158),(CE,2022,1,1),  (JD,2021,12,19),(EC,2022,1,1), (ID,2021,52,6),(AM,5782,10,28),(AH,1443,5,27), (JN,0,0,2459581))::DateTable,
((EN,0,0,738426),(CE,2022,9,26), (JD,2022,9,13), (EC,2022,9,26),(ID,2022,39,1),(AM,5783,7,1),  (AH,1444,2,29), (JN,0,0,2459849))::DateTable,
((EN,0,0,1146992),(CE, 3141,5,9),(JD,3141,4,17), (EC,3141,5,9), (ID,3141,19,5),(AM,6901,2,9),  (AH,2597,2,10), (JN,0,0,2868415))::DateTable
]
    @testset "DayNumbers" begin
    for D in TestDates
        println()
        PrintDateTable(D)
        num = D[1][4]

        for d in D[1:end]
            dn = DayNumberFromDate(d)
            if dn > 0
                @test dn == num
            end
        end
    end
    end
end

function TestDuration()
    @testset "Duration" begin
    @test 2 == Duration((CE, 2022, 1, 1), (ID, 2022, 1, 1), true)
    @test 2 == Duration((ID, 2022, 1, 1), (CE, 2022, 1, 1), true)
    @test 0 == Duration((CE, 2022, 1, 1), (CE, 2022, 1, 1), true)
    @test 0 == Duration((ID, 2022, 1, 1), (CE, 2022, 1, 3), true)
    # Nota bene, not 365! Duration is a paradise for off-by-one errors.
    @test 364 == Duration((CE, 2022, 1, 1), (CE, 2022, 12, 31), true)
    @test 365 == Duration((CE, 2022, 1, 1), (CE, 2023,  1,  1), true)
    end
end

function TestDayOfYear()
    DOY = [
        ((CE, 1000,02,28),  59),
        ((CE, 1000,02,29),   0),
        ((CE, 1111,11,11), 315),
        ((CE, 1756,01,27),  27),
        ((CE, 1900,02,28),  59),
        ((CE, 1900,02,29),   0),
        ((CE, 1901,02,28),  59),
        ((CE, 1901,02,29),   0),
        ((CE, 1949,11,29), 333),
        ((CE, 1990,10,20), 293),
        ((CE, 2000,02,28),  59),
        ((CE, 2000,02,29),  60),
        ((CE, 2020,02,28),  59),
        ((CE, 2020,02,29),  60),
        ((CE, 2020,07,30), 212),
        ((CE, 2020,08,20), 233),
        ((CE, 2021,09,07), 250),
        ((CE, 2022,01,01),   1),
        ((CE, 2022,09,26), 269),
        ((CE, 2023,07,19), 200),
        ((CE, 2040,02,28),  59),
        ((CE, 2040,02,29),  60),
        ((CE, 2040,09,08), 252),
        ((CE, 2222,02,22),  53),
        ((CE, 3141,05,09), 129)
    ]

    println("\n====================\n")

    @testset "DayOfYear" begin
    for x in DOY
        date, num = x
        doy = DayOfYear(date)
        print(CDateStr(date), " ::", lpad(num, 4, " "), lpad(doy, 4, " "))
        println(" ", doy == num)
        @test doy == num
    end
    end
end

function ShowDayOfYear()

    Dates = [ (22, 1, 1), (756, 1, 27), (949, 11, 29), (990, 10, 20),
    (2022, 1, 1), (1756, 1, 27), (1949, 11, 29), (1990, 10, 20)]

    @testset "DayOfYear2" begin
    for date in Dates
        ce = DayOfYear(CE, date)
        jd = DayOfYear(JD, date)
        ec = DayOfYear(EC, date)
        @test ce == jd
        @test ec == jd
        println(CDateStr(EC, date), " ::", lpad(DayOfYear(EC, date),  4, " "))
    end
    end
end

# Checked with https://www.arc.id.au/Calendar.html
function TestJulianDayNumbers1()

    test = [  #  EC,       EuroNum, JulianNum, ModJulianNum
    ((EC,800,12,25)::CDate, 292194, 2013617, -386384),  # Coronation of Carolus Magnus
    ((EC, 843,8,10)::CDate, 307762, 2029185, -370816),  # Treaty of Verdun
    ((EC,1204,4,12)::CDate, 439498, 2160921, -239080),  # Sack of Constantinople
    ((EC,1452,4,15)::CDate, 530083, 2251506, -148495),  # Birth of Leonardo da Vinci
    ((EC,1666,9, 2)::CDate, 608376, 2329799, -70202),   # Great Fire of London
    ((EC,1789,7,14)::CDate, 653251, 2374674, -25327),   # The Storming of the Bastille
    ((EC,1906,4,28)::CDate, 695906, 2417329, 17328)]    # Birth of Kurt Gödel

    @testset "JulianDayNumbers1" begin
    for t in test
        println()
        edn  = DayNumberFromDate(t[1])
        jdn  = ConvertOrdinalDate(edn, EN, JN)
        @test jdn  == t[3]
        println(edn, " ", jdn, " ")
        edn  = ConvertOrdinalDate(jdn, JN, EN)
        @test edn  == t[2]
        println(jdn, " ", edn)
    end

    en = ConvertOrdinalDate(2440423, JN, EN)
    jn = ConvertOrdinalDate(719000,  EN, JN)
    @test en  == 719000
    @test jn  == 2440423
    println(en, " ", jn)
    end
end

function TestSaveCalendars()
    SaveEuropeanMonth(1, 1)
    SaveEuropeanMonth(1111, 11)
    SaveEuropeanMonth(1582, 10)
    SaveEuropeanMonth(9999, 12)
    println("Note the jump in the European calendar!")
    PrintIsoWeek(1582, 41)
    PrintIsoWeek(2525, 25)
end

function DayOfLife(birthdate::CDate)
    if isValidDate(birthdate)
        y, m, d = yearmonthday(now())
        return Duration(birthdate, (CE, y, m, d)) + 1
    end
    @warn("Invalid Date: $birthdate")
    return InvalidDuration
end

function TestToday()
    println("\nMozart would be ", DayOfLife((EC, 1756, 1, 27)), " days old today.")
    println("\nToday on all calendars:\n")
    (year, month, day) = yearmonthday(now())
    date = (EC, year, month, day)
    CalendarDates(date, true)
    doy = DayOfYear(date)
    println("\nThis is the $doy-th day of a EC-year.")
end

# Test values from Dershowitz & Reingold.
function TestJulianDayNumbers()
# RD,      Weekday,  Julian Day,    Gregorian,      ISODate      Islamic
RD = [
[ 25469, wednesday, 1746893.5, ( 70,  9, 24),   (70, 39, 3),  ( -568,  4,  1)],
[ 49217, sunday,    1770641.5, ( 135, 10,  2),  (135, 39, 7), ( -501,  4,  6)],
[171307, wednesday, 1892731.5, ( 470,  1,  8),  (470,  2, 3), ( -157, 10, 17)],
[210155, monday,    1931579.5, ( 576,  5, 20),  (576, 21, 1), (  -47,  6,  3)],
[253427, saturday,  1974851.5, ( 694, 11, 10),  (694, 45, 6), (   75,  7, 13)],
[369740, sunday,    2091164.5, (1013,  4, 25), (1013, 16, 7), (  403, 10,  5)],
[400085, sunday,    2121509.5, (1096,  5, 24), (1096, 21, 7), (  489,  5, 22)],
[434355, friday,    2155779.5, (1190,  3, 23), (1190, 12, 5), (  586,  2,  7)],
[452605, saturday,  2174029.5, (1240,  3, 10), (1240, 10, 6), (  637,  8,  7)],
[470160, friday,    2191584.5, (1288,  4,  2), (1288, 14, 5), (  687,  2, 20)],
[473837, sunday,    2195261.5, (1298,  4, 27), (1298, 17, 7), (  697,  7,  7)],
[507850, sunday,    2229274.5, (1391,  6, 12), (1391, 23, 7), (  793,  7,  1)],
[524156, wednesday, 2245580.5, (1436,  2,  3), (1436,  5, 3), (  839,  7,  6)],
[544676, saturday,  2266100.5, (1492,  4,  9), (1492, 14, 6), (  897,  6,  1)],
[567118, saturday,  2288542.5, (1553,  9, 19), (1553, 38, 6), (  960,  9, 30)],
[569477, saturday,  2290901.5, (1560,  3,  5), (1560,  9, 6), (  967,  5, 27)],
[601716, wednesday, 2323140.5, (1648,  6, 10), (1648, 24, 3), ( 1058,  5, 18)],
[613424, sunday,    2334848.5, (1680,  6, 30), (1680, 26, 7), ( 1091,  6,  2)],
[626596, friday,    2348020.5, (1716,  7, 24), (1716, 30, 5), ( 1128,  8,  4)],
[645554, sunday,    2366978.5, (1768,  6, 19), (1768, 24, 7), ( 1182,  2,  3)],
[664224, monday,    2385648.5, (1819,  8,  2), (1819, 31, 1), ( 1234, 10, 10)],
[671401, wednesday, 2392825.5, (1839,  3, 27), (1839, 13, 3), ( 1255,  1, 11)],
[694799, sunday,    2416223.5, (1903,  4, 19), (1903, 16, 7), ( 1321,  1, 21)],
[704424, sunday,    2425848.5, (1929,  8, 25), (1929, 34, 7), ( 1348,  3, 19)],
[708842, monday,    2430266.5, (1941,  9, 29), (1941, 40, 1), ( 1360,  9,  8)],
[709409, monday,    2430833.5, (1943,  4, 19), (1943, 16, 1), ( 1362,  4, 13)],
[709580, thursday,  2431004.5, (1943, 10,  7), (1943, 40, 4), ( 1362, 10,  7)],
[727274, tuesday,   2448698.5, (1992,  3, 17), (1992, 12, 2), ( 1412,  9, 13)],
[728714, sunday,    2450138.5, (1996,  2, 25), (1996,  8, 7), ( 1416, 10,  5)],
[744313, wednesday, 2465737.5, (2038, 11, 10), (2038, 45, 3), ( 1460, 10, 12)],
[764652, sunday,    2486076.5, (2094,  7, 18), (2094, 28, 7), ( 1518,  3,  5)]
]

RDN = [ 25469, 49217, 171307, 210155, 253427, 369740, 400085, 434355, 452605, 470160,
    473837, 507850, 524156, 544676, 567118, 569477, 601716, 613424, 626596, 645554,
    664224, 671401, 694799, 704424, 708842, 709409, 709580, 727274, 728714, 744313, 764652]

@testset "JulianDayNumbers" begin
    for rd in RD
        rdn = rd[1]
        local jn = ConvertOrdinalDate(rdn + 2, EN, JN)
        CEdate = DateFromDayNumber(CE, rdn + 2)
        IDdate = DateFromDayNumber(ID, rdn + 2)
        @test Int64(rd[3] + 0.5) == jn  # Note: +0.5
        @test rd[4] == (CEdate[2], CEdate[3], CEdate[4])
        @test rd[5] == (IDdate[2], IDdate[3], IDdate[4])

        # We do rnot handle preleptic Ismlam dates.
        if rdn >= 253427
            AHdate = DateFromDayNumber(AH, rdn + 2)
            @test rd[6] == (AHdate[2], AHdate[3], AHdate[4])
            println(rdn, " ", jn, " -> ", CDateStr(AHdate))
        end
        println(rdn, " ", jn, " -> ", CDateStr(CEdate), " / ", CDateStr(IDdate))
        println(WeekDays[rd[2]], WeekDay(rdn))
    end
end
end

function ShowProfileYear()
    ProfileYearAsEuropean(EC, 2022, true)
    ProfileYearAsEuropean(CE, 2022, true)
    ProfileYearAsEuropean(JD, 2022, true)
    ProfileYearAsEuropean(AM, 5783, true)
    ProfileYearAsEuropean(AH, 1444, true)
    ProfileYearAsEuropean(ID, 2022, true)
end

function TestStartEndYear()

    RHstart = [  # Dates for Rosh Hashana
    ((CE, 1842, 09, 05), (AM, 5603, 7, 1)),
    ((CE, 2020, 09, 19), (AM, 5781, 7, 1)),
    ((CE, 2021, 09, 07), (AM, 5782, 7, 1)),
    ((CE, 2022, 09, 26), (AM, 5783, 7, 1)),
    ((CE, 2023, 09, 16), (AM, 5784, 7, 1)),
    ((CE, 2024, 10, 03), (AM, 5785, 7, 1)),
    ((CE, 2025, 09, 23), (AM, 5786, 7, 1)),
    ((CE, 2026, 09, 12), (AM, 5787, 7, 1)),
    ((CE, 2027, 10, 02), (AM, 5788, 7, 1)),
    ((CE, 2043, 10, 05), (AM, 5804, 7, 1))
    ]

    RHend = [
    ((CE, 1843, 09, 24), (AM, 5603, 6, 29)),
    ((CE, 2021, 09, 06), (AM, 5781, 6, 29)),
    ((CE, 2022, 09, 25), (AM, 5782, 6, 29)),
    ((CE, 2023, 09, 15), (AM, 5783, 6, 29)),
    ((CE, 2024, 10, 02), (AM, 5784, 6, 29)),
    ((CE, 2025, 09, 22), (AM, 5785, 6, 29)),
    ((CE, 2026, 09, 11), (AM, 5786, 6, 29)),
    ((CE, 2027, 10, 01), (AM, 5787, 6, 29)),
    ((CE, 2028, 09, 20), (AM, 5788, 6, 29)),
    ((CE, 2044, 09, 21), (AM, 5804, 6, 29))
    ]

    @testset "RoshHashana" begin
    for rh in RHstart
        d = ConvertDate(rh[1], AM)
        @test rh[2] == d
    end
    for rh in RHend
        d = ConvertDate(rh[1], AM)
        @test rh[2] == d
    end
    end
end

function TestAll()
    TestDateGregorian(); println()
    TestDateJulian();    println()
    TestDateHebrew();    println()
    TestDateIslamic();   println()
    TestDateIso();       println()
    TestIso();           println()
    TestConversions();   println()
    TestDateTables();    println()
    TestDuration();      println()
    TestJulianDayNumbers1(); println()
    TestJulianDayNumbers(); println()
    TestDayOfYear();     println()
    TestStartEndYear();  println()
    ShowDayOfYear();     println()
    ShowProfileYear();   println()
    TestToday();         println()
    # TestSaveCalendars()
end

TestAll()


text:
using Documenter
using ProbabilisticCircuits

# The DOCSARGS environment variable can be used to pass additional arguments to make.jl.
# This is useful on CI, if you need to change the behavior of the build slightly but you
# can not change the .travis.yml or make.jl scripts any more (e.g. for a tag build).
if haskey(ENV, "DOCSARGS")
    for arg in split(ENV["DOCSARGS"])
        (arg in ARGS) || push!(ARGS, arg)
    end
end

makedocs(
    sitename = "ProbabilisticCircuits.jl",
    format = Documenter.HTML(
        # Use clean URLs, unless built as a "local" build
        prettyurls = !("local" in ARGS),
        canonical = "https://juice-jl.github.io/ProbabilisticCircuits.jl/stable/",
        assets = ["assets/favicon.ico"],
        analytics = "UA-136089579-2",
        highlights = ["yaml"],
    ),
    doctest = true,
    modules = [ProbabilisticCircuits],
    linkcheck_ignore = [
        # We'll ignore links that point to GitHub's edit pages, as they redirect to the
        # login screen and cause a warning:
        r"https://github.com/([A-Za-z0-9_.-]+)/([A-Za-z0-9_.-]+)/edit(.*)"
    ],
    pages = [
        "Home" => "index.md",
        "Manual" => Any[
            "Installation" => "manual/installation.md",
            "Examples" => "manual/examples.md"
        ],
        "API" => Any[
            "Public" => "api/public.md",
            "Internals" => map(
                s -> "api/internals/$(s)",
                sort(readdir(joinpath(@__DIR__, "src/api/internals")))
            ),
        ],
    ],
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    target = "build",
    repo = "github.com/Juice-jl/ProbabilisticCircuits.jl.git",
    branch = "gh-pages",
    devbranch = "master",
    devurl = "dev",
    versions = ["stable" => "v^", "v#.#"],
)


text:
using NodeJS

cd(joinpath(dirname(@__FILE__), "..", "packages")) do
    node = NodeJS.nodejs_cmd()
    npm = NodeJS.npm_cmd()

    run(`$npm install --scripts-prepend-node-path=auto --unsafe-perm`)
    if haskey(ENV, "WEBIO_WEBPACK_ARGS")
        args = [ENV["WEBIO_WEBPACK_ARGS"]]
    else
        args =[]
    end
    run(`$npm run build-prod --scripts-prepend-node-path=auto --unsafe-perm -- $args`)
end