text:
### SAVE VERSION OF FUNCTIONS ###

"""
    node(net,rp)
    node(net,s,rp)
    node_f(net,rp)
    node_f(net,s,rp)

Determine node id from relative position and node state.

The node function provides a way to access nodes form a the set of nodes in certain states,
or from the set of all nodes in a simple way. The two-argument version returns the id of the
node at poition *rp* in network *net*. The three-argument version returns the id of the node at
poition *rp* within all nodes in state *s*.

All version of this function run in constant time, but fast (_f) verions sacrifice some safty
checks for better performance. See [basic concepts](concepts.md) for details.

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> randomgraph!(net)
Network of 1000 nodes and 2000 links

julia> nodestate!(net,123,2)

julia> nodestate!(net,345,2)

julia> node(net,2,1)
345

julia> node(net,2,2)
123

julia> node(net,1)
1

julia> destroynode!(net,1)

julia> node(net,1)
998
```
"""
function node(net::FastNet,cls::Int,rp::Int)
    task="Trying to retrieve id of node number $rp in state $cls"
    checknodestate(net,cls,task)
    if rp<1 || rp>net.nclen[cls]
        throw(ArgumentError("$task, but there is no relative position $rp in state $cls"))
    end
    node_f(net,cls,rp)
end

function node(net::FastNet,rp::Int)
    task="Trying to retrieve id of node number $rp"
    if rp<1 || rp>countnodes_f(net)
        throw(ArgumentError("$task, but there is no relative position $rp"))
    end
    node_f(net,rp)
end

"""
    nodecounts(net)
    nodecounts_f(net)

Return an Array containing the number of nodes in the vairous node states.

The time required for this function scales only with the number of node states
(it is independent of the number of nodes).

The alternative (_f) version of this function is identical to nodecounts and is provided
only for convenience.

See also [countnodes](#Fastnet.countnodes)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> randomgraph!(net)
Network of 1000 nodes and 2000 links

julia> for i=1:20
         n=node(net,1,1)
         nodestate!(net,n,2)
       end

julia> nodecounts(net)
2-element Vector{Int64}:
 980
  20
```
"""
function nodecounts(net::FastNet)
    nodecounts_f(net)
end

"""
    nodestate(net,nid)
    nodestate_f(net,nid)

Return the state of the node with id *nid* in network *net*.

All version of this function run in constant time, but fast (_f) verions sacrifice some safty
checks for better performance. See [basic concepts](concepts.md) for details.

See also [nodestate!](#Fastnet.nodestate!)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> randomgraph!(net)
Network of 1000 nodes and 2000 links

julia> nodestate(net,1)
1

julia> nodestate!(net,1,2)

julia> nodestate(net,1)
2
```
"""
function nodestate(net::FastNet,nid)
    task="Trying to determine the state of node $nid"
    checknodeid(net,nid,task)
    checknodeexists(net,nid,task)
    nodestate_f(net,nid)
end

"""
    nodestate!(net,nid,s)
    nodestate_f!(net,nid,s)

Set the node with id *nid* in net *net* to *s*.

Worst-case performance of both versions of this function is O(ks\\*k)+O(ns) where ks is the number of
tracked link states, k is the degree of the affected node and ns is the number of node states.

The fast (_f) verions sacrifice some safty checks for better performance.
See [basic concepts](concepts.md) for details.

See also [nodestate](#Fastnet.nodestate)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> randomgraph!(net)
Network of 1000 nodes and 2000 links

julia> nodestate(net,1)
1

julia> nodestate!(net,1,2)

julia> nodestate(net,1)
2
```
"""
function nodestate!(net::FastNet,nid,cls)
    task="Trying to change the state of node $nid to $cls"
    checknodestate(net,cls,task)
    checknodeexists(net,nid,task)
    nodestate_f!(net,nid,cls)
end

"""
    makenode!(net,s)
    makenode_f!(net,s)

Create a new node in state *s* in the network *net* and return it's id.

Worst-case performance of both versions of this function scales only with the number of node states.

The fast (_f) verions sacrifice some safty checks for better performance.
See [basic concepts](concepts.md) for details.

See also [destroynode!](#Fastnet.destroynode!), [makenodes!](#Fastnet.makenodes!)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> makenode!(net,2)
1

julia> net
Network of 1 nodes and 0 links

julia> nodestate(net,1)
2
```
"""
function makenode!(net::FastNet,cls)
    task="Trying to make a node of state $cls"
    checknodestate(net,cls,task)
    if net.nclen[net.C]<1
        throw(ArgumentError("Out of nodes -- Try creating a larger net"))
    end
    makenode_f!(net,cls)
end

"""
    destroynode!(net,nid)
    destroynode_f!(net,nid)

Destroy the node with id *nid* in network *net*.

Worst-case performance of both versions of this function is O(ks\\*k)+O(ns) where ks is the number of
tracked link states, k is the degree of the affected node and ns is the number of node states.

The fast (_f) verions sacrifice some safty checks for better performance.
See [basic concepts](concepts.md) for details.

See also [makenode!](#Fastnet.makenode!)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> n=makenode!(net,1);

julia> net
Network of 1 nodes and 0 links

julia> destroynode!(net,n)

julia> net
Network of 0 nodes and 0 links
```
"""
function destroynode!(net::FastNet,nid)
    task="Trying to destroy node $nid"
    checknodeid(net,nid,task)
    checknodeexists(net,nid,task)
    destroynode_f!(net,nid)
end

"""
    countnodes(net)
    countnodes(net,s)
    countnodes_f(net)
    countnodes_f(net,s)

Count the nodes in state *s*, or, if no state is provided, in the entire network.

Instead of the state *s* also an Array or Tuple of states can be passed.
In this case the total number of nodes in all of the listed states is returned.

All versions of this function run in constant time.
The fast (_f) verions sacrifice some safty checks for better performance.
See [basic concepts](concepts.md) for details.

If performance is critical use this function rather than nodecounts.

See also [nodecounts](#Fastnet.nodecounts)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> for i=1:20
         makenode!(net,1)
         end

julia> for i=1:10
         makenode!(net,2)
         end

julia> countnodes(net)
30

julia> countnodes(net,1)
20

julia> countnodes(net,(1,2))
30
```
"""
function countnodes(net::FastNet)
    countnodes_f(net)
end

function countnodes(net::FastNet,cls::Integer)
    task="Trying to count nodes in state $cls"
    checknodestate(net,cls,task)
    countnodes_f(net,cls)
end

function countnodes(net::FastNet,cls::Union{Vector,Tuple})
    ret=0
    for c in cls
        ret+=countnodes(net,c)
    end
    ret
end

"""
    firstlinkin(net,nid)
    firstlinkin_f(net,nid)

Return the link id of the first incoming link to the node with id *nid* in network *net*.

If there are no incoming links then the return value is 0

All versions of this function run in constant time.
The fast (_f) verions sacrifice some safty checks for better performance.
See [basic concepts](concepts.md) for details.

See also [firstlinkout](#Fastnet.firstlinkout), [nextlinkin](#Fastnet.nextlinkin)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> n1=makenode!(net,1);

julia> n2=makenode!(net,1);

julia> l1=makelink!(net,n1,n2);

julia> firstlink=firstlinkin(net,n2);

julia> firstlink==l1
true

julia> linksrc(net,firstlink)==n1
true

julia> linkdst(net,firstlink)==n2
true
```
"""
function firstlinkin(net::FastNet,nid)
    task="Trying to retrieve id of first incoming link to node $nid"
    checknodeid(net,nid,task)
    checknodeexists(net,nid,task)
    firstlinkin_f(net,nid)
end

"""
    firstlinkout(net,nid)
    firstlinkout_f(net,nid)

Return the link id of the first outgoing link from the node with id *nid* in network *net*.

If there are no outgoing links then the return value is 0

All versions of this function run in constant time.
The fast (_f) verions sacrifice some safty checks for better performance.
See [basic concepts](concepts.md) for details.

See also [firstlinkout](#Fastnet.firstlinkin), [nextlinkin](#Fastnet.nextlinkout)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> n1=makenode!(net,1);

julia> n2=makenode!(net,1);

julia> l1=makelink!(net,n1,n2);

julia> firstlink=firstlinkout(net,n1);

julia> firstlink==l1
true

julia> linksrc(net,firstlink)==n1
true

julia> linkdst(net,firstlink)==n2
true
```
"""
function firstlinkout(net::FastNet,nid)
    task="Trying to retrieve id of first outgoing link from node $nid"
    checknodeid(net,nid,task)
    checknodeexists(net,nid,task)
    firstlinkout_f(net,nid)
end

"""
    degree(net,nid)
    degree_f(net,nid)

Return the degree of the node with id *nid* in network *net*.

Here degree is interpreted as the number of times this node appears an an endpoint of a link,
hence self-loops contribute 2 to the degree of the node that they link to.

The worst case performance scales only with the degree of the affected node.
The fast (_f) verions sacrifice some safty checks for better performance.
See [basic concepts](concepts.md) for details.

See also [indegree](#Fastnet.indegree), [outdegree](#Fastnet.outdegree)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> n1=makenode!(net,1);

julia> n2=makenode!(net,1);

julia> n3=makenode!(net,1);

julia> makelink!(net,n1,n2);

julia> makelink!(net,n2,n3);

julia> degree(net,n1)
1

julia> degree(net,n2)
2
```
"""
function degree(net::FastNet,nid)
    task="Trying to calculate degree of node $nid"
    checknodeid(net,nid,task)
    checknodeexists(net,nid,task)
    degree_f(net,nid)
end

"""
    indegree(net,nid)
    indegree_f(net,nid)

Return the incoming degree of the node with id *nid* in network *net*.

The worst case performance scales only with the indegree of the affected node.
The fast (_f) verions sacrifice some safty checks for better performance.
See [basic concepts](concepts.md) for details.

See also [degree](#Fastnet.degree), [outdegree](#Fastnet.outdegree)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> n1=makenode!(net,1);

julia> n2=makenode!(net,1);

julia> n3=makenode!(net,1);

julia> makelink!(net,n1,n2);

julia> makelink!(net,n3,n2);

julia> indegree(net,n1)
0

julia> indegree(net,n2)
2
```
"""
function indegree(net::FastNet,nid)
    task="Trying to calculate in-degree of node $nid"
    checknodeid(net,nid,task)
    checknodeexists(net,nid,task)
    indegree_f(net,nid)
end

"""
    outdegree(net,nid)
    outdegree_f(net,nid)

Return the outgoing degree of the node with id *nid* in network *net*.

The worst case performance scales only with the outdegree of the affected node.
The fast (_f) verions sacrifice some safty checks for better performance.
See [basic concepts](concepts.md) for details.

See also [degree](#Fastnet.degree), [indegree](#Fastnet.indegree)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> n1=makenode!(net,1);

julia> n2=makenode!(net,1);

julia> n3=makenode!(net,1);

julia> makelink!(net,n1,n2);

julia> makelink!(net,n3,n2);

julia> outdegree(net,n1)
1

julia> outdegree(net,n2)
0
```
"""
function outdegree(net::FastNet,nid)
    task="Trying to calculate out-degree of node $nid"
    checknodeid(net,nid,task)
    checknodeexists(net,nid,task)
    outdegree_f(net,nid)
end

"""
    nodeexists(net,nid)
    nodeexists_f(net,nid)

Return if a node with id *nid* exists in *net*, false otherwise.

This function runs in constant time.
The fast (_f) verions sacrifice some safty checks for better performance.
See [basic concepts](concepts.md) for details.

See also [makenode!](#Fastnet.makenode!), [destroynode!](#Fastnet.destroynode!)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> n1=makenode!(net,1);

julia> nodeexists(net,n1)
true
```
"""
function nodeexists(net::FastNet,nid)
    if nid<1 || nid>net.N
        throw(ArgumentError("Trying to check if node $nid is in the network but that id does not refer to a node"))
    end
    nodeexists_f(net,nid)
end

"""
    makenodes!(net,N,s)
    makenodes_f!(net,N,s)

Create *N* nodes in state *s* in the network *net.

Worst case performance of this function scales only with the number of node states.
The fast (_f) verions sacrifice some safty checks for better performance.
See [basic concepts](concepts.md) for details.

See also [makenode!](#Fastnet.makenode!), [destroynode!](#Fastnet.destroynode!)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> makenodes!(net,100,1);

julia> net
Network of 100 nodes and 0 links
```
"""
function makenodes!(net::FastNet,N,cls)
    task="Trying to create $N nodes of state $cls"
    checknodestate(net::FastNet,cls,task)
    N>net.nclen[net.C] ? throw(ArgumentError("$task, but this is more nodes than allowed by net")) :
    makenodes_f!(net,N,cls)
end

"""
    randomnode(net)
    randomnode(net,s)
    randomnode_f(net)
    randomnode_f(net,s)

Return the id of a random node drawn from *net*.

If the second argument *s* is not provided the node will be drawn uniformly from
all nodes in the network. If *s* is an integer then the node will be drawn uniformly
from the nodes in state *s*. If *s* is an Array or Tuple of Ints then the node will be
drawn uniformly from the nodes in the states listed.

This function runs in constant time if *s* is integer or omitted. If *s* is an Array or Tuple the
worst case performance scales only with the number of node states.
The fast (_f) verions sacrifice some safty checks for better performance.
See [basic concepts](concepts.md) for details.

The safe versions of this function will throw an ArgumentError with an informative error message
when trying to pick a node from an empty set. With the fast (_f) version, trying to pick a node from
an empty set will also result in an ArgumentError being thrown, but in this case the message will be
something like "Range must be non-empty".

See also [node](#Fastnet.node)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> makenodes!(net,100,1);

julia> makenodes!(net,100,2);

julia> nodestate(net,randomnode(net,1))
1
julia> nodestate(net,randomnode(net,2))
2
```
"""
function randomnode(net::FastNet)
    task="Trying to select random node"
    checknonnull(net,task)
    randomnode_f(net)
end

function randomnode(net::FastNet,cls::Integer)
    task="Trying to select random node from state $cls"
    checknodestate(net,cls,task)
    checknonnull(net,cls,task)
    randomnode_f(net,cls)
end

function randomnode(net::FastNet,cls::Union{Array,Tuple})
    task="Trying to select a random node from a set of states"
    for c in cls
        checknodestate(net,c,task)
    end
    tot=countnodes(net,cls)
    if tot===0
        throw(ArgumentError("Trying to select a random node from a set of states, but there isn't any node in any of the states."))
    end
    randomnode_f(net,cls)
end

"""
    adjacent(net,a,b)
    adjacent_f(net,a,b)

Check if nodes with the id's *a* and *b* are adjacent in network *net*. If they are
return the id a of the link connecting them. Otherwise return 0.

If multiple links connect the nodes the function will return a link in the direction
a->b if such a link exists.

Calling adjacent(net,a,a) will return a self-loop on a if one exists.

The worst case performance of this function scales with the degree of node a.
The fast (_f) verions sacrifice some safty checks for better performance.
See [basic concepts](concepts.md) for details.

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> makenodes!(net,100,1);

julia> net
Network of 100 nodes and 0 links
```
"""
function adjacent(net::FastNet,nida,nidb)
    task="Checking if there is a link connecting nodes $nida and $nidb,"
    checknodeid(net,nida,task)
    checknodeid(net,nidb,task)
    adjacent_f(netmnida,nidb)
end


text:
"""
    pwm(x::Vector{<:Real},p::Int,r::Int,s::Int)

Compute the empirical probability weighted moments.

The probability weighted moments are defined by [Landwehr *et al. (1979)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/WR015i005p01055)
as follows:
```math
M_{p,r,s} = \\mathbb{E}\\left[ X^p F^r(X) \\left\\{ 1-F(X) \\right\\}^s  \\right].
```
The unbiased empirical estimates is computed using the formula given by [Landwehr *et al. (1979)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/WR015i005p01055).

# Reference

Landwehr, J. M., Matalas, N. C. and Wallis, J. R. (1979). Probability weighted moments compared with
    some traditional techniques in estimating Gumbel parameters and quantiles. *Water Resources Research*,
    15:1055–1064.

"""
function pwm(x::Vector{<:Real},p::Int,r::Int,s::Int)::Real

 @assert sign(p)>=0 "p should be a non-negative integer."
 @assert sign(r)>=0 "r should be a non-negative integer."
 @assert sign(s)>=0 "s should be a non-negative integer."

 y = sort(x)
 n = length(y)

 m = 1/n*sum( y[i]^p * binomial(i-1,r)/binomial(n-1,r) * binomial(n-i,s)/binomial(n-1,s) for i=1:n )

 return m

end

include(joinpath("probabilityweightedmoment", "probabilityweightedmoment_gev.jl"))
include(joinpath("probabilityweightedmoment", "probabilityweightedmoment_gp.jl"))
include(joinpath("probabilityweightedmoment", "probabilityweightedmoment_gumbel.jl"))

ext: jl


text:
@params struct SensFilter{_filtering, T, TV <: AbstractVector{T}} <: AbstractSensFilter
    filtering::Val{_filtering}
    elementinfo::ElementFEAInfo
    metadata::FilterMetadata
    rmin::T
    nodal_grad::TV
    last_grad::TV
    cell_weights::TV
end
Base.show(::IO, ::MIME{Symbol("text/plain")}, ::SensFilter) = println("TopOpt sensitivity filter")

SensFilter{true}(args...) = SensFilter(Val(true), args...)
SensFilter{false}(args...) = SensFilter(Val(false), args...)

function SensFilter(solver::AbstractFEASolver; rmin)
    return SensFilter(Val(true), solver, rmin)
end
function SensFilter(::Val{filtering}, solver::AbstractFEASolver, args...) where {filtering}
    return SensFilter(Val(filtering), whichdevice(solver), solver, args...)
end
function SensFilter(::Val{true}, ::CPU, solver::TS, rmin::T, ::Type{TI}=Int) where {T, TI<:Integer, TS<:AbstractFEASolver}
    metadata = FilterMetadata(solver, rmin, TI)
    TM = typeof(metadata)
    problem = solver.problem
    elementinfo = solver.elementinfo
    grid = problem.ch.dh.grid
    nnodes = getnnodes(grid)
    nodal_grad = zeros(T, nnodes)
    TV = typeof(nodal_grad)

    black = problem.black
    white = problem.white
    nel = length(black)
    nfc = sum(black) + sum(white)
    last_grad = zeros(T, nel-nfc)

    cell_weights = zeros(T, nnodes)

    return SensFilter(Val(true), elementinfo, metadata, rmin, nodal_grad, last_grad, cell_weights)
end

function SensFilter(::Val{false}, ::CPU, solver::TS, rmin::T, ::Type{TI}=Int) where {T, TS<:AbstractFEASolver, TI<:Integer}
    elementinfo = solver.elementinfo
    metadata = FilterMetadata(T, TI)
    nodal_grad = T[]
    last_grad = T[]
    cell_weights = T[]
    return SensFilter(Val(false), elementinfo, metadata, rmin, nodal_grad, last_grad, cell_weights)
end

(cf::SensFilter)(x) = x
function ChainRulesCore.rrule(cf::SensFilter{true}, x)
    x, Δ -> begin
        cf.rmin <= 0 && return (nothing, Δ)
        newΔ = copy(Δ)
        @unpack elementinfo, nodal_grad, cell_weights, metadata = cf
        @unpack black, white, varind, cellvolumes, cells = elementinfo
        @unpack cell_neighbouring_nodes, cell_node_weights = metadata
        node_cells = elementinfo.metadata.node_cells
        update_nodal_grad!(nodal_grad, node_cells, cell_weights, cells, cellvolumes, black, white, varind, Δ)
        normalize_grad!(nodal_grad, cell_weights)
        update_grad!(newΔ, black, white, varind, cell_neighbouring_nodes, cell_node_weights, nodal_grad)
        return (nothing, newΔ)
    end
end

function update_nodal_grad!(nodal_grad::AbstractVector, node_cells, cell_weights, cells, cellvolumes, black, white, varind, grad)
    T = eltype(nodal_grad)
    for n in 1:length(nodal_grad)
        nodal_grad[n] = zero(T)
        cell_weights[n] = zero(T)
        r = node_cells.offsets[n]:node_cells.offsets[n+1]-1
        for i in r
            c = node_cells.values[i][1]
            if black[c] || white[c]
                continue
            end
            ind = varind[c]
            w = cellvolumes[c]
            cell_weights[n] += w
            nodal_grad[n] += w * grad[ind]
        end
    end
    return nodal_grad
end

function normalize_grad!(nodal_grad::AbstractVector, cell_weights)
    for n in 1:length(nodal_grad)
        if cell_weights[n] > 0
            nodal_grad[n] /= cell_weights[n]
        end
    end
end

function update_grad!(grad::AbstractVector, black, white, varind, cell_neighbouring_nodes, cell_node_weights, nodal_grad)
    @inbounds for i in 1:length(black)
        if black[i] || white[i]
            continue
        end
        ind = varind[i]
        nodes = cell_neighbouring_nodes[ind]
        if length(nodes) == 0
            continue
        end
        weights = cell_node_weights[ind]
        grad[ind] = dot(view(nodal_grad, nodes), weights) / sum(weights)
    end

    return
end


text:
import Combinatorics: permutations

function parsedata(raw_data)
    data = []
    for line in raw_data
        line = split(line, "\t")
        line = parse.(Int, line)
        push!(data, line)
    end
    data
end

function solve_a(data)
    total = 0

    for line in data
        total += maximum(line) - minimum(line)
    end
    total
end

function solve_b(data)
    total = 0
    for line in data
        for (i, j) in permutations(line,2)
            if i % j == 0
                total += i ÷ j
            end
        end
    end
    total
end

raw_data = readlines("input.txt")
data = parsedata(raw_data)
sol_a = solve_a(data)
println("sol a: $(sol_a)")

sol_b = solve_b(data)
println("sol b: $(sol_b)")


text:
# Needed for defining the rule for `adjoint(A::ITensor)`
# which currently doesn't work by overloading `ChainRulesCore.rrule`
using ZygoteRules: @adjoint

@adjoint function Base.adjoint(x::Union{ITensor,MPS,MPO})
  y = prime(x)
  function adjoint_pullback(ȳ)
    x̄ = inv_op(prime, ȳ)
    return (x̄,)
  end
  return y, adjoint_pullback
end


text:
using SymmetricFormats, Test, LinearAlgebra

n=10
for elty in (Float32, Float64)
    @testset "spr! $elty" begin
        α = rand(elty)
        M = rand(elty, n, n)
        AL = Symmetric(M, :L)
        AU = Symmetric(M, :U)
        x = rand(elty, n)

        ALP_result_julia_lower = SymmetricFormats.pack(α*x*x' + AL, :L)
        ALP_result_blas_lower = SymmetricFormats.pack(AL, :L)
        SymmetricFormats.spr!('L', α, x, ALP_result_blas_lower)
        @test ALP_result_julia_lower ≈ ALP_result_blas_lower

        AUP_result_julia_upper = SymmetricFormats.pack(α*x*x' + AU, :U)
        AUP_result_blas_upper = SymmetricFormats.pack(AU, :U)
        SymmetricFormats.spr!('U', α, x, AUP_result_blas_upper)
        @test AUP_result_julia_upper ≈ AUP_result_blas_upper
    end
end


text:
using BitBlit2D, Colors

function checkered_board(rows=4, cols=4)
    w, h = 4, 4
    board = Image(cols*4, rows*4, colorant"yellow")
    square = Image(4, 4, colorant"brown")

    for r in 0:rows-1
        for c in 0:cols-1
            if r % 2 == c % 2
                draw!(board, square, c*w, r*h)
            end
        end
    end
    board
end


text:
struct DoubleQuadRule{P,Q}
  outer_quad_points::P
  inner_quad_points::Q
end

"""
    regularcellcellinteractions!(biop, tshs, bshs, tcell, bcell, interactions, strat)

Function for the computation of moment integrals using simple double quadrature.
"""
function momintegrals!(biop, tshs, bshs, tcell, bcell, z, strat::DoubleQuadRule)

    # memory allocation here is a result from the type instability on strat
    # which is on purpose, i.e. the momintegrals! method is chosen based
    # on dynamic polymorphism.
    womps = strat.outer_quad_points
    wimps = strat.inner_quad_points

    for womp in womps
        tgeo = womp.point
        tvals = womp.value
        M = length(tvals)
        jx = womp.weight

        for wimp in wimps
            bgeo = wimp.point
            bvals = wimp.value
            N = length(bvals)
            jy = wimp.weight

            j = jx * jy
            kernel = kernelvals(biop, tgeo, bgeo)

            for n in 1 : N
                bval = bvals[n]
                for m in 1 : M
                    tval = tvals[m]
                    igd = integrand(biop, kernel, tval, tgeo, bval, bgeo)
                    z[m,n] += j * igd
                end
            end
        end
    end

    return z
end


text:
#=
isrtdmok

reset crc

load pipe with data to write, any type
unload pipe to uart with crc

load pipe with bytes from uart, calculating running crc
load pipe with bytes from uart, no crc update
unload pipe to read data, any type
=#

# interface to rest of module is with iopipe not iopipe2

function sendpipewithcrc(i::RTDMinterface)  # always attach and send crc
  while ~eof(i.iopipe)
    byte = read(i.iopipe, UInt8)
    write(i.iopipe2, byte)
    i.crc.data = byte
    compute_crc(i.crc)
  end
  write(i.iopipe2, i.crc.crc)
  write(i.iouart, takebuf_array(i.iopipe2))
  resetcrc!(i)
end

function loadpipe(i::RTDMinterface, bytestoload::Integer, last2crc = bytestoload>1)
  # set last2crc = true to not include last two bytes read in crc calculation
  data = read(i.iouart, UInt8, bytestoread)
  for (index, element) in enumerate(data)
    if element < bytestoread-1
      i.crc.data = element
      compute_crc(i.crc)
    end
    write(i.iopipe,data)
  end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

module DSAModule

using LinearAlgebra, ForwardDiff, Tensors, NLsolve, Parameters

import ..AbstractMaterial, ..AbstractMaterialState
import ..Utilities: Symm2, Symm4, isotropic_elasticity_tensor, lame, debang
import ..integrate_material!  # for method extension

# parametrically polymorphic for any type representing ℝ
export GenericDSA, GenericDSADriverState, GenericDSAParameterState, GenericDSAVariableState

# specialization for Float64
export DSA, DSADriverState, DSAParameterState, DSAVariableState

@with_kw mutable struct GenericDSADriverState{T <: Real} <: AbstractMaterialState
    time::T = zero(T)
    strain::Symm2{T} = zero(Symm2{T})
end

"""Parameter state for DSA (dynamic strain aging) material.

This is similar to the Chaboche model, but with additional static recovery terms.

Parameters:

  - `E`: Young's modulus
  - `nu`: Poisson's ratio
  - `R0`: initial yield strength
  - `Kn`: plasticity multiplier divisor (drag stress)
  - `nn`: plasticity multiplier exponent
  - `C1`, `D1`: parameters governing behavior of backstress X1
  - `C2`, `D2`: parameters governing behavior of backstress X2
  - `Q`: shift parameter for yield strength evolution
  - `b`: multiplier for yield strength evolution
  - `w`: controls the average waiting time a dislocation is arrested at localized obstacles.

    It represents a strain increment produced when all arrested
    dislocations overcome localized obstacles, and move toward the next
    pinned configuration.

    In practice, this parameter controls how fast the effective aging time
    reacts to plastic flow: \$\\dot{t}_a = 1 - t_a \\dot{p} / w\$

  - `P1`, `P2`: controls the maximum hardening in the fully aged state.

    Has the units of stress.

  - `m`: controls the characteristic diffusion time. Depends on the type of diffusion.

    The value `1/3` is thought to represent pipe diffusion along dislocation lines.
    Another typical value is `2/3`.

  - `m1`, `m2`: The exponent of the power-law type static recovery of backstresses.

    The static recovery mechanism becomes activated at higher temperatures.

    This parameter controls the secondary creep and constant slope relaxation of
    stresses over a longer period of time. Higher values (>6..10) effectively
    deactivate static recovery, whereas lower values (<5) activate it.

  - `M1`, `M2`: The normalizer of the power-law type static recovery of backstresses.

    Has the units of stress. Can be used to activate/deactivate static recovery.
    Deactivation occurs with high values.

  - `ba`: Controls the rate of evolution of aging stress to its asymptotic value.

    Dimensionless. Similar to the isotropic hardening `b`.

  - `xi`: Controls the magnitude of the Marquis effect from the aging stress.

    The Marquis effect is that increased hardening due to aging shows as
    increased relaxation.

    Dimensionless. Support `[0,1]`.

    At `0`, the aging stress contributes solely to the size of the yield surface `R`
    (isotropic hardening).

    At `1`, the aging stress contributes solely to the viscoplastic drag stress `K`.
"""
@with_kw struct GenericDSAParameterState{T <: Real} <: AbstractMaterialState
    E::T = 0.0
    nu::T = 0.0
    R0::T = 0.0
    Kn::T = 0.0
    nn::T = 0.0
    C1::T = 0.0
    D1::T = 0.0
    C2::T = 0.0
    D2::T = 0.0
    Q::T = 0.0
    b::T = 0.0
    w::T = 0.0
    P1::T = 0.0
    P2::T = 0.0
    m::T = 0.0
    m1::T = 0.0
    m2::T = 0.0
    M1::T = 0.0
    M2::T = 0.0
    ba::T = 0.0
    xi::T = 0.0
end

"""Problem state for DSA material.

- `stress`: stress tensor
- `X1`: backstress 1
- `X2`: backstress 2
- `plastic_strain`: plastic part of strain tensor
- `cumeq`: cumulative equivalent plastic strain (scalar, ≥ 0)
- `R`: yield strength
- `ta`: effective aging time
- `Ra`: aging stress
- `jacobian`: ∂σij/∂εkl
"""
@with_kw struct GenericDSAVariableState{T <: Real} <: AbstractMaterialState
    stress::Symm2{T} = zero(Symm2{T})
    X1::Symm2{T} = zero(Symm2{T})
    X2::Symm2{T} = zero(Symm2{T})
    plastic_strain::Symm2{T} = zero(Symm2{T})
    cumeq::T = zero(T)
    R::T = zero(T)
    ta::T = zero(T)
    Ra::T = zero(T)
    jacobian::Symm4{T} = zero(Symm4{T})
end

# TODO: Does this eventually need a {T}?
@with_kw struct DSAOptions <: AbstractMaterialState
    nlsolve_method::Symbol = :trust_region
end

@with_kw mutable struct GenericDSA{T <: Real} <: AbstractMaterial
    drivers::GenericDSADriverState{T} = GenericDSADriverState{T}()
    ddrivers::GenericDSADriverState{T} = GenericDSADriverState{T}()
    variables::GenericDSAVariableState{T} = GenericDSAVariableState{T}()
    variables_new::GenericDSAVariableState{T} = GenericDSAVariableState{T}()
    parameters::GenericDSAParameterState{T} = GenericDSAParameterState{T}()
    dparameters::GenericDSAParameterState{T} = GenericDSAParameterState{T}()
    options::DSAOptions = DSAOptions()
end

DSADriverState = GenericDSADriverState{Float64}
DSAParameterState = GenericDSAParameterState{Float64}
DSAVariableState = GenericDSAVariableState{Float64}
DSA = GenericDSA{Float64}

"""
    state_to_vector(sigma::U, R::T, X1::U, X2::U, ta::T, Ra::T) where U <: Symm2{T} where T <: Real

Adaptor for `nlsolve`. Marshal the problem state into a `Vector`.
"""
function state_to_vector(sigma::U, R::T, X1::U, X2::U, ta::T, Ra::T) where U <: Symm2{T} where T <: Real
    return [tovoigt(sigma); R; tovoigt(X1); tovoigt(X2); ta; Ra]::Vector{T}
end

"""
    state_from_vector(x::AbstractVector{<:Real})

Adaptor for `nlsolve`. Unmarshal the problem state from a `Vector`.
"""
function state_from_vector(x::AbstractVector{T}) where T <: Real
    sigma::Symm2{T} = fromvoigt(Symm2{T}, @view x[1:6])
    R::T = x[7]
    X1::Symm2{T} = fromvoigt(Symm2{T}, @view x[8:13])
    X2::Symm2{T} = fromvoigt(Symm2{T}, @view x[14:19])
    ta::T = x[20]
    Ra::T = x[21]
    return sigma, R, X1, X2, ta, Ra
end

"""
    integrate_material!(material::GenericDSA{T}) where T <: Real

Material model with dynamic strain aging (DSA). This is similar to the Chaboche
material with two backstresses, with both kinematic and isotropic hardening, but
this model also features static recovery terms.

This model captures dynamic (and static) strain aging (DSA) induced hardening.
The related phenomena are:

  - Portevin le Chatelier effect. Serrated yield, plastic instabilities.
  - Discontinuous yielding
  - Inverse strain rate sensitivity (inverse SRS)
  - Secondary hardening in low cycle fatigue (LCF) tests

These typically occur in a certain temperature/strain rate regime, where the
dislocations are pinned due to the diffusion of solute atoms. In the most
effective conditions, the speed of diffusion is comparable to the applied
strain rate (speed of dislocations).

See:

    J.-L. Chaboche, A. Gaubert, P. Kanouté, A. Longuet, F. Azzouz, M. Mazière.
    Viscoplastic constitutive equations of combustion chamber materials including
    cyclic hardening and dynamic strain aging. International Journal of Plasticity
    46 (2013), 1--22. http://dx.doi.org/10.1016/j.ijplas.2012.09.011

Further reading:

    M. Mazière, H. Dierke. Investigations on the Portevin Le Chatelier critical
    strain in an aluminum alloy. Computational Materials Science 52(1) (2012),
    68--72. https://doi.org/10.1016/j.commatsci.2011.05.039
"""
function integrate_material!(material::GenericDSA{T}) where T <: Real
    p  = material.parameters
    v  = material.variables
    dd = material.ddrivers
    d  = material.drivers
    @unpack E, nu, R0, Kn, nn, C1, D1, C2, D2, Q, b, w, P1, P2, m, m1, m2, M1, M2, ba, xi = p
    lambda, mu = lame(E, nu)

    @unpack strain, time = d
    dstrain  = dd.strain
    dtime    = dd.time
    @unpack stress, X1, X2, plastic_strain, cumeq, R, jacobian, ta, Ra = v

    # elastic part
    jacobian = isotropic_elasticity_tensor(lambda, mu)
    stress  += dcontract(jacobian, dstrain)

    # resulting deviatoric plastic stress (accounting for backstresses Xm)
    seff_dev = dev(stress - X1 - X2)
    # von Mises yield function
    f = sqrt(1.5)*norm(seff_dev) - (R0 + R + (1 - xi) * Ra)  # using elastic trial problem state
    if f > 0.0
        g! = create_nonlinear_system_of_equations(material)
        x0 = state_to_vector(stress, R, X1, X2, ta + dtime, Ra)
        res = nlsolve(g!, x0; method=material.options.nlsolve_method, autodiff = :forward)
        converged(res) || error("Nonlinear system of equations did not converge!")
        x = res.zero
        stress, R, X1, X2, ta, Ra = state_from_vector(x)

        # using the new problem state
        seff_dev = dev(stress - X1 - X2)
        f = sqrt(1.5)*norm(seff_dev) - (R0 + R + (1 - xi) * Ra)

        dotp = ((f >= 0.0 ? f : 0.0) / (Kn + xi * Ra))^nn
        dp = dotp*dtime
        n = sqrt(1.5)*seff_dev/norm(seff_dev)

        plastic_strain += dp*n
        cumeq += dp

        # Compute the new Jacobian, accounting for the plastic contribution.
        drdx = ForwardDiff.jacobian(debang(g!), x)
        drde = zeros((length(x), 6))
        drde[1:6, 1:6] = -tovoigt(jacobian)  # elastic Jacobian. Follows from the defn. of g!.
        jacobian = fromvoigt(Symm4, (drdx\drde)[1:6, 1:6])
    else
        ta += dtime
    end
    variables_new = GenericDSAVariableState{T}(stress = stress,
                                               X1 = X1,
                                               X2 = X2,
                                               R = R,
                                               plastic_strain = plastic_strain,
                                               cumeq = cumeq,
                                               jacobian = jacobian,
                                               ta = ta,
                                               Ra = Ra)
    material.variables_new = variables_new
    return nothing
end

"""
    create_nonlinear_system_of_equations(material::GenericDSA{T}) where T <: Real

Create and return an instance of the equation system for the incremental form of
the evolution equations of the DSA material.

Used internally for computing the plastic contribution in `integrate_material!`.

The input `material` represents the problem state at the end of the previous
timestep. The created equation system will hold its own copy of that state.

The equation system is represented as a mutating function `g!` that computes the
residual:

```julia
    g!(F::V, x::V) where V <: AbstractVector{<:Real}
```

Both `F` (output) and `x` (input) are length-21 vectors containing
[sigma, R, X1, X2, ta, Ra], in that order. The tensor quantities
sigma, X1, X2 are encoded in Voigt format.

The function `g!` is intended to be handed over to `nlsolve`.
"""
function create_nonlinear_system_of_equations(material::GenericDSA{T}) where T <: Real
    p  = material.parameters
    v  = material.variables
    dd = material.ddrivers
    d  = material.drivers
    @unpack E, nu, R0, Kn, nn, C1, D1, C2, D2, Q, b, w, P1, P2, m, m1, m2, M1, M2, ba, xi = p
    lambda, mu = lame(E, nu)

    # Old problem state (i.e. the problem state at the time when this equation
    # system instance was created).
    #
    # Note this does not include the elastic trial; this is the state at the
    # end of the previous timestep.
    @unpack strain, time = d
    dstrain = dd.strain
    dtime = dd.time
    @unpack stress, X1, X2, plastic_strain, cumeq, R, ta, Ra = v

    jacobian = isotropic_elasticity_tensor(lambda, mu)

    # Compute the residual. F is output, x is filled by NLsolve.
    # The solution is x = x* such that g(x*) = 0.
    function g!(F::V, x::V) where V <: AbstractVector{<:Real}
        stress_new, R_new, X1_new, X2_new, ta_new, Ra_new = state_from_vector(x)  # tentative new values from nlsolve

        seff_dev = dev(stress_new - X1_new - X2_new)
        f = sqrt(1.5)*norm(seff_dev) - (R0 + R_new + (1 - xi) * Ra_new)

        dotp = ((f >= 0.0 ? f : 0.0) / (Kn + xi * Ra_new))^nn
        dp = dotp*dtime
        n = sqrt(1.5)*seff_dev/norm(seff_dev)

        # The equations are written in an incremental form.
        # TODO: multiply the equations by -1 to make them easier to understand in the context of the rest of the model.

        dstrain_plastic = dp*n
        dstrain_elastic = dstrain - dstrain_plastic
        tovoigt!(view(F, 1:6), stress - stress_new + dcontract(jacobian, dstrain_elastic))

        F[7] = R - R_new + b*(Q - R_new)*dp

        # HACK: The zero special case is needed here to make ForwardDiff happy.
        #
        # Otherwise, when ndX1_new = 0, the components 2:end of the automatic
        # derivative of JX1_new will be NaN, which causes the calculation of the
        # material jacobian to silently fail. This usually manifests itself as a
        # mysterious convergence failure, when this model is used in the strain
        # optimizer.
        ndX1_new = norm(dev(X1_new))
        if iszero(ndX1_new)
            JX1_new = 0.0
        else
            JX1_new = sqrt(1.5) * ndX1_new
        end
        sr1_new = (JX1_new^(m1 - 1) * X1_new) / (M1^m1)  # static recovery term
        tovoigt!(view(F, 8:13), X1 - X1_new + dp*(2.0/3.0*C1*n - D1*X1_new) - dtime*sr1_new)

        ndX2_new = norm(dev(X2_new))
        if iszero(ndX2_new)
            JX2_new = 0.0
        else
            JX2_new = sqrt(1.5) * ndX2_new
        end
        sr2_new = (JX2_new^(m2 - 1) * X2_new) / (M2^m2)  # static recovery term
        tovoigt!(view(F, 14:19), X2 - X2_new + dp*(2.0/3.0*C2*n - D2*X2_new) - dtime*sr2_new)

        Ras = P1 * (1.0 - exp(-P2 * ta_new^m))
        F[20] = ta - ta_new + dtime - (ta_new / w)*dp
        F[21] = Ra - Ra_new + ba*(Ras - Ra_new)*dp
        return nothing
    end
    return g!
end

end


text:
import PipeLayout: pwl_ineqs, reorient_fwdflow, pwl_inverse, pipesplit

struct SOC <: JunctionLocationSolver
    optimizer          # for underlying NLP model
    timelimit::Float64 # seconds

    function SOC(solver; timelimit=Inf)
        new(solver, timelimit)
    end
end

"Convex reformulation of NLP model with SOC"
function make_soc(inst::Instance, topo::Topology, solver::SOC)
    nodes, nnodes = topo.nodes, length(topo.nodes)
    arcs, narcs = topo.arcs, length(topo.arcs)
    terms, nterms = inst.nodes, length(inst.nodes)
    termidx = termindex(nodes, terms)
    ndiams = length(inst.diameters)
    xmin, xmax = extrema(node.x for node in terms)
    ymin, ymax = extrema(node.y for node in terms)

    # tail and head for arcs
    T = [a.tail for a in arcs]
    H = [a.head for a in arcs]

    # pressure bounds for all nodes
    πlb_min = minimum([b.lb for b in inst.pressure])^2
    πub_max = maximum([b.ub for b in inst.pressure])^2
    πlb = fill(πlb_min, nnodes)
    πub = fill(πub_max, nnodes)
    πlb[termidx] = [b.lb^2 for b in inst.pressure]
    πub[termidx] = [b.ub^2 for b in inst.pressure]

    # need correct orientation (= positive flow)
    q = uniq_flow(inst, topo)
    @assert all(q .> 0)

    c = [diam.cost for diam in inst.diameters]
    Dm5 = [diam.value^(-5) for diam in inst.diameters]
    C = inst.ploss_coeff .* q .* abs.(q)

    model = JuMP.Model(solver.optimizer)

    # node positions
    @variable(model, xmin ≤ x[1:nnodes] ≤ xmax)
    @variable(model, ymin ≤ y[1:nnodes] ≤ ymax)

    # differences
    @variable(model, xmin - xmax ≤ dx[1:narcs] ≤ xmax - xmin)
    @variable(model, ymin - ymax ≤ dy[1:narcs] ≤ ymax - ymin)

    # squared pressure on nodes
    @variable(model, πlb[v] ≤ π[v=1:nnodes] ≤ πub[v])

    # auxiliary variable for pipe cost
    @variable(model, t[1:narcs] ≥ 0)

    # auxiliary variable for second-order cones
    @variable(model, L[1:narcs] ≥ 0)

    # fix terminal positions
    @constraint(model, fix[t=1:nterms], x[termidx[t]] == terms[t].x)
    @constraint(model, fiy[t=1:nterms], y[termidx[t]] == terms[t].y)

    # compute differences
    @constraint(model, diffx[a=1:narcs], dx[a] == x[T[a]] - x[H[a]])
    @constraint(model, diffy[a=1:narcs], dy[a] == y[T[a]] - y[H[a]])

    # second order cone for length (dx[a]^2 + dy[a]^2 ≤ L[a]^2)
    @constraint(model, soc[a=1:narcs],
                [L[a], dx[a], dy[a]] in JuMP.SecondOrderCone())

    # inequalities for PWL function
    ie = pwl_ineqs(Dm5, c)
    @assert size(ie) == (ndiams - 1, 3)
    @constraint(model, pwl[a=1:narcs, s=1:ndiams-1],
                ie[s,1]/C[a]*(π[T[a]] - π[H[a]]) + ie[s,2]*t[a] ≥ ie[s,3]*L[a])

    # lower bound for diameter impact (y)
    Dmax = inst.diameters[end].value
    @constraint(model, ylb[a=1:narcs], Dmax^5/C[a]*(π[T[a]] - π[H[a]]) ≥ L[a])

    # lower bound for cost factor (instead of upper bound on y)
    cmin = inst.diameters[1].cost
    @constraint(model, yub[a=1:narcs], t[a] ≥ cmin * L[a])

    # minimize total pipe cost
    @objective(model, Min, sum(t[a] for a=1:narcs))

    model, x, y, t, π
end

function PipeLayout.optimize(inst::Instance, topo::Topology, solver::SOC)
    rotopo = reorient_fwdflow(inst, topo)

    model, x, y, t, π = make_soc(inst, topo, solver)
    JuMP.optimize!(model)
    status = JuMP.termination_status(model)
    objval = JuMP.objective_value(model)
    nodes = map(Node, zip(JuMP.value.(x), JuMP.value.(y)))

    # compute l from t
    tsol = JuMP.value.(t)
    L = pipelengths(Topology(nodes, topo.arcs))
    D, c = [d.value for d in inst.diameters], [d.cost for d in inst.diameters]
    zsol = tsol ./ L
    pwl_xs, pwl_ys = reverse(D.^(-5)), reverse(c) # must be sorted
    ysol = [pwl_inverse(pwl_xs, pwl_ys, z) for z in zsol]
    lsol = vcat([pipesplit(inst.diameters, y^(-1/5))' for y in ysol]...)
    @assert size(lsol) == (length(L), length(D))

    sol = Solution(nodes, lsol, JuMP.value.(π))
    Result(status, sol, objval)
end


text:
# # Analyses
#
# ## Histogram

using SplitApplyPlot, CairoMakie

df = (x=randn(1000), y=randn(1000), z=rand(["a", "b", "c"], 1000))
specs = data(df) * mapping(:x, layout=:z) * histogram(bins=range(-2, 2, length=15)) * visual()
draw(specs)
AbstractPlotting.save("hist.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](hist.svg)
#

data(df) * mapping(:x, :y, layout=:z) * histogram(bins=15) |> draw
AbstractPlotting.save("hist2D.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](hist2D.svg)
#

# ## Density

using SplitApplyPlot, CairoMakie

df = (x=randn(1000), y=randn(1000), z=rand(["a", "b", "c"], 1000))
data(df) * mapping(:x, layout=:z) * SplitApplyPlot.density() |> draw
AbstractPlotting.save("density.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](density.svg)
#

data(df) * mapping(:x, :y, layout=:z) * SplitApplyPlot.density(npoints=20) |> draw
AbstractPlotting.save("density2D.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](density2D.svg)
#

# ## Linear

using SplitApplyPlot, CairoMakie

df = (x=randn(30), y=randn(30), z=rand(["a", "b", "c"], 30))
specs = data(df) * mapping(:x, :y, color=:z) * (linear() + visual(Scatter))
draw(specs, axis=(; limits=(-2, 2, -2, 2)))
AbstractPlotting.save("linearanalysis.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](linearanalysis.svg)
#


text:
@diffusion_process FitzHughDiffusionConjugAuxLin begin
    :dimensions
    process --> 2
    wiener --> 1

    :parameters
    (ϵ, s, γ, β, σ, t, T) --> Float64
    (u, v) --> ℝ{2}

    :additional
    constdiff --> true
    linear --> true
end

function B(t, P::FitzHughDiffusionConjugAuxLin)
    @SMatrix [
        0.0  1.0;
        (P.ϵ-P.γ-3.0*P.ϵ*P.v[1]^2-6*P.ϵ*P.v[1]*P.v[2]) (P.ϵ-1.0-3.0*P.ϵ*P.v[1]^2)
    ]
end

function β(t, P::FitzHughDiffusionConjugAuxLin)
    ℝ{2}(
        0.0,
        2*P.ϵ*P.v[1]^3+P.s-P.β+6*P.ϵ*P.v[1]^2*P.v[2]
    )
end

σ(t, P::FitzHughDiffusionConjugAuxLin) = ℝ{2}(0.0, P.σ)


text:
push!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))

import ClimaCore.Geometry, LinearAlgebra, UnPack
import ClimaCore:
    Fields,
    Domains,
    Topologies,
    Meshes,
    DataLayouts,
    Operators,
    Geometry,
    Spaces

using OrdinaryDiffEq: OrdinaryDiffEq, ODEProblem, solve, SSPRK33

using Logging: global_logger
using TerminalLoggers: TerminalLogger
global_logger(TerminalLogger())

const FT = Float64

# https://github.com/CliMA/CLIMAParameters.jl/blob/master/src/Planet/planet_parameters.jl#L5
const MSLP = 1e5 # mean sea level pressure
const grav = 9.8 # gravitational constant
const R_d = 287.058 # R dry (gas constant / mol mass dry air)
const γ = 1.4 # heat capacity ratio
const C_p = R_d * γ / (γ - 1) # heat capacity at constant pressure
const C_v = R_d / (γ - 1) # heat capacit at constant volume
const R_m = R_d # moist R, assumed to be dry

const f = 5e-5
const ν = 0.01
const L = 2e2
const nelems = 30
const Cd = ν / (L / nelems)
const ug = 1.0
const vg = 0.0
const d = sqrt(2 * ν / f)
domain = Domains.IntervalDomain(
    Geometry.ZPoint{FT}(0.0),
    Geometry.ZPoint{FT}(L);
    boundary_tags = (:bottom, :top),
)
#mesh = Meshes.IntervalMesh(domain, Meshes.ExponentialStretching(7.5e3); nelems = 30)
mesh = Meshes.IntervalMesh(domain; nelems = nelems)

cspace = Spaces.CenterFiniteDifferenceSpace(mesh)
fspace = Spaces.FaceFiniteDifferenceSpace(cspace)

# https://github.com/CliMA/Thermodynamics.jl/blob/main/src/TemperatureProfiles.jl#L115-L155
# https://clima.github.io/Thermodynamics.jl/dev/TemperatureProfiles/#DecayingTemperatureProfile
function adiabatic_temperature_profile(z; T_surf = 300.0, T_min_ref = 230.0)
    u = FT(ug)
    v = FT(vg)
    return (u = u, v = v)
end

zc = Fields.coordinate_field(cspace)
Yc = adiabatic_temperature_profile.(zc.z)
w = Geometry.Cartesian3Vector.(zeros(Float64, fspace))

Y_init = copy(Yc)
w_init = copy(w)

# Y = (Yc, w)

function tendency!(dY, Y, _, t)

    (Yc, w) = Y.x
    (dYc, dw) = dY.x

    UnPack.@unpack u, v = Yc
    du = dYc.u
    dv = dYc.v

    # S 4.4.1: potential temperature density
    # Mass conservation

    u_1 = parent(u)[1]
    v_1 = parent(v)[1]
    u_wind = sqrt(u_1^2 + v_1^2)
    A = Operators.AdvectionC2C(
        bottom = Operators.SetValue(0.0),
        top = Operators.SetValue(0.0),
    )

    # u-momentum
    bcs_bottom =
        Operators.SetValue(Geometry.Cartesian3Vector(Cd * u_wind * u_1))  # Eq. 4.16
    bcs_top = Operators.SetValue(FT(ug))  # Eq. 4.18
    gradc2f = Operators.GradientC2F(top = bcs_top)
    divf2c = Operators.DivergenceF2C(bottom = bcs_bottom)
    @. du = divf2c(ν * gradc2f(u)) + f * (v - vg) - A(w, u)   # Eq. 4.8

    # v-momentum
    bcs_bottom =
        Operators.SetValue(Geometry.Cartesian3Vector(Cd * u_wind * v_1))  # Eq. 4.17
    bcs_top = Operators.SetValue(FT(vg))  # Eq. 4.19
    gradc2f = Operators.GradientC2F(top = bcs_top)
    divf2c = Operators.DivergenceF2C(bottom = bcs_bottom)
    @. dv = divf2c(ν * gradc2f(v)) - f * (u - ug) - A(w, v)   # Eq. 4.9

    return dY
end

using LinearAlgebra
using RecursiveArrayTools

Y = ArrayPartition(Yc, w)
dY = tendency!(similar(Y), Y, nothing, 0.0)

Δt = 2.0
ndays = 0
# Solve the ODE operator
prob = ODEProblem(tendency!, Y, (0.0, 60 * 60 * 50))
sol = solve(
    prob,
    SSPRK33(),
    dt = Δt,
    saveat = 600, # save 10 min
    progress = true,
    progress_message = (dt, u, p, t) -> t,
);

ENV["GKSwstype"] = "nul"
import Plots
Plots.GRBackend()

dirname = "ekman"
path = joinpath(@__DIR__, "output", dirname)
mkpath(path)

z_centers = parent(Fields.coordinate_field(cspace))
z_faces = parent(Fields.coordinate_field(fspace))

function ekman_plot(u; title = "", size = (1024, 600))
    u_ref =
        ug .-
        exp.(-z_centers / d) .*
        (ug * cos.(z_centers / d) + vg * sin.(z_centers / d))
    sub_plt1 = Plots.plot(
        u_ref,
        z_centers,
        marker = :circle,
        xlabel = "u",
        label = "Ref",
    )
    sub_plt1 =
        Plots.plot!(sub_plt1, parent(u.x[1].u), z_centers, label = "Comp")

    v_ref =
        vg .+
        exp.(-z_centers / d) .*
        (ug * sin.(z_centers / d) - vg * cos.(z_centers / d))
    sub_plt2 = Plots.plot(
        v_ref,
        z_centers,
        marker = :circle,
        xlabel = "v",
        label = "Ref",
    )
    sub_plt2 =
        Plots.plot!(sub_plt2, parent(u.x[1].v), z_centers, label = "Comp")

    return Plots.plot(
        sub_plt1,
        sub_plt2,
        title = title,
        layout = (1, 2),
        size = size,
    )
end

anim = Plots.@animate for (i, u) in enumerate(sol.u)
    ekman_plot(u, title = "Hour $(i)")
end
Plots.mp4(anim, joinpath(path, "ekman.mp4"), fps = 10)

Plots.png(ekman_plot(sol[end]), joinpath(path, "ekman_end.png"))

function linkfig(figpath, alt = "")
    # buildkite-agent upload figpath
    # link figure in logs if we are running on CI
    if get(ENV, "BUILDKITE", "") == "true"
        artifact_url = "artifact://$figpath"
        print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n")
    end
end

linkfig("output/$(dirname)/ekman_end.png", "ekman End")


text:
using Test
using Plots

import Distributions

using Scruff
using Scruff.Operators
using Scruff.MultiInterface
using Scruff.Models
using Scruff.SFuncs

import Scruff: make_initial, make_transition
import Scruff.Models: get_dt
import Scruff.Operators: Sample, sample, Logcpdf, logcpdf, Marginalize, marginalize

struct MyModel <: Model{Tuple{}, Tuple{Int}, Tuple{Int, Int}} end

function make_transition(m::MyModel, parenttimes, time)
    return Constant((parenttimes[1],time))
end

function make_initial(m::MyModel, time)
    return Cat([(1,0)], [1.0])
end

struct MyNormal <: SFunc{Tuple{}, Float32}
    mean::Float32
    var::Float32
end

@impl begin
    struct MyNormalExpectation end
    function expectation(sf::MyNormal, i::Tuple{})::Float32
        # Have access the sfunc as `sf` as well as `op_imp`
        return sf.mean
    end
end

@impl begin
    struct MyNormalVariance end
    function variance(sf::MyNormal, i::Tuple{})::Float32
        return sf.var
    end
end

@impl begin
    struct MyNormalSample end
    function sample(sf::MyNormal, i::Tuple{})::Float32
        return rand(Distributions.Normal(sf.mean, sqrt(sf.var)))
    end
end

@impl begin
    struct MyNormalLogcpdf end
    function logcpdf(sf::MyNormal, i::Tuple{}, o::Float32)::AbstractFloat
        return Distributions.logpdf(Distributions.Normal(sf.mean, sqrt(sf.var)), o)
    end
end
struct MyCondMuNormal <: SFunc{Tuple{Float32}, Float32}
    # A conditional distribution for a Normal conditioned on mu with fixed var.
    # CondMuNormal_var(mu) = N(mu,var), essentially.
    var::Float32
end

@impl begin
    struct MyCondMuNormalSample end
    function sample(sf::MyCondMuNormal, x::Tuple{Float32})::Float32
        return rand(Distributions.Normal(x[1], sqrt(sf.var)))
    end
end

@impl begin
    struct MyCondMuNormalLogcpdf end
    function logcpdf(sf::MyCondMuNormal, i::Tuple{Float32}, o::Float32)::AbstractFloat
        return Distributions.logpdf(Distributions.Normal(i[1], sqrt(sf.var)), o)
    end
end

@impl begin
    struct MyCondMuNormalMarginalize end
    function marginalize(sf::MyCondMuNormal, x::MyNormal)::MyNormal
        mu = expectation(x, tuple())
        var = variance(x, tuple()) + sf.var
        return MyNormal(mu, var)
    end
end

@impl begin
    struct MyCondMuNormalExpectation end
    function expectation(sf::MyCondMuNormal, x::Tuple{Float32, Tuple{}})::Float32
        return x[1]
    end
end

@impl begin
    struct MyCondMuNormalVariance end
    function variance(sf::MyCondMuNormal, i::Tuple{Float32, Tuple{}})::Float32
        return sf.var
    end
end

randomwalk = HomogeneousModel(MyNormal(0.0, 1.0), MyCondMuNormal(1.0), 2.0)

struct WienerProcess <: VariableTimeModel{Tuple{}, Tuple{Float32}, Float64}
    # A continuous limit of Gaussian random walk
    k::Float32 # "Rate" of random walk. Units of var/time
end
function make_transition(m::WienerProcess,parenttimes,time)
    var = m.k*abs(time-parenttimes[1])
    return MyCondMuNormal(var)
end
function make_initial(m::WienerProcess)
    return MyNormal(0,1)
end
wienerprocess = WienerProcess(0.1)

@testset "newcore" begin
    @testset "instantiate!" begin
        @testset "instant network" begin
            @testset "correctly instantiates variables and placeholders" begin
                v = MyModel()(:v)
                p = Placeholder{Tuple{Int,Int}}(:p)
                net = InstantNetwork(Variable[v], VariableGraph(v => [p]), Placeholder[p])
                run = Runtime(net)
                ensure_all!(run)
                inst1 = current_instance(run, p)
                inst2 = current_instance(run, v)
                @test inst1 isa PlaceholderInstance
                @test inst2 isa VariableInstance
                @test get_sfunc(inst2) isa Cat
            end
        end

        @testset "dynamic network" begin
            @testset "correctly calls make_initial or make_transition" begin
                v = MyModel()(:v)
                u = Constant(1)()(:u)
                net = DynamicNetwork(Variable[u,v], VariableGraph(), VariableGraph(v => [u]))
                run = Runtime(net)
                instantiate!(run, u, 1)
                instantiate!(run, v, 2)
                instantiate!(run, v, 4)
                inst2 = get_instance(run, v, 2)
                inst4 = get_instance(run, v, 4)
                @test get_sfunc(inst2) isa Cat
                @test get_sfunc(inst4) isa Constant
            end

            @testset "correctly instantiates variables and placeholders" begin
                v = MyModel()(:v)
                p = Placeholder{Int}(:p)
                net = DynamicNetwork(Variable[v], VariableGraph(),
                    VariableGraph(v => [p]), Placeholder[p])
                run = Runtime(net)
                inst1 = instantiate!(run, p, 1)
                inst2 = instantiate!(run, v, 2)
                @test isa(inst1, PlaceholderInstance)
                @test isa(inst2, VariableInstance)
            end

            @testset "passes the correct times to the model" begin
                v = MyModel()(:v)
                net = DynamicNetwork(Variable[], VariableGraph(), VariableGraph(v => [v]))
                run = Runtime(net)
                inst2 = instantiate!(run, v, 2)
                inst4 = instantiate!(run, v, 4)
                inst3 = instantiate!(run, v, 3)
                inst5 = instantiate!(run, v, 5)
                # Each instance should have a previous time of the previously existing latest instance, if any
                @test get_sfunc(inst4).x == (2,4)
                @test get_sfunc(inst3).x == (2,3)
                @test get_sfunc(inst5).x == (4,5)
            end
        end
    end

    @testset "models" begin
        @testset "HomogeneousModel" begin
            @testset "Initial model" begin
                sf = make_initial(randomwalk)
                @test sf isa MyNormal
                # Test some ops
                _sample = sample(sf, ())
                _logcpdf = logcpdf(sf, (), _sample)
            end

            @testset "Transition model" begin
                sf = make_transition(randomwalk,0.,2.)
                @test sf isa MyCondMuNormal
                # Test some ops
                _sample = sample(sf, (0.0f0,))
                _logcpdf = logcpdf(sf, (0.0f0,), _sample)
            end
        end

        @testset "VariableTimeModel" begin
            @testset "Initial model" begin
                sf = make_initial(wienerprocess)
                @test sf isa MyNormal
                # Test some ops
                _sample = sample(sf, ())
                _logcpdf = logcpdf(sf, (), _sample)
            end

            @testset "Transition model" begin
                sf = make_transition(wienerprocess,(0.,),2.)
                @test sf isa MyCondMuNormal
                # Test some ops
                _sample = sample(sf, (0.0f0,))
                _logcpdf = logcpdf(sf, (0.0f0,), _sample)
            end
        end

    end

end


text:
# this should be formatted
a = f(aaa, bbb, ccc)

# this should not be formatted
#! format: off
a = f(aaa,
    bbb,ccc)

c = 102000

d = @foo 10 20

e = "what the foocho"

# ok


text:
#__precompile__()
module KoalaEnsembles

# new:
export EnsembleRegressor, fit_weights!, weight_regularization_curve

# for use in this module:
import Koala: Regressor, BaseType, SupervisedMachine, type_parameters, params
import Koala: err, supports_ensembling, softwarn
import KoalaTrees: TreeRegressor
import StatsBase
import DataFrames: AbstractDataFrame

# to be extended (but not explicitly rexported):
import Koala: setup, fit, predict
import Koala: default_transformer_X, default_transformer_y, transform, inverse_transform

# development only:
# import ADBUtilities: @dbg, @colon

# constants:
const EPS = eps(Float64)

## Weighted ensembles of predictors

mutable struct WeightedEnsemble{P, Atom <: Regressor{P}} <: BaseType
    model::Atom
    ensemble::Vector{P}
    weights::Vector{Float64}
end

WeightedEnsemble(model::P, ensemble::Vector{P}, weights) where P =
    WeightedEnsemble{P, Atom}(model, ensemble, weights)

function Base.show(stream::IO,
                   object::WeightedEnsemble{P, Atom}) where {P, Atom<:Regressor{P}}
    abbreviated(n) = "..."*string(n)[end-2:end]
    type_string = string("{", Atom, "}")
    print(stream, string(typeof(object).name.name,
                         type_string,
                         "@", abbreviated(hash(object))))
end

function predict(wens::WeightedEnsemble, X)
    ensemble = wens.ensemble
    weights = wens.weights

    !isempty(weights) || error("Empty ensemble cannot make predictions.")

    n_regressors = length(ensemble)
    prediction = weights[1]*predict(wens.model, ensemble[1], X, false, false)
    i =2
    while i <= n_regressors
        prediction += weights[i]*predict(wens.model, ensemble[i], X, false, false)
        i += 1
    end
    return prediction
end

## The `EnsembleRegressor` type

mutable struct EnsembleRegressor{P, Atom <: Regressor{P}} <: Regressor{WeightedEnsemble{P, Atom}}

    atom::Atom
    weight_regularization::Float64
    bagging_fraction::Float64
    n::Int # number of predictors created (added) in call to fit (with add=true)

end

function clean!(model::EnsembleRegressor)
    message = ""
    if model.bagging_fraction > 1 || model.bagging_fraction <= 0
        message = message*"`bagging_fraction` should be "*
        "in the range (0,1]. Reset to 1. "
        model.bagging_fraction = 1.0
    end
    if model.weight_regularization > 1 || model.weight_regularization < 0
        message = message*"`weight_regularization` should be "*
        "in the range [0,1]. Reset to 1. "
        model.weight_regularization = 1.0
    end
    if !supports_ensembling(model.atom)
        message = message*"Ensembling not supported for given `atom`. "*
        "Resetting to `TreeRegressor()`"
        model.atom = TreeRegressor()
    end
    return message
end

# constructor to infer type automatically:
EnsembleRegressor(atom::Atom, weight_regularization,
                  bagging_fraction, n::Int) where {P, Atom<:Regressor{P}} =
                      EnsembleRegressor{P, Atom}(atom, weight_regularization, bagging_fraction, n)

# lazy keyword constructor:
function EnsembleRegressor(;atom=TreeRegressor(), weight_regularization=1,
    bagging_fraction=0.8, n::Int=100)

    model = EnsembleRegressor(atom, weight_regularization, bagging_fraction, n)
    softwarn(clean!(model))
    return model
end

function Base.show(stream::IO, object::EnsembleRegressor)
    abbreviated(n) = "..."*string(n)[end-2:end]
    type_params = type_parameters(object)
    if isempty(type_params)
        type_string = ""
    else
        type_string = string("{", type_params[2], "}")
    end
    print(stream, string(typeof(object).name.name,
                         type_string,
                         "@", abbreviated(hash(object))))
end

# transformers are inherited from atom:
default_transformer_X(model::EnsembleRegressor) = default_transformer_X(model.atom)
default_transformer_y(model::EnsembleRegressor) = default_transformer_y(model.atom)

function setup(model::EnsembleRegressor{P, Atom},
               Xt, yt, scheme_X, parallel, verbosity) where {P, Atom<:Regressor{P}}
    ensemble = Array{P}(0)
    return Xt, yt, scheme_X, ensemble
end

# Note: Whenceforth I use "X" and "y" instead of "Xt" and "yt"

function fit(model::EnsembleRegressor{P, Atom}, cache, add, parallel,
             verbosity; optimize_weights_only=false) where {P,
             Atom<:Regressor{P}}

    X, y, scheme_X, ensemble_so_far = cache
    n = model.n
    nbr_so_far = length(ensemble_so_far)
    n_patterns = length(y)
    n_train = round(Int, floor(model.bagging_fraction*n_patterns))

    # core ensemble building function:
    function get_ensemble(n, verbosity::Int)

        # initialize random number generator:
        srand((round(Int,time()*1000000)))

        ensemble = Vector{P}(n)

        for i in 1:n
            N = add ? i + nbr_so_far : i
            verbosity < 1 || print("\rComputing regressor number: $N          ")
            train_rows = StatsBase.sample(1:n_patterns, n_train, replace=false)
            atom_cache = setup(model.atom, X, y, scheme_X, false, verbosity - 1)
            atom_predictor, atom_report, atom_cache =
                fit(model.atom, atom_cache, false, false, verbosity - 1)
            ensemble[i] = atom_predictor
        end
        verbosity < 1 || println()

        return ensemble

    end

    ## Build or supplement ensemble if required

    if optimize_weights_only # build not required
        ensemble = ensemble_so_far
    else # build required
        if !parallel || nworkers() == 1 # build in serial
            ensemble = get_ensemble(n, verbosity)
        else # build in parallel
            if verbosity >= 1
                println("Ensemble-building in parallel on $(nworkers()) processors.")
            end
            chunk_size = div(n, nworkers())
            left_over = mod(n, nworkers())
            ensemble =  @parallel (vcat) for i = 1:nworkers()
                if i != nworkers()
                    get_ensemble(chunk_size, 0) # 0 means silent
                else
                    get_ensemble(chunk_size + left_over, 0)
                end
            end
        end
        # include an existing ensemble if required:
        if add
            ensemble = vcat(ensemble_so_far, ensemble)
        end
    end

    ## Optimize weights

    n = length(ensemble)

    if model.weight_regularization == 1
        weights = ones(n)/n
        verbosity < 1 || info("Weighting regressors uniformly.")
    else
        verbosity < 1 || print("\nOptimizing weights...")
        Y = Array{Float64}(n, n_patterns)
        for k in 1:n
            Y[k,:] = predict(model.atom, ensemble[k], X, false, false)
        end

        # If I rescale all predictions by the same amount it makes no
        # difference to the values of the optimal weights:
        ybar = mean(abs.(y))
        Y = Y/ybar

        A = Y*Y'
        b = Y*(y/ybar)

        scale = abs(det(A))^(1/n)

        if scale < EPS

            verbosity < 0 || warn("Weight optimization problem ill-conditioned. " *
                 "Using uniform weights.")
            weights = ones(n)/n

        else

            # need regularization, `gamma`, between 0 and infinity:
            if model.weight_regularization == 0
                gamma = 0
            else
                gamma = exp(atanh(2*model.weight_regularization - 1))
            end

            # add regularization and augment linear system for constraint
            # (weights sum to one)
            AA = hcat(A + scale*gamma*eye(n), ones(n))
            AA = vcat(AA, vcat(ones(n), [0.0])')
            bb = b + scale*gamma*ones(n)/n
            bb = vcat(bb, [1.0])

            weights = (AA \ bb)[1:n] # drop Lagrange multiplier
            verbosity < 1 || println("\r$n weights optimized.\n")

        end

    end

    predictor = WeightedEnsemble(model.atom, ensemble, weights)
    report = Dict{Symbol, Any}()
    report[:normalized_weights] = weights*length(weights)

    cache = (X, y, scheme_X, ensemble)

    return predictor, report, cache

end

predict(model::EnsembleRegressor, predictor, Xt, parallel, verbosity) =
    predict(predictor, Xt)

function fit_weights!(mach::SupervisedMachine{WeightedEnsemble{P, Atom},
                                              EnsembleRegressor{P, Atom}};
              verbosity=1, parallel=true) where {P, Atom <: Regressor{P}}

    mach.n_iter != 0 || error("Cannot fit weights to empty ensemble.")

    mach.predictor, report, mach.cache =
        fit(mach.model, mach.cache, false, parallel, verbosity;
            optimize_weights_only=true)
    merge!(mach.report, report)

    return mach
end

function weight_regularization_curve(mach::SupervisedMachine{WeightedEnsemble{P, Atom},
                                                           EnsembleRegressor{P, Atom}},
                                   test_rows;
                                   verbosity=1, parallel=true, range=linspace(0,1,101),
                                   raw=false) where {P, Atom <: Regressor{P}}

    mach.n_iter > 0 || error("No regressors in the ensemble. Run `fit!` first.")
    !raw || verbosity < 0 ||
        warn("Reporting errors for *transformed* target. Use `raw=false` "*
             " to report true errors.")

    if parallel && nworkers() > 1
        if verbosity >= 1
            println("Optimizing weights in parallel on $(nworkers()) processors.")
        end
        errors = pmap(range) do w
            verbosity < 2 || print("\rweight_regularization=$w       ")
            mach.model.weight_regularization = w
            fit_weights!(mach; parallel=false, verbosity=verbosity - 1)
            err(mach, test_rows, raw=raw)
        end
    else
        errors = Float64[]
        for w in range
            verbosity < 1 || print("\rweight_regularization=$w       ")
            mach.model.weight_regularization = w
            fit_weights!(mach; parallel= parallel, verbosity=verbosity - 1)
            push!(errors, err(mach, test_rows, raw=raw))
        end
        verbosity < 1 || println()

        mach.report[:weight_regularization_curve] = (range, errors)
    end

    return range, errors
end

end # module


text:
using WeibullCount

begin
    x = 2
    c = 0.36787944117144233
    λ = 20.08553692

    d = WeibullCountModel(λ, c)

    @test pdf(d, x) ≥ 0
    @test all(pdf(d, [x]) .≥ 0)

    @test cdf(d, x) ≥ 0
    @test all(cdf(d, [x]) .≥ 0)
end


text:
## Let's solve the bidding problem
using SigmoidalProgramming
using Random

# problem data
Random.seed!(4)

function bidding(nvar=500; problemtype=:paper)
    nconstr = 1
    l = fill(-10, nvar)
    u = fill(10, nvar)
    A = ones(nconstr, nvar)

    # boring bidding
    if problemtype==:boring
    B = [.4*nvar]
    z = zeros(nvar)
    fs = Function[x -> i/(1 + exp(-x)) for i=1:nvar]
    dfs = Function[x-> i*exp(-x)/(1 + exp(-x))^2 for i=1:nvar]

    # bidding example from the paper
    elseif problemtype==:paper
        v = 4*rand(nvar)
        beta = -3*v
        alpha = 10
        B = [.2*sum(v)]
        z = -beta/alpha
        fs = Function[x->logistic(alpha*x + beta[i]) - logistic(beta[i]) for i=1:nvar]
        dfs = Function[x->alpha*logistic_prime(alpha*x + beta[i]) for i=1:nvar]
    else
        error("did not recognize problem type $problemtype")
    end

    problem = LinearSP(fs, dfs, z, A, B)

    # branch and bound
    pq, bestnodes, lbs, ubs = @time solve_sp(l, u, problem, verbose=false;
                                             TOL=nvar*1e-3, maxiters=1000)
    println("$(length(lbs)) iterations")
    node = bestnodes[end] #dequeue!(pq)
    # println("best node has node.ub = $(node.ub) and solution $(node.x)")
    println("lb: ",lbs[end])
    println("ub: ",ubs[end])
    println("ratio: ",(ubs[end]-lbs[end])/lbs[end])

    return pq, bestnodes, lbs, ubs
end

# for size=[1,10,20,50,100,500,1000,2000,5000,10000]
#     println("bidding $size")
#     bidding(size)
# end

pq, bestnodes, lbs, ubs = bidding(20)


text:
using Revise, TensND, LinearAlgebra, SymPy, Tensors, OMEinsum, Rotations, Test

# Spheroidal
ϕ = symbols("ϕ", real = true)
p = symbols("p", real = true)
p̄ = √(1 - p^2)
q = symbols("q", positive = true)
q̄ = √(q^2 - 1)
c = symbols("c", positive = true)
coords = (ϕ, p, q)
OM = Tens(c * [p̄ * q̄ * cos(ϕ), p̄ * q̄ * sin(ϕ), p * q])
rules = Dict(
    sqrt(1 - p^2) * sqrt(q^2 - 1) => sqrt(-(p^2 - 1) * (q^2 - 1)),
    sqrt((p^2 - q^2) / (p^2 - 1)) * sqrt(1 - p^2) => sqrt(q^2 - p^2),
)
rules = Dict(
    sqrt(-(p^2 - 1) * (q^2 - 1)) => sqrt(1 - p^2) * sqrt(q^2 - 1),
    sqrt((p^2 - q^2) / (p^2 - 1)) * sqrt(1 - p^2) => sqrt(q^2 - p^2),
)
Spheroidal = CoorSystemSym(OM, coords; rules = rules)

ϕ, p = symbols("ϕ p", real = true);
p̄, q, q̄, c = symbols("p̄ q q̄ c", positive = true);
coords = (ϕ, p, q);
tmp_coords = (p̄, q̄);
params = (c,);
OM = Tens(c * [p̄ * q̄ * cos(ϕ), p̄ * q̄ * sin(ϕ), p * q]);
Spheroidal = CoorSystemSym(
    OM,
    coords,
    tmp_coords,
    params;
    tmp_var = Dict(1 - p^2 => p̄^2, q^2 - 1 => q̄^2),
    to_coords = Dict(p̄ => √(1 - p^2), q̄ => √(q^2 - 1)),
);
simplify(LAPLACE(OM[1]^2, Spheroidal))
m = 2;
n = 5;
P = sympy.assoc_legendre;
T = P(n, m, p) * P(n, m, q) * cos(m * ϕ);
simplify(LAPLACE(T, Spheroidal))


text:
module MEngProject
include("./LaminartFunc.jl")
include("./LaminartInitFunc.jl")
include("./LaminartEqImfilter.jl")
include("./LaminartEqImfilterGPU_FFT.jl")
include("./LaminartEqImfilterGPU_FIR.jl")
include("./LaminartEqImfilterGPU_IIR.jl")
include("./LaminartEqConv.jl")
include("./LaminartKernels.jl")
include("./Utils.jl")
include("./Parameters.jl")
export LaminartFunc, LaminartInitFunc, LaminartEqImfilter, LaminartEqConv, LaminartKernels, Utils, Parameters, LaminartEqImfilterGPU_FFT, LaminartEqImfilterGPU_FIR, LaminartEqImfilterGPU_IIR

end


text:
module MechaJulia
import Requests
import JSON

const auth_token = readall("../src/token")
const METADATA   = "IainNZ/METADATA.jl"

const MECHAJULIA_HEAD = """
*This is an automatic message from [MechaJulia](https://github.com/MechaJulia/MechaJulia),
a Julia bot that aims to help maintain Julia and its package ecosystem. If you think there
is something wrong with this message, or have ideas for a new feature, please
[file an issue](https://github.com/MechaJulia/MechaJulia/issues/new)*

---
"""

const MSG_VER000 = """
* It looks like you are adding a version 0.0.0, but this is discouraged.
  Consider tagging a "real" first version like 0.0.1 or, if you feel the
  package is not ready, encourage testers to obtain it with `Pkg.clone`"""

function is_there_ver000(file_list)
    """
    Given a list of files parsed from the JSON response from
      GET /repos/:owner/:repo/pulls/:number/files
      https://developer.github.com/v3/pulls/#list-pull-requests-files
    look for a new 0.0.0 version being added (which is discouraged).
    """
    for file in file_list
        # Only look at added files, anything else probably is either
        # tidy up or a fix to a wrong SHA
        file["status"] != "added" && continue

        if contains(file["filename"], "versions/0.0.0")
            return true
        end
    end
    return false
end

function does_require_match(file_list)
    """
    Given a list of files parsed from the JSON response from
      GET /repos/:owner/:repo/pulls/:number/files
      https://developer.github.com/v3/pulls/#list-pull-requests-files
    look for `requires` files, and cross-reference them with the REQUIRE
    file in the package itself
    """
    for file in file_list
        file["status"] == "deleted" && continue
        !contains(file["filename"])

function check_pr_ok()
    # Pull all pull requests
    # GET /repos/:owner/:repo/pulls
    # https://developer.github.com/v3/pulls/#list-pull-requests
    url = "https://api.github.com/repos/$METADATA/pulls?access_token=$(auth_token)"
    raw = Requests.get(url; json = {"sort" => "updated"})
    prs = raw.data |> JSON.parse

    # Walk through them all looking for things to comment on that we haven't already
    for pr in prs
        println(pr["title"])
        println(pr["number"])
        println(pr["_links"]["issue"]["href"])

        message = MECHAJULIA_HEAD

        # Look at files being committed
        number = int(pr["number"])
        url    = "https://api.github.com/repos/$METADATA/pulls/$number/files?access_token=$(auth_token)"
        raw    = Requests.get(url)
        files  = raw.data |> JSON.parse
        if is_there_ver000(files)
            message *= MSG_VER000
        end

        println(message)
    end
end

end

max_stars_repo_name: MechaJulia/MechaJulia

max_issues_repo_name: MechaJulia/MechaJulia

max_forks_repo_name: MechaJulia/MechaJulia


text:
module HCBE
using SimpleHypergraphs, Random,  ProgressBars, LightGraphs, ProgressMeter, GraphPlot, StatsBase, JSON, StatsPlots

include("./algorithms/HypergraphClustering.jl")
include("./algorithms/Weighted.jl")
include("./algorithms/ClusteringUtil.jl")
include("./IO.jl")
end # module


text:
# TODO default analysis for metadata here

abstract type MetadataAnalysis <: AbstractAnalysis end


text:
#=
The following file contains various function to help smoke test and benchmark
arbitrary AbstractRBMs and configurations.
=#

using Test
using SparseArrays

"""
The TestReporter uses the ratio test or D'Alembert's criterion
to monitor convergence. This is helpful in automating tests that
confirm that a given RBM is learning something.

NOTE: current we just store the calculated ratios from each epoch,
but we could probably just do an online mean calculation.
"""
mutable struct TestReporter
    ratios::Vec{Float64}
    prev::Float64
    log::Bool

    function TestReporter(log=false)
        new(Array{Float64, 1}(undef, 0), NaN, log)
    end
end

"""
Computes and adds a new ratio to the reporter. If log is set to
true in the reporter a message containing the score, time and mean ratio
will be printed.
"""
function report(r::TestReporter, rbm::AbstractRBM,
                epoch::Int, epoch_time::Float64, score::Float64)
    if epoch == 1
        r.ratios = Array{Float64, 1}(undef, 0)
        r.prev = NaN
    end

    if !isnan(r.prev)
        # We're ratio test, also known as D'Alembert's criterion
        ratio = r.prev == 0.0 ? abs(score / nextfloat(r.prev)) : abs(score / r.prev)

        push!(r.ratios, ratio)
    end
    r.prev = score

    if r.log
        println("[Epoch $epoch] Score: $score [$(epoch_time)s] (Mean Ratio: $(mean(r.ratios)) )")
    end
end

"""
Simply return true or false as to whether the mean ratio is less than 1
"""
converged(r::TestReporter) = mean(r.ratios) < 1.0

"""
Simply returns a DefaultContext Dict.
"""
function DefaultContext()
    return Dict{Any,Any}(
        :weight_decay_kind => :l2,
        :weight_decay_rate => 0.001,
        :sparsity_cost => 0.001,
        :sparsity_target => 0.02,
        :lr => 0.1,
        :momentum => 0.9,
        :batch_size => 100,
        :n_epochs => 5,
        :n_gibbs => 1,
        :reporter => TestReporter()
    )
end

"""
Generates synthetic random datasets with several modifiable properties.

Args:
* T::Type - the type (or precision) of the resulting matrix
* n_features::Int - the number of features in each observation

Optional:
* n_classes::Int - the number of unique classes or categories in the resulting dataset. (default=10)
* n_obs::Int - total number of observations in resulting dataset. (default=1000)
* sparsity::Float64 - specifies the density if a sparse matrix is desired. (default=-1.0)
                      if less than 0.0 a dense matrix is created.
* binary::Bool - whether or not to round the result to 0.0 or 1.0

Returns:
* Mat{T}(n_features, n_obs) - with the various properties specified.
"""
function generate_dataset(T::Type, n_features; n_classes=10, n_obs=1000, sparsity=-1.0, binary=true)
    function rand_fill!(X::Mat, prototypes::Mat)
        for i in 1:size(X, 2)
            j = rand(1:size(prototypes, 2))
            X[:, i] = prototypes[:, j]
        end

        if binary
            return round.(X)
        else
            return X
        end
    end

    if sparsity > 0
        prototypes = map(T, sprand(n_features, n_classes, sparsity))
        X = map(T, spzeros(n_features, n_obs))
        return rand_fill!(X, prototypes)
    else
        prototypes = rand(T, n_features, n_classes)
        X = zeros(n_features, n_obs)
        return rand_fill!(X, prototypes)
    end
end

"""
Runs simple smoke tests on the provided RBM.

Args:
* rbm::AbstractRBM - the rbm to test

Optional:
* opts::Dict - an alternate context to use when calling `fit`. (default=DEFAULT_CONTEXT)
* n_obs::Int - the number of observations to generate for the synthetic datasets. (default=1000)
* debug::Bool - whether or not to print each epoch. (default=false)
                Only applies if the reporter in `ctx` is `TestReporter`

NOTE: Only using dense arrays for the dataset cause the conditional rbm doesn't support
sparse ones yet.
"""
function test(rbm::AbstractRBM{T,V,H}; opts::Dict=DefaultContext(), input_size=-1, n_obs=1000, debug=false) where {T,V,H}
    n_hid, n_vis = size(rbm.W)

    if input_size > 0
        n_vis = input_size
    end

    ctx = deepcopy(opts)
    DX = generate_dataset(T, n_vis; n_obs=n_obs)
    SX = generate_dataset(T, n_vis; n_obs=n_obs, sparsity=0.3)

    if debug && isa(ctx[:reporter], TestReporter)
        ctx[:reporter].log = true
    end

    for X in (DX, SX)
        @info("Testing against $(typeof(X).name)")
        @info("Testing fit")
        fit(rbm, X, ctx)

        @info("Testing for convergence")
        @test converged(ctx[:reporter])

        # Check that none of the array fields contain Infs or NaNs
	for name in fieldnames(typeof(rbm))
            if fieldtype(typeof(rbm), name) <: AbstractArray
                field = getfield(rbm, name)

                @test !any(isnan, field) && !any(isinf, field)
            end
        end

        @info("Testing transform")
        Y = transform(rbm, X)
        @test !any(isnan, Y) && !any(isinf, Y)

        @info("Testing generate")
        Y = generate(rbm, X)
        @test !any(isnan, Y) && !any(isinf, Y)
    end
end


text:
using Test
using MathOptInterface
const MOI = MathOptInterface

using MathOptSetDistances
const MOD = MathOptSetDistances

using FiniteDifferences
using LinearAlgebra

const bfdm = FiniteDifferences.backward_fdm(5, 1)
const ffdm = FiniteDifferences.forward_fdm(5, 1)
const cfdm = FiniteDifferences.central_fdm(5,1)

import ChainRulesCore
const CRC = ChainRulesCore

"""
A multivariate Gaussian generator without points too close to 0
"""
function safe_randn(n)
    v = 2 * randn(n)
    for i in eachindex(v)
        while v[i] ≈ 0
            v[i] = 2 * randn()
        end
    end
    return v
end

@testset "Test gradients with finite differences" begin
    Ntrials = 10
    @testset "Dimension $n" for n in (1, 3, 10)
        vector_sets = (
            MOI.Zeros(n),
            MOI.Reals(n),
            MOI.Nonnegatives(n),
            MOI.Nonpositives(n),
        )
        vs = [safe_randn(n) for _ in 1:Ntrials]
        for _ in 1:Ntrials
            push!(vs, rand(-5:2:5, n))
        end
        @testset "Vector set $s" for s in vector_sets
            @testset "Vector $v" for v in vs
                dΠ = MOD.projection_gradient_on_set(MOD.DefaultDistance(), v, s)
                grad_fdm1 = FiniteDifferences.jacobian(bfdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
                grad_fdm2 = FiniteDifferences.jacobian(ffdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
                @test size(grad_fdm1) == size(grad_fdm2) == size(dΠ)
                # finite diff can mess up for low v
                @test all(eachindex(v)) do idx
                    dΠ[idx,idx] ≈ grad_fdm1[idx,idx] || dΠ[idx,idx] ≈ grad_fdm2[idx,idx]
                end
            end
        end
        @testset "PSD cone" begin
            s = MOI.PositiveSemidefiniteConeTriangle(n)
            for _ in 1:Ntrials
                L = 3 * tril(rand(n, n))
                M = L * L'
                @testset "Positive definite" begin
                    v = MOD.vec_symm(M)
                    dΠ = MOD.projection_gradient_on_set(MOD.DefaultDistance(), v, s)
                    grad_fdm1 = FiniteDifferences.jacobian(ffdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
                    grad_fdm2 = FiniteDifferences.jacobian(bfdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
                    @test size(grad_fdm1) == size(grad_fdm2) == size(dΠ)
                    @test dΠ ≈ I
                end
                @testset "Negative definite" begin
                    v = MOD.vec_symm(-M)
                    dΠ = MOD.projection_gradient_on_set(MOD.DefaultDistance(), v, s)
                    grad_fdm1 = FiniteDifferences.jacobian(ffdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
                    grad_fdm2 = FiniteDifferences.jacobian(bfdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
                    @test size(grad_fdm1) == size(grad_fdm2) == size(dΠ)
                    if !isapprox(det(M), 0, atol=10e-6)
                        @test all(dΠ .≈ 0)
                    end
                end
            end
        end
        @testset "SOC" begin
            s = MOI.SecondOrderCone(n+1)
            for _ in 1:Ntrials
                x = safe_randn(n)
                @testset "SOC interior and negative bound" begin
                    t = LinearAlgebra.norm2(x) + 2 * rand()
                    v = vcat(t, x)
                    dΠ = MOD.projection_gradient_on_set(MOD.DefaultDistance(), v, s)
                    grad_fdm1 = FiniteDifferences.jacobian(ffdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
                    grad_fdm2 = FiniteDifferences.jacobian(bfdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
                    @test size(grad_fdm1) == size(grad_fdm2) == size(dΠ)
                    @test dΠ ≈ grad_fdm1 || dΠ ≈ grad_fdm2
                    v = vcat(-t, x)
                    dΠ = MOD.projection_gradient_on_set(MOD.DefaultDistance(), v, s)
                    grad_fdm1 = FiniteDifferences.jacobian(ffdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
                    grad_fdm2 = FiniteDifferences.jacobian(bfdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
                    @test size(grad_fdm1) == size(grad_fdm2) == size(dΠ)
                    @test dΠ ≈ grad_fdm1 || dΠ ≈ grad_fdm2
                end
                @testset "Out of cone point" begin
                    for tscale in (0.1, 0.5, 0.9)
                        t = tscale * LinearAlgebra.norm2(x)
                        v = vcat(t, x)
                        dΠ = MOD.projection_gradient_on_set(MOD.DefaultDistance(), v, s)
                        grad_fdm1 = FiniteDifferences.jacobian(ffdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
                        grad_fdm2 = FiniteDifferences.jacobian(bfdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
                        @test size(grad_fdm1) == size(grad_fdm2) == size(dΠ)
                        @test dΠ ≈ grad_fdm1 || dΠ ≈ grad_fdm2
                        t = tscale
                        xr = x / norm(x)
                        v = vcat(t, xr)
                        dΠ = MOD.projection_gradient_on_set(MOD.DefaultDistance(), v, s)
                        # theoretical expression for unit vector
                        @test 2dΠ ≈ [
                            1 xr'
                            xr ((t + 1) * I - t * xr * xr')
                        ]
                    end
                end
            end
        end
    end
    @testset "Indefinite matrix" begin
        s = MOI.PositiveSemidefiniteConeTriangle(2)
        Q = [
            1 0
            0 -1
        ]
        Qi = Q'
        B = [
            0 1/2
            1/2 1
        ]
        for _ in 1:5
            # scale factor
            f = 20 * rand() + 5
            A = [
                -f 0
                0 f
            ]
            Λ = Diagonal([-f, f])
            Λp = Diagonal([0, f])
            @test A ≈ Q * Λ * Qi
            v = MOD.vec_symm(A)
            Πv = MOD.projection_on_set(MOD.DefaultDistance(), v, s)
            Π = MOD.unvec_symm(Πv, 2)
            @test Π ≈ Q * Λp * Qi
            DΠ = MOD.projection_gradient_on_set(MOD.DefaultDistance(), v, s)
            # directional derivative
            for _ in 1:Ntrials
                Xd = randn(2,2)
                xd = MOD.vec_symm(Xd)
                @test DΠ * xd ≈ MOD.vec_symm(Q * (B .* (Q' * Xd * Q)) * Q')
            end
            grad_fdm1 = FiniteDifferences.jacobian(ffdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
            grad_fdm2 = FiniteDifferences.jacobian(bfdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
            @test grad_fdm1 ≈ DΠ || grad_fdm2 ≈ DΠ
        end
    end
    @testset "Scalar $ST" for ST in (MOI.LessThan, MOI.GreaterThan, MOI.EqualTo)
        s = ST(10 * randn())
        for v in 1:Ntrials
            v = randn()
            while v ≈ MOI.constant(s)
                v += 2 * randn()
            end
            dfor = ffdm(x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)
            dback = bfdm(x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)
            dΠ = MOD.projection_gradient_on_set(MOD.DefaultDistance(), v, s)
            @test ≈(dfor, dΠ, atol=1e-5) || ≈(dback, dΠ, atol=1e-5)
        end
    end

    @testset "Exp Cone" begin
        function det_case_exp_cone(v; dual=false)
            v = dual ? -v : v
            if MOD._in_exp_cone(v; dual=false)
                return 1
            elseif MOD._in_exp_cone(-v; dual=true)
                return 2
            elseif v[1] <= 0 && v[2] <= 0 #TODO: threshold here??
                return 3
            else
                return 4
            end
        end

        Random.seed!(0)
        s = MOI.ExponentialCone()
        sd = MOI.DualExponentialCone()
        case_p = zeros(4)
        case_d = zeros(4)
        # Adjust tolerance down because a 1-2 errors when projection ends up
        #   very close to the z axis
        # For intuition, see Fig 5.1 https://docs.mosek.com/modeling-cookbook/expo.html
        #   Note that their order is reversed: (x, y, z) = (x3, x2, x1) [theirs]
        tol = 1e-6
        for ii in 1:100
            v = 5*randn(3)
            @testset "Primal Cone" begin
                case_p[det_case_exp_cone(v; dual=false)] += 1
                dΠ = MOD.projection_gradient_on_set(MOD.DefaultDistance(), v, s)
                grad_fdm1 = FiniteDifferences.jacobian(ffdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
                grad_fdm2 = FiniteDifferences.jacobian(bfdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
                @test size(grad_fdm1) == size(grad_fdm2) == size(dΠ)
                @test ≈(dΠ, grad_fdm1,atol=tol) || ≈(dΠ, grad_fdm2, atol=tol)
            end

            @testset "Dual Cone" begin
                case_d[det_case_exp_cone(v; dual=true)] += 1
                dΠ = MOD.projection_gradient_on_set(MOD.DefaultDistance(), v, sd)
                grad_fdm1 = FiniteDifferences.jacobian(ffdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, sd), v)[1]'
                grad_fdm2 = FiniteDifferences.jacobian(bfdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, sd), v)[1]'
                @test size(grad_fdm1) == size(grad_fdm2) == size(dΠ)
                @test ≈(dΠ, grad_fdm1,atol=tol) || ≈(dΠ, grad_fdm2, atol=tol)
            end
        end
        @test all(case_p .> 0) && all(case_d .> 0)
    end

    @testset "Power Cone" begin
        function det_case_pow_cone(x, α; dual=false)
            v = dual ? -x : x
            s = MOI.PowerCone(α)
            if MOD._in_pow_cone(v, s)
                return 1
            elseif MOD._in_pow_cone(-v, MOI.dual_set(s))
                return 2
            elseif abs(v[3]) <= 1e-8
                return 3
            else
                return 4
            end
        end

        case_p = zeros(4)
        case_d = zeros(4)
        Random.seed!(0)
        tol = 1e-5
        for ii in 1:100
            v = 5*randn(3)
            for α in [0.5; rand(0.05:0.05:0.95)]
                if ii % 10 == 1
                    v[3] = 0.0
                end
                s = MOI.PowerCone(α)
                sd = MOI.dual_set(s)
                @testset "Primal Cone" begin
                    case = det_case_pow_cone(v, α; dual=false)
                    case_p[case] += 1
                    dΠ = MOD.projection_gradient_on_set(MOD.DefaultDistance(), v, s)
                    grad_fdm1 = FiniteDifferences.jacobian(ffdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
                    grad_fdm2 = FiniteDifferences.jacobian(bfdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
                    grad_fdm3 = FiniteDifferences.jacobian(cfdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
                    @test size(grad_fdm1) == size(grad_fdm2) == size(dΠ)
                    @test ≈(dΠ, grad_fdm1,atol=tol) || ≈(dΠ, grad_fdm2, atol=tol) || ≈(dΠ, grad_fdm3, atol=tol)
                    if !(≈(dΠ, grad_fdm1,atol=tol) || ≈(dΠ, grad_fdm2, atol=tol) || ≈(dΠ, grad_fdm3, atol=tol))
                        error("v=$v\ndΠ = $dΠ\ncase=$case\nFD=$grad_fdm3")
                    end
                end

                @testset "Dual Cone" begin
                    case = det_case_pow_cone(v, α; dual=true)
                    case_d[case] += 1
                    dΠ = MOD.projection_gradient_on_set(MOD.DefaultDistance(), v, sd)
                    grad_fdm1 = FiniteDifferences.jacobian(ffdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, sd), v)[1]'
                    grad_fdm2 = FiniteDifferences.jacobian(bfdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, sd), v)[1]'
                    grad_fdm3 = FiniteDifferences.jacobian(cfdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, sd), v)[1]'
                    @test size(grad_fdm1) == size(grad_fdm2) == size(dΠ)
                    @test ≈(dΠ, grad_fdm1,atol=tol) || ≈(dΠ, grad_fdm2, atol=tol) || ≈(dΠ, grad_fdm3, atol=tol)
                    if !(≈(dΠ, grad_fdm1,atol=tol) || ≈(dΠ, grad_fdm2, atol=tol) || ≈(dΠ, grad_fdm3, atol=tol))
                        error("v=$v\ndΠ = $dΠ\ncase=$case\nFD=$grad_fdm3")
                    end
                end
            end
        end
        @test all(case_p .> 0) && all(case_d .> 0)
    end
end