text:
# # Rabbit, Fox, Hawk
# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../rabbit_fox_hawk.mp4" type="video/mp4">
# </video>
# ```
# This model is a variation on the [Predator-prey dynamics](@ref) example. It uses a 3-dimensional
# [`ContinuousSpace`](@ref), a realistic terrain for the agents, and pathfinding (with multiple
# pathfinders).
#
# Agents in this model are one of three species of animals: rabbits, foxes and hawks. Rabbits
# eat grass, and are hunted by foxes and hawks. While rabbits and foxes are restricted to
# walk on suitable portions of the map, hawks are capable of flight and can fly over a much
# larger region of the map.
#
# Similar to the [Predator-prey dynamics](@ref) example, agent types are distinguished using a
# `type` field. Agents also have an additional `energy` field, which is consumed to move around
# and reproduce. Eating food (grass or rabbits) replenishes `energy` by a fixed amount.
using Agents, Agents.Pathfinding
using Random
using FileIO # To load images you also need ImageMagick available to your project

mutable struct Animal <: AbstractAgent
    id::Int
    pos::NTuple{3,Float64}
    type::Symbol ## one of :rabbit, :fox or :hawk
    energy::Float64
end

## Some utility functions to create specific types of agents, and find the norm of a vector
Rabbit(id, pos, energy) = Animal(id, pos, :rabbit, energy)
Fox(id, pos, energy) = Animal(id, pos, :fox, energy)
Hawk(id, pos, energy) = Animal(id, pos, :hawk, energy)
norm(vec) = √sum(vec .^ 2)

# The environment is generated from a heightmap: a 2D matrix, where each value denotes the
# height of the terrain at that point. We segregate the model into 4 regions based on the
# height:
# - Anything below `water_level` is waster and cannot be walked on
# - The region between `water_level` and `grass_level` is flatland, that can be walked on
# - The part of the map between `grass_level` and `mountain_level` is too high for animals to
#   walk over, but it can be flown over
# - The terrain above `mountain_level` is completely inaccessible
#
# Grass is the food source for rabbits. It can grow anywhere from `water_level` to `grass_level`.
# The spread of grass across the terrain is specified using a BitArray. A value of
# 1 at a location indicates the presence of grass there, which can be consumed when it is eaten
# by a rabbit. The probability of grass growing is proportional to how close it is to the water.
#
# The `initialize_model` function takes in the URL to our heightmap, the thresholds for the 4
# regions, and some additional parameters for the model. It then creates and returns a model
# with the specified heightmap and containing the specified number of rabbits, foxes and hawks.

function initialize_model(
    heightmap_url,
    water_level = 8,
    grass_level = 20,
    mountain_level = 35;
    n_rabbits = 160,  ## initial number of rabbits
    n_foxes = 30,  ## initial number of foxes
    n_hawks = 30,  ## initial number of hawks
    Δe_grass = 25,  ## energy gained from eating grass
    Δe_rabbit = 30,  ## energy gained from eating one rabbit
    rabbit_repr = 0.06,  ## probability for a rabbit to (asexually) reproduce at any step
    fox_repr = 0.03,  ## probability for a fox to (asexually) reproduce at any step
    hawk_repr = 0.02, ## probability for a hawk to (asexually) reproduce at any step
    rabbit_vision = 6,  ## how far rabbits can see grass and spot predators
    fox_vision = 10,  ## how far foxes can see rabbits to hunt
    hawk_vision = 15,  ## how far hawks can see rabbits to hunt
    rabbit_speed = 1.3, ## movement speed of rabbits
    fox_speed = 1.1,  ## movement speed of foxes
    hawk_speed = 1.2, ## movement speed of hawks
    regrowth_chance = 0.03,  ## probability that a patch of grass regrows at any step
    dt = 0.1,   ## discrete timestep each iteration of the model
    seed = 42,  ## seed for random number generator
)

    ## Download and load the heightmap. The grayscale value is converted to `Float64` and
    ## scaled from 1 to 40
    heightmap = floor.(Int, convert.(Float64, load(download(heightmap_url))) * 39) .+ 1
    ## The x and y dimensions of the pathfinder are that of the heightmap
    dims = (size(heightmap)..., 50)
    ## The region of the map that is accessible to each type of animal (land-based or flying)
    ## is defined using `BitArrays`
    land_walkmap = BitArray(falses(dims...))
    air_walkmap = BitArray(falses(dims...))
    for i in 1:dims[1], j in 1:dims[2]
        ## land animals can only walk on top of the terrain between water_level and grass_level
        if water_level < heightmap[i, j] < grass_level
            land_walkmap[i, j, heightmap[i, j]+1] = true
        end
        ## air animals can fly at any height upto mountain_level
        if heightmap[i, j] < mountain_level
            air_walkmap[i, j, (heightmap[i, j]+1):mountain_level] .= true
        end
    end

    ## Generate the RNG for the model
    rng = MersenneTwister(seed)

    ## Note that the dimensions of the space do not have to correspond to the dimensions
    ## of the pathfinder. Discretisation is handled by the pathfinding methods
    space = ContinuousSpace((100., 100., 50.); periodic = false)

    ## Generate an array of random numbers, and threshold it by the probability of grass growing
    ## at that location. Although this causes grass to grow below `water_level`, it is
    ## effectively ignored by `land_walkmap`
    grass = BitArray(
        rand(rng, dims[1:2]...) .< ((grass_level .- heightmap) ./ (grass_level - water_level)),
    )
    properties = (
        ## The pathfinder for rabbits and foxes
        landfinder = AStar(space; walkmap = land_walkmap),
        ## The pathfinder for hawks
        airfinder = AStar(space; walkmap = air_walkmap, cost_metric = MaxDistance{3}()),
        Δe_grass = Δe_grass,
        Δe_rabbit = Δe_rabbit,
        rabbit_repr = rabbit_repr,
        fox_repr = fox_repr,
        hawk_repr = hawk_repr,
        rabbit_vision = rabbit_vision,
        fox_vision = fox_vision,
        hawk_vision = hawk_vision,
        rabbit_speed = rabbit_speed,
        fox_speed = fox_speed,
        hawk_speed = hawk_speed,
        heightmap = heightmap,
        grass = grass,
        regrowth_chance = regrowth_chance,
        water_level = water_level,
        grass_level = grass_level,
        dt = dt,
    )

    model = ABM(Animal, space; rng, properties)

    ## spawn each animal at a random walkable position according to its pathfinder
    for _ in 1:n_rabbits
        add_agent_pos!(
            Rabbit(
                nextid(model), ## Using `nextid` prevents us from having to manually keep track
                               ## of animal IDs
                random_walkable(model, model.landfinder),
                rand(model.rng, Δe_grass:2Δe_grass),
            ),
            model,
        )
    end
    for _ in 1:n_foxes
        add_agent_pos!(
            Fox(
                nextid(model),
                random_walkable(model, model.landfinder),
                rand(model.rng, Δe_rabbit:2Δe_rabbit),
            ),
            model,
        )
    end
    for _ in 1:n_hawks
        add_agent_pos!(
            Hawk(
                nextid(model),
                random_walkable(model, model.airfinder),
                rand(model.rng, Δe_rabbit:2Δe_rabbit),
            ),
            model,
        )
    end

    return model
end

# The `animal_step!` function dispatches to the proper function depending on the type of agent.
# The stepping functions for each type of agent are similar: They lose energy per step, and
# die if their energy ever reaches 0. They also have a random probability to reproduce at an
# iteration. Agents all move towards their food. In the case of rabbits, they also move away
# from any nearby predators.

function animal_step!(animal, model)
    if animal.type == :rabbit
        rabbit_step!(animal, model)
    elseif animal.type == :fox
        fox_step!(animal, model)
    else
        hawk_step!(animal, model)
    end
end

# Rabbits eat grass at their position, if it exists. If they see a predator, they run away.
# The direction in which they flee is dependent on all predators in their vision, with closer
# ones contributing more to the chosen direction. If there are no predators to flee from,
# rabbits walk around randomly.

function rabbit_step!(rabbit, model)
    ## Eat grass at this position, if any
    if get_spatial_property(rabbit.pos, model.grass, model) == 1
        model.grass[get_spatial_index(rabbit.pos, model.grass, model)] = 0
        rabbit.energy += model.Δe_grass
    end

    ## The energy cost at each step corresponds to the amount of time that has passed
    ## since the last step
    rabbit.energy -= model.dt
    ## All animals die if their energy reaches 0
    if rabbit.energy <= 0
        kill_agent!(rabbit, model, model.landfinder)
        return
    end

    ## Get a list of positions of all nearby predators
    predators = [
        x.pos for x in nearby_agents(rabbit, model, model.rabbit_vision) if
            x.type == :fox || x.type == :hawk
            ]
    ## If the rabbit sees a predator and isn't already moving somewhere
    if !isempty(predators) && is_stationary(rabbit, model.landfinder)
        ## Try and get an ideal direction away from predators
        direction = (0., 0., 0.)
        for predator in predators
            ## Get the direction away from the predator
            away_direction = (rabbit.pos .- predator)
            ## In case there is already a predator at our location, moving anywhere is
            ## moving away from it, so it doesn't contribute to `direction`
            all(away_direction .≈ 0.) && continue
            ## Add this to the overall direction, scaling inversely with distance.
            ## As a result, closer predators contribute more to the direction to move in
            direction = direction .+ away_direction ./ norm(away_direction) ^ 2
        end
        ## If the only predator is right on top of the rabbit
        if all(direction .≈ 0.)
            ## Move anywhere
            chosen_position = random_walkable(rabbit.pos, model, model.landfinder, model.rabbit_vision)
        else
            ## Normalize the resultant direction, and get the ideal position to move it
            direction = direction ./ norm(direction)
            ## Move to a random position in the general direction of away from predators
            position = rabbit.pos .+ direction .* (model.rabbit_vision / 2.)
            chosen_position = random_walkable(position, model, model.landfinder, model.rabbit_vision / 2.)
        end
        plan_route!(rabbit, chosen_position, model.landfinder)
    end

    ## Reproduce with a random probability, scaling according to the time passed each
    ## step
    rand(model.rng) <= model.rabbit_repr * model.dt && reproduce!(rabbit, model)

    ## If the rabbit isn't already moving somewhere, move to a random spot
    if is_stationary(rabbit, model.landfinder)
        plan_route!(
            rabbit,
            random_walkable(rabbit.pos, model, model.landfinder, model.rabbit_vision),
            model.landfinder
        )
    end

    ## Move along the route planned above
    move_along_route!(rabbit, model, model.landfinder, model.rabbit_speed, model.dt)
end

# Foxes hunt for rabbits, and eat rabbits within a unit radius of its position.

function fox_step!(fox, model)
    ## Look for nearby rabbits that can be eaten
    food = [x for x in nearby_agents(fox, model) if x.type == :rabbit]
    if !isempty(food)
        kill_agent!(rand(model.rng, food), model, model.landfinder)
        fox.energy += model.Δe_rabbit
    end

    ## The energy cost at each step corresponds to the amount of time that has passed
    ## since the last step
    fox.energy -= model.dt
    ## All animals die once their energy reaches 0
    if fox.energy <= 0
        kill_agent!(fox, model, model.landfinder)
        return
    end

    ## Random chance to reproduce every step
    rand(model.rng) <= model.fox_repr * model.dt && reproduce!(fox, model)

    ## If the fox isn't alreadu moving somewhere
    if is_stationary(fox, model.landfinder)
        ## Look for any nearby rabbits
        prey = [x for x in nearby_agents(fox, model, model.fox_vision) if x.type == :rabbit]
        if isempty(prey)
            ## Move anywhere if no rabbits were found
            plan_route!(
                fox,
                random_walkable(fox.pos, model, model.landfinder, model.fox_vision),
                model.landfinder,
            )
            return
        end
        ## Move toward a random rabbit
        plan_route!(fox, rand(model.rng, map(x -> x.pos, prey)), model.landfinder)
    end

    move_along_route!(fox, model, model.landfinder, model.fox_speed, model.dt)
end

# Hawks function similarly to foxes, except they can also fly. They dive down for prey and
# fly back up after eating it.

function hawk_step!(hawk, model)
    ## Look for rabbits nearby
    food = [x for x in nearby_agents(hawk, model) if x.type == :rabbit]
    if !isempty(food)
        ## Eat (kill) the rabbit
        kill_agent!(rand(model.rng, food), model, model.airfinder)
        hawk.energy += model.Δe_rabbit
        ## Fly back up
        plan_route!(hawk, hawk.pos .+ (0., 0., 7.), model.airfinder)
    end

    ## The rest of the stepping function is similar to that of foxes, except hawks use a
    ## different pathfinder
    hawk.energy -= model.dt
    if hawk.energy <= 0
        kill_agent!(hawk, model, model.airfinder)
        return
    end

    rand(model.rng) <= model.hawk_repr * model.dt && reproduce!(hawk, model)

    if is_stationary(hawk, model.airfinder)
        prey = [x for x in nearby_agents(hawk, model, model.hawk_vision) if x.type == :rabbit]
        if isempty(prey)
            plan_route!(
                hawk,
                random_walkable(hawk.pos, model, model.airfinder, model.hawk_vision),
                model.airfinder,
            )
        else
            plan_route!(hawk, rand(model.rng, map(x -> x.pos, prey)), model.airfinder)
        end
    end

    move_along_route!(hawk, model, model.airfinder, model.hawk_speed, model.dt)
end

# This function is called when an animal reproduces. The animal loses half its energy, and
# a copy of it is created and added to the model.

function reproduce!(animal, model)
    animal.energy = ceil(Int, animal.energy / 2)
    add_agent_pos!(Animal(nextid(model), animal.pos, animal.type, animal.energy), model)
end

# The model stepping function simulates the growth of grass

function model_step!(model)
    ## To prevent copying of data, obtain a view of the part of the grass matrix that
    ## doesn't have any grass, and grass can grow there
    growable = view(
        model.grass,
        model.grass .== 0 .& model.water_level .< model.heightmap .<= model.grass_level,
    )
    ## Grass regrows with a random probability, scaling with the amount of time passing
    ## each step of the model
    growable .= rand(model.rng, length(growable)) .< model.regrowth_chance * model.dt
end

# ## Visualization
# Now we use `InteractiveDynamics` to create a visualization of the model running in 3D space
#
# The agents are color-coded according to their `type`, to make them easily identifiable in
# the visualization.

using InteractiveDynamics
using GLMakie
GLMakie.activate!() # hide

animalcolor(a) =
    if a.type == :rabbit
        :brown
    elseif a.type == :fox
        :orange
    else
        :blue
    end

# We use `surface!` to plot the terrain as a mesh, and colour it using the `:terrain`
# colormap. Since the heightmap dimensions don't correspond to the dimensions of the space,
# we explicitly provide ranges to specify where the heightmap should be plotted.
function static_preplot!(ax, model)
    surface!(
        ax,
        (100/205):(100/205):100,
        (100/205):(100/205):100,
        model.heightmap;
        colormap = :terrain
    )
end

# Passing in a sample heightmap to the `initialize_model` function we created returns the generated
# model.
heightmap_url =
    "https://raw.githubusercontent.com/JuliaDynamics/" *
    "JuliaDynamics/master/videos/agents/rabbit_fox_hawk_heightmap.png"
model = initialize_model(heightmap_url)

abm_video(
    "rabbit_fox_hawk.mp4",
    model,
    animal_step!,
    model_step!;
    resolution = (700, 700),
    frames = 300,
    framerate = 20,
    ac = animalcolor,
    as = 1.0,
    static_preplot!
)

nothing # hide
# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../rabbit_fox_hawk.mp4" type="video/mp4">
# </video>
# ```


text:
L = 20
labels = rand(Bool, L);
scores = add_noise.(labels, 0.6)

curve = PR(labels, scores)
B = BinaryMetrics(labels, scores)
curve = PR(B)
println(curve)

OP = op(B)
@test op(curve) == OP
@test precision(curve) == precision(B,OP)
@test recall(curve) == tpr(B,OP)

# comparing with sklearn
data = CSV.read(joinpath(@__DIR__, "data", "data.csv"))
pr_sklearn = CSV.read(joinpath(@__DIR__, "data", "pr_sklearn.csv"))

precision_sklearn = pr_sklearn[:,:precision]
recall_sklearn = pr_sklearn[:,:recall]

curve = PR(data.labels, data.predictions)
@test norm(recall(curve)-reverse(recall_sklearn)[2:end]) <= 1e-8
@test norm(precision(curve)-reverse(precision_sklearn)[2:end]) <= 1e-8
@test auc(curve)-0.7814245832778571 < 1e-8

ext: jl


text:
#########################################################
## Randomized singular value and eigen- decompositions ##
########################################################

# This file provides a rudimentary implementation of the randomized singular
# value decomposition and spectral (eigen-) decomposition as described in
# \cite{Halko2011}.
#
# Reference:
# @article{Halko2011,
#    author = {Halko, N and Martinsson, P G and Tropp, J A},
#    doi = {10.1137/090771806},
#    journal = {SIAM Review},
#    month = jan,
#    number = {2},
#    pages = {217--288},
#    title = {Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions},
#    volume = {53},
#    year = {2011}
# }

import Base.LinAlg: Eigen, SVD

export rsvd, reig

@doc doc"""
Computes the partial singular value decomposition of `A` using a randomized
algorithm.

Inputs:

    `A`: input matrix
    `n`: Number of singular value/vector pairs to find
    `p`: Number of extra vectors to include in computation

Output:

    `F`: An `SVD` `Factorization` object

Warning:

    This variant of the randomized singular value decomposition is the most
    commonly found implementation but is not recommended for accurate
    computations, as it often has trouble finding the `n` largest singular pairs,
    but rather finds `n` large singular pairs which may not necessarily be the
    largest.

Implementation note:

    This function calls `rrange()`, which uses naive randomized rangefinding to
    compute a basis for a subspace of dimension `n` (Algorithm 4.1 of
    \cite{Halko2011}), followed by `svdfact_restricted()`, which computes the
    exact SVD factorization on the restriction of `A` to this randomly selected
    subspace (Algorithm 5.1 of \cite{Halko2011}).

    Alternatively, you can mix and match your own randomized algorithm using
    any of the randomized range finding algorithms to find a suitable subspace
    and feeding the result to one of the routines that computes the SVD
    restricted to that subspace.
""" ->
function rsvd(A, n::Int, p::Int=0)
    Q = rrange(A, n, p=p)
    svdfact_restricted(A, Q)
end

@doc doc"""
Like `rsvd`, but returns only the singular values.

Inputs:

    as for `rsvd`.

Output:

    A vector containing the estimated singular values of `A`

""" ->
function rsvdvals(A, n::Int, p::Int=0)
    Q = rrange(A, n, p=p)
    svdvals_restricted(A, Q)
end

@doc doc"""
Computes an orthonormal basis for a subspace of `A` of dimension `l` using
naive randomized rangefinding.

Inputs:

    `A`: Input matrix. Must support `size(A)` and premultiply
    `l`: The number of basis vectors to compute
    `p`: The oversampling parameter. The number of extra basis vectors to use
       in the computation, which get discarded at the end.

Output:

    `Q`: A dense matrix of dimension `size(A,1)` x l containing the basis
    vectors of the computed subspace of `A`

Reference:

    Algorithm 4.1 of \cite{Halko2011}

Warning:

    The Reference explicitly discourages using this algorithm.

Implementation note:

    Whereas \cite{Halko2011} recommends classical Gram-Schmidt with double
    reorthogonalization, we instead compute the basis with qrfact(), which
    for dense A computes the QR factorization using Householder reflectors.
""" ->
function rrange(A, l::Int; p::Int=0)
    p≥0 || error()
    m, n = size(A)
    if l > m
	warn("Cannot find $l linearly independent vectors of $m x $n matrix")
	warn("Truncating to l=$m, p=0")
        l=m
	p=0
    end
    Ω = randn(n, l+p)
    Y = A*Ω
    Q = full(qrfact!(Y)[:Q])
    Q = p==0 ? Q : Q[:,1:l]
    @assert l==size(Q, 2)
    Q
end

@doc doc"""
Computes an orthogonal basis for a subspace of `A` of dimension `l` using adaptive
randomized rangefinding.

Similar to `rrange()`, but determines the oversampling parameter adaptively given
a threshold `ϵ`.

Inputs:

    `A`: Input matrix. Must support `size(A)` and premultiply
    `l`: The number of basis vectors to compute
    `ϵ`: Threshold to determine adaptive fitting
    `maxiter`: Maximum number of iterations to run. (Default: 10.)

Output:

    `Q`: A dense matrix of dimension `size(A,1)` x l containing the basis
         vectors of the computed subspace of `A`

Reference:

    Algorithm 4.2 of \cite{Halko2011}

""" ->
function rrange_adaptive(A, r::Integer, ϵ::Real=eps(); maxiter::Int=10)
    m, n = size(A)

    Ω = randn(n,r)
    Y = A*Ω
    Q = zeros(m,0)

    const tol=ϵ/(10*√(2/π))
    for j=1:maxiter
        Tol = maximum([norm(Y[:,i]) for i=j:(j+r-1)])
        Tol > tol || break

        y = sub(Y,:,j)
        y = Y[:,j] = y - Q*(Q'y)
        q = y/norm(y)
        Q = [Q q]

        ω = randn(n)
        y = A*ω
        y = y - Q*(Q'y)
        Y = [Y y]
        Yb = sub(Y, :, (j+1):(j+r-1))
        Yb = Y[:, (j+1):(j+r-1)] = Yb - q * q'Yb

        j==maxiter && warn("Maximum number of iterations reached with norm $Tol > $tol")
    end
    Q
end

@doc doc"""
Computes an orthonormal basis for a subspace of `A` of dimension `l` using
randomized rangefinding by subspace iteration.

Inputs:

    `A` : Input matrix. Must support `size(A)` and premultiply
    `l` : The number of basis vectors to compute
    `At`: (Optional) the transpose of `A`. Defaults to `A'`.
    `p` : The oversampling parameter. The number of extra basis vectors to use
          in the computation, which get discarded at the end. (Default: 0)
    `q` : Number of subspace iterations. (Default: 0)

Note:
    Running with `q=0` is functionally equivalent to `rrange()`

Output:
    `Q`: A dense matrix of dimension `size(A,1)` x `l` containing the basis
         vectors of the computed subspace of `A`

Reference:

    Algorithm 4.4 of \cite{Halko2011}

Implementation note:

    Whereas the Reference recommends classical Gram-Schmidt with double
    reorthogonalization, we instead compute the basis with `qrfact()`, which
    for dense A computes the QR factorization using Householder reflectors.
""" ->
function rrange_si(A, l::Int; At=A', p::Int=0, q::Int=0)
    const basis=_->full(qrfact(_)[:Q])
    p>=0 || error()
    n = size(A, 2)
    Ω = randn(n,l+p)
    Y = A*Ω
    Q = basis(Y)
    for iter=1:q #Do some power iterations
        Ỹ=At*Q
        Q̃=basis(Ỹ)
        Y=A*Q̃
        Q=basis(Y)
    end
    p==0 ? Q : Q[:,1:l]
end

@doc doc"""
Computes an orthonormal basis for a subspace of `A` of dimension `l` using
naive randomized rangefinding using stochastic randomized Fourier transforms.

Inputs:

    `A` : Input matrix. Must support `size(A)` and premultiply
    `l` : The number of basis vectors to compute
    `p` : The oversampling parameter. The number of extra basis vectors to use
          in the computation, which get discarded at the end. (Default: 0)

Output:
    `Q`: A dense matrix of dimension `size(A,1)` x `l` containing the basis
         vectors of the computed subspace of `A`

Note:

    Similar to `rrange()`, but does not use Gaussian random matrices.

Reference:

    Algorithm 4.5 of \cite{Halko2011}

Implementation note:

    Whereas the Reference recommends classical Gram-Schmidt with double
    reorthogonalization, we instead compute the basis with `qrfact()`, which
    for dense A computes the QR factorization using Householder reflectors.

""" ->
function rrange_f(A, l::Int; p::Integer=0)
    p ≥ 0 || error()
    n = size(A, 2)
    Ω = srft(l+p)
    Y = A*Ω
    Q = full(qrfact!(Y)[:Q])
    p==0 ? Q : Q[:,1:l]
end

@doc doc"""
Computes the SVD factorization of `A` restricted to the subspace spanned by `Q`
using exact projection.

Inputs:

    `A`: Input matrix. Must support postmultiply
    `Q`: Matrix containing basis vectors of the subspace whose restriction to is
       desired

Output:

    `F`: An `SVD` `Factorization` object

Reference:

    Algorithm 5.1 of \cite{Halko2011}
""" ->
function svdfact_restricted(A, Q)
    B=Q'A
    S=svdfact!(B)
    SVD(Q*S[:U], S[:S], S[:Vt])
end

@doc doc"""
Like `svdfact_restricted`, but returns only the singular values.

Inputs:

    as for `svdfact_restricted`.

Output:

    `v`: A vector containing the estimated singular values of `A`
""" ->
function svdvals_restricted(A, Q)
    B=Q'A
    S=svdvals!(B)
end

@doc doc"""
Computes the SVD factorization of `A` restricted to the subspace spanned by `Q`
using row extraction.

Inputs:

    `A`: Input matrix. Must support postmultiply
    `Q`: Matrix containing basis vectors of the subspace whose restriction to is
         desired. Need not be orthogonal or normalized.

Note:

    \cite[Remark 5.2]{Halko2011} recommends input of `Q` of the form `Q=A*Ω`
    where `Ω` is a sample computed by `randn(n,l)` or even `srft(l)`.

Output:

    `F`: An `SVD` `Factorization` object

Note:

    A faster but less accurate variant of `svdfact_restricted()` which uses the
    interpolative decomposition `idfact()`.

Reference:

    Algorithm 5.2 of \cite{Halko2011}
""" ->
function svdfact_re(A, Q)
    F = idfact(Q)
    X, J = F[:B], F[:P]
    R′, W′= qr(A[J, :])
    Z = X*R′
    S=svdfact(Z)
    SVD(S[:U], S[:S], S[:Vt]*W′)
end

@doc doc"""
Computes the spectral (`Eigen`) factorization of `A` restricted to the subspace
spanned by `Q` using row extraction.

Inputs:

    `A`: Input matrix. Must be `Hermitian` and support pre- and post-multiply
    `Q`: Orthonormal matrix containing basis vectors of the subspace whose
         restriction to is desired.

Output:

    `F`: An `Eigen` `Factorization` object

Reference:

    Algorithm 5.3 of \cite{Halko2011}
""" ->
function eigfact_restricted(A::Hermitian, Q)
    B = Q'A*Q
    E = eigfact!(B)
    Eigen(E[:values], Q*E[:vectors])
end

@doc doc"""
Computes the spectral (`Eigen`) factorization of `A` restricted to the subspace
spanned by `Q` using row extraction.

Inputs:

    `A`: Input matrix. Must be `Hermitian` and support pre- and post-multiply
    `Q`: Matrix containing basis vectors of the subspace whose restriction to is
         desired. Need not be orthogonal or normalized.

Note:

    \cite[Remark 5.2]{Halko2011} recommends input of `Q` of the form `Q=A*Ω`
    where `Ω` is a sample computed by `randn(n,l)` or even `srft(l)`.

Output:

    `F`: An `Eigen` `Factorization` object

Note:

    A faster but less accurate variant of `eigfact_restricted()` which uses the
    interpolative decomposition `idfact()`.

Reference:

    Algorithm 5.4 of \cite{Halko2011}
""" ->
function eigfact_re(A::Hermitian, Q)
    X, J = idfact(Q)
    F = qrfact!(X)
    V, R = F[:Q], F[:R]
    Z=R*A[J, J]*R'
    E=eigfact(Z)
    Eigen(E[:values], V*E[:vectors])
end

@doc doc"""
Computes the spectral (`Eigen`) factorization of `A` restricted to the subspace
spanned by `Q` using the Nyström method.

Inputs:

    `A`: Input matrix. Must be positive semidefinite.
    `Q`: Orthonormal matrix containing basis vectors of the subspace whose
         restriction to is desired.

Output:

    `F`: An `Eigen` `Factorization` object

Note:

    More accurate than `eigfact_restricted()` but is restricted to matrices
    that can be Cholesky decomposed.

Reference:

    Algorithm 5.5 of \cite{Halko2011}
""" ->
function eigfact_nystrom(A, Q)
    B₁=A*Q
    B₂=Q'*B₁
    C=cholfact!(B₂)
    F=B₁*inv(C)
    S=svdfact!(F)
    Eigen(S[:S].^2, S[:U])
end

@doc doc"""
Computes the spectral (`Eigen`) factorization of `A` using only one matrix
product involving `A`.

Inputs:

    `A`: Input matrix.
    `Ω`: Sample matrix for the column space, e.g. `randn(n, l)` or `srft(l)`
    `Ω̃;: Sample matrix for the row space. Not neeeded for `Hermitian` matrices
    `At`: Computes transpose of input matrix. Default: `A'`

Output:

    `F`: An `Eigen` `Factorization` object

Reference:

    Algorithm 5.6 of \cite{Halko2011}
""" ->
function eigfact_onepass(A::Hermitian, Ω)
    Y=A*Ω; Q = full(qrfact!(Y)[:Q])
    B=(Q'Y)\(Q'Ω)
    E=eigfact!(B)
    Eigen(E[:values], Q*E[:vectors])
end

function eigfact_onepass(A, Ω, Ω̃; At=A')
    Y=A *Ω; Q = full(qrfact!(Y)[:Q])
    Ỹ=At*Ω; Q̃ = full(qrfact!(Ỹ)[:Q])
    #Want least-squares solution to (5.14 and 5.15)
    B=(Q'Y)\(Q̃'Ω)
    B̃=(Q̃'Ỹ)\(Q'Ω̃)
    #Here is a very very very hacky way to solve the problem
    B=0.5(B + B̃')
    E=eigfact!(B)
    Eigen(E[:values], Q*E[:vectors])
end

@doc doc"""
Computes the spectral (`Eigen`) decomposition of `A` using a randomized
algorithm.

Inputs:

    `A`: input matrix
    `n`: Number of eigenpairs to find

Output:

    `F`: An `Eigen` `Factorization` object

Implementation Note:

    This is a wrapper around `eigfact_onepass()` which uses the randomized
    samples found using `srft(l)`.
""" ->
reig(A::Hermitian, l::Int) = eigfact_onepass(A, srft(l))
reig(A, l::Int) = eigfact_onepass(A, srft(l), srft(l))


text:
using TestExtension
using Test

@testset "@test" begin
    @test true 37+5
    @test false 37+5 end

end


text:
# %% Modules

include("../../Modules/Old_net/DASTDP.jl")
include("../../Modules/Old_net/NeuronModel.jl")

using .DASTDP
using .NeuronModel
using Statistics

# %% Constants

const T = 800
randarray = rand(1:797,50)
const A = randarray
const B = randarray .+ 1
const S = randarray .+ 2

syn_SA_tmp = []
syn_SB_tmp = []

for pre_neuron in pre[A]
    indexes_A = [cartind for cartind in pre_neuron if cartind[1] in S]
    append!(syn_SA_tmp,indexes_A)
end

for pre_neuron in pre[B]
    indexes_B = [cartind for cartind in pre_neuron if cartind[1] in S]
    append!(syn_SB_tmp,indexes_B)
end

const syn_SA = syn_SA_tmp
const syn_SB = syn_SB_tmp

# learning rate related
const DAinc = 0.5
const STDPinc = 0.1
const sm = 4
const LTPinc = 1.0 # 1.35
const LTDinc = 1.8 # 2.15
const I_inc = 18 # 150

# %% Network structure

mutable struct NeuralNet
    v::Array{Float64}
    u::Array{Float64}
    s::Array{Float64,2}
    sd::Array{Float64,2}
    STDP::Array{Float64,2}
    firings::Array{Int64,2}
    DA::Float64
    rew::Array{Int64}
    card::Array{Int64}                       # card[A,B]
    target::Int64                            # target = 1 : group A , target = 2 : group B
    I::Array{Float64}
    shist::Array{Float64,2}
    phist::Array{Float64,2}

    function NeuralNet()
        v = -65.0*ones(N)
        u = 0.2*v
        s = vcat(1.0 .* ones(Ne,M),-1.0 .* ones(Ni,M))
        sd = 0.0 .* zeros(N,M)
        STDP = 0.0 .* zeros(N,1001+D)
        firings = [-D 0]
        DA = 0.0
        rew = []
        card = [0,0]
        target = 1
        I = Float64[]
        shist = zeros(10*T,6)
        phist = zeros(div(T,10), 2)
        new(v,u,s,sd,STDP,firings,DA,rew,card,target,I,shist,phist)
    end
end

# %% Functions

function stimuli_fire(I::Array{Float64})
    I[S] = I[S] .+ I_inc
    return I
end

function cardAB(fired::Array{Int64},card::Array{Int64})
    for neuron in fired
        if neuron in A
            card[1] += 1
        elseif neuron in B
            card[2] += 1
        end
    end
    return card
end

function reward(rew::Array{Int64},card::Array{Int64},target::Int64,time::Int64)
    if card[target]>card[3-target]
        delay = 1000*(card[3-target]/card[target])
        append!(rew,time + 1 + round(delay))
    end
    return rew
end

function probAB(card::Array{Int64},phist::Array{Float64,2},trial::Int64)
    if sum(card)!=0
        pA_new = card[1]/(sum(card))
        pB_new = card[2]/(sum(card))
    else
        pA_new = 0
        pB_new = 0
    end
    if trial>=1
        ind = min(trial,50)
        meanwindA = mean(phist[trial+1-ind:trial,1])
        meanwindB = mean(phist[trial+1-ind:trial,2])
        pA_new = mean([meanwindA,pA_new])
        pB_new = mean([meanwindB,pB_new])
    end
    return pA_new,pB_new
end

function means(s::Array{Float64,2},sd::Array{Float64,2})
    @inbounds mean_SA = mean(s[syn_SA])
    @inbounds mean_SB = mean(s[syn_SB])
    @inbounds max_SA = maximum(s[syn_SA])
    @inbounds min_SA = minimum(s[syn_SA])
    @inbounds max_SB = maximum(s[syn_SB])
    @inbounds min_SB = minimum(s[syn_SB])
    return mean_SA,mean_SB,max_SA,min_SA,max_SB,min_SB
end

# %% Main loop

net = NeuralNet()

@inbounds for sec in 0:T-1
    @time @inbounds for msec in 1:1000
        net.I = 13*(rand(N).-0.5)
        time = 1000*sec+msec
        trial = div(sec,10)+1
        fired = findall(x->x>=thresh,net.v)
        net.v,net.u = izhikevicmodel_fire(net.v,net.u,fired)
        net.STDP = STDP_fire(net.STDP,fired,msec)
        net.sd = LTP(net.STDP,net.sd,fired,msec,LTPinc)
        net.firings = vcat(net.firings,hcat(msec.*ones(length(fired)),fired))
        net.I,net.sd = LTD(net.STDP,net.sd,net.s,net.firings,net.I,msec,LTDinc)

        if msec==1
            net.I = stimuli_fire(net.I)
            net.card = [0,0]
        end
        if msec<=20
            net.card = cardAB(fired,net.card)
        end
        if msec==20
            net.rew = reward(net.rew,net.card,net.target,time)
            net.phist[trial,:] .= probAB(net.card,net.phist,trial)
        end

        net.v,net.u = izhikevicmodel_step(net.v,net.u,net.I)
        net.STDP,net.DA = DA_STDP_step(net.STDP,net.DA,msec)
        if msec%10==0
            net.s,net.sd = synweight_step(net.sd,net.s,net.DA,sm)
            if msec%100==0
                net.shist[div(time,100),:] .= means(net.s,net.sd)
            end
        end
        net.DA = DA_inc(net.rew,net.DA,time)

    end
    if sec==400
        net.target = 2
    end
    if sec%100==0
        println("spike frequency = $(length(net.firings))")
        print("\rsec = $sec")
    end
    net.STDP,net.firings = time_reset(net.STDP,net.firings)
end

# %% Plot probability
using Plots
gr()
x1 = collect(1:length(net.phist[:,1]))
y1 = net.phist[:,1]
x2 = x1
y2 = net.phist[:,2]
fig = plot()
plot!(x1,y1,color="blue",label="Group A", legend = true)
plot!(x2,y2,color="green",label="Group B", legend = true)
xlabel!("Trial")
ylabel!("Probability of response")

# %% Plot synaptic weight in A and B
using Plots
gr()
x1 = 0.1.*collect(1:length(net.shist[:,1]))
y1 = net.shist[:,1]
x2 = x1
y2 = net.shist[:,2]
fig = plot()
plot!(x1,y1,color="blue",label="S->A", legend = true)
plot!(x2,y2,color="green",label="S->B", legend = true)
xlabel!("Time (sec)")
ylabel!("mean synaptic weigth (mV)")

# %% Plot enveloppe s for A and B
using Plots
gr()
x1 = 0.1.*collect(1:length(net.shist[:,1]))
y1 = net.shist[:,3]
x2 = x1
y2 = net.shist[:,5]
x3 = x1
y3 = net.shist[:,4]
x4 = x1
y4 = net.shist[:,6]
fig = plot()
plot!(x1,y1,color="blue",label="S->A", legend = true)
plot!(x2,y2,color="green",label="S->B", legend = true)
plot!(x3,y3,color="blue")
plot!(x4,y4,color="green")
xlabel!("Time (sec)")
ylabel!("mean synaptic weigth derivative")

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
using LinearAlgebra

mutable struct RadialBasis{F,Q,X,Y,L,U,C,S} <: AbstractSurrogate
    phi::F
    dim_poly::Q
    x::X
    y::Y
    lb::L
    ub::U
    coeff::C
    scale_factor::S
end

mutable struct RadialFunction{Q,P}
    q::Q
    phi::P
end

linearRadial = RadialFunction(0,z->norm(z))
cubicRadial = RadialFunction(1,z->norm(z)^3)
multiquadricRadial = RadialFunction(1,z->sqrt(norm(z)^2+1))
thinplateRadial = RadialFunction(2,z->norm(z)^2*log(norm(z)))

"""
    RadialBasis(x,y,lb::Number,ub::Number; rad::RadialFunction = linearRadial,scale::Real=1.0)

Constructor for RadialBasis surrogate.
"""
function RadialBasis(x, y, lb::Number, ub::Number; rad::RadialFunction=linearRadial, scale_factor::Real=1.0)
    q = rad.q
    phi = rad.phi
    coeff = _calc_coeffs(x, y, lb, ub, phi, q,scale_factor)
    return RadialBasis(phi, q, x, y, lb, ub, coeff,scale_factor)
end

"""
RadialBasis(x,y,lb,ub,rad::RadialFunction, scale_factor::Float = 1.0)

Constructor for RadialBasis surrogate
"""
function RadialBasis(x, y, lb, ub; rad::RadialFunction = linearRadial, scale_factor::Real=1.0)
    q = rad.q
    phi = rad.phi
    coeff = _calc_coeffs(x, y, lb, ub, phi, q, scale_factor)
    return RadialBasis(phi, q, x, y, lb, ub, coeff,scale_factor)
end

function _calc_coeffs(x, y, lb, ub, phi, q, scale_factor)
    D = _construct_rbf_interp_matrix(x, first(x), lb, ub, phi, q, scale_factor)
    Y = _construct_rbf_y_matrix(y, first(y), length(y) + q)
    coeff = D \ Y
    return coeff
end

function _construct_rbf_interp_matrix(x, x_el::Number, lb, ub, phi, q, scale_factor)
    n = length(x)
    m = n + q
    D = zeros(eltype(x_el), m, m)

    @inbounds for i = 1:n
        for j = 1:n
            D[i,j] = phi( (x[i] .- x[j]) ./ scale_factor )
        end
        if i <= n
            for k = 1:q
                D[i,n+k] = _scaled_chebyshev(x[i], k-1, lb, ub)
            end
        end
    end
    D_sym = Symmetric(D, :U)
    return D_sym
end

function _construct_rbf_interp_matrix(x, x_el, lb, ub, phi, q, scale_factor)
    n = length(x)
    nd = length(x_el)
    central_point = _center_bounds(x_el, lb, ub)
    sum_half_diameter = sum((ub[k]-lb[k])/2 for k = 1:nd)
    mean_half_diameter = sum_half_diameter / nd
    m = n+q
    D = zeros(eltype(x_el), m, m)

    @inbounds for i = 1:n
        for j = 1:n
            D[i,j] = phi( (x[i] .- x[j]) ./ scale_factor)
        end
        if i < n + 1
            for k = 1:q
                D[i,n+k] = centralized_monomial(x[i], k-1, mean_half_diameter, central_point)
            end
        end
    end
    D_sym = Symmetric(D, :U)
    return D_sym
end

_construct_rbf_y_matrix(y, y_el::Number, m) = [i <= length(y) ? y[i] : zero(y_el) for i = 1:m]
_construct_rbf_y_matrix(y, y_el, m) = [i <= length(y) ? y[i][j] : zero(first(y_el)) for i=1:m, j=1:length(y_el)]

"""
    centralized_monomial(vect,alpha,mean_half_diameter,central_point)

Returns the value at point vect[] of the alpha degree monomial centralized.

#Arguments:
-'vect': vector of points i.e [x,y,...,w]
-'alpha': degree
-'mean_half_diameter': half diameter of the domain
-'central_point': central point in the domain
"""
function centralized_monomial(vect,alpha,mean_half_diameter,central_point)
    if iszero(alpha) return one(eltype(vect)) end
    centralized_product = prod(vect .- central_point)
    return (centralized_product / mean_half_diameter)^alpha
end

"""
Calculates current estimate of value 'val' with respect to the RadialBasis object.
"""
function (rad::RadialBasis)(val)
    approx = _approx_rbf(val, rad)
    return _match_container(approx, first(rad.y))
end

function _approx_rbf(val::Number, rad)
    n = length(rad.x)
    q = rad.dim_poly
    lb = rad.lb
    ub = rad.ub
    approx = zero(rad.coeff[1, :])
    for i = 1:n
        approx += rad.coeff[i, :] * rad.phi( (val .- rad.x[i]) / rad.scale_factor)
    end
    for k = 1:q
        approx += rad.coeff[n+k, :] * _scaled_chebyshev(val, k-1, lb, ub)
    end
    return approx
end
function _approx_rbf(val, rad)
    n = length(rad.x)
    d = length(rad.x[1])
    q = rad.dim_poly
    lb = rad.lb
    ub = rad.ub
    sum_half_diameter = sum((ub[k]-lb[k])/2 for k = 1:d)
    mean_half_diameter = sum_half_diameter/d
    central_point = _center_bounds(first(rad.x), lb, ub)

    approx = zero(rad.coeff[1, :])
    @views approx += sum( rad.coeff[i, :] * rad.phi( (val .- rad.x[i]) ./rad.scale_factor) for i = 1:n)
    for k = 1:q
        @views approx += rad.coeff[n+k, :] .* centralized_monomial(val, k-1, mean_half_diameter, central_point)
    end
    return approx
end

_scaled_chebyshev(x, k, lb, ub) = cos(k*acos(-1 + 2*(x-lb)/(ub-lb)))
_center_bounds(x::Tuple, lb, ub) = ntuple(i -> (ub[i] - lb[i])/2, length(x))
_center_bounds(x, lb, ub) = (ub .- lb) ./ 2

"""
    add_point!(rad::RadialBasis,new_x,new_y)

Add new samples x and y and update the coefficients. Return the new object radial.
"""
function add_point!(rad::RadialBasis,new_x,new_y)
    if (length(new_x) == 1 && length(new_x[1]) == 1) || ( length(new_x) > 1 && length(new_x[1]) == 1 && length(rad.lb)>1)
        push!(rad.x,new_x)
        push!(rad.y,new_y)
    else
        append!(rad.x,new_x)
        append!(rad.y,new_y)
    end
    rad.coeff = _calc_coeffs(rad.x,rad.y,rad.lb,rad.ub,rad.phi,rad.dim_poly, rad.scale_factor)
    nothing
end


text:
### TAKE TWO

using Primes: primes, factor, isprime

### Helpers

nroot(x::Number, n::Number) = x^(1/n)

function factors(n::Integer)
    f = [one(n)]

    for (p,e) in factor(n)
        f = reduce(vcat, [f*p^j for j in 1:e], init=f)
    end

    return ifelse(isone(length(f)), [one(n), n], sort!(f))
end

"""
    ord(q::Integer, n::Integer) -> Integer

Computes the multaplicative order, ord_q n, of a number n mod q.
"""
function ord(q::Integer, n::Integer)
    isone(gcd(n, q)) || throw(error("$n and $q are not coprime"))
    res = one(q)

    for (p, e) in factor(q)
        m = p^e
        t = div(q, p) * (p-1)
        for f in factors(t)
            if isone(powermod(n, f, q))
                res = lcm(res, f)
                break
            end
        end
    end

    return res
end

### Main code
#=
\begin{enumerate} \item For each prime power $q$ such that $2^q\leq n$, write down a positive integer $r_q$ such that if $n$ is a $q$\th power, then $n=r^q_q$.\item Find a finite coprime set $P$ of integers larger than 1 such that each of $n,r_2,r_3,r_4, r_5,r_7,\ldots$ is a product of powers of elements of $P$.  (In this paper, ``coprime'' means ``pairwise coprime.'')\item Factor $n$ as $\Pi_{p\in P}p^{n_p}$, and compute $k=\text{gcd}n_p:p\in P$. \end{enumerate}=#

#=
From \texttt{Powers2.pdf}:\begin{quotation}
	Define $r_q$ as an integer 2-adically close to $n^{1/q}$, as explained in \texttt{Powers.pdf}, section 21.\end{quotation}

From \texttt{Powers.pdf}, section 21:\begin{quotation}
	In this section, I describe a \emph{2-adic variant} of Algorithm X.  With this variant, I can work with integers rather than floating point numbers; I no longer need to guard bits; I can jump directly into Newton's method without a preliminary binary search; and a proper error analysis takes a few lines rather than several pages.\par
	It will be convenient to restrict attention to odd $n$.  See section 22 to handle even $n$.\par
	\textbf{Motivation.}  To check if $n$ is a $k$\th power, I compute a tentative $k$\th root of $n$\emdash an integer x such that no integer other than $x$ can possibly be the $k$\th root of $n$.  Then I test whether $x^k=n$.\par
	To find $x$ in Part II and Part III, I constructed a number that was \emph{close} to a $k$\th root of $n$ in the usual metric.  I used the same metric again to check whether $x^k=n$: I computed $x^k$ in low precision to see whether it was close to $n$.\par
	Nothing in the original problem suggests this metric.  The 2-adic variant uses a different metric, where $i$ and $j$ are close if $i-j$ is divisible by a high power of 2.\par
	\textbf{Notation.}  This section deviates from the notation of parts I, II, and III: $r$, $y$, and $z$ are odd integers rather than positive floating-point numbers.\par
	\textbf{Lemma 21.1.}  \emph{If $2i\equiv 2j \mod 2^{b+1}$, and $b\geq 1$, then $i^2\equiv j^2\mod w^{b+1}.$}\par
	\textbf{2-adoc approximate powers.}  Fix positive integers $k$ and $b$.  For any integer $m$ define $\text{pow}_{2,b}(m,k)=m^k\mod 2^b$.  See section 6 for methods of computing $\text{pow}_{2, b}$ without many multiplications.  As I compute $\text{pow}_{2,b}(m, k)$, I keep track of an ``overflow bit'' to figure out whether $m^k\mod 2^b=m^k$.
\end{quotation}

From \texttt{Powers.pdf}, section 6:\begin{quotation}
	Let $r$ be a positive floating-point number, and let $k$ be positive integers.  Then $\text{pow}_{b}(r, k)$, the \textbf{$\mathbf{b}$-bit approximate $\mathbf{k}th$ power of $\mathbf{r}$} is a floating-point approximation to $r^k$.  In this section, I show how to compute $\text{pow}_b(r,k)$ in $M$-time at most $P(k)M(b), where $P(k)\leq 2\lfloor\text{lg} k\rfloor$.\par
	Define $P(k)$ for $k\geq 1$ as follows: $P(1)=0;P(2k)=P(k)+1;P(2k+1)=P(2k)+1$.\par
	\textbf{Lemma 6.1.}  $P(k)\leq 2\lfloor\text{lg} k\rfloor$.\par
	\emph{Proof.} $\ldots$\par
	\textbf{Lemma 6.2} $\text{pow}_b(r,k)\leq r^k<\text{pow}_b(r,k)(1+2^{1-b})^{2k-1}.$\par
	\emph{Proof.} $\ldots$\par
	\emph{Algorithm P.} $\ldots$\par
	\textbf{Lemma 6.3.}  Algorithm P computes $\text{pow}_b(r, k)$ in $M$-time at most $P(k)M(b)$.\par
	\emph{Proof.} $\ldots$\par
\end{quotation}

From \texttt{Powers.pdf}, section 22:\begin{quotation}

\end{quotation}
=#

lg(x) = log2(x)

#=
From Knuths The Art of Computer Programming, semi-numerical algorithms, 2nd edition, 4.3.1, exercise 16.

Design an algorithm that divides a nonnegative n-place integer $(u_1u_2\ldots u_n)_b$ by $v$, where $v$ is a single-precision number (i.e., $0 < v < b$), producing the quotient $(w_1w_2\ldots 2_n)_b$ and remainder $r$.
=#

"""
	divrem(b::Integer, u::Integer, v::Real)

This algorithm divides a non-negative n-place integer (u1 u2 ... un) by v, where v is a single-precision number (i.e., 0 < v < b), producing the quotient (w1 w2 ... wn) and remainder r.

This algorithm was written by Donald Knuth, in The Art of Computer Programming, semi-numerical algorithms, 2nd edition, exercise 4.3.1, 16.
"""
function Base.divrem(b::Integer, u::Integer, v::Real)
	r, j = 0, 1
	n = ndigits(u) #; base = b
	u̲ = round.(reverse(digits(u)), base = b)
	w̲ = Vector{Int}(undef, n)

	while j ≤ n
		w̲[j] = floor(Int, (r*b + u̲[j]) / v)
		r = mod(r*b + u̲[j], v)
		j += 1
	end

	# w = sum(w̲[i] * 10 ^ (length(w̲) - i) for i in eachindex(w̲))
	w = parse(Int, join(string.(w̲)))
	# w = foldr((i, j) -> 10 * j + i, reverse(w̲))
	#
	# n = 0
	# reduce(w̲) do i, j
	# 	10 * j + i
	# 	# global n = n + ndigits(w̲[i])
	# end
	# n = 0
	# w = sum(eachindex(w̲)) do i
	# 	w̲[i] * 10 ^ (length(w̲) - i + n)
	# 	global n = n + ndigits(w̲[i])
	# end
	# n = 0

	return w, r
end

"""
	div(b::Integer, u::Integer, v::Real)

This algorithm divides a non-negative n-place integer (u1 u2 ... un) by v, where v is a single-precision number (i.e., 0 < v < b), producing the quotient (w1 w2 ... wn) and ommiting the remained.  To return the quotient and the remainder, see `divrem`.
"""
Base.div(b::Integer, u::Integer, v::Real) = first(divrem(b, u, v))

"""
	rem(b::Integer, u::Integer, v::Real)

This algorithm divides a non-negative n-place integer (u1 u2 ... un) by v, where v is a single-precision number (i.e., 0 < v < b), producing a remainder r.  To return the quotient and the remainder, see `divrem`.
"""
Base.rem(b::Integer, u::Integer, v::Real) = last(divrem(b, u, v))

#=
In this section, I define truncation to $b$ bits, written $\text{trunc}_b$, and show that $\frac{r}{\text{trunc}_br}$ is between $1$ and $1+2^{1-b}$.  More generally, for any positive integer $k$, I define $\text{div}_b(r,k)$ as a floating point approximation to $\frac{r}{k}$, so that $\frac{r}{k\text{div}_b(r,k)}$ is between $1$ and $1+2^{1-b}$.\\ Fix $b\geq 1$.  Set $\text{div}_b(a,n,k)=(a+f-\left\lceil\text{lg} k\right\rceil - b, \left\lfloor\frac{n}{2^{f-\left\lceil\text{lg} k\right\rceil-b}k}\right\rfloor)$, where $2^{f-1}\leq n<w^f$.  (Note that $f-\left\lceil \text{lg} k\right\rceil - b$ may be negative.)  This map induces a map, also denoted $\text{div}_b$, upon positive floating-point numbers: \begin{equation}\text{div}_b(2^an,k)=2^{a+f-\left\lceil\text{lg} k\right\rceil - b}\left\lfloor\frac{n}{2^{f-\left\lceil\text{lg}k\right\rceil - b}k}\right\rfloor\quad\quad\quad\text{if }2^{f-1}\leq n<2^f.\end{equation}\\ To compute $\text{div}_b(r,k)$ I also use an algorithm designed for dividing by small integers; time spent computing $\text{div}_b$ is not $M$-time.
=#

"""
	trunc(b::Integer, n::Real)

Truncates a number `n` to `b` bits.
"""
function Base.trunc(b::Integer, n::Real)
	return n & (1 << b - 1)
	# return div(b, n, 1)
end

#=
Given a positive floating-point number r and two positive integers b, k, to prime pow_b(r, k).  Let $\text{pow}_b(r,k)$ be the $b$-bit approximate $k$th power of $r$.
\begin{equation}\begin{aligned}
\text{pow}_b(r,1) &= \text{trunc}_br; \\
\text{pow}_b(r,2k) &= \text{trunc}_b(\text{pow}_b(r,k)^2); \\
\text{pow}_b(r,2k+1) &= \text{trunc}_b(\text{pow}_b(r,2k)\text{pow}_b(r,1))
\end{aligned}\end{equation}
=#

"""
	pow(b::Integer, r::Real, k::Integer)

Given a positive floating-point number r and two positive integers b and k, pow(b, r, k) computes the b-bit approximate kth power of r.  pow_b(r, k) ≤ r^k < pow_b(r, k)(1 + 2^(1 - b))^(2k - 1).
"""
function pow(b::Integer, r::Real, k::Integer)
	if r ≤ 0 || b ≤ 0 || k ≤ 0
		throw(error("r, b, and k must be positive numbers to compute pow_b(r, k)."))
	end
	isone(k) && return trunc(b, r)
	iszero(mod(k, 2)) && return trunc(b, pow(b, r, k ÷ 2)^2) # use \div because k is even
	return trunc(b, pow(b, r, k - 1) * trunc(b, r))
end

"""
	pow2(b::Integer, r::Real, k::Integer)

Given a positive floating-point number r and two positive integers b and k, pow2(b, r, k) computes the 2-adic b-bit approximate kth power of r.
"""
function pow2(b::Integer, m::Integer, k::Integer)
	if k ≤ 0 || b ≤ 0
		throw(error("k, and b must be positive integers to find a 2-adic approximate power."))
	end

	return mod(m^k, 2^b)
end

"""
	tentative_check_kth_roots(n::Integer, x::Integer, k::Integer)

A straight forward algorithm for checking whether x^k = n.
"""
function tentative_check_kth_roots(n::Integer, x::Integer, k::Integer)
	if n ≤ 0 || x ≤ 0 || k ≤ 0
		throw(error("n, x, and k must be positive integers to compute check whether n=x^k."))
	end

	f = floor(Integer, lg(2*n))

	isone(x) && return ifelse(isone(n), true, false) # return ifelse(isone(n), 0, 2)

	b = 1
	while true
		# r = pow(x, b, k) # 2, b
		r = pow2(b, x, k)
		isequal(mod(n, 2^b), r) || return false # 2
		b ≥ f && return ifelse(isequal(r, x^k), true, false) # ifelse(isequal(r, x^k), 0, 2)

		b = min(2*b, f)
	end
end

#=
\textbf{Lemma 21.2}.  Algorithm C2 prints 0 if and only if n=x^k.\par
\emph{Proof.} $\ldots$\par
\textbf{2-adic approximate multiplication and division}.  Fix $b\geq 1$.  For $m$ an integer and $k$ a positive integer write $\text{mul}_{2, b}(m, k) = km \mod 2^b.$\par
If $k$ is odd, write $\text{div}_{2, b}(m, k)$ for the unique integer between 0 inclusive and $2^b$ exclusive such that $m\equiv k \text{div}_{2, b}(m, k) \mod 2^b$.
=#

#=
Write a MIX program that multiplies $(u_1u_2\ldots u_n)_b$ by $v$, where $v$ is a single-precision number (i.e., $0 \leq v < b$), producing the answer $(w_0w_1\ldots w_n)_b$.  How much running time is required?

		ENTX	N		1	ADD		CARRY	N
		JOV		OFLO	1	JNOV	*+2		N
		ENTX	0		1	INCX	1		K
	2H	STX		CARRY	N	STA		W,1		N
		LDA		U,1		N	DEC1	1		N
		MUL		V		N	J1P		2B		N
		SLC		5		N	STX		W		1

The running time is 23N + K + 5 cycles, and K is rougly (1/2)N.
=#

"""
	mul2(b::Integer, m::Integer, k::Integer)

The notation mul(r, k) is used in the Bernstein paper to denote k * r.  The notation mul_{2, b}(m, k) represents 2-adic multiplication.
"""
function mul2(b::Integer, m::Integer, k::Integer)
	if k ≤ 0
		throw(error("Cannot compute mul(m, k) when k is negative using the 2-adic method."))
	end

	return mod(k*m, 2^b)
end

"""
	div2(b::Integer, m::Integer, k::Integer)
div_b is defined as a floating-point approximation to r/k, so that r/k div_b(r, k) is between 1 and 1 + 2^(1 - b).  div2 is the 2-adic approximation of this.
"""
function div2(b::Integer, m::Integer, k::Integer)
	if k ≤ 0
		throw(error("Cannot compute mul(m, k) when k is negative using the 2-adic method."))
	end

	for m in 0:2^b
		if isequal(m, mod(div(b, m, k), 2^b))
			return div(b, m, k)
		end
	end
end

"""
	odd_nroot2(b::Integer, y::Integer, k::Integer)

Fix an odd integer y and a positive odd integer k.  odd_nroot2 finds an approximate negative kth root of y by Newton's method.
"""
function odd_nroot2(b::Integer, y::Integer, k::Integer)
	if k ≤ 0 || b ≤ 0
		throw(error("Cannot use this method for negative number k and b."))
	end

	if iseven(y) || iseven(k)
		throw(error("Cannot use this method for even y and k."))
	end

	b′ = ceil(Integer, b / 2)

	isone(b) && return 1
	z = odd_nroot2(b′, y, k)
	r₂ = mul2(b, z, k + 1)
	r₃ = mod(y * pow2(b, z, k + 1), 2^b)
	return r₄ = div2(b, r₂ - r₃, k)
end

#=
\textbf{Lemma 21.3.}  If $k$ is odd and $r=\text{nroot}_{2,b}(y, k)$, then $r^k y\mod 2^b = 1$.\par
\emph{Proof.} $\ldots$\par
=#

"""
	sqrt_nroot2(b::Integer, y::Integer)

Fix odd integer y.  For each b ≥ 1, define and construct nroot_(2, b)(y, 2).
"""
function sqrt_nroot2(b::Integer, y::Integer)
	# if b ≤ 0
	# 	throw(error("Cannot use this method for negative number b."))
	# end
	#
	# if iseven(y)
	# 	throw(error("Cannot use this method for even y."))
	# end

	b′ = ceil(Integer, (b + 1) / 2)

	isone(b) && return ifelse(isone(mod(y, 4)), 1, 0)
	isequal(b, 2) && return ifelse(isone(mod(y, 8)), 1, 0)
	iszero(z) && return 0

	r₂ = mul2(b + 1, z, 3)
	r₃ = mod(y * pow2(b + 1, z, 3), 2^(b + 1))
	return r₄ = mod(r₂ - r₃, 2^b)
end

#=
\textbf{Lemma 21.4.}  Set $r = \text{nroot}_{2, b}(y, 2)$.  If $i^2y\mod 2^{b + 1}=1$ for some odd integer $i$, then $r\neq 0$.  If $r\neq 0$ then $r^2y\mod 2^{b+1} = 1$.\par
\emph{Proof.} $\ldots$\par
=#

"""
	perfect_power_decomp_odd(n::Integer, k::Integer, y::Integer)

Given a positive odd integer n, an integer k ≥ 2 such that either k = 2 or k is odd, and an odd integer y, `perfect_power_decomp_odd` checks if n is a kth power.
"""
function perfect_power_decomp_odd(n::Integer, k::Integer, y::Integer)
	if iseven(n) && !isequal(k, 2)
		throw(error("Cannot use this algorithm on even n."))
	end

	if k < 2
		throw(error("k must be ≥ 2."))
	else
		if iseven(k) && !isequal(k, 2)
			throw(error("k must be 2 or odd."))
		end
	end

	f = floor(Integer, lg(2*n))
	b = ceil(Integer, f / k)

	r = odd_nroot2(b, y, k)
	(isequal(k, 2) && iszero(r)) && return 0
	tentative_check_kth_roots(n, r, k) && return r
	(isequal(k, 2) && tentative_check_kth_roots(n, 2^b - r, k)) && return 2^b - r
	return 0
end

#=
\textbf{Lemma 21.5.}  Set $f=\lfloor\text{lg}2n\rfloor$ and $b=\lceilf/k\rceil$.  Assume that $yn\mod 2^{b+1}=1$.  If $n$ is a $k$\th power, Algorithm K2 (above) prints $n^{1/k}$.  If $n$ is not a $k$\th power, Algorithm K2 prints 0.\par
\emph{Proof.} $\ldots$\par
=#

"""
	odd_decomp_perf_pow(n::Integer)

Given an odd integer n ≥ 2, `odd_decomp_perf_pow` decomposes n as a perfect power if possible.
"""
function perfectpower(n::Integer)
    # isprime(n) && return n, 1

	# if iseven(n) && !isequal(n, 2)
	# 	throw(error("This method must have odd n."))
	# end

	if n < 2
		throw(error("n must be ≥ 2."))
	end

    isprime(n) && return n, 1

    # sqrt_val = sqrt(n)
    # isinteger(sqrt_val) && return round(Integer, sqrt_val), round(Integer, sqrt_val)

    if iseven(n)
        if ispow2(n)
            return 2, round(Integer, log2(n))
        else
            return n, 1
        end
    end

	f = floor(Integer, lg(2*n))

	y = odd_nroot2(ceil(Integer, f / 2) + 1, n, 1)

	for p in primes(f - 1)
		x = perfect_power_decomp_odd(n, p, y)
		x > 0 && return x, p
	end

	return n, 1
end

#=
\textbf{Lemma 21.6.}  If $n$ is a perfect power, the above algorithm (Algorithm X2) prints a prime number $p$ and a positive integer $x$ such that $x^p=n$.  If $n$ is not a perfect power, Algorithm X2 print (n, 1).
\emph{Proof.} $\ldots$\par
=#

#=
As usual, fix $n\geq 2$.  In this section I discuss several tricks based on computing $n$ mod $q$ for on or more primes $q$.\par
\textbf{If $n$ has no small prime divisors, lower the exponent bound.}  If $n$ is odd and $n=x^k$, then $x$ is also odd.  So $x\geq 3$ and $k\leq \log_3 n$.  More generally, one may compute $n\mod q$ for all primes $p < T$, for some bound $T$.  If $n\mod q$ is always nonzero, one need not check exponents past $\log_t n$.\par
\textbf{If n has a prime divisor, find its order.}  What if $n\mmod 1 = 0$?  First compute the number $\text{ord}_q n$ of factors $q$ in $n$, together with $n/q^{\text{ord}_qn}$.  Then check, for each $p$ dividing $\text{ord}_qn$ whether $n/q^{\text{ord}_qn}$ is a $p$\th power.  Otherwise, $n$ cannot be a perfect power.  (Note that $n/q^{\text{ord}_qn}$ may be 1, in which case no trsting is necessary.)\par Recall that the 2-adic method in section 21 requires $n$ to be odd.  This is not a serious restriction.  If $n$ is even, the method here ends up checking whether $n/2^{\text{ord}_2n}$ is a $p$th power, for various primes $p$; and $n/2^{\text{ord}_2n}$ is odd.\par There are several plausible ways to compute the number of factors $q$ in $n$.\par If $q=2$ then $\text{ord}_qn$ is the number of 0 buts at the bottom of $n$'s binary expansion.  \par If $q>2$, I do a binary search upon the number of factors.  The idea is to compute $n\mod q^c$ and $\lfloor n/q^c\rfloor$ for some integers $c\approx (\log_q n)/2$.  If $n\mod q^c\neq 0$ then $\text{ord}_qn=\text{ord}_q(n\mod q^c)$; if $n\mod q^c=0$ then $\text{ord}_qn=c+\text{ord}_q(n/q^c)$.  Chop $c$ in half and repeat.  This method takes essentially linear time given fast multiplication and the algorithms from section 8.\par I could instead do a linear search; this amounts to always taking $c=1$ in the above description.  This will be faster than a binary search on average.  I could compromise with a sequence of $c$ that is at first optimistic but backs off quickly if necessary.  For example: begin with $c=1$, double $c$ if $n\mod q^c = 0$, and chop $c$ in half if $n\mod q^c\neq 0$.\par
\textbf{Check the character of residues of n}.  If $n$ is a $k$\th power, and $q$ is a prime with $q\mod k=1$, then $n^{(q - 1) / k}\mod q$ is either 0 or 1.  A non-$k$\th power has roughly a $1/k$ chance of passing this test.\par
\textbf{Check the residues of tentative roots.}  If $n=x^k$ then $n\mod q=x^k\mod q$.  So, given $n\mod q$, I can try to weed out a tentative root $x$ by calculating the $k\th$ power modulo $q$ of $x\mod q$.  In practice this test is quite powerful: fi $n\neq x^k$ then very few primes $q$ divide $n-x^k$.\par In this test $q$ need not be prime.  It might be convenient to check whether $n$ agrees with $x^k$ modulo $2^{32}$, for example, though this is redundant if $x$ was constructed by 2-adic methods.\par One could develop a fast randomized power-testing algorithm along these lines.  Start from a tentative root $x$.  First check if $x^k\leq 2n$.  Then check if $n\mod q$ equals $x^k\mod q$ for a set of ``random'' primes $q$ with product larger than $n$.  This test will succeed if and only if $n=x^k$.  If $n\neq x^k$ then one will, on average, test very few $q$'s.\par
\textbf{Check for small divisors of tentative roots}.  If $n$ is not divisible by any primes $q<T$, and $n=x^k$, then $x$ is not divisible by any primes $q<T$.  So I can throw away any tentative root $x$ that has prime factors smaller than $T$.  This is much weaker than testing whether $n\mod d=x^k\mod q$ for each $q<T$, but it is also much faster.
=#

function lower_exp_bound(n::Integer, x::Integer, k::Integer)
    # if isodd(n)
        # continue
    # end
end


text:
function init(init_state, energy)
    init_energy = energy(init_state)
    init_trial = Trial(init_state, init_energy)
    return init_trial, init_trial
end

function anneal_one(now_trial, move, energy)
    prev_trial = deepcopy(now_trial)
    now_state = move(now_trial.state)
    now_energy = energy(now_state)
    now_trial = Trial(now_state, now_energy)
    return now_trial, prev_trial
end

function check_state_update(now_trial, prev_trial, t)
    # diff_energy > 0 は劣化
    diff_energy = prev_trial.energy - now_trial.energy
    uniform = Distributions.Uniform(0, 1)
    if diff_energy > 0.0 && exp(- diff_energy / t) < rand(uniform)
        return false
    else
        return true
    end
end

function update(now_trial, prev_trial, best_trial)
    if best_trial.energy <= now_trial.energy
        best_trial = deepcopy(now_trial)
    end
    return now_trial, best_trial
end

function move(state)
    # need original function
    return  state
end

function energy(state)
    # need original function
    return  sum(state)
end

max_stars_repo_name: ariririri/heuristicopt

max_issues_repo_name: ariririri/heuristicopt


text:
function minimize(objective::ObjWrapper, x0,
                   scheme::QuasiNewton,
                   options=OptOptions())

    minimize(objective, (x0, nothing), (scheme, Backtracking()), options)
end
function minimize(objective::ObjWrapper, x0,
                   approach::Tuple, options=OptOptions())

    minimize(objective, (x0, nothing), approach, options)
end
function minimize(objective::ObjWrapper, state0::Tuple,
                  scheme::QuasiNewton,
                  options=OptOptions())
    minimize(objective, state0, (scheme, Backtracking()), options)
end
function minimize(objective::T1, state0::Tuple, approach::Tuple{<:Any, <:LineSearch},
                  options::OptOptions=OptOptions(),
                  linesearch::T2 = Backtracking()
                  ) where {T1<:ObjWrapper, T2}

    x0, B0 = state0
    T = eltype(x0)
    x, fx, ∇fx, z, fz, ∇fz, B = prepare_variables(objective, approach, x0, copy(x0), B0)
    # first iteration
    x, z, fz, ∇fz, B, is_converged = iterate(x, fx, ∇fx, z, fz, ∇fz, B, approach, objective, options)

    iter = 0
    while iter <= options.maxiter && !is_converged
        iter += 1

        # take a step and update approximation
        x, z, fz, ∇fz, B, is_converged = iterate(x, fx, ∇fx, z, fz, ∇fz, B, approach, objective, options, false)
        if norm(x.-z, Inf) ≤ T(1e-16)
            break
        end
    end

    return z, fz, ∇fz, iter
end

function iterate(x, fx::Tf, ∇fx, z, fz, ∇fz, B, approach, objective, options, is_first=nothing) where Tf
    # split up the approach into the hessian approximation scheme and line search
    scheme, linesearch = approach

    # Move nexts into currs
    fx = fz
    x = copy(z)
    ∇fx = copy(∇fz)

    # Update current gradient and calculate the search direction
    d = find_direction(B, ∇fx, scheme) # solve Bd = -∇fx
    φ = LineObjective(objective, x, d, fx, dot(∇fx, d))

    # # Perform line search along d
    α, f_α, ls_success = find_steplength(linesearch, φ, Tf(1))

    # # Calculate final step vector and update the state
    s = @. α * d
    z = @. x + s

    # Update approximation
    fz, ∇fz, B = update_obj(objective, d, s, ∇fx, z, ∇fz, B, scheme, is_first)

    # Check for gradient convergence
    is_converged = converged(z, ∇fz, options.g_tol)

    return x, z, fz, ∇fz, B, is_converged
end


text:
module ContactLQR

using RigidBodyDynamics
using ForwardDiff

"""
`care(A, B, Q, R)`

Compute 'X', the solution to the continuous-time algebraic Riccati equation,
defined as A'X + XA - (XB)R^-1(B'X) + Q = 0, where R is non-singular.

Algorithm taken from:
Laub, "A Schur Method for Solving Algebraic Riccati Equations."
http://dspace.mit.edu/bitstream/handle/1721.1/1301/R-0859-05666488.pdf

Implementation from https://github.com/JuliaControl/ControlSystems.jl/blob/master/src/matrix_comps.jl
by Jim Crist and other contributors.
"""
function care(A, B, Q, R)
    G = try
        B*inv(R)*B'
    catch
        error("R must be non-singular.")
    end

    Z = [A  -G;
        -Q  -A']

    S = schurfact(Z)
    S = ordschur(S, real(S.values).<0)
    U = S.Z

    (m, n) = size(U)
    U11 = U[1:div(m, 2), 1:div(n,2)]
    U21 = U[div(m,2)+1:m, 1:div(n,2)]
    return U21/U11
end

"""
`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
A = [0 1; 0 0]
B = [0;1]
C = [1 0]
sys = ss(A,B,C,0)
Q = eye(2)
R = eye(1)
L = lqr(sys,Q,R)

u(t,x) = -L*x # Form control law,
t=0:0.1:5
x0 = [1,0]
y, t, x, uout = lsim(sys,u,t,x0)
plot(t,x, lab=["Position", "Velocity"]', xlabel="Time [s]")
```

Implementation from https://github.com/JuliaControl/ControlSystems.jl/blob/master/src/synthesis.jl
by Jim Crist and other contributors.
"""
function lqr(A, B, Q, R)
    S = care(A, B, Q, R)
    K = R\B'*S
    return K, S
end

"""
    `dare(A, B, Q, R)`

Compute `X`, the solution to the discrete-time algebraic Riccati equation,
defined as A'XA - X - (A'XB)(B'XB + R)^-1(B'XA) + Q = 0, where A and R
are non-singular.

Algorithm taken from:
Laub, "A Schur Method for Solving Algebraic Riccati Equations."
http://dspace.mit.edu/bitstream/handle/1721.1/1301/R-0859-05666488.pdf

Implementation from https://github.com/JuliaControl/ControlSystems.jl/blob/master/src/matrix_comps.jl
by Jim Crist and other contributors.
"""
function dare(A, B, Q, R)
    G = try
        B*inv(R)*B'
    catch
        error("R must be non-singular.")
    end

    Ait = try
        inv(A)'
    catch
        error("A must be non-singular.")
    end

    Z = [A + G*Ait*Q   -G*Ait;
         -Ait*Q        Ait]

    S = schurfact(Z)
    S = ordschur(S, abs.(S.values).<=1)
    U = S.Z

    (m, n) = size(U)
    U11 = U[1:div(m, 2), 1:div(n,2)]
    U21 = U[div(m,2)+1:m, 1:div(n,2)]
    return U21/U11
end

"""
    `dlqr(A, B, 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 discrete time model `x[k+1] = Ax[k] + Bu[k]`.

Implementation from https://github.com/JuliaControl/ControlSystems.jl/blob/master/src/synthesis.jl
by Jim Crist and other contributors.
```
"""
function dlqr(A, B, Q, R)
    S = dare(A, B, Q, R)
    K = (B'*S*B + R)\(B'S*A)
    return K, S
end

function contact_jacobian(state::MechanismState{T}, contacts) where T
    mechanism = state.mechanism
    if isempty(contacts)
        return zeros(0, num_velocities(state))
    end
    world = root_body(mechanism)
    contact_jacobians = map(contacts) do contact
        body = RigidBodyDynamics.body_fixed_frame_to_body(mechanism, contact.frame)
        p = transform(state, contact, default_frame(world))
        Jp = point_jacobian(state, path(mechanism, root_body(mechanism), body), p)
        @assert Jp.frame == default_frame(world)
        Array(Jp)
    end
    Jc = vcat(contact_jacobians...)
    Jc = Jc[[i for i in 1:size(Jc, 1) if !(norm(Jc[i, :]) < 1e-9)], :]
    Jc
end

function dynamics_with_contact_constraint(state::MechanismState, input::AbstractVector, Jc::AbstractMatrix)
    v = velocity(state)
    M = full(mass_matrix(state))
    C_plus_g = dynamics_bias(state)
    St = eye(num_velocities(state))
    for joint in joints(state.mechanism)
        bounds = effort_bounds(joint)
        for (i, b) in enumerate(bounds)
            if b.upper == b.lower == 0
                j = velocity_range(state, joint)[i]
                St[j, j] = 0
            end
        end
    end
    Jct = Jc'
    Jcbar = inv(M) * Jct * inv(Jc * inv(M) * Jct)
    Φ = (I - Jcbar * Jc) * inv(M)
    ϕ = -Φ * (C_plus_g)

    v̇ = Φ * St * input + ϕ
    q̇ = configuration_derivative(state)
    vcat(q̇, v̇)
end

"""
From "Full Dynamics LQR Control of a Humanoid Robot: An Experimental Study on
Balancing and Squatting" by Sean Mason et al.
"""
function contact_linearize(state0, input0, Jc)
    if norm(velocity(state0)) > 0
        error("Only static postures supported")
    end
    mechanism = state0.mechanism

    function dynamics(x, u)
        q = x[1:num_positions(state0)]
        v = x[num_positions(state0) + (1:num_velocities(state0))]
        state = MechanismState(mechanism, q, v)
        dynamics_with_contact_constraint(state, u, Jc)
    end

    A = ForwardDiff.jacobian(vcat(configuration(state0), velocity(state0))) do x
        dynamics(x, input0)
    end
    B = ForwardDiff.jacobian(input0) do u
        dynamics(Vector(state0), u)
    end
    A, B, dynamics(Vector(state0), input0)
end

"""
    A_d, B_d, c_d = zero_order_hold(A, B, c, Δt)

Discretize the continuous-time linear system:

    ẋ = Ax + Bu + c

into the discrete-time linear system:

    x[i+1] = Ax[i] + Bu[i] + c

using a zero-order (piecewise constant) hold on u(t) with discretization
interval Δt.
"""
function zero_order_hold(A, B, c, Δt)
    nx = size(A, 1)
    nu = size(B, 2)
    matc = zeros(nx + nu + 1, nx + nu + 1)
    matc[1:nx, :] = hcat(A, B, c)
    matd = expm(matc * Δt)
    A_d = matd[1:nx, 1:nx]
    B_d = matd[1:nx, nx + (1:nu)]
    c_d = matd[1:nx, nx + nu + 1]
    A_d, B_d, c_d
end

"""
From "Balancing and Walking Using Full Dynamics LQR Control With Contact
Constraints" by Sean Mason et al.
"""
function contact_lqr(state::MechanismState, input::AbstractVector, Q::AbstractMatrix, R::AbstractMatrix, contacts::AbstractVector{<:Point3D})
    Jc = contact_jacobian(state, contacts)
    A, B, c = contact_linearize(state, input, Jc)
    A[abs.(A) .< 1e-6] .= 0
    N = nullspace([Jc zeros(Jc); zeros(Jc) Jc])
    Am = N' * A * N
    Bm = N' * B
    Rm = R
    Qm = N' * Q * N
    Km, Sm = lqr(Am, Bm, Qm, Rm)
    K = Km * N'
    S = N * Sm * N'
    return A, B, c, K, S
end

function contact_dlqr(state::MechanismState, input::AbstractVector, Q::AbstractMatrix, R::AbstractMatrix, Δt, contacts::AbstractVector{<:Point3D}=[])
    Jc = contact_jacobian(state, contacts)
    A, B, c = contact_linearize(state, input, Jc)
    A[abs.(A) .< 1e-6] .= 0
    J_qdot_to_v = RigidBodyDynamics.configuration_derivative_to_velocity_jacobian(state)
    Jcv = Jc
    Jcq = Jc * J_qdot_to_v
    M = [Jcq                              zeros(size(Jcq, 1), size(Jc, 2));
         zeros(size(Jc, 1), size(Jcq, 2)) Jc]
    N = nullspace(M)
    Am = N' * A * N
    Bm = N' * B
    Rm = R
    Qm = N' * Q * N
    Am_d, Bm_d, cm_d = zero_order_hold(Am, Bm, zeros(size(Am, 1)), Δt)
    Km_d, Sm_d = dlqr(Am_d, Bm_d, Qm, Rm)
    K_d = Km_d * N'
    S_d = N * Sm_d * N'
    return zero_order_hold(A, B, c, Δt)..., K_d, S_d
end

end


text:
using Documenter
using DemoCards
using JSON

root = joinpath(@__DIR__, "..")

# these are tracked by git
raw_problems_dir = joinpath(root, "src", "problems")
raw_theme_dir = joinpath(root, "docs", "leetcode_theme")

# where the problems files are generated -- this won't be tracked by git
page_root = abspath(root, "docs", "problems")
page_problems = abspath(page_root, "problems")

if isdir(page_root)
    rm(page_root; force=true, recursive=true)
end
mkpath(page_problems)

# 1. copy the whole problems file dir
cp(raw_problems_dir, page_problems; force=true)
rm(joinpath(page_problems, "problems.jl")) # no need to render this

# 2. provide a configuration file to tell DemoCards that we want to use number order instead of string order

open(joinpath(page_problems, "config.json"), "w") do io
    files = filter(x->x!="config.json", readdir(page_problems))
    sort!(files; by=x->parse(Int, split(x, ".")[1]))
    JSON.print(io, Dict("title"=>"Problems", "order"=>files))
end

## build docs

# 1. generate demo files
demopage, postprocess_cb = makedemos("problems")

# 2. normal Documenter usage
format = Documenter.HTML(prettyurls=get(ENV, "CI", nothing) == "true")
makedocs(format = format,
         pages = [
            "Home" => "index.md",
            demopage,
         ],
         sitename = "LeetCode")

# 3. hide the code
js_dir = abspath(root, "docs", "build", "assets")
css_dir = abspath(js_dir, "themes")

# 3.1 modify the css file to seting the style of code
for css_file in readdir(css_dir; join=true)
    # change the permissions mode of files
    chmod(css_file, 0o777)

    # add code in the end of line
    css_data = read(joinpath(raw_theme_dir, "leetcode.css"))
    open(css_file, "a") do io
        write(io, css_data)
    end
end

# 3.2 modify the documenter.js to setting the mouse actions
documenter_file = joinpath(js_dir, "documenter.js")

# add code in the end of line
js_data = read(joinpath(raw_theme_dir, "leetcode.js"))
open(documenter_file, "a") do io
    write(io, js_data)
end

# 4. postprocess after makedocs
postprocess_cb()

deploydocs(repo = "github.com/JuliaCN/LeetCode.jl.git")


text:
module IntegrationTests

using Test

@testset "Quadratures" begin include("QuadraturesTests.jl") end

@testset "DuffyQuadratures" begin include("DuffyQuadraturesTests.jl") end

@testset "QuadraturesFactories" begin include("QuadratureFactoriesTests.jl") end

@testset "CellQuadratures" begin include("CellQuadraturesTests.jl") end

@testset "Triangulations" begin include("TriangulationsTests.jl") end

@testset "CellIntegration" begin include("CellIntegrationTests.jl") end

@testset "BoundaryTriangulations" begin include("BoundaryTriangulationsTests.jl") end

@testset "NormalVectors" begin include("NormalVectorsTests.jl") end

@testset "SkeletonTriangulations" begin include("SkeletonTriangulationsTests.jl") end

@testset "BoundaryCellFields" begin include("BoundaryCellFieldsTests.jl") end

@testset "SkeletonCellFields" begin include("SkeletonCellFieldsTests.jl") end

end # module


text:
export assemble, symboltable

"Thrown if assembler doesn't know the mnemonic parsed"
struct InvalidMnemonicError <: Exception
    mnemonic::String
end

"Remove comment from source code line"
remove_comment(line) = split(line, "//") |> first |> strip

"""
    symboltable(filename::AbstractString)
Parse `filename` and find all labels used in assembly program.
"""
function symboltable(filename::AbstractString)
    lines = readlines(filename)
    symboltable(lines)
end

"""
    symboltable(lines::AbstractArray{<:AbstractString})
Record address of labels used in code in a symbol table.
Return this as a dictionary.

    symtable = symboltable("foobar.lmc")

    # Get address of LOOP label
    address = symtable["LOOP"]
"""
function symboltable(lines::AbstractArray{<:AbstractString})
    labels = Dict{String, Int}()
    possible_used_labels = Set{String}()

    address = 0

    for line in lines
        codeline = remove_comment(line)
        words = split(codeline)

        if isempty(words) continue end

        if length(words) > 1
            if haskey(opcodes, words[2]) || words[2] == "DAT"
                label = words[1]
                labels[label] = address
            end
            if isletter(words[end][1])
                push!(possible_used_labels, words[end])
            end
        end

        address += 1
    end

    # We don't need labels which are never used
    filter(labels) do (k, _)
        k in possible_used_labels
    end
end

"""
    assemble_mnemonic(words, labels) -> Union{Int, Nothing}
Assemble a single mnemonic split into multiple parts. Assume `words` is non-empty.
Words could contain only a label in which case `nothing` would be returned.

```julia-repl
julia> assemble_mnemonic(["somelabel"], Dict("somelabel" => 4))

julia> assemble_mnemonic(["ADD", "42"])
142
```
"""
function assemble_mnemonic(words::Vector{<:AbstractString}, labels=Dict{String, Int}())
    i = findfirst(words) do word
        in(word, keys(opcodes)) || word == "DAT"
    end

    if i == nothing
        for word in words
            if word ∉ keys(labels)
                throw(InvalidMnemonicError(word))
            end
        end
    end

    # Check if we have a data directive and handle it
    if words[i] == "DAT"
		if length(words) > 2
            return parse(Int, words[i+1])
		else
			return 0 # If you value is set default is 0
		end
    end

    # Deal with regular assembly code
    mnemonic = words[i]
    opcode = opcodes[mnemonic]
    if opcode != 0 && rem(opcode, 100) == 0
        arg = words[i+1]
        operand = if haskey(labels, arg)
            labels[arg]
        else
            parse(Int, arg)
        end

        return opcode + operand
    else
        return opcode
    end
end

"""
     assemble(filename)
Turn assembly code found in `filename` into numeric codes (machine code).
Returned as an array of 3-digit decimal numbers.
"""
function assemble(filename::AbstractString)
    lines = readlines(filename)
    program = Int[]

    labels = symboltable(lines)
    @debug "found labels" labels

    for (i, line) in enumerate(lines)
        codeline = remove_comment(line)
        words = split(codeline)

        if isempty(words)
            continue
        end

        try
            instruction = assemble_mnemonic(words, labels)
            if instruction == nothing
                continue
            else
                @debug "parsed line $i:" line words instruction
                push!(program, instruction)
            end
        catch ex
            if isa(ex, InvalidMnemonicError)
                @error "Line $i: Encountered invalid mnemonic '$(ex.mnemonic)'"
            else
                rethrow()
            end
        end
    end
    program
end

"""
     assemble(filename, outfile)
Turn assembly code found in `filename` and store result in `outfile`
"""
function assemble(filename::AbstractString, outfile::AbstractString)
    program = assemble(filename)
    open(outfile, "w") do io
        for instruction in program
            println(io, instruction)
        end
    end
    program
end


text:
struct ParticleFilterSubState{U,I,L}
    traces::SubArray{U,1,Vector{U},I,L}
    new_traces::SubArray{U,1,Vector{U},I,L}
    log_weights::SubArray{Float64,1,Vector{Float64},I,L}
    parents::SubArray{Int,1,Vector{Int},I,L}
end

Gen.get_traces(state::ParticleFilterSubState) = state.traces
Gen.get_log_weights(state::ParticleFilterSubState) = state.log_weights

const ParticleFilterView{U} =
    Union{ParticleFilterState{U}, ParticleFilterSubState{U}} where {U}

function Base.view(state::ParticleFilterState{U},
                   indices::AbstractVector) where {U}
    L = Base.viewindexing((indices,)) == IndexLinear()
    return ParticleFilterSubState{U,Tuple{typeof(indices)},L}(
        view(state.traces, indices),
        view(state.new_traces, indices),
        view(state.log_weights, indices),
        view(state.parents, indices)
    )
end

Base.getindex(state::ParticleFilterState, indices) =
    Base.view(state, indices)

Base.firstindex(state::ParticleFilterState) = 1
Base.lastindex(state::ParticleFilterState) = length(state.traces)


text:
#=

  Set covering in Julia.

  Example from Steven Skiena, The Stony Brook Algorithm Repository
  http://www.cs.sunysb.edu/~algorith/files/set-cover.shtml
  """
  Input Description: A set of subsets S_1, ..., S_m of the
  universal set U = {1,...,n}.

  Problem: What is the smallest subset of subsets T subset S such
  that \cup_{t_i in T} t_i = U?
  """
  Data is from the pictures INPUT/OUTPUT.

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function set_covering_skiena(problem,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   # "all_solutions"=> all_solutions,
                                                            "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    n, m = size(problem)

    @variable(model, x[1:n], Bin)
    @variable(model, 1 <= z <= n, Int)

    for b in eachrow(transpose(problem))
        @constraint(model, sum(x.*b) >= 1)
    end
    @constraint(model,z == sum(x))

    @objective(model,Min,z)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                z_val = convert.(Integer,JuMP.value.(z; result=sol))
                println("x:$x_val")
                println("z:$z_val")
                println()

            end
        end
    else
        println("status:$status")
    end

    return status
end

set_covering_skiena_problems = Dict(
:1 => resize_matrix([[1,1,0,0,0,0,0,0,0,0,0,0],
                    [0,1,0,0,0,0,0,1,0,0,0,0],
                    [0,0,0,0,1,1,0,0,0,0,0,0],
                    [0,0,0,0,0,1,1,0,0,1,1,0],
                    [0,0,0,0,0,0,0,0,1,1,0,0],
                    [1,1,1,0,1,0,0,0,1,1,1,0],
                    [0,0,1,1,0,0,1,1,0,0,1,1]])
)

@time set_covering_skiena(set_covering_skiena_problems[:1])

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using DocumenterGitHubActionsExample
using Documenter

makedocs(;
    modules=[DocumenterGitHubActionsExample],
    authors="Morten Piibeleht <morten.piibeleht@gmail.com> and contributors",
    repo="https://github.com/mortenpi/DocumenterGitHubActionsExample.jl/blob/{commit}{path}#L{line}",
    sitename="DocumenterGitHubActionsExample.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://mortenpi.github.io/DocumenterGitHubActionsExample.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/mortenpi/DocumenterGitHubActionsExample.jl",
)

max_stars_repo_name: mortenpi/DocumenterGitHubActionsExample

max_issues_repo_name: mortenpi/DocumenterGitHubActionsExample

max_forks_repo_name: mortenpi/DocumenterGitHubActionsExample


text:
using DataFrames: DataFrame

export issmooth,
    isdecomposable,
    isstruct_decomposable,
    isdeterministic,
    iscanonical,
    implied_literals,
    sat_prob,
    issatisfiable,
    istautology,
    model_count,
    prob_equiv_signature,
    prob_equiv,
    infer_vtree

#####################
# circuit evaluation infrastructure
#####################

import Base: foreach # extend

function foreach(node::LogicCircuit, f_con::Function, f_lit::Function,
                f_a::Function, f_o::Function, seen=nothing)
    f_leaf(n) = isliteralgate(n) ? f_lit(n) : f_con(n)
    f_inner(n) = is⋀gate(n) ? f_a(n) : f_o(n)
    foreach(node, f_leaf, f_inner, seen)
    nothing # returning nothing helps save some allocations and time
end

import DirectedAcyclicGraphs: foldup # extend

"""
    foldup(node::LogicCircuit,
        f_con::Function,
        f_lit::Function,
        f_a::Function,
        f_o::Function)::T where {T}

Compute a function bottom-up on the circuit.
`f_con` is called on constant gates, `f_lit` is called on literal gates,
`f_a` is called on conjunctions, and `f_o` is called on disjunctions.
Values of type `T` are passed up the circuit and given to `f_a` and `f_o` through a callback from the children.
"""
function foldup(node::LogicCircuit, f_con::Function, f_lit::Function,
                f_a::Function, f_o::Function, ::Type{T}, cache=nothing)::T where {T}
    f_leaf(n) = isliteralgate(n) ? f_lit(n)::T : f_con(n)::T
    f_inner(n, call) = is⋀gate(n) ? f_a(n, call)::T : f_o(n, call)
    foldup(node, f_leaf, f_inner, T, cache)::T
end

import DirectedAcyclicGraphs: foldup_aggregate # extend

"""
    foldup_aggregate(node::LogicCircuit,
        f_con::Function,
        f_lit::Function,
        f_a::Function,
        f_o::Function,
        ::Type{T})::T where T

Compute a function bottom-up on the circuit.
`f_con` is called on constant gates, `f_lit` is called on literal gates,
`f_a` is called on conjunctions, and `f_o` is called on disjunctions.
Values of type `T` are passed up the circuit and given to `f_a` and `f_o` in an aggregate vector from the children.
"""
function foldup_aggregate(node::LogicCircuit, f_con::Function, f_lit::Function,
                          f_a::Function, f_o::Function, ::Type{T}, cache=nothing) where T
    function f_leaf(n)
        isliteralgate(n) ? f_lit(n)::T : f_con(n)::T
    end
    function f_inner(n, cs)
        is⋀gate(n) ? f_a(n, cs)::T : f_o(n, cs)::T
    end
    foldup_aggregate(node, f_leaf::Function, f_inner::Function, T, cache)::T
end

#####################
# variable-scope-based queries
#####################

import .Utils.variables # extend

"""
    variables(root::LogicCircuit)::BitSet

Get a bitset of variables mentioned in the circuit `root`.
"""
function variables(root::LogicCircuit, cache=nothing)::BitSet
    f_con(_) = BitSet()
    f_lit(n) = BitSet(variable(n))
    f_inner(n, call) = mapreduce(call, union, children(n))
    foldup(root, f_con, f_lit, f_inner, f_inner, BitSet, cache)
end

"Get the variable in the circuit with the largest index"
function max_variable(root::LogicCircuit, cache=nothing)::Var
    f_con(n) = Var(0)
    f_lit(n) = variable(n)
    f_inner(n, call) = mapreduce(call, max, children(n))
    foldup(root, f_con, f_lit, f_inner, f_inner, Var, cache)
end

"""
    issmooth(root::LogicCircuit)::Bool

Is the circuit smooth?
"""
function issmooth(root::LogicCircuit, cache=nothing)::Bool
    result::Bool = true
    f_con(_) = BitSet()
    f_lit(n) = BitSet(variable(n))
    f_a(_, cs) = reduce(union, cs)
    f_o(_, cs) = begin
        scope = reduce(union, cs)
        result = result && all(c -> c == scope, cs)
        scope
    end
    foldup_aggregate(root, f_con, f_lit, f_a, f_o, BitSet, cache)
    result
end

"""
    isdecomposable(root::LogicCircuit)::Bool

Is the circuit decomposable?
"""
function isdecomposable(root::LogicCircuit, cache=nothing)::Bool
    result::Bool = true
    f_con(_) = BitSet()
    f_lit(n) = BitSet(variable(n))
    f_a(_, cs) = begin
        result = result && isdisjoint(cs...)
        reduce(union, cs)
    end
    f_o(_, cs) = reduce(union, cs)
    foldup_aggregate(root, f_con, f_lit, f_a, f_o, BitSet, cache)
    result
end

"""
    isstruct_decomposable(root::LogicCircuit)::Bool

Is the circuit structured-decomposable?
"""
function isstruct_decomposable(root::LogicCircuit)::Bool
    infer_vtree(root) !== nothing
end

"""
    infer_vtree(root::LogicCircuit)::Vtree
Infer circuits vtree if the circuit is struct-decomposable.
Otherwise return `nothing`.
"""
function infer_vtree(root::LogicCircuit, cache=nothing)::Union{Vtree, Nothing}
    vtr_dict = Dict{Var, PlainVtreeLeafNode}()

    f_con(_) = (nothing, true) # give constants nothing
    f_lit(n) = begin
        if variable(n) ∉ keys(vtr_dict)
            vtr_dict[variable(n)] = PlainVtree(variable(n))
        end
        (vtr_dict[variable(n)], true)
    end
    f_a(n, call) = begin
        if num_children(n) != 2
            return (call(children(n)[1])[1], false)
        end

        (left, ld) = call(children(n)[1])
        (right, rd) = call(children(n)[2])
        if ((left === nothing) | (right === nothing)) & (left !== right)
            (left === nothing ? right : left, true)
        elseif !isdisjoint(variables(left), variables(right))
            (left, false)
        elseif !has_parent(left) & !has_parent(right)
            (PlainVtreeInnerNode(left, right), ld & rd)
        elseif has_parent(left) & has_parent(right) & (parent(left) == parent(right))
            (parent(left), ld & rd)
        else
            (left, false)
        end
    end
    f_o(n, call) = begin
        @assert num_children(n) > 0 "Or node has no children"
        ccalls = map(x -> call(x), children(n))
        (ccalls[1][1], all(x -> x[1] == ccalls[1][1], ccalls) & all(x -> x[2], ccalls))
    end
    res = foldup(root, f_con, f_lit, f_a, f_o, Tuple{Union{Vtree, Nothing}, Bool})
    res[2] ? res[1] : nothing
end

#####################
# structural properties
#####################

"""
    iscanonical(circuit::LogicCircuit, k::Int; verbose = false)

Does the given circuit have canonical Or gates, as determined by a probabilistic equivalence check?
"""
function iscanonical(circuit::LogicCircuit, k::Int; verbose = false)
   signatures = prob_equiv_signature(circuit, k)
   decision_nodes_by_signature = DirectedAcyclicGraphs.groupby(n -> signatures[n], ⋁_nodes(circuit))
   for (signature, nodes) in decision_nodes_by_signature
      if length(nodes) > 1
         if verbose
            println("Equivalent Nodes:")
            for node in nodes
               println("  - Node: $node MC: $(model_count(node))")
            end
         end
         return false
      end
   end
   return true
end

"""
    isdeterministic(root::LogicCircuit)::Bool

Is the circuit determinstic?
Note: this function is generally intractable for large circuits.
"""
function isdeterministic(root::LogicCircuit, cache=nothing)::Bool
    mgr = sdd_mgr_for(root)
    result::Bool = true
    f_con(c) = compile(mgr,constant(c))
    f_lit(n) = compile(mgr,literal(n))
    f_a(_, cs) = reduce(&, cs)
    f_o(_, cs) = begin
        for i = 1:length(cs)
            for j = i+1:length(cs)
               result &= isfalse(cs[i] & cs[j])
            end
        end
        reduce(|, cs)
    end
    foldup_aggregate(root, f_con, f_lit, f_a, f_o, Sdd, cache)
    result
end

"""
    implied_literals(root::LogicCircuit)::Union{BitSet, Nothing}

Compute at each node literals that are implied by the formula.
`nothing` at a node means all literals are implied (i.e. the node's formula is false)

This algorithm is sound but not complete - all literals returned are correct, but some true implied literals may be missing.
"""
function implied_literals(root::LogicCircuit, cache=nothing)
    f_con(c) = constant(c) ? BitSet() : nothing
    f_lit(n) = BitSet([literal(n)])
    f_a(_, cs) = if any(isnothing, cs)
        nothing
    else
        reduce(union, cs)
    end
    f_o(_, cs) = begin
        reduce(intersect, filter(issomething, cs))
    end
    foldup_aggregate(root, f_con, f_lit, f_a, f_o, Union{BitSet, Nothing}, cache)
end

#####################
# algebraic model counting queries
#####################

"""
    sat_prob(root::LogicCircuit; varprob::Function)::Rational{BigInt}

Get the probability that a random world satisties the circuit.
Probability of each variable is given by `varprob` Function which defauls to 1/2 for every variable.
"""
function sat_prob(root::LogicCircuit, cache=nothing;
                  varprob::Function = v -> BigInt(1) // BigInt(2))::Rational{BigInt}
    f_con(n) = istrue(n) ? one(Rational{BigInt}) : zero(Rational{BigInt})
    f_lit(n) = ispositive(n) ? varprob(variable(n)) : one(Rational{BigInt}) - varprob(variable(n))
    f_a(n, call) = mapreduce(call, *, children(n))
    f_o(n, call) = mapreduce(call, +, children(n))
    foldup(root, f_con, f_lit, f_a, f_o, Rational{BigInt}, cache)
end

"""
    issatisfiable(root::LogicCircuit)::Bool

Determine whether the logical circuit is satisfiable (has a satisfying assignment).
Requires decomposability of the circuit."""
issatisfiable(root::LogicCircuit)::Bool =
    !iszero(sat_prob(root))

"""
istautology(root::LogicCircuit)::Bool

Determine whether the logical circuit is a tautology (every assignment satisfies it; the sentence is valid).
Requires decomposability and determinism of the circuit."""
istautology(root::LogicCircuit)::Bool =
    isone(sat_prob(root))

"""
    model_count(root::LogicCircuit, num_vars_in_scope::Int = num_variables(root))::BigInt

Get the model count of the circuit.
The `num_vars_in_scope` is set to number of variables in the circuit, but sometimes need to set different values,
for example, if not every variable is mentioned in the circuit.
"""
function model_count(root::LogicCircuit, num_vars_in_scope::Int = num_variables(root))::BigInt
    # note that num_vars_in_scope for computing the model count can be more than num_variables(circuit); in particular when variables in the application are missing from the circuit
    BigInt(sat_prob(root) * BigInt(2)^num_vars_in_scope)
end

const Signature = Vector{Rational{BigInt}}

"""
    prob_equiv_signature(circuit::LogicCircuit, k::Int)

Get a signature for each node using probabilistic equivalence checking.
Note that this implentation may not have any formal guarantees as such.
"""
function prob_equiv_signature(circuit::LogicCircuit, k::Int, signs=Dict{Union{Var,LogicCircuit},Signature}())
    # uses probability instead of integers to circumvent smoothing, no mod though
    if signs === nothing
        signs = Dict{Union{Var,LogicCircuit},Signature}()
    end
    prime::Int = 7919 #TODO set as smallest prime larger than num_variables
    randprob() = BigInt(1) .// rand(1:prime,k)
    do_signs(v::Var) = get!(randprob, signs, v)
    f_con(n) = (signs[n] = (istrue(n) ? ones(Rational{BigInt}, k) : zeros(Rational{BigInt}, k)))
    f_lit(n) = (signs[n] = (ispositive(n) ? do_signs(variable(n)) : BigInt(1) .- do_signs(variable(n))))
    f_a(n, call) = (signs[n] = (mapreduce(c -> call(c), (x,y) -> (x .* y), children(n))))
    f_o(n, call) = (signs[n] = (mapreduce(c -> call(c), (x,y) -> (x .+ y), children(n))))
    foldup(circuit, f_con, f_lit, f_a, f_o, Signature)
    signs
end

"""
Check equivalence using probabilistic equivalence checking.
Note that this implentation may not have any formal guarantees as such.
"""
function prob_equiv(circuit1::LogicCircuit, circuit2::LogicCircuit, k::Int)
    signature1 = prob_equiv_signature(circuit1, k)
    signature2 = prob_equiv_signature(circuit2, k, signature1)
    signature1[circuit1] == signature2[circuit2]
end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
module Traits

export traits, trait, Trait

struct Trait{T} end

"Produces type which will only match those which have trait"
@inline trait(x::Type{X}) where X = Trait{T} where {T >: X}
@inline trait(x1, x2) = trait(Union{x1, x2})
@inline trait(x1, x2, x3) = trait(Union{x1, x2, x3})
@inline trait(x1, x2, x3, x4) = trait(Union{x1, x2, x3, x4})
@inline trait(x1, x2, x3, x4, x5) = trait(Union{x1, x2, x3, x4, x5})
@inline trait(x1, x2, x3, x4, x5, xs...) = trait(Union{x1, x2, x3, x4, x5, xs...})

"`traits(::Type{T}` traits of `T`"
function traits end

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
struct RobustCovariance <: CovarianceEstimator end

robust() = RobustCovariance()

function S_hat(x::RegressionModel, ::RobustCovariance)
    m = modelmatrix(x)
    r = residuals(x)
    X2 = zeros(size(m, 1), size(m, 2) * size(r, 2))
    index = 0
    for k in 1:size(r, 2)
        for j in 1:size(m, 2)
            index += 1
            for i in 1:size(m, 1)
                X2[i, index] = m[i, j] * r[i, k]
            end
        end
    end
    S2 = X2' * X2
    Symmetric(rmul!(S2, size(m, 1) / dof_residual(x)))
end

function StatsBase.vcov(x::RegressionModel, v::RobustCovariance)
    xtx = inv(crossmodelmatrix(x))
    pinvertible(Symmetric(xtx * S_hat(x, v) * xtx))
end