text:
module YaoTensorNetwork

using BitBasis
using Yao, LuxurySparse
using DelimitedFiles
using Random
using Requires
using YaoBlocks.Optimise

const control_seperated = Ref(false)

export seperatecontrol
"""
    seperatecontrol([val::Bool])

Global switch, seperate control into rank tree tensors if true.
set value if `val` provided.
"""
function seperatecontrol(val::Bool)
    control_seperated[] = val
end

seperatecontrol() = control_seperated[]

include("utils.jl")
include("EinGraph.jl")
include("GraphBuilder.jl")
include("convert.jl")
include("tobasic_ext.jl")

@init @require OMEinsum = "ebe7aa44-baf0-506c-a96f-8464559b3922" include("contract.jl")

end # module


text:
export singleLayer,getSingleLayer

mutable struct singleLayer{T, TK <: AbstractConvKernel{T}, TN <: Union{batchNormNN{T}, normLayer{T}}} <: AbstractMeganetElement{T}
        activation  :: Function # activation function
        activation! :: Function # in place activation function
        K           :: TK # transformation type
        nLayer      :: TN # normalization layer
        Bin         :: Array{T} # bias inside nonlinearity
        Bout        :: Array{T} # bias outside nonlinearity

end

function getSingleLayer(TYPE::Type, K,nLayer;Bin=zeros(TYPE,nFeatOut(K),0),Bout=zeros(TYPE,nFeatOut(K),0),
                        activation=tanhActivation,activation_inplace=tanhActivation!)
	singleLayer(activation,activation_inplace,K,nLayer,Bin,Bout);
end

function splitWeights(this::singleLayer{T},theta::Array{T}) where {T <: Number}
    th1 = theta[1:nTheta(this.K)]
    cnt = length(th1)
    th2 = theta[cnt+(1:size(this.Bin,2))]
    cnt += length(th2)
    th3 = theta[cnt+(1:size(this.Bout,2))]
    cnt += length(th3)
    th4 = theta[cnt+(1:nTheta(this.nLayer))]
    cnt += length(th4)
    if cnt!=length(theta); error("splitWeights: length does not match"); end
    return th1, th2, th3, th4
end

function apply(this::singleLayer{T},theta::Array{T},Yin::Array{T},tmp,doDerivative=false) where {T <: Number}

    if isempty(tmp)
        tmp = Array{Any}(2)
        tmp[1] = Array{Any}(0)
        tmp[2] = Array{Any}(0)
    end
    nex = div(length(Yin),nFeatIn(this))
    Y   = reshape(Yin,:,nex)
    th1,th2,th3,th4 = splitWeights(this,theta)

    Yout::Array{T,2}     =  getOp(this.K,th1)*Y
    Yout .+= this.Bin * th2
    Yout,dummy,tmp[1] = apply(this.nLayer,th4,Yout,tmp[1],doDerivative)

    Yout,tmp[2]  = this.activation!(Yout,tmp[2],doDerivative)

    Yout .+= this.Bout*th3
    Ydata = Yout
    return Ydata, Yout, tmp
end

function nTheta(this::singleLayer)
    return nTheta(this.K)+size(this.Bin,2) + size(this.Bout,2) + nTheta(this.nLayer)
end

function nFeatIn(this::singleLayer)
    return nFeatIn(this.K)
end

function nFeatOut(this::singleLayer)
    return nFeatOut(this.K)
end

function nDataOut(this::singleLayer)
    return nFeatOut(this.K)
end

function initTheta(this::singleLayer{T}) where {T <: Number}
    return [vec(initTheta(this.K)); convert(T,0.01)*ones(T,size(this.Bin,2),1) ; convert(T,0.01)*ones(T,size(this.Bout,2),1); initTheta(this.nLayer) ]
end

function Jthetamv(this::singleLayer{T},dtheta::Array{T},theta::Array{T},Yin::Array{T},tmp) where {T <: Number}
    dA::Array{T,2}             = tmp[2]
    nex            = div(length(Yin),nFeatIn(this))
    Y              = reshape(Yin,:,nex)

    th1,th2,th3,th4     = splitWeights(this,theta)
    dth1,dth2,dth3,dth4 = splitWeights(this,dtheta)
    dZ::Array{T,2} = Jthetamv(this.K,dth1,th1,Y) .+ this.Bin*dth2
    Kop = getOp(this.K,th1)
    dZ  = Jmv(this.nLayer,dth4,dZ,th4,Kop*Y.+this.Bin*th2,tmp[1])[2]
    dZ .*= dA
    dZ .+= this.Bout*dth3
    return dZ, dZ
end

function JYmv(this::singleLayer{T},dYin::Array{T},theta::Array{T},Y::Array{T},tmp) where {T <: Number}
    dA  = tmp[2]
    nex = div(length(dYin),nFeatIn(this))
    th1,th2,th3,th4 = splitWeights(this,theta)
    Kop         = getOp(this.K,th1)
    dY          = reshape(dYin,:,nex)
    dZ          = Kop*dY
    dZ          = JYmv(this.nLayer,dZ,th4,Kop*Y.+this.Bin*th2,tmp[1])[2]
    # dZ          = dA.*dZ
    dZ .*= dA
    return dZ,dZ
end

function Jmv(this::singleLayer{T},dtheta::Array{T},dYin::Array{T},theta::Array{T},Yin::Array{T},tmp) where {T <: Number}
    dA::Array{T,2}  = tmp[2]
    nex = div(length(Yin),nFeatIn(this))
    th1,th2,th3,th4 = splitWeights(this,theta)
    dth1,dth2,dth3,dth4 = splitWeights(this,dtheta)

    dY = reshape(dYin,:,nex);
    Kop = getOp(this.K,th1)
    dZ::Array{T, 2} = Kop*dY;

    Y   = reshape(Yin,:,nex);
    dZ += Jthetamv(this.K,dth1,th1,Y) .+ this.Bin*dth2
    dZ = Jmv(this.nLayer,dth4,dZ,th4,Kop*Y.+this.Bin*th2,tmp[1])[2]

    dZ .*= dA
    dZ .+= this.Bout*dth3
    return dZ,dZ
end

function JTmv(this::singleLayer{T},Zin::Array{T},dummy::Array{T},theta::Array{T},Y::Array{T},tmp) where {T <: Number}
    dA::Array{T,2}   = tmp[2]
    nex  = div(length(Y),nFeatIn(this))
    Z    = reshape(Zin,:,nex)
    th1,th2,th3,th4  = splitWeights(this,theta)
    Kop = getOp(this.K,th1)

    dth3      = vec(sum(this.Bout'*Z,2))
    dAZ       = dA.*Z
    dth4,dAZ  = JTmv(this.nLayer,dAZ,zeros(T,0),th4,Kop*Y.+this.Bin*th2,tmp[1]) # this not type stable
    dth2      = vec(sum(this.Bin'*reshape(dAZ,:,nex),2))
    dth1      = JthetaTmv(this.K, dAZ,theta,Y) # this not type stable

    dY   = Kop'*reshape(dAZ,:,nex)
    dtheta = [vec(dth1); vec(dth2); vec(dth3); vec(dth4)]

    return dtheta, dY

end

function JthetaTmv(this::singleLayer{T},Zin::Array{T},dummy::Array{T},theta::Array{T},Y::Array{T},tmp) where {T <: Number}
    dA        = tmp[2]
    nex       = div(length(Zin),nFeatOut(this))
    th1,th2,th3,th4  = splitWeights(this,theta)

    Z         = reshape(Zin,:,nex);
    dAZ       = dA.*Z;
    dth3      = vec(sum(this.Bout'*Z,2));
    Kop       = getOp(this.K,th1)
    dth4,dAZ  = JTmv(this.nLayer,dAZ,(T)[],th4,Kop*Y.+this.Bin*th2,tmp[1])
    dth1      = JthetaTmv(this.K,dAZ,theta,Y);
    dth2      = vec(sum(this.Bin'*reshape(dAZ,:,nex),2));
    return [vec(dth1); vec(dth2); vec(dth3); vec(dth4)];
end

function JYTmv(this::singleLayer{T},Zin::Array{T},dummy::Array{T},theta::Array{T},Y::Array{T},tmp) where {T <: Number}
    dA::Array{T,2}   = tmp[2]
    nex  = div(length(Y),nFeatIn(this))
    th1,th2,th3,th4 = splitWeights(this,theta)
    Kop = getOp(this.K,th1)
    Z    = reshape(Zin,:,nex)
    dAZ::Array{T,2}  = dA.*Z
    dAZ = JYTmv(this.nLayer,dAZ,(T)[],th4,Kop*Y.+this.Bin*th2,tmp[1])
    ret::Array{T,2} = Kop'*reshape(dAZ,:,nex)
    return ret #TODO: @lars or eldad rename this variable as I'm not sure what to call it
end

ext: jl


text:
#####
##### Convinience functions
#####

unpack_grid(grid) = grid.Nx, grid.Ny, grid.Nz, grid.Lx, grid.Ly, grid.Lz

#####
##### Input validation
#####

"""Validate that an argument tuple is the right length and has elements of type `argtype`."""
function validate_tupled_argument(arg, argtype, argname)
    length(arg) == 3        || throw(ArgumentError("length($argname) must be 3."))
    all(isa.(arg, argtype)) || throw(ArgumentError("$argname=$arg must contain $argtype s."))
    all(arg .> 0)           || throw(ArgumentError("Elements of $argname=$arg must be > 0!"))
    return nothing
end

function validate_grid_size_and_length(sz, len, halo, x, y, z)
    validate_tupled_argument(sz, Integer, "size")
    validate_tupled_argument(halo, Integer, "halo")

    # Find domain endpoints or domain length, depending on user input:
    if !isnothing(len) # the user has specified a length!
        (!isnothing(x) || !isnothing(y) || !isnothing(z)) &&
            throw(ArgumentError("Cannot specify both length and x, y, z keyword arguments."))

        validate_tupled_argument(len, Number, "length")

        Lx, Ly, Lz = len

        # An "oceanic" default domain
        x = (0, Lx)
        y = (0, Ly)
        z = (-Lz, 0)

    else # isnothing(length) === true implies that user has not specified a length
        (isnothing(x) || isnothing(y) || isnothing(z)) &&
            throw(ArgumentError("Must supply length or x, y, z keyword arguments."))

        function coord2xyz(c)
            c == 1 && return "x"
            c == 2 && return "y"
            c == 3 && return "z"
        end

        for (i, c) in enumerate((x, y, z))
            name = coord2xyz(i)
            length(c) == 2       || throw(ArgumentError("$name length($c) must be 2."))
            all(isa.(c, Number)) || throw(ArgumentError("$name=$c should contain numbers."))
            c[2] >= c[1]         || throw(ArgumentError("$name=$c should be an increasing interval."))
        end

        Lx = x[2] - x[1]
        Ly = y[2] - y[1]
        Lz = z[2] - z[1]
    end

    return Lx, Ly, Lz, x, y, z
end

@inline get_grid_spacing(z::Function, k) = z(k)
@inline get_grid_spacing(z::AbstractArray{T,1}, k) where T = z[k]

function generate_variable_grid_spacings_from_zF(FT, zF_source, Nz)
    zF  = zeros(FT, Nz+1)
    zC  = zeros(FT, Nz)
    ΔzF = zeros(FT, Nz)
    ΔzC = zeros(FT, Nz-1)

    for k in 1:Nz+1
        zF[k] = get_grid_spacing(zF_source, k)
    end

    for k in 1:Nz
        zC[k] = (zF[k] + zF[k+1]) / 2
        ΔzF[k] = zF[k+1] - zF[k]
    end

    for k in 1:Nz-1
        ΔzC[k] = zC[k+1] - zC[k]
    end

    return zF, zC, ΔzF, ΔzC
end

function validate_and_generate_variable_grid_spacing(FT, zF_source, Nz, z₁, z₂)
    isnothing(zF_source) && throw(ArgumentError("Must pass zF to VerticallyStretchedCartesianGrid"))

    zF, zC, ΔzF, ΔzC = generate_variable_grid_spacings_from_zF(FT, zF_source, Nz)

    !isapprox(zF[1],   z₁) && throw(ArgumentError("Bottom face zF[1]=$(zF[1]) must equal bottom endpoint z₁=$z₁"))
    !isapprox(zF[end], z₂) && throw(ArgumentError("Top face zF[end]=$(zF[end]) must equal top endpoint z₂=$z₂"))

    return zF, zC, ΔzF, ΔzC
end


text:
using Images

using PCG.Types
using PCG.Geometry
using PCG.Universes
using PCG.Topologies.HexGreedTopologies
using PCG.Neighborhoods
using PCG.Operations
using PCG.Spaces
using PCG.Recorders.GreedImages
using PCG.Recorders.HexGreedImages
using PCG.Recorders.TurnsLogger

const RADIUS = 40
const CELL_SIZE = Size(5, 5)
const TURNS = 100
const DEBUG = false

struct State <: Checkable
    value::Int64
end

function Operations.check(element::Element, parameters::State)
    return element.properties.state == parameters
end

const DEAD = State(1)
const ALIVE = State(2)

struct Properties <: AbstractProperties
    state::State
end

const SPRITE_ALIVE_CELL = HexSprite(RGBA(1, 1, 1, 1), CELL_SIZE)
const SPRITE_DEAD_CELL = HexSprite(RGBA(0, 0, 0, 1), CELL_SIZE)

function GreedImages.choose_sprite(recorder::GreedImageRecorder, element)
    if element |> ALIVE |> exists
        return SPRITE_ALIVE_CELL
    end

    # in case of dead cell
    return SPRITE_DEAD_CELL
end

function process(universe::Universe, turns::Int64)

    ring = Neighborhood(universe.topology, ring_distance)

    universe() do element
        if element |> Fraction(0.2) |> exists
            element << (state=ALIVE,)
        end
    end

    universe(turns=turns) do element

        if element |> ALIVE |> ring() |> ALIVE |> count ∉ 2:3
            element << (state=DEAD,)
        end

        if element |> DEAD |> ring() |> ALIVE |> count == 3
            element << (state=ALIVE,)
        end

    end

end

topology = HexGreedTopology(RADIUS)

universe = initialize(topology,
                      Properties(DEAD),
                      DEBUG ? Recorder[] : [TurnsLoggerRecorder(),
                                            GreedImageRecorder(CELL_SIZE, 100, "output.gif")])

@time process(universe, 1)

@time process(universe, TURNS)

finish_recording!(universe)

max_stars_repo_name: Tiendil/pcg-julia

max_issues_repo_name: Tiendil/pcg-julia

max_forks_repo_name: Tiendil/pcg-julia


text:
"""
The edge-wise objective function to optimise. These values are natively supported:

* `Load`: the load of each edge is considered.
* `KleinrockLoad`: Kleinrock function is applied on the load.
* `FortzThorupLoad`: Fortz-Thorup linearisation of Kleinrock function is applied on the load.

Kleinrock function is defined as:

``K(x) = \\frac{x}{1-x}``
"""
abstract type EdgeWiseObjectiveFunction end

supports_min(::EdgeWiseObjectiveFunction) = false
supports_max(::EdgeWiseObjectiveFunction) = false

struct Load <: EdgeWiseObjectiveFunction end
struct FortzThorupLoad <: EdgeWiseObjectiveFunction end

struct KleinrockLoad <: EdgeWiseObjectiveFunction
    use_nonconvex_bilinear_formulation_for_equality::Bool

    KleinrockLoad(b::Bool=false) = new(b)

    # Solver is supposed to be able to handle bilinear constraints!
    # Outside generic global/MINLP solvers, there is just Gurobi to handle this:
    #     optimizer_with_attributes(Gurobi.Optimizer, "NonConvex" => 2)
end

supports_min(::Union{Load, KleinrockLoad, FortzThorupLoad}) = true
supports_max(::Load) = true
supports_max(k::KleinrockLoad) = k.use_nonconvex_bilinear_formulation_for_equality
supports_max(::FortzThorupLoad) = false

struct AlphaFairness <: EdgeWiseObjectiveFunction
    α::Float64
    force_power_cone::Bool # When there is a choice between a SOCP-based and a power-cone-based formulation, use the latter.

    AlphaFairness() = new(1.0, false)
    function AlphaFairness(α::Float64, force_power_cone::Bool=false)
        if force_power_cone
            if α == 0.0
                @warn "The power-cone formulation is forced to be used, " *
                      "but it is not available for α = $(α). " *
                      "The parameter `force_power_cone` is therefore ignored."
            elseif ! (α in [0.5, 1.5, 2.0]) # Easily SOCP-representable.
                @warn "The power-cone formulation is forced to be used, " *
                      "but it is the only available one for α = $(α). " *
                      "The parameter `force_power_cone` is therefore ignored."
            end
        end
        return new(α, force_power_cone)
    end
end

supports_min(af::AlphaFairness) = af.α == 0.0 # Only for LP model.
supports_max(::AlphaFairness) = true

"""
The aggregation-function objective to optimise, giving one objective function aggregating
all the components for each edge (indicated by `EdgeWiseObjectiveFunction`). These values are natively supported:

* `MinimumTotal`: the total considered value is minimised.
* `MinimumMaximum`: the maximum considered value is minimised.
* `MinMaxFair`: a min-max-fair solution is computed with respect to the considered value.
  By definition, this means that each user has their own value; not a single user can reduce
  this value without others seeing an increase.

In theory, a max-min fair solution is kept when any nondecreasing function is applied on the criteria.
It becomes a min-max fair solution if the function is nonincreasing. In practice, it means that
using `MinMaxFair` with `KleinrockLoad` or `Load` should give the same solution. The solver will
not forbid such combinations, though.
"""
abstract type AggregationObjectiveFunction end

abstract type MinimisingAggregationObjectiveFunction <: AggregationObjectiveFunction end
abstract type MaximisingAggregationObjectiveFunction <: AggregationObjectiveFunction end

struct MinimumTotal <: MinimisingAggregationObjectiveFunction end
struct MinimumMaximum <: MinimisingAggregationObjectiveFunction end
struct MinMaxFair <: MinimisingAggregationObjectiveFunction
    ε::Float64 # Numerical accuracy when adding MMF constraints.

    MinMaxFair(ε::Float64) = new(ε)
    MinMaxFair() = new(1.0e-5)
end

struct MaximumTotal <: MaximisingAggregationObjectiveFunction end
struct MaximumMinimum <: MaximisingAggregationObjectiveFunction end
struct MaxMinFair <: MaximisingAggregationObjectiveFunction
    ε::Float64 # Numerical accuracy when adding MMF constraints.

    MaxMinFair(ε::Float64) = new(ε)
    MaxMinFair() = new(1.0e-5)
end

"""
The type of formulation that is used to solve the routing problem. This package natively supports two main families of formulations:

* `FlowFormulation`: each edge is explicitly present in the model
* `PathFormulation`: the model uses path variables. The path-based formulations can be solved in two ways:
  * precomputing paths: only the set of precomputed paths is considered; if not enough paths are precomputed,
    the solution is only approximate, and no approximation factor is guaranteed
  * column generation: new paths are added on the fly, when they might improve the solution
"""
abstract type FormulationType end

struct FlowFormulation <: FormulationType end
struct PathFormulation <: FormulationType end

"""
For robust and oblivious uncertainties, this package natively offers multiple algorithm choices:

* `Automatic`: let the package decide which algorith to use
* `CuttingPlane`: use an iterative cutting-plane algorithm (solves many small problems)
* `DualReformulation`: use a dual-reformulation of the uncertainty set (solves one large problem)

For other types of uncertainty, use the value `Automatic`.
"""
abstract type AlgorithmChoice end

struct Automatic <: AlgorithmChoice end
struct CuttingPlane <: AlgorithmChoice end
struct DualReformulation <: AlgorithmChoice end
# TODO: Different algorithms for MMF?

"""
The way uncertainty is implemented. These values are natively supported:

* `NoUncertaintyHandling`: the uncertainty is not taken into account, only nominal values are used
* `StochasticUncertainty`: the uncertainty is implemented through stochastic optimisation,
  i.e. with several (potentially weighted) values of the unknown parameters
* `RobustUncertainty`: the uncertainty is implemented through robust optimisation,
  i.e. the unknown parameters belong to a given (bounded) uncertainty set
* `ObliviousUncertainty`: the uncertainty is implemented through oblivious optimisation,
  i.e. the unknown parameters belong to an unbouded uncertainty set
"""
abstract type UncertaintyHandling end

struct NoUncertaintyHandling <: UncertaintyHandling end
struct StochasticUncertainty <: UncertaintyHandling end
struct RobustUncertainty <: UncertaintyHandling end
struct ObliviousUncertainty <: UncertaintyHandling end

"""
The parameters of the problem that can see uncertainty. These values are natively supported:

* `NoUncertainty`: there is no uncertainty in the problem
* `UncertainDemand`: the demand to route are not fully known, but the topology has no uncertainty
* `UncertainCapacity`: the links have unknown capacities, but the demands to route are fully known
"""
abstract type UncertainParameters end

struct NoUncertainty <: UncertainParameters end
struct UncertainDemand <: UncertainParameters end
struct UncertainCapacity <: UncertainParameters end

"""
Parameterising the way a routing-optimisation problem must be solved. Seven parameters must be chosen (independently
of each other):

* an edge-wise objective function, which assigns a value for each edge, of type `EdgeWiseObjectiveFunction`
* an aggregation technique to build the complete objective function from the edge-wise pieces, of type `AggregationObjectiveFunction`
* a formulation type, of type `FormulationType` (for path-based formulations, column generation can be separately enabled)
* whether column generation is enabled, only for path-based formulations
* a way to heandle the uncertainty, of type `UncertaintyHandling`
* a precise algorithm to solve the instance, of type `AlgorithmChoice`
* a definition of which parameters can be afflicted by uncertainty, of type `UncertainParameters`

The exact algorithm used to solve the instance is based on those seven parameters.
"""
struct ModelType
    # TODO: integer vs. continuous flows?
    edge_obj::EdgeWiseObjectiveFunction
    agg_obj::AggregationObjectiveFunction
    type::FormulationType
    cg::Bool
    algo::AlgorithmChoice
    unc::UncertaintyHandling
    uncparams::UncertainParameters

    function ModelType(edge_obj::EdgeWiseObjectiveFunction,
                       agg_obj::AggregationObjectiveFunction,
                       type::FormulationType=FlowFormulation(),
                       cg::Bool=false,
                       algo::AlgorithmChoice=Automatic(),
                       unc::UncertaintyHandling=NoUncertaintyHandling(),
                       uncparams::UncertainParameters=NoUncertainty())
        is_obj_load = edge_obj in [Load(), KleinrockLoad(), FortzThorupLoad()]
        if is_obj_load && uncparams == UncertainCapacity()
            error("Load-related edge-wise objective functions like $(edge_obj) are not supposed to be used with capacity uncertainty.")
        end

        if type == FlowFormulation() && cg
            error("Column generation is not defined for flow formulations")
        end

        if (unc != ObliviousUncertainty() && unc != RobustUncertainty()) && algo != Automatic()
            # TODO: MMF?
            error("Different algorithms are only available for robust- and oblivious-uncertainty models.")
        end
        if (unc == ObliviousUncertainty() || unc == RobustUncertainty()) && algo == Automatic()
            algo = DualReformulation() # Usually the most efficient one. Could implement heuristics based on numerical experiments.
        end

        if (unc == NoUncertaintyHandling() && uncparams != NoUncertainty()) || (unc == NoUncertaintyHandling() && uncparams != NoUncertainty())
            text = "When there is no uncertainty, the parameters UncertaintyHandling"
            if unc != NoUncertaintyHandling() && uncparams == NoUncertainty()
                text *= " ($unc != $(NoUncertaintyHandling()))"
            end
            text *= " and UncertainParameters"
            if unc == NoUncertaintyHandling() && uncparams != NoUncertainty()
                text *= " ($uncparams != $(NoUncertainty()))"
            end
            text *= " must match."
            error(text)
        end

        if uncparams == UncertainCapacity()
            error("No model implemented uncertainty on link capacity")
        end

        # Check that the parameter detection worked.
        if unc == ObliviousUncertainty() || unc == RobustUncertainty()
            @assert algo != Automatic()
        end

        return new(edge_obj, agg_obj, type, cg, algo, unc, uncparams)
    end
end

function Base.copy(mt::ModelType)
    return ModelType(mt.edge_obj, mt.agg_obj, mt.type, mt.cg, mt.algo, mt.unc, mt.uncparams)
end


text:
using MOCNeutronTransport
@testset "Quadrilateral_3D" begin
    for T in [Float32, Float64]
        @testset "Constructors" begin
            p₁ = Point_3D(T, 0)
            p₂ = Point_3D(T, 1)
            p₃ = Point_3D(T, 1, 1)
            p₄ = Point_3D(T, 0, 1)
            quad = Quadrilateral_3D((p₁, p₂, p₃, p₄))
            @test quad.points == (p₁, p₂, p₃, p₄)

            # single constructor
            quad = Quadrilateral_3D(p₁, p₂, p₃, p₄)
            @test quad.points == (p₁, p₂, p₃, p₄)
        end

        @testset "Methods" begin
            p₁ = Point_3D(T, 0)
            p₂ = Point_3D(T, 1)
            p₃ = Point_3D(T, 1, 1)
            p₄ = Point_3D(T, 0, 1)
            quad = Quadrilateral_3D((p₁, p₂, p₃, p₄))

            # interpolation
            quad(0, 0) ≈ p₁
            quad(1, 0) ≈ p₂
            quad(0, 1) ≈ p₃
            quad(1, 1) ≈ p₄
            quad(1//2, 1//2) ≈ Point_3D(T, 1//2, 1//2)

            # area
            a = area(quad)
            @test typeof(a) == typeof(T(1))
            @test a == T(1)

            # intersect
            # line is not coplanar with Quadrilateral_3D
#            p₄ = Point_3D(T, 9//10, 1//10, -5)
#            p₅ = Point_3D(T, 9//10, 1//10,  5)
#            l = LineSegment_3D(p₄, p₅)
#            bool, point = intersect(l, quad)
#            @test bool
#            @test point ≈ Point_3D(T, 9//10, 1//10, 0)
#
#            # line is coplanar with Quadrilateral_3D
#            p₄ = Point_3D(T, 1//2, -1)
#            p₅ = Point_3D(T, 1//2,  2)
#            l = LineSegment_3D(p₄, p₅)
#            bool, point = intersect(l, quad)
#            @test !bool
#
#            # no intersection non-coplanar
#            p₄ = Point_3D(T, 2, 1//10, -5)
#            p₅ = Point_3D(T, 2, 1//10,  5)
#            l = LineSegment_3D(p₄, p₅)
#            bool, point = intersect(l, quad)
#            @test !bool
#
#            # no intersection coplanar
#            p₄ = Point_3D(T, 2, -1)
#            p₅ = Point_3D(T, 2,  2)
#            l = LineSegment_3D(p₄, p₅)
#            bool, point = intersect(l, quad)
#            @test !bool
        end
    end
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
export HestonParams, evaluateCharFunc, evaluateLogCharFunc, computeCumulants
export BlackParams, evaluateCharFuncAndDerivative
export CVCharFunc, HestonCVCharFunc

struct HestonParams{T}
    v0::T
    κ::T
    θ::T
    ρ::T
    σ::T
end

struct BlackParams{T}
    σ::T
end
DefaultCharFunc(params::HestonParams{Float64}) =
    DefaultCharFunc{HestonParams{Float64},Complex}(params)

function evaluateLogCharFunc(
    cf::CharFunc{BlackParams{T},CR},
    z::CT,
    τ::T,
)::CR where {T,CR,CT}
    cc1 = oneim(cf)
    p = cf.model
    return -p.σ^2 * τ / 2 * z * (z + cc1)
end

function evaluateLogCharFuncAndDerivative(
    cf::CharFunc{BlackParams{T},CR},
    z::CT,
    τ::T,
)::Tuple{CR,CR} where {T,CR,CT}
    cc1 = oneim(cf)
    p = cf.model
    phi = -p.σ^2 * τ / 2 * z * (z + cc1)
    phi_d = -p.σ^2 * τ / 2 * (2 * z + cc1)
    return phi, phi_d
end

struct CVCharFunc{MAINT,CONTROLT,CR} <: CharFunc{MAINT,CR}
    main::CharFunc{MAINT,CR}
    control::CharFunc{CONTROLT,CR}
end

model(cf::CVCharFunc) = model(cf.main)
oneim(cf::CVCharFunc) = oneim(cf.main)

HestonCVCharFunc(heston::CharFunc{HestonParams{T},CR}) where {T,CR} =
    CVCharFunc{HestonParams{T},BlackParams{T},CR}(
        heston,
        DefaultCharFunc{BlackParams{T},CR}(
            BlackParams{T}(sqrt(model(heston).v0))
        ),
    )

function evaluateCharFuncAndDerivative(
    p::CVCharFunc{MAINT,CONTROLT,CR},
    z::CT,
    τ::T,
)::Tuple{CR,CR} where {T,CR,CT,MAINT,CONTROLT}
    phi, phi_d = evaluateCharFuncAndDerivative(p.main, z, τ)
    phiB, phiB_d = evaluateCharFuncAndDerivative(p.control, z, τ)
    return (phi - phiB, phi_d - phiB_d)
end

#cinf(params::HestonCVParams{T}, τ::T) where {T} = cinf(params.heston, τ)

cinf(params::HestonParams{T}, τ::T) where {T} =
    (params.v0 + params.κ * params.θ * τ) / params.σ * sqrt(1 - params.ρ^2)

#CC the Nemo arbField or Complex type. z: Nemo complex number or Complex number.
@inline function evaluateCharFunc(
    p::CharFunc{MT,CR},
    z::CT,
    τ::T,
)::CR where {T,CR,CT,MT}
    return exp(evaluateLogCharFunc(p, z, τ))
end

@inline evaluateLogCharFunc(
    p::CharFunc{HestonParams{T},CR},
    z::CT,
    τ::T,
) where {T,CR,CT} = evaluateLogCharFuncAL(p, z, τ)

@inline function evaluateLogCharFuncAL(
    cf::CharFunc{HestonParams{T},CR},
    z::CT,
    τ::T,
)::CR where {T,CR,CT}
    #follows Andersen-Lake fast implementation
    p = model(cf)
    v0 = p.v0
    κ = p.κ
    θ = p.θ
    ρ = p.ρ
    σ = p.σ
    cc1 = oneim(cf)
    α = -(z * (z + cc1)) * σ^2
    β = κ - σ * ρ * z * cc1
    D = sqrt(β^2 - α)
    rr = (real(β) * real(D) + imag(β) * imag(D)) > 0 ? α / (β + D) : (β - D)
    y = D == 0 ? -τ / 2 : expm1(-D * τ) / (2 * D)
    l = log1p(-rr * y)
    A = (rr * τ - 2 * l) * (κ * θ / σ^2)
    B = z * (z + cc1) * y / (1 - rr * y)
    return A + B * v0
end

#Faster evaluation of the characteristic function phi(ix) where x is real. Useful for optimal alpha.
@inline function evaluateLogCharFuncAtImag(
    cf::CharFunc{HestonParams{T},CR},
    imz::T,
    τ::T,
)::T where {T,CR}
    #follows Andersen-Lake fast implementation
    p = model(cf)
    v0 = p.v0
    κ = p.κ
    θ = p.θ
    ρ = p.ρ
    σ = p.σ
    cc1 = oneim(cf)
    α = (imz * (imz + one(T))) * σ^2
    β = κ + σ * ρ * imz
    D2 = β^2 - α
    local ch, sh
    if D2 >= 0
        D = sqrt(D2)
        ch = cosh(D * τ / 2)
        sh = (D == 0) ? τ / 2 : sinh(D * τ / 2) / D
    else
        D = sqrt(-D2)
        sh, ch = sincos(D * τ / 2)
        sh /= D
    end
    A = κ * θ / σ^2 * (β * τ - log((ch + β * sh)^2))
    B = imz * (imz + 1) * sh / (ch + β * sh)
    return A + B * v0
end

@inline function evaluateCharFuncAndDerivative(
    p::CharFunc{MAINT,CR},
    z::CT,
    τ::T,
)::Tuple{CR,CR} where {T,CR,CT,MAINT}
    arg, arg_d = evaluateLogCharFuncAndDerivative(p, z, τ)
    phi = exp(arg)
    phi_d = arg_d * phi
    return phi, phi_d
end

@inline function evaluateLogCharFuncAndDerivative(
    cf::CharFunc{HestonParams{T},CR},
    z::CT,
    τ::T,
)::Tuple{CR,CR} where {T,CR,CT}
    #derivative towards real part of z
    p = model(cf)
    v0 = p.v0
    κ = p.κ
    θ = p.θ
    ρ = p.ρ
    σ = p.σ
    cc1 = oneim(cf)
    α = -(z * (z + cc1)) * σ^2
    α_d = -(2 * z + cc1) * σ^2
    β_d = -σ * ρ * cc1
    β = κ + β_d * z

    D = sqrt(β^2 - α)
    D_d = (2 * β * β_d - α_d) / (2 * D)
    local rr
    local rr_d
    if (real(β) * real(D) + imag(β) * imag(D)) > 0
        rr = α / (β + D)
        rr_d = α_d / (β + D) - rr * (β_d + D_d) / (β + D)
    else
        rr = (β - D)
        rr_d = β_d - D_d
    end
    local y
    local y_d
    if D != 0
        em = expm1(-D * τ)
        y = em / (2 * D)
        y_d = (y * (-2 * D_d * (D * τ + 1)) + (-D_d * τ)) / (2 * D) # ((em+1)*(-D_d * τ)*(2*D) - em*(2*D_d))/(2*D)^2
    else
        y = -τ / 2
        y_d = zero(T)
    end
    l = log1p(-rr * y)
    rry_d = (rr_d * y + rr * y_d)
    l_d = -rry_d / (1 - rr * y)
    A = (rr * τ - 2 * l) * (κ * θ / σ^2)
    A_d = (rr_d * τ - 2 * l_d) * (κ * θ / σ^2)
    Bnum = z * (z + cc1) * y
    Bdenom = (1 - rr * y)
    B = Bnum / Bdenom
    B_d = ((2 * z + cc1) * y + y_d * z * (z + cc1)) / Bdenom + rry_d * B / Bdenom
    return A + B * v0, A_d + B_d * v0
end

struct CuiCharFunc{MT,CR} <: CharFunc{MT,CR} #model type, return type (e.g. Complex or acb),
    model::MT
end

CuiCharFunc(params::HestonParams{Float64}) =
    CuiCharFunc{HestonParams{Float64},Complex,Type}(params, Complex)
@inline model(cf::CuiCharFunc) = cf.model
@inline field(cf::CuiCharFunc) = cf.field
@inline evaluateLogCharFunc(
    p::CuiCharFunc{HestonParams{T},CR},
    z::CT,
    τ::T,
) where {T,CR,CF,CT} = evaluateLogCharFuncAL(p, z, τ)

function evaluateLogCharFuncCui(
    cf::CharFunc{HestonParams{T},CR},
    z::CT,
    τ::T,
)::CR where {T,CR,CT}
    p = model(cf)

    v0 = p.v0
    κ = p.κ
    θ = p.θ
    ρ = p.ρ
    σ = p.σ
    cc1 = oneim(cf)
    iu = cc1 * z
    ξ = κ - σ * ρ * iu
    α = iu * (-iu + 1)
    d = sqrt(ξ^2 + σ^2 * α)
    dht = d * (τ / 2)
    sh1 = sinh(dht)
    ch1 = cosh(dht)
    A2v = (d * ch1 + ξ * sh1)
    A1 = α * sh1
    # edt = ch1 - sh1
    # logB = (κ * τ / 2 - dht) - log1p( (ξ - d) * sh1/ d * edt )

    edt = (ch1 + sh1)
    logB = (κ * τ / 2 - dht) - log(A2v / (d * edt))
    return -κ * θ * ρ * τ / σ * iu - A1 / A2v * v0 + 2 * κ * θ / σ^2 * logB
end

function computeCumulants(p::HestonParams{T}, τ::T) where {T}
    lambda = p.κ
    ubar = p.θ
    u0 = p.v0
    eta = p.σ
    rho = p.ρ
    term = τ
    if lambda == 0 && eta == 0
        c1 = -u0 * term / 2
        c2 = u0 * term
        c4 = 0
        return c1, c2, c4
    elseif lambda == 0
        c1 = -u0 * term / 2
        c2 = u0 * term * (1 + term * term * eta * (eta * term / 12 - rho / 2))
        c4 =
            u0 *
            eta^2 *
            term^3 *
            (
                2 * rho^2 - 2 * eta * term * rho +
                eta^4 * term^4 * 17 / 1680 +
                eta^4 * term^5 * rho^2 * 11 / 20 - eta * term * rho^3 / 2 +
                eta^2 * term^2 * 3 / 10 - eta^3 * term^3 * rho * 17 / 120 + 1
            )
        return c1, c2, c4
    end
    elt = exp(-lambda * term)
    c1 = (1 - elt) * (ubar - u0) / (2 * lambda) - ubar * term / 2
    lambda2 = lambda * lambda
    lambda3 = lambda * lambda2
    eta2 = eta * eta
    elt2 = elt * elt
    # correct c2
    c2 =
        1 / (8 * lambda3) * (
            2 *
            u0 *
            (
                lambda2 * 4 * (elt * term * rho * eta + 1 - elt) +
                lambda * (4 * eta * rho * (elt - 1) - 2 * elt * term * eta2) +
                eta2 * (1 - elt2)
            ) +
            ubar * (
                8 * lambda3 * term +
                lambda2 * 8 * (-eta * rho * term * (1 + elt) + elt - 1) +
                lambda * (2 * eta2 * term * (1 + 2 * elt) + eta * rho * 16 * (1 - elt)) +
                eta2 * (-5 + 4 * elt + elt2)
            )
        )
    # paper c2
    # c2 = 1 / (8 * lambda3) * (eta*term*lambda*elt*(u0-ubar)*(8*lambda*rho-4*eta) + lambda*rho*eta*(1-elt)*(16*ubar-8*u0) + 2*ubar*lambda*term*(-4*lambda*rho*eta+eta2+4*lambda2) + eta2*((ubar-2*u0)*elt2+ubar*(6*elt-7)+2*u0) + 8*lambda2*(u0-ubar)*(1-elt))

    lambda4 = lambda2 * lambda2
    lambda5 = lambda4 * lambda
    lambda6 = lambda4 * lambda2
    lambda7 = lambda3 * lambda4
    eta3 = eta * eta2
    term2 = term * term
    term3 = term2 * term
    rho2 = rho * rho
    eta4 = eta2 * eta2
    elt3 = elt2 * elt
    elt4 = elt2 * elt2
    rho3 = rho2 * rho
    c4 =
        -1 / lambda7 *
        2 *
        ubar *
        eta2 *
        (
            (
                (term3 * rho3 * eta - 3 * term2 * rho2) * lambda6 -
                3 *
                term *
                (term2 * rho2 * eta2 - 4 * term * rho * (rho2 + 1) * eta + 8 * rho2 + 2) *
                lambda5 / 2 +
                (
                    3 * term3 * rho * eta3 / 4 - 21 * term2 * (rho2 + 3 / 14) * eta2 / 2 +
                    (18 * term * rho3 + 24 * term * rho) * eta - 18 * rho2 - 3
                ) * lambda4 -
                (
                    term3 * eta3 - 42 * term2 * rho * eta2 +
                    (240 * term * rho2 + 54 * term) * eta - 192 * rho3 - 192 * rho
                ) *
                eta *
                lambda3 / 8 -
                3 *
                eta2 *
                (term2 * eta2 - 35 * term / 2 * rho * eta + 40 * rho2 + 15 / 2) *
                lambda2 / 4 - 27 * eta3 * (term * eta - 20 * rho / 3) * lambda / 16 -
                21 * eta4 / 16
            ) * elt +
            (
                (-3 / 4 + 3 * term * rho * eta - 3 * term2 * rho2 * eta2 / 2) * lambda4 +
                3 *
                eta *
                (term2 * rho * eta2 + (-4 * term * rho2 - 3 * term / 2) * eta + 4 * rho) *
                lambda3 / 2 -
                3 *
                eta2 *
                (term2 * eta2 - 14 * term * rho * eta + 20 * rho2 + 6) *
                lambda2 / 8 + (-15 / 16 * term * eta4 + 9 * eta3 * rho / 2) * lambda -
                21 * eta4 / 32
            ) * elt2 +
            3 *
            eta2 *
            (
                (term * rho * eta - 1) * lambda2 +
                (-term / 2 * eta2 + 2 * rho * eta) * lambda - eta2 / 2
            ) *
            elt3 / 8 - 3 * eta4 * elt4 / 128 +
            (-6 * term * rho2 - 3 * term / 2) * lambda5 +
            ((6 * term * rho3 + 9 * term * rho) * eta + 18 * rho2 + 15 / 4) * lambda4 -
            9 * eta * ((rho2 + 0.25) * term * eta + 8 * rho3 / 3 + 10 * rho / 3) * lambda3 +
            15 * eta2 * (term * rho * eta + 10 * rho2 + 11 / 5) * lambda2 / 4 +
            (-33 / 2 * eta3 * rho - 15 / 32 * term * eta4) * lambda +
            279 * eta4 / 128
        )
    c4 +=
        u0 / lambda7 * (
            2 *
            eta2 *
            (
                (
                    (term3 * rho3 * eta - 3 * term2 * rho2) * lambda6 -
                    3 *
                    term *
                    (
                        term2 * rho2 * eta2 - 2 * term * rho * (rho2 + 2) * eta +
                        4 * rho2 +
                        2
                    ) *
                    lambda5 / 2 +
                    (
                        3 * term3 * rho * eta3 / 4 - 6 * (rho2 + 3 / 8) * term2 * eta2 +
                        6 * term * rho * (rho2 + 2) * eta - 6 * rho2
                    ) * lambda4 -
                    eta *
                    (
                        term3 * eta3 - 24 * term2 * rho * eta2 +
                        (72 * term * rho2 + 18 * term) * eta - 48 * rho3
                    ) *
                    lambda3 / 8 -
                    3 * eta2 * (term2 * eta2 - 7 * term * rho * eta - 3) * lambda2 / 8 -
                    3 * eta3 * (term * eta + 10 * rho) * lambda / 16 + 3 * eta4 / 8
                ) * elt +
                (
                    (6 * term * rho * eta - 3 * term2 * rho2 * eta2 - 3 / 2) * lambda4 +
                    3 *
                    (
                        term2 * rho * eta2 +
                        (-3 * term * rho2 - 3 * term / 2) * eta +
                        3 * rho
                    ) *
                    eta *
                    lambda3 -
                    3 *
                    eta2 *
                    (term2 * eta2 - 10 * term * rho * eta + 12 * rho2 + 3) *
                    lambda2 / 4 - 9 * eta3 * (term * eta - 10 * rho / 3) * lambda / 8 -
                    3 * eta4 / 8
                ) * elt2 +
                9 *
                eta2 *
                (
                    (term * rho * eta - 1) * lambda2 +
                    (-term / 2 * eta2 + 5 / 3 * rho * eta) * lambda - eta2 / 3
                ) *
                elt3 / 8 - 3 * eta4 * elt4 / 32 -
                6 *
                ((rho2 + 1 / 4) * lambda2 - 5 * lambda * rho * eta / 4 + 5 * eta2 / 16) *
                (lambda * rho * eta - eta2 / 4 - lambda2)
            )
        )
    return c1, c2, c4
end


text:
using Nemo, GibbsTypePriors, JLD

grid_k = [collect(1:25:1000); 1000]

Vnk_numerical_accuracy = GibbsTypePriors.Vnk_NGG.(1000, grid_k, 1.2, 0.6) |> x -> accuracy_bits.(x)
save("test/graphical_tests/saves_for_graphical_tests/accuracy_Vnk_1000.jld", "Vnk_numerical_accuracy", Vnk_numerical_accuracy)

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors


text:
"""`lqr(A, B, Q, R)`

Calculate the optimal gain matrix `K` for the state-feedback law `u = K*x` that
minimizes the cost function:

J = integral(x'Qx + u'Ru, 0, inf).

For the continuous time model `dx = Ax + Bu`.

`lqr(sys, Q, R)`

Solve the LQR problem for state-space system `sys`. Works for both discrete
and continuous time systems.

See also `LQG`

Usage example:
```julia
using LinearAlgebra # For identity matrix I
A = [0 1; 0 0]
B = [0;1]
C = [1 0]
sys = ss(A,B,C,0)
Q = I
R = I
L = lqr(sys,Q,R)

u(x,t) = -L*x # Form control law,
t=0:0.1:5
x0 = [1,0]
y, t, x, uout = lsim(sys,u,t,x0=x0)
plot(t,x, lab=["Position" "Velocity"], xlabel="Time [s]")
```
"""
function lqr(A, B, Q, R)
    S = care(A, B, Q, R)
    K = R\B'*S
    return K
end

"""`kalman(A, C, R1, R2)` kalman(sys, R1, R2)`

Calculate the optimal Kalman gain

See also `LQG`

"""
kalman(A, C, R1,R2) = Matrix(lqr(A',C',R1,R2)')

function lqr(sys::StateSpace, Q, R)
    if iscontinuous(sys)
        return lqr(sys.A, sys.B, Q, R)
    else
        return dlqr(sys.A, sys.B, Q, R)
    end
end

function kalman(sys::StateSpace, R1,R2)
    if iscontinuous(sys)
        return Matrix(lqr(sys.A', sys.C', R1,R2)')
    else
        return Matrix(dlqr(sys.A', sys.C', R1,R2)')
    end
end

"""`dlqr(A, B, Q, R)`, `dlqr(sys, Q, R)`

Calculate the optimal gain matrix `K` for the state-feedback law `u[k] = K*x[k]` that
minimizes the cost function:

J = sum(x'Qx + u'Ru, 0, inf).

For the discrte time model `x[k+1] = Ax[k] + Bu[k]`.

See also `lqg`

Usage example:
```julia
using LinearAlgebra # For identity matrix I
h = 0.1
A = [1 h; 0 1]
B = [0;1]
C = [1 0]
sys = ss(A,B,C,0, h)
Q = I
R = I
L = dlqr(A,B,Q,R) # lqr(sys,Q,R) can also be used

u(x,t) = -L*x # Form control law,
t=0:h:5
x0 = [1,0]
y, t, x, uout = lsim(sys,u,t,x0=x0)
plot(t,x, lab=["Position"  "Velocity"], xlabel="Time [s]")
```
"""
function dlqr(A, B, Q, R)
    S = dare(A, B, Q, R)
    K = (B'*S*B + R)\(B'S*A)
    return K
end

"""`dkalman(A, C, R1, R2)` kalman(sys, R1, R2)`

Calculate the optimal Kalman gain for discrete time systems

"""
dkalman(A, C, R1,R2) = Matrix(dlqr(A',C',R1,R2)')

"""`place(A, B, p)`, `place(sys::StateSpace, p)`

Calculate gain matrix `K` such that
the poles of `(A-BK)` in are in `p`.

Uses Ackermann's formula."""
function place(A, B, p)
    n = length(p)
    n != size(A,1) && error("Must define as many poles as states")
    n != size(B,1) && error("A and B must have same number of rows")
    if size(B,2) == 1
        acker(A,B,p)
    else
        error("place only implemented for SISO systems")
    end
end

function place(sys::StateSpace, p)
    return place(sys.A, sys.B, p)
end

#Implements Ackermann's formula for placing poles of (A-BK) in p
function acker(A,B,P)
    n = length(P)
    #Calculate characteristic polynomial
    poly = mapreduce(p -> Poly([1, -p]), *, P, init=Poly(one(eltype(P))))
    q = zero(Array{promote_type(eltype(A),Float64),2}(undef, n,n))
    for i = n:-1:0
        q += A^(n-i)*poly[i]
    end
    S = Array{promote_type(eltype(A),eltype(B),Float64),2}(undef, n,n)
    for i = 0:(n-1)
        S[:,i+1] = A^i*B
    end
    return [zeros(1,n-1) 1]*(S\q)
end

"""
`feedback(L)` Returns L/(1+L)
`feedback(P1,P2)` Returns P1/(1+P1*P2)
"""
feedback(L::TransferFunction) = L/(1+L)
feedback(P1::TransferFunction, P2::TransferFunction) = P1/(1+P1*P2)

#Efficient implementations
function feedback(L::TransferFunction{T}) where T<:SisoRational
    if size(L) != (1,1)
        error("MIMO TransferFunction feedback isn't implemented, use L/(1+L)")
    end
    P = numpoly(L)
    Q = denpoly(L)
    tf(P, P+Q, L.Ts)
end

function feedback(L::TransferFunction{T}) where {T<:SisoZpk}
    if size(L) != (1,1)
        error("MIMO TransferFunction feedback isn't implemented, use L/(1+L)")
    end
    #Extract polynomials and create P/(P+Q)
    k = L.matrix[1].k
    denpol = numpoly(L)[1]+denpoly(L)[1]
    kden = denpol[end] # Get coeff of s^n
    # Create siso system
    sisozpk = T(L.matrix[1].z, roots(denpol), k/kden)
    return TransferFunction{T}(fill(sisozpk,1,1), L.Ts)
end

"""
`feedback(sys)`

`feedback(sys1,sys2)`

Forms the negative feedback interconnection
```julia
>-+ sys1 +-->
  |      |
 (-)sys2 +
```
If no second system is given, negative identity feedback is assumed
"""
function feedback(sys::Union{StateSpace, DelayLtiSystem})
    ninputs(sys) != noutputs(sys) && error("Use feedback(sys1, sys2) if number of inputs != outputs")
    feedback(sys,ss(Matrix{numeric_type(sys)}(I,size(sys)...)))
end

function feedback(sys1::StateSpace,sys2::StateSpace)
    sum(abs.(sys1.D)) != 0 && sum(abs.(sys2.D)) != 0 && error("There can not be a direct term (D) in both sys1 and sys2")
    A = [sys1.A+sys1.B*(-sys2.D)*sys1.C sys1.B*(-sys2.C); sys2.B*sys1.C  sys2.A+sys2.B*sys1.D*(-sys2.C)]
    B = [sys1.B; sys2.B*sys1.D]
    C = [sys1.C  sys1.D*(-sys2.C)]
    ss(A,B,C,sys1.D)
end

"""
`feedback2dof(P,R,S,T)` Return `BT/(AR+ST)` where B and A are the numerator and denomenator polynomials of `P` respectively
`feedback2dof(B,A,R,S,T)` Return `BT/(AR+ST)`
"""
function feedback2dof(P::TransferFunction,R,S,T)
    !issiso(P) && error("Feedback not implemented for MIMO systems")
    tf(conv(poly2vec(numpoly(P)[1]),T),zpconv(poly2vec(denpoly(P)[1]),R,poly2vec(numpoly(P)[1]),S))
end

feedback2dof(B,A,R,S,T) = tf(conv(B,T),zpconv(A,R,B,S))


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

using Markdown
using InteractiveUtils

# ╔═╡ 4fc6ce8c-bd32-4a84-9b77-b8dc2faa81ce
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
    ])
    using Plots, PlutoUI

	Plots.gr(fmt=:png, lw=2)

	PlutoUI.TableOfContents()
end

# ╔═╡ 2e1253e8-2aca-4d2b-950d-31e0ac0befac
md"""
# Figure 4.11

Surface plots.

Reference: [surface plots @ PlotsGallery.jl](https://goropikari.github.io/PlotsGallery.jl/src/surface.html)
"""

# ╔═╡ 0073304f-63e2-4e30-8358-e056362536e2
x1 = y1 = range(-1.0, 1.0, length=51)

# ╔═╡ 6ad81170-fe26-4bda-a326-494699938661
z1(x, y) = x^2 + 0.5y^2

# ╔═╡ 4f28ccba-f9c7-4068-88e8-37848e4061a1
x2 = range(-2.75, 2.75, length=80)

# ╔═╡ 89346c28-6c84-4204-8dff-6bac59de8bd5
y2 = range(-0.75, 0.75, length=80)

# ╔═╡ 6bbb0bf2-ae1f-4790-8847-dac28a74dd0d
z2(x, y) = (.2x^2-1)^2 + y^2

# ╔═╡ 646f3b77-a4bf-4a3c-b971-85a2841b4543
begin
	p1 = surface(x1, y1, z1, title="Single-well potential")
	p2= contourf(x1, y1, z1)
	p3 = surface(x2, y2, z2, title="Double-well potential")
	p4 = contourf(x2, y2, z2)

	plot(p1, p2, p3, p4,size=(1000, 800))
end

# ╔═╡ 460b4f90-bdec-11eb-2a85-0d7050badf3b
md"""
## Appendix
Running environment and some auxillary functions.

**You can use [this helper tool](https://fonsp.com/article-test-3/pkghelper.html) to generate these commands!**
"""

# ╔═╡ 375b4858-0248-49dd-bc35-4da8e461f0b2
begin
	# Convenience functions
    hill(x, k) = x / (x + k)
    hill(x, k, n) = hill(x^n, k^n)
    exprel(x) = ifelse(x≈zero(x), one(x), x / expm1(x))
end

# ╔═╡ Cell order:
# ╠═2e1253e8-2aca-4d2b-950d-31e0ac0befac
# ╠═0073304f-63e2-4e30-8358-e056362536e2
# ╠═6ad81170-fe26-4bda-a326-494699938661
# ╠═4f28ccba-f9c7-4068-88e8-37848e4061a1
# ╠═89346c28-6c84-4204-8dff-6bac59de8bd5
# ╠═6bbb0bf2-ae1f-4790-8847-dac28a74dd0d
# ╠═646f3b77-a4bf-4a3c-b971-85a2841b4543
# ╠═460b4f90-bdec-11eb-2a85-0d7050badf3b
# ╠═375b4858-0248-49dd-bc35-4da8e461f0b2
# ╠═4fc6ce8c-bd32-4a84-9b77-b8dc2faa81ce

max_stars_repo_name: stevengogogo/pluto-notebooks

max_issues_repo_name: stevengogogo/pluto-notebooks

max_forks_repo_name: stevengogogo/pluto-notebooks


text:
using RiskAdjustedLinearizations, ForwardDiff, Random

Random.seed!(1793)

aR = rand(2)
bR = rand(3)
aD = rand(2)
bD = rand(3)

AR = rand(ForwardDiff.Dual, 2, 2, 3)
BR = rand(ForwardDiff.Dual, 3, 2, 3, 2)
AD = rand(ForwardDiff.Dual, 2, 2, 3)
BD = rand(ForwardDiff.Dual, 3, 2, 3, 2)
@testset "dualarray and dualvector" begin
    @test RiskAdjustedLinearizations.dualarray(AR, BR) == similar(AR)
    @test RiskAdjustedLinearizations.dualarray(BR, AR) == similar(BR)
    @test RiskAdjustedLinearizations.dualarray(AR, BD) == similar(AR, eltype(BD))
    @test RiskAdjustedLinearizations.dualarray(BD, AR) == similar(BD, eltype(AR))
    @test RiskAdjustedLinearizations.dualarray(AD, BD) == similar(AD)
    @test RiskAdjustedLinearizations.dualarray(BD, AD) == similar(BD)
    @test RiskAdjustedLinearizations.dualarray(AD, BR) == similar(AD)
    @test RiskAdjustedLinearizations.dualarray(BR, AD) == similar(BR, eltype(AD))

    @test RiskAdjustedLinearizations.dualvector(aR, bR) == similar(aR)
    @test RiskAdjustedLinearizations.dualvector(bR, aR) == similar(bR)
    @test RiskAdjustedLinearizations.dualvector(aR, bD) == similar(aR, eltype(bD))
    @test RiskAdjustedLinearizations.dualvector(bD, aR) == similar(bD, eltype(aR))
    @test RiskAdjustedLinearizations.dualvector(aD, bD) == similar(aD)
    @test RiskAdjustedLinearizations.dualvector(bD, aD) == similar(bD)
    @test RiskAdjustedLinearizations.dualvector(aD, bR) == similar(aD)
    @test RiskAdjustedLinearizations.dualvector(bR, aD) == similar(bR, eltype(aD))
end

nothing


text:
# Appends curves to an existing curve
function append_to_curve!(target::Curve, newstuff::Curve)
    s_end = target[end].s
    for c in newstuff
        push!(target, CurvePt(c.pos, c.s+s_end, c.k, c.kd))
    end
    return target
end

# Get the lane object by index from the roadway
roadway_lane(roadway::Roadway, laneid::Int) = roadway.segments[laneid].lanes[1]

# Get the frenet position for provided lane id
function AutomotiveSimulator.Frenet(roadway::Roadway, laneid::Int, s::Float64, t::Float64 = 0., ϕ::Float64 = 0.)
    Frenet(roadway_lane(roadway, laneid), s, t, ϕ)
end

## Straight roadway with a crosswalk
struct Crosswalk
    crosswalk::Lane
end

# Render the crosswalk
function AutomotiveVisualization.add_renderable!(rendermodel::RenderModel, env::Crosswalk)
    curve = env.crosswalk.curve
    n = length(curve)
    pts = Array{Float64}(undef, 2, n)
    for (i,pt) in enumerate(curve)
        pts[:,i] .= pt.pos[1:2]
    end

    add_instruction!(
        rendermodel, render_dashed_line,
        (pts, colorant"white", env.crosswalk.width, 1.0, 1.0, 0.0, 0)
    )
    rendermodel
end

# The crosswalk is lane 2 in the function below
const DEFAULT_CROSSWALK_LANE = 2

function road_with_crosswalk(; roadway_length = 50., crosswalk_length = 20., crosswalk_width = 6.0, crosswalk_pos = roadway_length/2, )
    # Generate a straight 2-lane roadway and a crosswalk lane
    roadway = gen_straight_roadway(2, roadway_length)
    crosswalk_start = VecE2(crosswalk_pos, -crosswalk_length/2)
    crosswalk_end = VecE2(crosswalk_pos, crosswalk_length/2)
    crosswalk_lane = gen_straight_curve(crosswalk_start, crosswalk_end, 2)
    crosswalk = Lane(LaneTag(2,1), crosswalk_lane, width = crosswalk_width)
    cw_segment = RoadSegment(2, [crosswalk])
    push!(roadway.segments, cw_segment) # append it to the roadway

    # Describes which lanes each lane should yield to (i.e. lane 4 yields to 1 and 2)
    yields_way = Dict(1=>[2], 2=>[])

    # The entry point of the intersection for each lane
    intersection_enter_loc = Dict(
        1 => VecSE2(crosswalk_pos - crosswalk_width /2., 0., 0),
        2 => VecSE2(crosswalk_pos, -DEFAULT_LANE_WIDTH / 2., 0)
        )

    # the exit point of the intersection for each lane
    intersection_exit_loc = Dict(
        1 => VecSE2(crosswalk_pos + crosswalk_width /2., 0., 0),
        2 => VecSE2(crosswalk_pos, DEFAULT_LANE_WIDTH / 2., 0)
        )

    goals = Dict{Int64, Array{Int}}(
        1 => [1],
        2 => [2],
        )

    should_blink = Dict{Int64, Bool}(
        1 => false,
        2 => false,
        )

    roadway, Crosswalk(crosswalk), yields_way, intersection_enter_loc, intersection_exit_loc, goals, should_blink
end

ped_roadway, crosswalk, ped_yields_way, ped_intersection_enter_loc, ped_intersection_exit_loc, ped_goals, ped_should_blink = road_with_crosswalk()

## T-Intersection Roadway

# LP is the left-most point in between upper and lower lanes
# RP is the Right-most point in between upper and lower lanes
# BP is the Bottom-most point in between left and right lanes
# C is the center of the intersection
# dx is half lane shift in the x direction
# dy is half lane shift in the y direction
# r is the turn radius for the intersection
function T_intersection(; LP = VecE2(-50,0), RP = VecE2(50,0), BP = VecE2(0,-50), C = VecE2(0,0), dx = VecE2(DEFAULT_LANE_WIDTH / 2, 0), dy = VecE2(0, DEFAULT_LANE_WIDTH / 2), r = 4.)

    # Setup the roadway
    roadway = Roadway()

    # Append right turn coming from the left
    curve = gen_straight_curve(LP - dy, C - dy - r*dx, 2)
    append_to_curve!(curve, gen_bezier_curve(VecSE2(C - dy - r*dx, 0), VecSE2(C - r*dy - dx, -π/2), 0.6r, 0.6r, 51)[2:end])
    append_to_curve!(curve, gen_straight_curve(C - r*dy - dx, BP - dx, 2))
    lane1 = Lane(LaneTag(length(roadway.segments)+1,1), curve)
    push!(roadway.segments, RoadSegment(lane1.tag.segment, [lane1]))

    # Append straight right
    curve = gen_straight_curve(LP - dy, RP - dy, 2)
    lane2 = Lane(LaneTag(length(roadway.segments)+1,1), curve)
    push!(roadway.segments, RoadSegment(lane2.tag.segment, [lane2]))

    # Append straight left
    curve = gen_straight_curve(RP + dy, LP + dy, 2)
    lane3 = Lane(LaneTag(length(roadway.segments)+1,1), curve)
    push!(roadway.segments, RoadSegment(lane3.tag.segment, [lane3]))

    # Append left turn coming from the right
    curve = gen_straight_curve(RP + dy, C + r*dx + dy, 2)
    append_to_curve!(curve, gen_bezier_curve(VecSE2(C + r*dx + dy,-π), VecSE2(C - dx - r*dy, -π/2), r, r, 51)[2:end])
    append_to_curve!(curve, gen_straight_curve(C - dx - r*dy, BP - dx, 2))
    lane4 = Lane(LaneTag(length(roadway.segments)+1,1), curve)
    push!(roadway.segments, RoadSegment(lane4.tag.segment, [lane4]))

    # Append left turn coming from below
    curve = gen_straight_curve(BP + dx, C + dx - r*dy, 2)
    append_to_curve!(curve, gen_bezier_curve(VecSE2(C + dx - r*dy, π/2), VecSE2(C - r*dx + dy, -π), r, r, 51)[2:end])
    append_to_curve!(curve, gen_straight_curve(C - r*dx + dy, LP + dy, 2))
    lane5 = Lane(LaneTag(length(roadway.segments)+1,1), curve)
    push!(roadway.segments, RoadSegment(lane5.tag.segment, [lane5]))

    # Append right turn coming from below
    curve = gen_straight_curve(BP + dx, C + dx - r*dy, 2)
    append_to_curve!(curve, gen_bezier_curve(VecSE2(C + dx - r*dy, π/2), VecSE2(C + r*dx - dy, 0), 0.6r, 0.6r, 51)[2:end])
    append_to_curve!(curve, gen_straight_curve(C + r*dx - dy, RP - dy, 2))
    lane6 = Lane(LaneTag(length(roadway.segments)+1,1), curve)
    push!(roadway.segments, RoadSegment(lane6.tag.segment, [lane6]))

    # Setup dictionaries that describe priorities, goals and intersection information
    ################################################################################

    # Describes which lanes each lane should yield to (i.e. lane 4 yields to 1 and 2)
    yields_way = Dict(1=>[], 2=>[], 3=>[], 4=>[1,2], 5=>[2,3,4], 6=>[1,2,])

    # The entry point of the intersection for each lane
    intersection_enter_loc = Dict(
        1=> VecSE2(C - dy - r*dx, 0),
        2=> VecSE2(C - dy - r*dx, 0),
        3=> VecSE2(C + dy + r*dx, π),
        4=> VecSE2(C + dy + r*dx, π),
        5=> VecSE2(C + dx - r*dy, π/2),
        6=> VecSE2(C + dx - r*dy, π/2),
        )

    # the exit point of the intersection for each lane
    intersection_exit_loc = Dict(
        1=> VecSE2(C - r*dy - dx, -π/2),
        2=> VecSE2(C - dy + r*dx, 0),
        3=> VecSE2(C + dy - r*dx, π),
        4=> VecSE2(C - r*dy - dx, -π/2),
        5=> VecSE2(C + dy - r*dx, π),
        6=> VecSE2(C - dy + r*dx, 0),
        )

    goals = Dict{Int64, Array{Int}}(
        1 => [1,2],
        2 => [1,2],
        3 => [3,4],
        4 => [3,4],
        5 => [5,6],
        6 => [5,6]
    )

    # Whether or not the blinker should be on for each lane
    should_blink = Dict{Int64, Bool}(
        1 => true,
        2 => false,
        3 => false,
        4 => true,
        5 => true,
        6 => true,
        )

    return roadway, yields_way, intersection_enter_loc, intersection_exit_loc, goals, should_blink
end

# Strings describing the goal of each lane
Tint_goal_label = Dict(1 => "turn right", 2 =>"straight", 3=>"straight", 4=>"turn left", 5=>"turn left", 6=> "turn right")

# whether the signal is on the right side
Tint_signal_right = Dict(
    1 => true,
    2 => true,
    3 => false,
    4 => false,
    5 => false,
    6 => true
    )

# Construct a global T-intersection to be used
Tint_roadway, Tint_yields_way, Tint_intersection_enter_loc, Tint_intersection_exit_loc, Tint_goals, Tint_should_blink = T_intersection()

# Construct template TIDM models (one for each roadway)

Tint_TIDM_template = TIDM(yields_way = Tint_yields_way,
                    intersection_enter_loc = Tint_intersection_enter_loc,
                    intersection_exit_loc = Tint_intersection_exit_loc,
                    goals = Tint_goals,
                    should_blink = Tint_should_blink
                    )

ped_TIDM_template = TIDM(yields_way = ped_yields_way,
                        intersection_enter_loc = ped_intersection_enter_loc,
                        intersection_exit_loc = ped_intersection_exit_loc,
                        goals = ped_goals,
                        should_blink = ped_should_blink
                        )


text:
function spline2d(xr, yr, mx, my)

    out = zeros(Float64, mx*my)

    for i=1:mx
        for j=1:my
            k = i + (j-1)*mx
            a = xr^(i-1) / factorial(i-1)
            b = yr^(j-1) / factorial(j-1)
            out[k] = a * b
        end
    end

    return out
end

function stencil2d(mx, my)

    # (i,j) indices for center point
    ic = Int(ceil(mx/2.0))
    jc = Int(ceil(my/2.0))

    M = zeros(Float64, (mx*my, mx*my))

    for i=1:mx
        for j=1:my
            k = i + (j-1)*mx
            M[k,:] = spline2d(i-ic, j-jc, mx, my)
        end
    end

    return Matrix(transpose(inv(M)))
end

function interpolation_row(x, y, mx, my, MinvT, grd::Grid)

    x0, y0 = grd.r0
    dx = grd.dx
    dy = grd.dy
    Nx = grd._Nx
    Ny = grd._Ny

    # (i,j) indices for center point of stencil
    ics = Int(ceil(mx/2.0))
    jcs = Int(ceil(my/2.0))

    # indices for center point in grid
    # NOTE: the addition of 1e-5 before calling floor is to compensate for annoying floating point imprecision.
    icg = Int(floor((x - x0)/dx + 1e-5)) + 1 + grd._Nbuffer
    jcg = Int(floor((y - y0)/dy + 1e-5)) + 1 + grd._Nbuffer

    # values for row
    xr = (x - x0)/dx - (icg - 1)
    yr = (y - y0)/dy - (jcg - 1)
    spline = spline2d(xr, yr, mx, my)
    row_dat = MinvT * spline

    # collective index for first nonzero element of output
    k0 = (icg - ics + 1) + (jcg - jcs)*Nx

    # indices for nonzero elements of row
    row_js = zeros(Int, mx*my)
    for i=1:mx
        for j=1:my
            row_js[i+(j-1)*mx] = k0 + (i-1) + (j-1)*Nx
        end
    end

    return row_dat, row_js
end

function interpolation_function(V::Vector{Float64}, mx, my, MinvT, grd::Grid)

    f(x, y) = begin
        row_dat, row_js = interpolation_row(x, y, mx, my, MinvT, grd)
        row = sparse(ones(Int64, mx*my), row_js, row_dat, 1, grd._Nx*grd._Ny) * transpose(grd.Proj)

        return (row * V)[1]
    end

    return f
end


text:
function ssaargmap(f, stmt)
    urs = userefs(stmt)
    urs === () && return stmt
    for op in urs
        val = op[]
        if isa(val, Union{SSAValue, Argument})
            op[] = f(val)
        end
    end
    urs[]
end

function replace_code!(ci::CodeInfo, code::IRCode, nargs)
    cfg = compute_basic_blocks(code.stmts)
    if !isempty(code.new_nodes)
        code = compact!(code)
    end
    # All but the first `nargs` slots will now be unused
    resize!(ci.slottypes, nargs+1)
    resize!(ci.slotnames, nargs+1)
    resize!(ci.slotflags, nargs+1)
    # For every used SSAValues, we register one base format ssa value
    used = Set{Int}()
    foreach(stmt->scan_ssa_use!(used, stmt), code.stmts)
    mapping = Dict{Int, Int}()
    n = 0
    resize!(ci.ssavaluetypes, length(used))
    for ssa in sort(collect(used))
        mapping[ssa] = n
        n += 1
        ci.ssavaluetypes[n] = code.types[ssa]
    end
    # Find all jump targets (we need to insert LabelNodes for them) and
    # jump origins (we insert a label node on the statement after, to
    # make sure we can track them)
    dest_blocks = Set{Int}()
    jump_origins = Set{Int}()
    for stmt in code.stmts
        if isa(stmt, GotoNode)
            push!(dest_blocks, stmt.label)
        elseif isa(stmt, GotoIfNot)
            push!(dest_blocks, stmt.dest)
        elseif isa(stmt, PhiNode)
            for edge in stmt.edges
                push!(jump_origins, edge)
            end
        end
    end
    block_start = Dict(first(code.cfg.blocks[x].stmts)=>x for x in dest_blocks)
    comefrom_labels = Set{Int}(last(code.cfg.blocks[x].stmts)+1 for x in jump_origins)
    block_terminators = Dict(last(block.stmts)=>i for (i,block) in pairs(code.cfg.blocks))
    function rename(val)
        if isa(val, SSAValue)
            if haskey(mapping, val.id)
                return SSAValue(mapping[val.id])
            end
        elseif isa(val, Argument)
            return SlotNumber(val.n)
        end
        return val
    end
    # Now translate the code
    new_code = Vector{Any}()
    label_mapping = Dict{Int, Int}()
    terminator_mapping = Dict{Int, Int}()
    fixup = Int[]
    for (idx, stmt) in pairs(code.stmts)
        if haskey(block_start, idx)
            push!(new_code, LabelNode(length(new_code) + 1))
            label_mapping[block_start[idx]] = length(new_code)
        elseif idx in comefrom_labels
            push!(new_code, LabelNode(length(new_code) + 1))
        end
        if isa(stmt, GotoIfNot)
            new_stmt = Expr(:gotoifnot, rename(stmt.cond), stmt.dest)
            push!(fixup, length(new_code)+1)
        elseif isa(stmt, ReturnNode)
            if isdefined(stmt, :val)
                new_stmt = Expr(:return, rename(stmt.val))
            else
                # Unreachable, so no issue with this
                new_stmt = nothing
            end
        elseif isa(stmt, SSAValue)
            new_stmt = rename(stmt)
        elseif isa(stmt, PhiNode)
            new_stmt = ssaargmap(rename, stmt)
            push!(fixup, length(new_code)+1)
        elseif isa(stmt, GotoNode)
            push!(fixup, length(new_code)+1)
            new_stmt = stmt
        else
            new_stmt = ssaargmap(rename, stmt)
        end
        if haskey(mapping, idx)
            new_stmt = Expr(:(=), SSAValue(mapping[idx]), new_stmt)
        end
        if haskey(block_terminators, idx)
            terminator_mapping[block_terminators[idx]] = length(new_code)+1
        end
        push!(new_code, new_stmt)
    end
    for i in fixup
        val = new_code[i]
        if isexpr(val, :(=))
            val = val.args[2]
        end
        val = if isa(val, PhiNode)
            # Translate from BB edges to statement edges
            edges = map(val.edges) do edge
                terminator_mapping[edge]
            end
            PhiNode(convert(Vector{Any}, edges), val.values)
        elseif isa(val, GotoNode)
            GotoNode(label_mapping[val.label])
        elseif isexpr(val, :gotoifnot)
            Expr(:gotoifnot, val.args[1], label_mapping[val.args[2]])
        else
            @show val
            error()
        end
        if isexpr(new_code[i], :(=))
            new_code[i].args[2] = val
        else
            new_code[i] = val
        end
    end
    ci.code = new_code
    ci
end


text:
module Pylustrator

using PyCall, PyPlot, MacroTools

const pylustrator = PyNULL()

function start()
    copy!(pylustrator, pyimport("pylustrator"))
    pylustrator.start()
end

macro pylustrator(expr)
    # use rmlines so stacktrace returns expected line number
    return rmlines(quote
        # get filename and line position
        stack_position = stacktrace()[1]
        filename = abspath(string(stack_position.file))
        lineno = stack_position.line
        # @assert filename == @__FILE__
        ct = pylustrator.change_tracker
        ct.custom_stack_position = ct.CustomStackPosition(filename, lineno)
        # wrap generated code in a py"""...""" string
        ct.custom_prepend = "py\"\"\"from matplotlib import pyplot as plt;"
        ct.custom_append = "\"\"\""
        $expr
    end)
end

export @pylustrator, @py_str

end # module


text:
# ---
# title: 700. Search in a Binary Search Tree
# id: problem700
# author: zhwang
# date: 2022-02-22
# difficulty: Easy
# categories: Tree
# link: <https://leetcode.com/problems/search-in-a-binary-search-tree/description/>
# hidden: true
# ---
#
# Given the root node of a binary search tree (BST) and a value. You need to
# find the node in the BST that the node's value equals the given value. Return
# the subtree rooted with that node. If such node doesn't exist, you should
# return NULL.
#
# For example,
#
#
#
#     Given the tree:
#             4
#            / \
#           2   7
#          / \
#         1   3
#
#     And the value to search: 2
#
#
# You should return this subtree:
#
#
#
#           2
#          / \
#         1   3
#
#
# In the example above, if we want to search the value `5`, since there is no
# node with value `5`, we should return `NULL`.
#
# Note that an empty tree is represented by `NULL`, therefore you would see the
# expected output (serialized tree format) as `[]`, not `null`.
#
#
## @lc code=start
using LeetCode

function search_bst(root::TreeNode, val::Int)::Union{TreeNode, Nothing}
    while !isnothing(root)
        root.val == val && return root
        root = (root.val > val) ? root.left : root.right
    end
end
## @lc code=end
## @lc test=start
@testset "700.search-in-a-binary-search-tree.jl" begin
    @test search_bst(TreeNode{Int}([4,2,7,1,3]), 2) == TreeNode{Int}([2,1,3])
    @test search_bst(TreeNode{Int}([4,2,7,1,3]), 5) === nothing
end
## @lc test=end

## @lc info=start
# link: [solution to 700](https://leetcode-cn.com/problems/search-in-a-binary-search-tree/solution/700-er-cha-sou-suo-shu-zhong-de-sou-suo-5su23/)
## @lc info=end


text:
# open Demo/UIControls/UIControls.xcworkspace
# Run

using Swifter

vc = initial("http://localhost:8080")

@query vc.label.text = "Hello Swift"
sleep(0.5)

@query vc.textField.text = "Hello Julia"
sleep(0.5)

@query vc.button.tap()
sleep(0.5)

@query vc.switch_.tap()
sleep(0.5)

@query vc.segmentedControl.tap(title: "Second")
sleep(0.5)

@query vc.segmentedControl.tap(index: 2)
sleep(0.5)

@query vc.slider.tap(value: 1)
sleep(0.5)

@query vc.progressView.setProgress(1, animated: true)

max_stars_repo_name: wookay/Console

max_issues_repo_name: wookay/Console

max_forks_repo_name: wookay/Console


text:
# This is the HMM from the particle filter chapter (ch5).

using Stochy, DataStructures

@pp function hmm(states, observations)
    if observations == list()
        reverse(states)
    else
        local prevstate = first(states),
              state = flip(prevstate ? 0.9 : 0.1)
        factor(state == first(observations) ? 0 : -1)
        hmm(state..states, tail(observations))
    end
end

dist = @pp enum() do
    hmm(list(false), list(true,true,true,true))
end

println("Exact:\n", dist)

dist = @pp smc(100) do
    hmm(list(false), list(true,true,true,true))
end

println("SMC:\n", dist)


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -23.072345451695163 <= q <= 18.38163265505858)

                     add_NL_constraint(m, :((-0.6895881478464707 + 0.42824131807966515*(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4]))/(1 + exp(-(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4])))) + -0.8758217593243995*(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4]))/(1 + exp(-(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4])))) + -0.3819126296826152*(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4]))/(1 + exp(-(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4])))) + -0.7523924735752416*(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4]))/(1 + exp(-(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4])))))/(1 + exp(-(-0.6895881478464707 + 0.42824131807966515*(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4]))/(1 + exp(-(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4])))) + -0.8758217593243995*(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4]))/(1 + exp(-(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4])))) + -0.3819126296826152*(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4]))/(1 + exp(-(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4])))) + -0.7523924735752416*(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4]))/(1 + exp(-(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4]))))))) + (0.5903696565764314 + 0.6045010784464115*(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4]))/(1 + exp(-(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4])))) + 0.19992590084675577*(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4]))/(1 + exp(-(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4])))) + -0.13493555263327828*(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4]))/(1 + exp(-(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4])))) + 0.0007730145672497635*(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4]))/(1 + exp(-(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4])))))/(1 + exp(-(0.5903696565764314 + 0.6045010784464115*(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4]))/(1 + exp(-(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4])))) + 0.19992590084675577*(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4]))/(1 + exp(-(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4])))) + -0.13493555263327828*(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4]))/(1 + exp(-(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4])))) + 0.0007730145672497635*(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4]))/(1 + exp(-(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4]))))))) + (-0.987583954230804 + 0.6297203029880079*(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4]))/(1 + exp(-(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4])))) + -0.939181196020844*(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4]))/(1 + exp(-(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4])))) + 0.23747324715066886*(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4]))/(1 + exp(-(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4])))) + 0.2556344288953589*(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4]))/(1 + exp(-(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4])))))/(1 + exp(-(-0.987583954230804 + 0.6297203029880079*(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4]))/(1 + exp(-(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4])))) + -0.939181196020844*(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4]))/(1 + exp(-(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4])))) + 0.23747324715066886*(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4]))/(1 + exp(-(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4])))) + 0.2556344288953589*(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4]))/(1 + exp(-(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4]))))))) + (-0.44105426128556546 + -0.39452594471110514*(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4]))/(1 + exp(-(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4])))) + 0.16620534874676274*(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4]))/(1 + exp(-(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4])))) + 0.6245292013722641*(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4]))/(1 + exp(-(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4])))) + 0.9983459126930239*(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4]))/(1 + exp(-(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4])))))/(1 + exp(-(-0.44105426128556546 + -0.39452594471110514*(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4]))/(1 + exp(-(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4])))) + 0.16620534874676274*(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4]))/(1 + exp(-(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4])))) + 0.6245292013722641*(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4]))/(1 + exp(-(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4])))) + 0.9983459126930239*(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4]))/(1 + exp(-(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using GivensBoson
using Test
using Random
using LinearAlgebra
include("./testutil.jl")
include("./function_test.jl")

@testset "One Zero Mode Cases" begin
    A = [1 1;1 1]
    ϵ = 1E-12
    S,V = given_eigen_solver(copy(A),hamiltonian_type="Symmetry",zeromode=true)
    # do not deal zero modes. Canonical and Right
    test_canonicality(V,ϵ)
    test_right(S,V,A,ϵ)
    test_full_right(S,V,A,ϵ)

    deal_zeromodes!(S,V,type=:abnormal)

    # deal zero modes. Not canonical but diagonal.
    test_diag(S,ϵ)
    test_right(S,V,A,ϵ)
    #TODO: zeromode canonicality

    #TODO: if type=normal test
end

@testset "GivensBoson.jl - random case with partile-hole symmetry" begin
    for N = 4:2:40
        Random.seed!(42)
        ϵ = 1E-6
        A = test_random_particle_hole_positive_definite_hamiltonian(N)
        S,V = given_eigen_solver(copy(A),hamiltonian_type="Symmetry",zeromode=true)
        deal_zeromodes!(S,V,type=:abnormal)
        test_canonicality(V,ϵ)
        test_diag(S,ϵ)
        test_right(S,V,A,ϵ)
    end
end

@testset "GivensBoson -Schwinger Boson PBC- AntiSymmetry" begin
    for N = 4:2:60
        ϵ = 1E-4
        A = schwinger_boson_hamiltonian(N,"PBC")
        S,V = given_eigen_solver(copy(A),hamiltonian_type="AntiSymmetry",zeromode=true)
        test_canonicality_zeromode(V,ϵ)
        test_right(S,V,A,ϵ)
        # test_diag(S,ϵ)
    end
end

@testset "GivensBoson -Schwinger Boson OBC- AntiSymmetry" begin
    for N = 4:2:64
        ϵ = 1E-9
        A = schwinger_boson_hamiltonian(N,"OBC")
        origin_A = copy(A)
        S,V = given_eigen_solver(A,hamiltonian_type="AntiSymmetry")
        test_canonicality_zeromode(V,ϵ)
        test_right(S,V,A,ϵ)
        # test_diag(S,ϵ)
    end
end

# @testset "ReCanonicalize - GivensBoson Equal OBC test" begin
#     include("../src/ReCanonicalize.jl")
#     for N = 2:2:40
#         A = 1
#         λ = 1
#         bc = "OBC"
#         begin
#             Q = zeros(N,N)
#             for id in CartesianIndices((1:N,1:N))
#             if id[1]-id[2]==1
#                 Q[id] = 1.0*A
#             elseif id[2] - id[1] ==1
#                 Q[id] = -1.0*A
#             end
#             end

#             if bc == "PBC"
#             Q[N,1] = -1.0*A
#             Q[1,N] =  1.0*A
#             end
#             Q = Q/2.0
#             Q
#         end
#         Lambda = diagm([λ for i = 1:N])
#         M = [Lambda Q;transpose(Q) Lambda]
#         A = M
#         origin_A = copy(A)
#         # S,V = given_eigen_solver(A)
#         ϵ = 1E-8

#         S1,V1 = given_eigen_solver(origin_A)
#         S2,V2 = recanonicalize(A,"AntiSymmetry")
#         @test norm(S1-diagm(S2))<ϵ
#         @test norm(V1-V2)<ϵ
#     end
# end


text:
using Pkg; Pkg.activate(".")
using HypergraphModularity
using StatsBase
using DataFrames
using SparseArrays
using Clustering
using CSV

function recoveryExperiment(; data, n_rounds = 2, core = 0, Γ₀ = 100.0, γ₀ = 1.0, kmax = 20, save_as = "recovery_throughput/experiment.csv")

    # collector
    DF = DataFrame()

    H, Z = read_hypergraph_data(data, kmax)

    H, Z = kcore(H, Z, core)
    n = length(H.D)

    kmin = minimum([k for k in keys(H.E) if length(H.E[k]) > 0])

    println("$(core)-core of $(data) has $n nodes")

    for k ∈ reverse(kmin:kmax)
        # remove edges of size larger than k
        removeEdges!(H; remove = collect((k+1):(kmax+1)))

        m = sum(length(H.E[k]) for k in keys(H.E))

        # polyadic warmstart
        D = big(sum(H.D))

        ωᵢ = 2.0
        ωₒ = 0.6
        function ω(p, α)
            num = p[1] == 1 ? ωᵢ : ωₒ
            denom = D^sum(p[2])
            return Γ₀*num / denom
        end

        Ω̂ = allOrNothingIntensityFunction(ω, maximum(keys(H.E)))
        Ẑ_w = SuperNode_PPLouvain(H, Ω̂; α = 0, verbose = false, Z0 = collect(1:n));

        # polyadic experiments
        for i ∈ 1:n_rounds
            Ẑ, t = @timed SuperNode_PPLouvain(H, Ω̂; α = 0, verbose = false, scan_order = "random")
            Ω̄ = estimateΩEmpirically(H, Ẑ; aggregator = p -> [length(p) == 1, sum(p)])

            Q = modularity(H, Ẑ, Ω̂; α = 0)
            ℓ = length(unique(Ẑ))
            ARI = randindex(Z, Ẑ)[1]
            df = DataFrame(data = data, kmax = k, ℓ = ℓ, Q = Q, t = t, ARI = ARI, method = "Polyadic", n = n, core = core, round = i, thread = Threads.threadid(), m = m)
            DF = vcat(DF, df)
        end

        # construct dyadic graph
        H̄ = projectedGraph(H)

        # dyadic warmstart
        D̄ = big(sum(H̄.D))

#         function ω_d(p, α)
#             num = p[1] == 1 ? ωᵢ : ωₒ
#             denom = D̄^sum(p[2])
#             return γ₀*num / denom
#         end

#         Ω̄ = allOrNothingIntensityFunction(ω_d, 2)
#         Z̄_w = SuperNode_PPLouvain(H, Ω̄; α = 0, verbose = false, Z0 = collect(1:n));

        best_Z̄ = zero(Z)
        best_Q = -Inf

        Ω̄ = estimateΩEmpirically(H̄, Ẑ_w; aggregator = p -> [length(p) == 1, sum(p)])
        for i ∈ 1:n_rounds
            Z̄, t = @timed SuperNode_PPLouvain(H̄, Ω̄; α = 0, verbose = false)
            Ω̄ = estimateΩEmpirically(H̄, Z̄; aggregator = p -> [length(p) == 1, sum(p)])
            Q = modularity(H̄, Z̄, Ω̄; α = 0)
            ℓ = length(unique(Z̄))
            ARI = randindex(Z, Z̄)[1]
            df = DataFrame(data = data, kmax = k, ℓ = ℓ, Q = Q, method = "Dyadic", t = t, ARI = ARI, n = n, core = core, round = i, thread = Threads.threadid(), m = m)
            DF = vcat(DF, df)

            if Q > best_Q
                best_Q = Q
                best_Z̄ = Z̄
            end
        end

        # dyadic with weighted projection

        γ = computeDyadicResolutionParameter(H, Ẑ_w; mode = "γ", weighted = true)

        for i ∈ 1:n_rounds
            Z̄, t = @timed CliqueExpansionModularity(H, γ; weighted = true)
            γ = computeDyadicResolutionParameter(H, Z̄; mode = "γ", weighted = true)
            Ω̄ = estimateΩEmpirically(H̄, Z̄; aggregator = p -> [length(p) == 1, sum(p)])
            Q = modularity(H̄, Z̄, Ω̄; α = 0)
            ℓ = length(unique(Z̄))
            ARI = randindex(Z, Z̄)[1]
            df = DataFrame(data = data, kmax = k, ℓ = ℓ, Q = Q, method = "Dyadic (weighted)", t = t, ARI = ARI, n = n, core = core, round = i, thread = Threads.threadid(), m = m)
            DF = vcat(DF, df)
        end

        # refinement
        Ẑ = best_Z̄
        for i ∈ 1:n_rounds
            Ω̂ = estimateΩEmpirically(H, Ẑ; aggregator = p -> [length(p) == 1, sum(p)])
            Ẑ, t = @timed SuperNode_PPLouvain(H, Ω̂; Z0 = Ẑ, α = 0, scan_order = "random", verbose = false)
            Q = modularity(H, Ẑ, Ω̂; α = 0)
            ℓ = length(unique(Ẑ))
            ARI = randindex(Z, Ẑ)[1]
            df = DataFrame(data = data, kmax = k, ℓ = ℓ, Q = Q, t = t, ARI = ARI, method = "Polyadic Refinement", n = n, core = core, round = i, thread = Threads.threadid(), m = m)
            DF = vcat(DF, df)
        end
        path = "fig/recovery_throughput/$(k)_$(save_as).csv"
        CSV.write(path, DF)
    end

    return DF
end

## Main experiment
n_rounds = 20
kmax = 25

control = [
    Dict(:data => "TrivagoClickout-fix",   :core => 2,  :n_rounds => n_rounds, :Γ₀ => 1000.0, :γ₀ => 100000.0, :kmax => kmax, :save_as => "TrivagoClickout-fix-2-core"),
    Dict(:data => "TrivagoClickout",   :core => 2,  :n_rounds => n_rounds, :Γ₀ => 1000.0, :γ₀ => 100000.0, :kmax => kmax, :save_as => "TrivagoClickout-2-core")
];

Threads.@threads for d in control
    recoveryExperiment(;d...)
end

max_stars_repo_name: PhilChodrow/HypergraphModularity

max_issues_repo_name: PhilChodrow/HypergraphModularity

max_forks_repo_name: PhilChodrow/HypergraphModularity


text:
include("../src/TrueSkill.jl")
using .TrueSkill
global const ttt = TrueSkill
using CSV; using Dates

println("Code 1")
mu = 0.0; sigma = 6.0; beta = 1.0; gamma = 0.03; draw = 0.0

println("Code 2")
a1 = ttt.Player(ttt.Gaussian(mu, sigma), beta, gamma); a2 = ttt.Player()
a3 = ttt.Player(); a4 = ttt.Player()

println("Code 3")
team_a = [ a1, a2 ]
team_b = [ a3, a4 ]
teams = [team_a, team_b]
g = ttt.Game(teams)

@time ttt.posteriors(ttt.Game(teams))

println("Code 4")
lhs = g.likelihoods
ev = g.evidence
ev = round(ev, digits=3)
println(ev)

println("Code 5")
pos = ttt.posteriors(g)
print(pos[1][1])
print(lhs[1][1] * a1.prior)

println("Code 6")
ta = [a1]
tb = [a2, a3]
tc = [a4]
teams_3 = [ta, tb, tc]
result = [1., 0., 0.]
g = ttt.Game(teams_3, result, p_draw=0.25)

@time ttt.posteriors(ttt.Game(teams_3, result, p_draw=0.25))

println("Code 7")
c1 = [["a"],["b"]]
c2 = [["b"],["c"]]
c3 = [["c"],["a"]]
composition = [c1, c2, c3]
h = ttt.History(composition, gamma = 0.0)

@time ttt.History(composition)

println("Code 8")
lc = ttt.learning_curves(h)
print(lc["a"])
print(lc["b"])

println("Code 9")
ttt.convergence(h)
lc = ttt.learning_curves(h)
print(lc["a"])
print(lc["b"])

@time ttt.convergence(h, iterations=1)

println("Code 10")
using Random; Random.seed!(999); N = 1000
function skill(experience, middle, maximum, slope)
    return maximum/(1+exp(slope*(-experience+middle)))
end
target = skill.(1:N, 500, 2, 0.0075)
opponents = Random.randn.(1000)*0.5 .+ target

println("Code 11")
composition = [[["a"], [string(i)]] for i in 1:N]
results = [ r ? [1.,0.] : [0.,1.] for r in (Random.randn(N).+target .> Random.randn(N).+opponents) ]
times = [i for i in 1:N]
priors = Dict{String,ttt.Player}()
for i in 1:N  priors[string(i)] = ttt.Player(ttt.Gaussian(opponents[i], 0.2))  end
@time h = ttt.History(composition, results, times, priors, gamma=0.015)
ttt.convergence(h)
@time ttt.convergence(h, iterations=1)

mu = [tp[2].mu for tp in ttt.learning_curves(h)["a"]]

println("Code 12")
data = CSV.read("input/history.csv")

times = Dates.value.(data[:,"time_start"] .- Date("1900-1-1"))
composition = [ r.double == "t" ? [[r.w1_id,r.w2_id],[r.l1_id,r.l2_id]] :
[[r.w1_id],[r.l1_id]] for r in eachrow(data) ]

h = ttt.History(composition=composition, times = times, sigma = 1.6, gamma = 0.036)
ttt.convergence(h,epsilon=0.01, iterations=10)

println("Code 13")
players = Set(vcat((composition...)...))
priors = Dict([(p, ttt.Player(ttt.Gaussian(0., 1.6), 1.0, 0.036) ) for p in players])

composition_ground = [ r.double == "t" ? [[r.w1_id, r.w1_id*r.ground, r.w2_id, r.w2_id*r.ground],[r.l1_id, r.l1_id*r.ground, r.l2_id, r.l2_id*r.ground]] : [[r.w1_id, r.w1_id*r.ground],[r.l1_id, r.l1_id*r.ground]] for r in eachrow(data) ]

h_ground = ttt.History(composition=composition_ground, times = times, sigma = 1.0, gamma = 0.01, beta = 0.0, priors = priors)
ttt.convergence(h_ground,epsilon=0.01, iterations=10)

println("Code 14")
N1 = ttt.Gaussian(mu = 1.0, sigma = 1.0); N2 = ttt.Gaussian(1.0, 2.0)

println("Code 15")
p1 = ttt.performance(a1)
p2 = ttt.performance(a2)
p3 = ttt.performance(a3)
p4 = ttt.performance(a4)

println("Code 16")
ta = p1 + p2; tb = p3 + p4

println("Code 17")
d = ta - tb

println("Code 18")
e = 1.0 - ttt.cdf(d, 0.0)

println("Code 19")
na = length(team_a)
nb = length(team_b)
sd = sqrt(na + nb)*beta
p_draw = 0.25
margin = ttt.compute_margin(p_draw, sd)

println("Code 20")
g = ttt.Game(teams, p_draw = 0.25)
post = ttt.posteriors(g)

println("Code 21")
d_approx = ttt.approx(d, margin, false)

println("Code 22")
approx_lh_d = d_approx / d

println("Code 23")
mu = a1.prior.mu
sigma2 = a1.prior.sigma^2
phi = d.mu
v2 = d.sigma^2
phi_div = approx_lh_d.mu
v2_div = approx_lh_d.sigma^2
prior = a1.prior
posterior = post[1][1]
println( prior * ttt.Gaussian(mu-phi+phi_div, sqrt(v2 + v2_div - sigma2)) )

println(posterior)


text:
module TestECOS

using Test
using MathOptInterface
import ECOS

const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

function test_runtests()
    solver = MOI.OptimizerWithAttributes(ECOS.Optimizer, MOI.Silent() => true)
    model = MOI.instantiate(solver, with_bridge_type = Float64)
    MOI.Test.runtests(
        model,
        MOI.Test.Config(
            atol = 1e-3,
            rtol = 1e-3,
            exclude = Any[
                MOI.ConstraintBasisStatus,
                MOI.VariableBasisStatus,
                MOI.ConstraintName,
                MOI.VariableName,
                MOI.ObjectiveBound,
            ],
        ),
        exclude = String[
            # Expected test failures:
            #   Problem is a nonconvex QP
            "test_basic_ScalarQuadraticFunction_EqualTo",
            "test_basic_ScalarQuadraticFunction_GreaterThan",
            "test_basic_ScalarQuadraticFunction_Interval",
            "test_basic_VectorQuadraticFunction_",
            "test_quadratic_SecondOrderCone_basic",
            "test_quadratic_nonconvex_",
            #   MathOptInterface.jl issue #1431
            "test_model_LowerBoundAlreadySet",
            "test_model_UpperBoundAlreadySet",
        ],
    )
    return
end

function test_RawOptimizerAttribute()
    model = ECOS.Optimizer()
    MOI.set(model, MOI.RawOptimizerAttribute("abstol"), 1e-5)
    @test MOI.get(model, MOI.RawOptimizerAttribute("abstol")) ≈ 1e-5
    MOI.set(model, MOI.RawOptimizerAttribute("abstol"), 2e-5)
    @test MOI.get(model, MOI.RawOptimizerAttribute("abstol")) == 2e-5
end

function test_iteration_limit()
    # Problem data
    v = [5.0, 3.0, 1.0]
    w = [2.0, 1.5, 0.3]

    solver = MOI.OptimizerWithAttributes(ECOS.Optimizer, MOI.Silent() => true)
    model = MOI.instantiate(solver, with_bridge_type = Float64)

    maxit = 1
    MOI.set(model, MOI.RawOptimizerAttribute("maxit"), maxit)
    MOI.set(model, MOI.Silent(), true)

    x = MOI.add_variables(model, 3)
    for xj in x
        MOI.add_constraint(
            model,
            MOI.SingleVariable(xj),
            MOI.Interval(0.0, 1.0),
        )
    end

    # Constraint
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(w, x), 0.0),
        MOI.LessThan(3.0),
    )

    # Set the objective
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(v, x), 0.0),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(model)

    @test MOI.get(model, MOI.TerminationStatus()) == MOI.ITERATION_LIMIT
    @test MOI.get(model, MOI.BarrierIterations()) == maxit
end

end  # module

TestECOS.runtests()


text:
##  Copyright 2022 Myrl Marmarelis
##
##  Licensed under the Apache License, Version 2.0 (the "License");
##  you may not use this file except in compliance with the License.
##  You may obtain a copy of the License at
##
##    http://www.apache.org/licenses/LICENSE-2.0
##
##  Unless required by applicable law or agreed to in writing, software
##  distributed under the License is distributed on an "AS IS" BASIS,
##  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
##  See the License for the specific language governing permissions and
##  limitations under the License.

using Base.Threads
using ProgressMeter

# integrate with IPW to attain unbiased overall dose-response curves?
function bound_potential_outcome(operate::Function, # monotonic
    covariates::AbstractVector{SVector{CIn,T}},
    predictors::Vector{P}, predictor_param::PredictionParam{T},
    propensity_models::Vector{PP}, propensity_param::PredictionParam{T},
    weight_scheme::WeightScheme{T}, trust_scheme::TrustScheme{T},
    treatment_index::Int; n_resamples::Int, n_draws::Int, order::Int, verbose::Bool=true,
    sensitivity::T, treatment::T, generator::Distribution, static_draws::AbstractVector{T}=T[],
    weight_threshold::T=T(1e-8) )::NTuple{2,Array{T,3}} where {
      In, CIn, POut, PPIn, PPOut, T, P <: Predictor{In,POut,T}, PP <: Predictor{PPIn,PPOut,T} } # In-1 is invalid syntactically
  @assert In == CIn+1
  n_predictors = length(predictors)
  n_propensities = length(propensity_models)
  n_instances = length(covariates)
  treatment_index = wrap_treatment_index(treatment_index, In)
  instances = [ Instance{In, 0, T}(i,
    vcat(
        covariate[1:treatment_index-1],
        [treatment],
        covariate[treatment_index:end]
      ) |> SVector{In,T},
    SVector{0,T}([]) ) for (i, covariate) in enumerate(covariates) ]
  n_outputs = n_instance_outputs(predictors[1], predictor_param)
  progress = Progress(n_propensities * n_instances, desc="Evaluating.. ", enabled=verbose)
  denominator_bounds = map(Iterators.product(propensity_models, instances)) do (model, instance)
    next!(progress)
    propensity_instance = make_propensity_instance(instance, treatment_index)
    parameters = map_prediction(
      identity, model, propensity_param, propensity_instance.input, 1)
    bounds = bound_ignorance_weights(trust_scheme, propensity_param,
      parameters, order, sensitivity, treatment)
  end
  denominators_up = [b[2] for b in denominator_bounds]
  denominators_lo = [b[1] for b in denominator_bounds]
  # replicate `bootstrap_prediction_statistic` here?
  bootstrapped_upper = fill(T(NaN), n_outputs, n_resamples, n_instances)
  bootstrapped_lower = fill(T(NaN), n_outputs, n_resamples, n_instances)
  progress = Progress(n_resamples*n_instances, desc="Bootstrapping.. ", enabled=verbose) # better hope this is thread-safe
  @threads for i in 1:n_resamples
    # enable e.g. Bernoulli draws to simply cover the entire support {0,1} once
    draws = length(static_draws) > 0 ? static_draws : convert.(T, rand(generator, n_draws))
    functions = operate.(draws)
    generator_densities = convert.(T, pdf.(generator, draws)) # a la importance sampling
    # this and `denominator_bounds` are (models x instances) matrices
    numerators = map(Iterators.product(predictors, instances)) do (model, instance)
      next!(progress)
      densities = map_prediction(model, predictor_param, instance.input) do params
        # apply the same sample to every output dimension
        [ density(predictor_param, params, draw) / scale
          for (draw, scale) in zip(draws, generator_densities) ]
      end
      reduce(hcat, densities) # draws x outputs
    end # matrix of matrices
    # the block above, up to the loop beginning, were heretofore not in the resampler.
    # but this hopefully averages out some of the monte carlo biases
    predictor_weights = draw_weights(weight_scheme, n_predictors)
    propensity_weights = draw_weights(weight_scheme, n_propensities)
    # I realized that it's permitted to do vector math on a matrix of matrices! (as I do in the following line)
    numerator = (permutedims(numerators) * predictor_weights) ./ sum(predictor_weights)
    denominator_up = (denominators_up' * propensity_weights) ./ sum(propensity_weights)
    denominator_lo = (denominators_lo' * propensity_weights) ./ sum(propensity_weights)
    for k in 1:n_outputs, j in 1:n_instances
      ratios = @view numerator[j][:, k]
      expectee = functions
      if denominator_lo[j] > weight_threshold
        upper = 1/denominator_lo[j]
        lower = 1/denominator_up[j]
        result = bound_expectation(
          ratio -> ratio .* (lower, upper), ratios, expectee)
        bootstrapped_upper[k, i, j] = result[2]
        bootstrapped_lower[k, i, j] = result[1]
      end
      if k == n_outputs
        report = () -> []
        next!(progress, showvalues=report)
      end
    end
  end
  (bootstrapped_lower, bootstrapped_upper)
end

# for comparison to the binary-treatment solutions of Kallus or Jesson,
# repeat the SAME procedure just with binary propensities on roughly discretized treatments.
# the above need not be modified at all?

max_stars_repo_name: marmarelis/TreatmentCurves

max_issues_repo_name: marmarelis/TreatmentCurves

max_forks_repo_name: marmarelis/TreatmentCurves


text:
export wPCA, sPCA, rPCA, wpca, spca, rpca

using MultivariateStats: preprocess_mean, pcacov, PCA, fit, transform
using StatsBase: AnalyticWeights
using Statistics: covm

"""
    f = wPCA(x::Weighted, d=2)
    f(x)::Weighted
Works out the `MultivariateStats.PCA` to project `x` into `d` dimensions.
To ignore weights, use `wPCA(array(x), d)` instead.
Returned object can be used as a function `f(y) ≈ transform(f,array(y))`, preserving `weights(y)`.

    sPCA(x, d=2)
    rPCA(y)
Version with `s` saves the function above to a global variable, and also applies it, returning `f(x)`.
Then you can recall it, and apply it to some other data, with `rPCA(y)`.
For example:
```
julia> xx = wrandn(7,50); yy = xx[:, 1:10] |> normalise
julia> plot(xx, sPCA); plot!(yy, rPCA)
```
Now equivalent to using PCA-plot function: `pplot(xx); pplot!(yy)`.
"""
function wPCA(x::Union{Matrix, Weighted}, outdim::Int=2) ## from  fit{T<:AbstractFloat}(::Type{PCA}, X::DenseMatrix{T};
    d, n = size(x)
    mv = preprocess_mean(array(x), nothing)
    aw = AnalyticWeights(weights(x))
    C = covm(array(x), isempty(mv) ? 0 : mv, aw, 2; corrected=true) ## corrected=true from depwarn?
    M = pcacov(C, mv; maxoutdim=outdim, pratio=0.999)
end

(pca::PCA)(x::Weighted) = Weighted(transform(pca, array(x)), weights(x), unclamp(addname(x.opt, "-PC")) )

@doc @doc(wPCA)
function sPCA(x::Union{Matrix, Weighted}, outdim::Int=2)
    global SAVED_PCA = wPCA(x, outdim)
    return rPCA(x)
end

SAVED_PCA = identity ## initially ::Function, later ::MultivariateStats.PCA

@doc @doc(wPCA)
rPCA(x...) = SAVED_PCA(x...) ## rPCA is always of type Function

## lower-case versions:
const wpca = wPCA
const spca = sPCA
const rpca = rPCA


text:
"""

$(TYPEDEF)

Structure to contain a list of all the distances smaller than
the cutoff.

$(TYPEDFIELDS)

"""
struct CutoffDistances

  # To store all distances smaller than the cutoff
  nd::Vector{Int} # Number of small distances (vector size = 1), just to be mutable
  d::Vector{Float64} # All distances smaller than the cutoff
  iat::Vector{Int} # atom of the solute
  jat::Vector{Int} # atom of the solvent
  imol::Vector{Int} # molecule of the solute
  jmol::Vector{Int} # molecule of the solvent

  # Size of the arrays
  maxdim::Vector{Int}

end

# Generator
CutoffDistances(natoms::Int) =
  CutoffDistances(zeros(Int,1), # nd
                  zeros(Float64,natoms), # d
                  zeros(Int,natoms), # iat
                  zeros(Int,natoms), # jat
                  zeros(Int,natoms), # imol
                  zeros(Int,natoms), # jmol
                  [ natoms ]) # maxdim

# Function that zeroes all the values in this structure
function reset!(c::CutoffDistances)
  c.nd[1] = 0
  @. c.d = 0.
  @. c.iat = 0
  @. c.jat = 0
  @. c.imol = 0
  @. c.jmol = 0
  return nothing
end

function reset!(c::Vector{CutoffDistances})
  for i in 1:size(c,1)
    reset!(c[i])
  end
  return nothing
end