text:
using Symata

@sym begin
  ClearAll(f, y0, a, l, x)
  y0(x_) := (e*x)*(l^3 - 2*l*x^2 + x^3)

  # Stodola-Vianello method, e.g. see Den Hartog, Advanced Strength of Materials, pg 268
  f(maxiters_, y0_) := Module([i, yder, res=ConstantArray(0.0,[maxiters,6])],
    (
      # Initialize the first row of the res matrix
      res[1,2] = y0(x),
      For(i=1, i<=4, Increment(i),
        [
          yder(x_) = Simplify(D(y0, x, i))
          res[1, 2+i] = yder(x)
        ]
      ),
      Println("N[res[1, 1]]: ", res[1,1]),
      # Start the loop to improve accuracy stored in subsequent rows,
      # res[i, 1] contains the Euler ratio.
      For(i=2, i<=maxiters, Increment(i),
        begin
          res[i, 6] = Simplify(w/(EI)*((x-l)*res[i-1,4] + res[i-1,3]))
          # Above the basic Euler equation for a pinned-pinned flagpole
          ClearAll(K1, K11, tmpy5, tmpy4, tmp)
          tmpy5(x_) = Simplify(Integrate(res[i,6], x) + K1)
          tmpy4(x_) = Simplify(Integrate(tmpy5(x), x))
          # Use first pinned constraint at x=l to resolve K1
          K11 = Solve(-tmpy4(l), K1)[1]
          res[i,5] = Simplify( tmpy5(x) ./ K11 )
          res[i,4] = Simplify( tmpy4(x) ./ K11 )
          ClearAll(K2, K22, tmpy3, tmpy2, tmp)
          tmpy3(x_) = Simplify(Integrate(res[i,4], x) + K2)
          tmpy2(x_) = Simplify(Integrate(tmpy3(x), x))
          # Use second pinned constraint at x=l to resolve K2
          K22 = Solve(-tmpy2(l), K2)[1]
          res[i,3] = Simplify( tmpy3(x) ./ K22 )
          res[i,2] = Simplify( tmpy2(x) ./ K22 )
          tmp = Simplify( (((l^3)*w)/EI)*(res[i-1,3]/res[i, 3]) )
          res[i,1] = N(Simplify(ReplaceAll(tmp, x => l)))
          Println("N[res[$i, 1]]: ", N(res[i,1], 4))
        end
      ),
      Println(),
      Return(res)
    )
  )

  r = f(8, y0(x))
  Println(N(Transpose(r)[1], 4))
end


text:
using JuMP
using Ipopt

# Define sets
# Streams
S = 2:7

# Components (chemical species)
C = ["A", "B"]

# Units
U = ["U1", "U2"]

# Define stream/unit connectivity using dictionaries
inlets = Dict{String,Integer}("U1"=>2, "U2"=>5)
outletV = Dict{String,Integer}("U1"=>3, "U2"=>6)
outletL = Dict{String,Integer}("U1"=>4, "U2"=>7)

# steams
feeds = [2, 5]

# Equilibrium coefficients

K = Dict{Tuple{String, String}, Float64}(
	("U1","A")=>1.008,		("U1","B")=>0.9,
	("U2","A")=>1.099,		("U2","B")=>0.9 )

# Feed flowrate specification
# (for components, same for all feed streams)
feedflow = Dict{String, Float64}("A"=>0.55, "B"=>0.45)

# Define model
m = Model(solver=IpoptSolver())

##### General flowsheet unit model #####

# Total molar flowrate
@variable(m, f[S] >= 0)

# Component molar flowrate
@variable(m, fc[S,C] >= 0)

# Component mole fraction
@variable(m, 0 <= x[S,C] <= 1)

# Unit model
for u in U

	i = inlets[u]
	v = outletV[u]
	l = outletL[u]

	# Overall mass balance
	@constraint(m, f[i] == f[v] + f[l])

	for c in C
		# Component mass balance
		@constraint(m, fc[i,c] == fc[v,c] + fc[l,c])

		# vapor-liquid equilibrium
		@constraint(m, x[v, c] == K[(u,c)]*x[l,c])

	end

	# Unit summation (Rachford-Rice equation)
	@constraint(m, sum(x[v,c] - x[l,c] for c in C) == 0)

end

# Stream model: mole fraction specification
for s in S
	for c in C
		@constraint(m, fc[s,c] == f[s]*x[s,c])
	end
end

##### Economics and recovery/purity specifications

@variable(m, 0.6 <= purityA <= 1.0)
@variable(m, 0.551 <= recoveryA <= 1.0)

# Equipment cost
ecost = Dict{String,Float64}("U1"=>1.5, "U2"=>1.0)

# Set total feed to one
@constraint(m, sum(f[s] for s in feeds) == 1)

# Calculate purity of A
@constraint(m, purityA*sum(sum(f[s] for s in outletV[u]) for u in U) == sum(sum(fc[s,"A"] for s in outletV[u]) for u in U))

# Calculate recovery of A
@constraint(m, recoveryA == sum(sum(fc[s,"A"] for s in outletV[u]) for u in U) / feedflow["A"])

# Objective
@objective(m, Min, sum(sum(f[s] for s in inlets[u])*ecost[u] for u in U) - 100*purityA)

for s in feeds
	for c in C
		setupperbound(x[s,c], feedflow[c])
		setlowerbound(x[s,c], feedflow[c])
	end
end

##### Initialization

# Guess all molar flowrates are one
# And molar fractions match feed
for s in S

	setvalue(f[s], 1.0)

	for c in C
		setvalue(fc[s,c], feedflow[c])
		setvalue(x[s,c], feedflow[c])
	end
end

print(m)

solve(m)

ext: jl


text:
using JSServe, WGLMakie, AbstractPlotting

function handler(session, request)
    return scatter(1:4, color=1:4)
end

dir = joinpath(@__DIR__, "exported")
isdir(dir) || mkdir(dir)
JSServe.export_standalone(handler, dir)
# Then serve it with e.g. LiveServer
using LiveServer

LiveServer.serve(dir=dir)


text:
## =====================================================================================
# Primary supported functionality
## =====================================================================================

export plotBelief
export plotPose
export plotVariable2D
export plotSLAM2D, plotSLAM2D_KeyAndRef, plotSLAM2DLandmarks, plotSLAM2DPoses
export getColorsByLength

# FIXME what to do about this?
export  plot

# needs better testing (might be broken following upgrade to Manifolds.jl)
export plotPose2Vels
export plotFactor
export plotFactorValues, plotFactorMeasurements, plotVariableGivenFactor
export plotPairVariables, plotPairPose2, plotPose3Pairs
export localProduct, plotLocalProductCylinder, plotLocalProduct
export saveplot
export plotMarginalContour, accumulateMarginalContours

# plot tree functionality (needs testing, possible broken with recent upgrades)
export plotUpMsgsAtCliq
export plotTreeProductUp, plotTreeProductDown, plotCliqUpMsgs, plotCliqDownMsgs

## =====================================================================================
## Known FIXME issues
## =====================================================================================

# export reportFactor
# export plotPriorsAtCliq

## =====================================================================================
## CURRENTLY UNMAINTAINED, utility functions that must be refactored, consolidated, standardized, or deprecated
## =====================================================================================

export  investigateMultidimKDE,
  plotKDEofnc,
  plotKDEresiduals,
  investigateMultidimKDE,
  plotLbl,
  animateVertexBelief,
  # Associated with RoME
  plotTrajectoryArrayPose2,
  saveImgSeq,
  plotProductVsKDE

## =====================================================================================
## Old Victoria Park dataset features that were good but unmaintained
## =====================================================================================

export plotLsrScanFeats, progressExamplePlot, plotTrckStep

#


text: # demoing defining the circuit first, then compiling to fix norms.

max_stars_repo_name: anandijain/MyQuantum

max_issues_repo_name: anandijain/MyQuantum

max_forks_repo_name: anandijain/MyQuantum


text:
import SpecialFunctions: erfc

"""
Ewald term: electrostatic energy per unit cell of the array of point
charges defined by `model.atoms` in a uniform background of
compensating charge yielding net neutrality.
"""
struct Ewald end
(::Ewald)(basis) = TermEwald(basis)

struct TermEwald{T} <: Term
    energy::T  # precomputed energy
end
function TermEwald(basis::PlaneWaveBasis{T}) where {T}
    TermEwald(T(energy_ewald(basis.model)))
end

function ene_ops(term::TermEwald, basis::PlaneWaveBasis, ψ, occ; kwargs...)
    (E=term.energy, ops=[NoopOperator(basis, kpt) for kpt in basis.kpoints])
end

@timing "forces: Ewald" function compute_forces(term::TermEwald, basis::PlaneWaveBasis{T},
                                                ψ, occ; kwargs...) where {T}
    # TODO this could be precomputed
    forces = zero(basis.model.positions)
    energy_ewald(basis.model; forces)
    forces
end

function energy_ewald(model::Model{T}; kwargs...) where {T}
    isempty(model.atoms) && return zero(T)

    # DFTK currently assumes that the compensating charge in the electronic and nuclear
    # terms is equal and of opposite sign. See also the PSP correction term, where
    # n_electrons is used synonymously for sum of charges
    charges = T.(charge_ionic.(model.atoms))
    @assert sum(charges) == model.n_electrons
    energy_ewald(model.lattice, charges, model.positions; kwargs...)
end

"""
Compute the electrostatic interaction energy per unit cell between point
charges in a uniform background of compensating charge to yield net
neutrality. The `lattice` and `recip_lattice` should contain the
lattice and reciprocal lattice vectors as columns. `charges` and
`positions` are the point charges and their positions (as an array of
arrays) in fractional coordinates. If `forces` is not nothing, minus the derivatives
of the energy with respect to `positions` is computed.
"""
function energy_ewald(lattice, charges, positions; η=nothing, forces=nothing)
    T = eltype(lattice)
    for i=1:3
        if iszero(lattice[:, i])
            # TODO should something more clever be done here? For now
            # we assume that we are not interested in the Ewald
            # energy of non-3D systems
            return zero(T)
        end
    end
    energy_ewald(lattice, compute_recip_lattice(lattice), charges, positions; η, forces)
end

# This could be factorised with Pairwise, but its use of `atom_types` would slow down this
# computationally intensive Ewald sums. So we leave it as it for now.
function energy_ewald(lattice, recip_lattice, charges, positions; η=nothing, forces=nothing)
    T = eltype(lattice)
    @assert T == eltype(recip_lattice)
    @assert length(charges) == length(positions)
    if η === nothing
        # Balance between reciprocal summation and real-space summation
        # with a slight bias towards reciprocal summation
        η = sqrt(sqrt(T(1.69) * norm(recip_lattice ./ 2T(π)) / norm(lattice))) / 2
    end
    if forces !== nothing
        @assert size(forces) == size(positions)
        forces_real = copy(forces)
        forces_recip = copy(forces)
    end

    # Numerical cutoffs to obtain meaningful contributions. These are very conservative.
    # The largest argument to the exp(-x) function
    max_exp_arg = -log(eps(T)) + 5  # add some wiggle room
    max_erfc_arg = sqrt(max_exp_arg)  # erfc(x) ~= exp(-x^2)/(sqrt(π)x) for large x

    # Precomputing summation bounds from cutoffs.
    # In the reciprocal-space term we have exp(-||B G||^2 / 4η^2),
    # where B is the reciprocal-space lattice, and
    # thus use the bound  ||B G|| / 2η ≤ sqrt(max_exp_arg)
    Glims = estimate_integer_lattice_bounds(recip_lattice, sqrt(max_exp_arg) * 2η)

    # In the real-space term we have erfc(η ||A(rj - rk - R)||),
    # where A is the real-space lattice, rj and rk are atomic positions and
    # thus use the bound  ||A(rj - rk - R)|| * η ≤ max_erfc_arg
    poslims = [maximum(rj[i] - rk[i] for rj in positions for rk in positions) for i in 1:3]
    Rlims = estimate_integer_lattice_bounds(lattice, max_erfc_arg / η, poslims)

    #
    # Reciprocal space sum
    #
    # Initialize reciprocal sum with correction term for charge neutrality
    sum_recip::T = - (sum(charges)^2 / 4η^2)

    for G1 in -Glims[1]:Glims[1], G2 in -Glims[2]:Glims[2], G3 in -Glims[3]:Glims[3]
        G = Vec3(G1, G2, G3)
        (G == zero(G)) && continue
        Gsq = sum(abs2, recip_lattice * G)
        cos_strucfac = sum(Z * cos2pi(dot(r, G)) for (r, Z) in zip(positions, charges))
        sin_strucfac = sum(Z * sin2pi(dot(r, G)) for (r, Z) in zip(positions, charges))
        sum_strucfac = cos_strucfac^2 + sin_strucfac^2
        sum_recip += sum_strucfac * exp(-Gsq / 4η^2) / Gsq
        if forces !== nothing
            for (ir, r) in enumerate(positions)
                Z = charges[ir]
                dc = -Z*2T(π)*G*sin2pi(dot(r, G))
                ds = +Z*2T(π)*G*cos2pi(dot(r, G))
                dsum = 2cos_strucfac*dc + 2sin_strucfac*ds
                forces_recip[ir] -= dsum * exp(-Gsq / 4η^2)/Gsq
            end
        end
    end

    # Amend sum_recip by proper scaling factors:
    sum_recip *= 4T(π) / compute_unit_cell_volume(lattice)
    if forces !== nothing
        forces_recip .*= 4T(π) / compute_unit_cell_volume(lattice)
    end

    #
    # Real-space sum
    #
    # Initialize real-space sum with correction term for uniform background
    sum_real::T = -2η / sqrt(T(π)) * sum(Z -> Z^2, charges)

    for R1 in -Rlims[1]:Rlims[1], R2 in -Rlims[2]:Rlims[2], R3 in -Rlims[3]:Rlims[3]
        R = Vec3(R1, R2, R3)
        for i = 1:length(positions), j = 1:length(positions)
            # Avoid self-interaction
            R == zero(R) && i == j && continue
            Zi = charges[i]
            Zj = charges[j]
            Δr = lattice * (positions[i] - positions[j] - R)
            dist = norm(Δr)
            energy_contribution = Zi * Zj * erfc(η * dist) / dist
            sum_real += energy_contribution
            if forces !== nothing
                # `dE_ddist` is the derivative of `energy_contribution` w.r.t. `dist`
                dE_ddist = Zi * Zj * η * (-2exp(-(η * dist)^2) / sqrt(T(π)))
                dE_ddist -= energy_contribution
                dE_ddist /= dist
                dE_dti = lattice' * ((dE_ddist / dist) * Δr)
                forces_real[i] -= dE_dti
                forces_real[j] += dE_dti
            end
        end
    end
    energy = (sum_recip + sum_real) / 2  # Divide by 2 (because of double counting)
    if forces !== nothing
        forces .= (forces_recip .+ forces_real) ./ 2
    end
    energy
end


text:
using PyCall
using Conda

function onnxruntime_infer(f, inputs...)
	if "onnxruntime" ∉ Conda._installed_packages() || Conda._installed_packages_dict()["onnxruntime"][1] !== v"1.8.1"
		# TODO: Add some kind of warning if an incompatible default python installation is used due to ENV["PYTHON"] not being set to ""
		Conda.pip_interop(true)
		Conda.pip("install --no-warn-script-location", "onnxruntime==1.8.1")
	end

	modfile = "tmpmodel.onnx"
	try
		save(modfile, f, size.(inputs)...)

		ort = pyimport("onnxruntime")
		sess = ort.InferenceSession(modfile);
		ins = Dict(map(ii -> ii.name, sess.get_inputs()) .=> PyReverseDims.(inputs))
		return JuReverseDims.(Tuple(sess.run(nothing, ins)))
	finally
		rm(modfile; force=true)
	end
end

JuReverseDims(a::PyObject) = JuReverseDims(PyArray(a))
JuReverseDims(a::PyArray) = JuReverseDims(Array(a))
JuReverseDims(a::AbstractArray{T,N}) where {T, N} = permutedims(a, N:-1:1)


text:
@testset "Puzzle is DATORSPLE; User guesses DATORPLES; Guess is not a word" begin
    # Arrange
    dictionary = SwedishDictionary(["DATORSPEL"])
    game = NiancatGame(dictionary, TestingGameService())
    team = Team(1, "defaultteam", "")
    user = User(1, "name", "displayname", team)
    gamecommand(game, user, SetPuzzle("DATORSPLE"))

    # Act
    response = gamecommand(game, user, Guess("DATORPLES"))

    # Assert
    @test response == Incorrect("DATORPLES")
end

@testset "Puzzle is DATORSPLE; User guesses DATORSPEL; Guess is correct" begin
    # Arrange
    dictionary = SwedishDictionary(["DATORSPEL"])
    game = NiancatGame(dictionary, TestingGameService())
    team = Team(1, "defaultteam", "")
    user = User(1, "name", "displayname", team)
    gamecommand(game, user, SetPuzzle("DATORSPLE"))

    # Act
    response = gamecommand(game, user, Guess("DATORSPEL"))

    # Assert
    @test response == Correct("DATORSPEL")
end

@testset "Puzzle is DATORSPLE; User guesses other solution LEDARPOST; Guess is correct" begin
    # Arrange
    dictionary = SwedishDictionary(["DATORSPEL", "LEDARPOST"])
    game = NiancatGame(dictionary, TestingGameService())
    team = Team(1, "defaultteam", "")
    user = User(1, "name", "displayname", team)
    gamecommand(game, user, SetPuzzle("DATORSPLE"))

    # Act
    response = gamecommand(game, user, Guess("LEDARPOST"))

    # Assert
    @test response == Correct("LEDARPOST")
end

@testset "Puzzle is PUSSGRUKA; User guesses PUSSGURKA; Guess is correct" begin
    # Arrange
    dictionary = SwedishDictionary(["PUSSGURKA"])
    game = NiancatGame(dictionary, TestingGameService())
    team = Team(1, "defaultteam", "")
    user = User(1, "name", "displayname", team)
    gamecommand(game, user, SetPuzzle("PUSSGRUKA"))

    # Act
    response = gamecommand(game, user, Guess("PUSSGURKA"))

    # Assert
    @test response == Correct("PUSSGURKA")
end

@testset "Puzzle is PUSSGRUKA; User guesses DATORSPEL; Guess is incorrect" begin
    # Arrange
    dictionary = SwedishDictionary(["DATORSPEL", "PUSSGURKA"])
    game = NiancatGame(dictionary, TestingGameService())
    team = Team(1, "defaultteam", "")
    user = User(1, "name", "displayname", team)
    gamecommand(game, user, SetPuzzle("PUSSGRUKA"))

    # Act
    response = gamecommand(game, user, Guess("DATORSPEL"))

    # Assert
    @test response isa Incorrect
    @test response.word == Word("DATORSPEL")
end

@testset "Puzzle is PUSSGRUKA; User guesses PUSSGURKA; Notification that the user has guessed correctly" begin
    # Arrange
    gameservice = TestingGameService()
    dictionary = SwedishDictionary(["PUSSGURKA"])
    game = NiancatGame(dictionary, gameservice)
    team = Team(1, "defaultteam", "")
    user = User(1, "name", "displayname", team)
    gamecommand(game, user, SetPuzzle("PUSSGRUKA"))

    # Act
    gamecommand(game, user, Guess("PUSSGURKA"))

    # Assert
    @test length(gameservice.notifications) != 0
    @test lastnotification(gameservice) == CorrectNotification(user)
end

@testset "Puzzle is PUSSGRUKA; User guesses PUSSGURKA; Send a game event that the user solved it" begin
    # Arrange
    gameservice = TestingGameService()
    dictionary = SwedishDictionary(["PUSSGURKA"])
    game = NiancatGame(dictionary, gameservice)
    team = Team(1, "defaultteam", "")
    user = User(1, "name", "displayname", team)
    gamecommand(game, user, SetPuzzle("PUSSGRUKA"))

    # Act
    gamecommand(game, user, Guess("PUSSGURKA"))

    # Assert
    @test length(gameservice.userevents) != 0
    event = lastuserevent(gameservice)
    @test event.eventtype == Gameface.UserEvent_Solution
    @test event.user == user
    @test event.round == gameround(game)
    @test event.data == normalizedword(dictionary, "PUSSGURKA")
end


text:
## Spatial
function inversefourier(nx::Int,ny::Int,u::Array{ComplexF64,2})
    umn = zeros(ComplexF64,2*ny-1,2*nx-1)
    uxy = zeros(Float64,2*ny-1,2*nx-1)
    for m1 = 0:1:nx-1
        n1min = m1 == 0 ? 1 : -ny + 1
        for n1 = n1min:1:ny-1
            umn[n1 + ny,m1+nx] = u[n1+ny,m1+1]
            umn[-n1 + ny,-m1+nx] = conj(u[n1+ny,m1+1])
        end
    end
    real(ifft(ifftshift(umn)))
end

function inversefourier(nx::Int,ny::Int,u::Array{Array{ComplexF64,2},1})
    umn = zeros(ComplexF64,2*ny-1,2*nx-1,length(u))
    uxy = zeros(Float64,2*ny-1,2*nx-1,length(u))
    for i in eachindex(u)
        for m1 = 0:1:nx-1
            n1min = m1 == 0 ? 1 : -ny + 1
            for n1 = n1min:1:ny-1
                umn[n1 + ny,m1+nx,i] = u[i][n1+ny,m1+1]
                umn[-n1 + ny,-m1+nx,i] = conj(u[i][n1+ny,m1+1])
            end
        end
        uxy[:,:,i] = real(ifft(ifftshift(umn[:,:,i])))
    end
    uxy
end

function inversefourier(nx::Int,ny::Int,Λ::Int,u::Array{ArrayPartition{Complex{Float64},Tuple{Array{Complex{Float64},2},Array{Complex{Float64},4}}},1})
    umn = zeros(ComplexF64,2*ny-1,2*nx-1,length(u))
    uxy = zeros(Float64,2*ny-1,2*nx-1,length(u))
    for i in eachindex(u)
        for m1 = 0:1:Λ
            n1min = m1 == 0 ? 1 : -ny + 1
            for n1 = n1min:1:ny-1
                umn[n1 + ny,m1+nx,i] = u[i].x[1][n1+ny,m1+1]
                umn[-n1 + ny,-m1+nx,i] = conj(u[i].x[1][n1+ny,m1+1])
            end
        end
        uxy[:,:,i] = real(ifft(ifftshift(umn[:,:,i])))
    end
    uxy
end

## Velocity
function velocity(lx::Float64,ly::Float64,nx::Int,ny::Int,u::Array{Array{ComplexF64,2},1})
    uk = zeros(ComplexF64,2*ny-1,2*nx-1,length(u))
    ux = zeros(Float64,2*ny-1,2*nx-1,length(u))
    for i in eachindex(u)
        for m1 = 0:nx-1
            n1min = m1 == 0 ? 1 : -(ny-1)
            for n1 = n1min:ny-1
                kx = 2.0*Float64(pi)/lx*m1
                ky = 2.0*Float64(pi)/ly*n1
                uk[n1 + ny,m1 + nx,i] = -ky*u[i][n1 + ny,m1 + 1]/(kx^2+ky^2)
                uk[-n1 + ny,-m1 + nx,i] = conj(uk[n1 + ny,m1 + nx,i])
            end
        end
        ux[:,:,i] = real(ifft(ifftshift(uk[:,:,i])))
    end
    ux
end

function zonalvelocity(lx::Float64,ly::Float64,nx::Int,ny::Int,Λ::Int,u::Array{Array{ComplexF64,2},1})
    uk = zeros(ComplexF64,2*ny-1,2*nx-1,length(u))
    ux = zeros(Float64,2*ny-1,2*nx-1,length(u))
    for i in eachindex(u)
        for m1 = 0:Λ
            n1min = m1 == 0 ? 1 : -(ny-1)
            for n1 = n1min:ny-1
                kx = 2.0*Float64(pi)/lx*m1
                ky = 2.0*Float64(pi)/ly*n1
                uk[n1 + ny,m1 + nx,i] = -ky*u[i][n1 + ny,m1 + 1]/(kx^2+ky^2)
                uk[-n1 + ny,-m1 + nx,i] = conj(uk[n1 + ny,m1 + nx,i])
            end
        end
        ux[:,:,i] = real(ifft(ifftshift(uk[:,:,i])))
    end
    ux
end

function zonalvelocity(lx::Float64,ly::Float64,nx::Int,ny::Int,Λ::Int,u::Array{ArrayPartition{Complex{Float64},Tuple{Array{Complex{Float64},2},Array{Complex{Float64},4}}},1})
    uk = zeros(ComplexF64,2*ny-1,2*nx-1,length(u))
    ux = zeros(Float64,2*ny-1,2*nx-1,length(u))
    for i in eachindex(u)
        for m1 = 0:Λ
            n1min = m1 == 0 ? 1 : -(ny-1)
            for n1 = n1min:ny-1
                kx = 2.0*Float64(pi)/lx*m1
                ky = 2.0*Float64(pi)/ly*n1
                uk[n1 + ny,m1 + nx,i] = -ky*u[i].x[1][n1 + ny,m1 + 1]/(kx^2+ky^2)
                uk[-n1 + ny,-m1 + nx,i] = conj(uk[n1 + ny,m1 + nx,i])
            end
        end
        ux[:,:,i] = real(ifft(ifftshift(uk[:,:,i])))
    end
    ux
end

## Energy for NL/GQL
function e_lohi(lx::Float64,ly::Float64,nx::Int,ny::Int,Λ::Int,u::Array{Array{ComplexF64,2},1})
    e_lo = zeros(Float64,length(u))
    e_hi = fill!(similar(e_lo),0)
    for i in eachindex(u)
        for m1 = 0:1:Λ
            n1min = m1 == 0 ? 1 : -(ny-1)
            for n1 = n1min:1:ny-1

                kx = 2.0*Float64(pi)/lx*m1
                ky = 2.0*Float64(pi)/ly*n1
                e_lo[i] += abs(u[i][n1 + ny,m1 + 1])^2/(kx^2 + ky^2)
            end
        end
        for m1 = Λ+1:1:nx-1
            for n1 = -(ny-1):1:ny-1
                kx = 2.0*Float64(pi)/lx*m1
                ky = 2.0*Float64(pi)/ly*n1
                e_hi[i] += abs(u[i][n1 + ny,m1 + 1])^2/(kx^2 + ky^2)
            end
        end
    end
    e_lo,e_hi
end

function e_lohi(lx::Float64,ly::Float64,nx::Int,ny::Int,Λ::Int,u::Array{ArrayPartition{Complex{Float64},Tuple{Array{Complex{Float64},2},Array{Complex{Float64},4}}},1})
    e_lo = zeros(Float64,length(u))
    e_hi = fill!(similar(e_lo),0)
    for i in eachindex(u)
        for m1 = 0:1:Λ
            n1min = m1 == 0 ? 1 : -(ny-1)
            for n1 = n1min:1:ny-1

                kx = 2.0*Float64(pi)/lx*m1
                ky = 2.0*Float64(pi)/ly*n1
                e_lo[i] += abs(u[i].x[1][n1 + ny,m1 + 1])^2/(kx^2 + ky^2)
            end
        end
        for m1 = Λ+1:1:nx-1
            for n1 = -(ny-1):1:ny-1
                kx = 2.0*Float64(pi)/lx*m1
                ky = 2.0*Float64(pi)/ly*n1
                e_hi[i] += abs(u[i].x[2][n1 + ny,m1 - Λ,n1 + ny,m1 - Λ])/(kx^2 + ky^2)
            end
        end
    end
    e_lo,e_hi
end

function fourierenergy(lx::Float64,ly::Float64,nx::Int,ny::Int,u::Array{Array{ComplexF64,2},1})
    E = zeros(Float64,2*ny-1,2*nx-1,length(u))
    for i in eachindex(u)
        for m1 = 0:1:nx-1
            n1min = m1 == 0 ? 1 : -ny + 1
            for n1 = n1min:1:ny-1
                kx = 2.0*Float64(pi)/lx*m1
                ky = 2.0*Float64(pi)/ly*n1
                E[n1 + ny,m1+nx,i] = abs(u[i][n1+ny,m1+1])^2/(kx^2 + ky^2)
                E[-n1 + ny,-m1+nx,i] = E[n1 + ny,m1+nx,i]
            end
        end
    end
    E
end

function fourierenergy(lx::Float64,ly::Float64,nx::Int,ny::Int,Λ::Int,u::Array{ArrayPartition{Complex{Float64},Tuple{Array{Complex{Float64},2},Array{Complex{Float64},4}}},1})
    E = zeros(Float64,2*ny-1,2*nx-1,length(u))
    for i in eachindex(u)
        for m1 = 0:1:Λ
            n1min = m1 == 0 ? 1 : -(ny-1)
            for n1 = n1min:1:ny-1

                kx = 2.0*Float64(pi)/lx*m1
                ky = 2.0*Float64(pi)/ly*n1
                E[n1 + ny,m1+nx,i] += abs(u[i].x[1][n1 + ny,m1 + 1])^2/(kx^2 + ky^2)
                E[-n1 + ny,-m1+nx,i] = E[n1 + ny,m1+nx,i]
            end
        end
        for m1 = Λ+1:1:nx-1
            for n1 = -(ny-1):1:ny-1
                kx = 2.0*Float64(pi)/lx*m1
                ky = 2.0*Float64(pi)/ly*n1
                E[n1 + ny,m1+nx,i] += abs(u[i].x[2][n1 + ny,m1 - Λ,n1 + ny,m1 - Λ])/(kx^2 + ky^2)
                E[-n1 + ny,-m1+nx,i] = E[n1 + ny,m1+nx,i]
            end
        end
    end
    E
end

function energy(lx::Float64,ly::Float64,nx::Int,ny::Int,u::Array{Array{ComplexF64,2},1})
    E = zeros(Float64,length(u))
    Z = fill!(similar(E),0)
    for i in eachindex(u)

        for m1 = 0:nx-1
            n1min = m1 == 0 ? 1 : -(ny-1)
            for n1 = n1min:ny-1

                kx = 2.0*Float64(pi)/lx*m1
                ky = 2.0*Float64(pi)/ly*n1

                E[i] += abs(u[i][n1 + ny,m1 + 1])^2/(kx^2 + ky^2)
                Z[i] += abs(u[i][n1 + ny,m1 + 1])^2

            end
        end
    end
    E,Z
end

# energy for GCE2
function energy(lx::Float64,ly::Float64,nx::Int,ny::Int,Λ::Int,u::Array{ArrayPartition{Complex{Float64},Tuple{Array{Complex{Float64},2},Array{Complex{Float64},4}}},1})

    E = zeros(Float64,length(u))
    Z = fill!(similar(E),0)

    for i in eachindex(u)

        for m1 = 0:1:Λ
            n1min = m1 == 0 ? 1 : -(ny-1)
            for n1 = n1min:1:ny-1

                kx = 2.0*Float64(pi)/lx*m1
                ky = 2.0*Float64(pi)/ly*n1

                E[i] += abs(u[i].x[1][n1 + ny,m1 + 1])^2/(kx^2 + ky^2)
                Z[i] += abs(u[i].x[1][n1 + ny,m1 + 1])^2

            end
        end

        for m1 = Λ+1:1:nx-1
            for n1 = -(ny-1):1:ny-1

                kx = 2.0*Float64(pi)/lx*m1
                ky = 2.0*Float64(pi)/ly*n1

                E[i] += abs(u[i].x[2][n1 + ny,m1 - Λ,n1 + ny,m1 - Λ])/(kx^2 + ky^2)
                Z[i] += abs(u[i].x[2][n1 + ny,m1 - Λ,n1 + ny,m1 - Λ])

            end
        end
    end
    E,Z
end

## modal strength
function modalstrength(lx::Float64,ly::Float64,nx::Int,ny::Int,u::Array{Array{ComplexF64,2},1})

    E = zeros(Float64,length(u),(2*ny-1)*nx)

    for i in eachindex(u)

        for m1 = 0:1:nx-1
            n1min = m1 == 0 ? 1 : -(ny-1)
            for n1 = n1min:1:ny-1

                kx = 2.0*Float64(pi)/lx*m1
                ky = 2.0*Float64(pi)/ly*n1

                E[i,m1*(2*ny-1) + n1+ny] = abs(u[i][n1 + ny,m1 + 1])

            end
        end

    end
    E
end

function modalstrength(lx::Float64,ly::Float64,nx::Int,ny::Int,Λ::Int,u::Array{ArrayPartition{Complex{Float64},Tuple{Array{Complex{Float64},2},Array{Complex{Float64},4}}},1})

    E = zeros(Float64,length(u),(2*ny-1)*nx)

    for i in eachindex(u)

        for m1 = 0:1:Λ

            n1min = m1 == 0 ? 1 : -(ny-1)
            for n1 = n1min:1:ny-1

                kx = 2.0*Float64(pi)/lx*m1
                ky = 2.0*Float64(pi)/ly*n1

                E[i,m1*(2*ny-1) + n1+ny] = abs(u[i].x[1][n1 + ny,m1 + 1])

            end
        end

        for m1 = Λ+1:1:nx-1
            for n1 = -(ny-1):1:ny-1

                kx = 2.0*Float64(pi)/lx*m1
                ky = 2.0*Float64(pi)/ly*n1

                E[i,m1*(2*ny-1) + n1+ny] = sqrt(abs(u[i].x[2][n1 + ny,m1 - Λ,n1 + ny,m1 - Λ]))

            end
        end

    end
    E
end

## zonal energy
function zonalenergy(lx::Float64,ly::Float64,nx::Int,ny::Int,u::Array{Array{ComplexF64,2},1})

    P = zeros(Float64,length(u),nx)
    O = fill!(similar(P),0)

    for i in eachindex(u)

        for m1 = 0:1:nx-1
            n1min = m1 == 0 ? 1 : -(ny-1)
            for n1 = n1min:1:ny-1

                kx = 2.0*Float64(pi)/lx*m1
                ky = 2.0*Float64(pi)/ly*n1

                P[i,m1+1] += abs(u[i][n1 + ny,m1 + 1])^2/(kx^2 + ky^2)
                O[i,m1+1] += abs(u[i][n1 + ny,m1 + 1])^2

            end
        end

    end
    P,O
end

function zonalenergy(lx::Float64,ly::Float64,nx::Int,ny::Int,Λ::Int,u::Array{ArrayPartition{Complex{Float64},Tuple{Array{Complex{Float64},2},Array{Complex{Float64},4}}},1})

    P = zeros(Float64,length(u),nx)
    O = fill!(similar(P),0)

    for i in eachindex(u)

        for m1 = 0:1:Λ

            n1min = m1 == 0 ? 1 : -(ny-1)
            for n1 = n1min:1:ny-1

                kx = 2.0*Float64(pi)/lx*m1
                ky = 2.0*Float64(pi)/ly*n1

                P[i,m1+1] += abs(u[i].x[1][n1 + ny,m1 + 1])^2/(kx^2 + ky^2)
                O[i,m1+1] += abs(u[i].x[1][n1 + ny,m1 + 1])^2

            end
        end

        for m1 = Λ+1:1:nx-1
            for n1 = -(ny-1):1:ny-1

                kx = 2.0*Float64(pi)/lx*m1
                ky = 2.0*Float64(pi)/ly*n1

                P[i,m1+1] += abs(u[i].x[2][n1 + ny,m1 - Λ,n1 + ny,m1 - Λ])/(kx^2 + ky^2)
                O[i,m1+1] += abs(u[i].x[2][n1 + ny,m1 - Λ,n1 + ny,m1 - Λ])

            end
        end

    end
    P,O
end

## mean vorticity NL/GQL
function meanvorticity(lx::Float64,ly::Float64,nx::Int,ny::Int,u::Array{Array{ComplexF64,2},1})

    ζf = zeros(ComplexF64,length(u),2*ny-1)
    ζy = zeros(Float64,length(u),2*ny-1)

    for i in eachindex(u)

        for n1 = 1:1:ny-1

            ζf[i,n1+ny] = u[i][n1+ny,1]
            ζf[i,-n1+ny] = conj(u[i][n1+ny,1])

        end

        ζf[i,ny] = 0.0 + 0.0im
        ζy[i,:] .= real(ifft(ifftshift(ζf[i,:])))

    end

    return ζy

end

# mean vorticity GCE2
function meanvorticity(lx::Float64,ly::Float64,nx::Int,ny::Int,u::Array{ArrayPartition{Complex{Float64},Tuple{Array{Complex{Float64},2},Array{Complex{Float64},4}}},1})

    ζf = zeros(ComplexF64,length(u),2*ny-1)
    ζy = zeros(Float64,length(u),2*ny-1)

    for i in eachindex(u)

        for n1 = 1:1:ny-1

            ζf[i,n1+ny] = u[i].x[1][n1+ny,1]
            ζf[i,-n1+ny] = conj(u[i].x[1][n1+ny,1])

        end

        ζf[i,ny] = 0.0 + 0.0im
        ζy[i,:] .= real(ifft(ifftshift(ζf[i,:])))

    end

    return ζy

end

## time averaged meanvorticity
function meanvorticity(lx::Float64,ly::Float64,nx::Int,ny::Int,t::Array{Float64,1},u::Array{Array{ComplexF64,2},1};t_begin::Float64=100.0)

    ζf = zeros(ComplexF64,length(u),2*ny-1)
    ζy = zeros(Float64,length(u),2*ny-1)

    i_begin = findall(x->x>t_begin,t)[1]

    for i in eachindex(u)

        for n1 = 1:1:ny-1
            if i > i_begin
                T = 0.0
                temp1 = 0.0
                temp2 = 0.0
                for j = i_begin+1:i
                    dt = t[j] - t[j-1]
                    T += dt
                    temp1 += u[j][n1+ny,1]*dt
                    temp2 += conj(u[j][n1+ny,1])*dt
                end
                ζf[i,n1+ny] = temp1/T
                ζf[i,-n1+ny] = temp2/T
            else
                ζf[i,n1+ny] = u[i][n1+ny,1]
                ζf[i,-n1+ny] = conj(u[i][n1+ny,1])
            end
        end

        ζf[i,ny] = 0.0 + 0.0im
        ζy[i,:] .= real(ifft(ifftshift(ζf[i,:])))

    end

    return ζy

end

function meanvorticity(lx::Float64,ly::Float64,nx::Int,ny::Int,t::Array{Float64,1},u::Array{ArrayPartition{Complex{Float64},Tuple{Array{Complex{Float64},2},Array{Complex{Float64},4}}},1};t_begin::Float64=100.0)

    ζf = zeros(ComplexF64,length(u),2*ny-1)
    ζy = zeros(Float64,length(u),2*ny-1)

    i_begin = findall(x->x>t_begin,t)[1]

    for i in eachindex(u)
        for n1 = 1:1:ny-1
            if i > i_begin
                T = 0.0
                temp1 = 0.0
                temp2 = 0.0
                for j = i_begin+1:i
                    dt = t[j] - t[j-1]
                    T += dt
                    temp1 += u[j].x[1][n1+ny,1]*dt
                    temp2 += conj(u[j].x[1][n1+ny,1])*dt
                end
                ζf[i,n1+ny] = temp1/T
                ζf[i,-n1+ny] = temp2/T
            else
                ζf[i,n1+ny] = u[i].x[1][n1+ny,1]
                ζf[i,-n1+ny] = conj(u[i].x[1][n1+ny,1])
            end
        end

        ζf[i,ny] = 0.0 + 0.0im
        ζy[i,:] .= real(ifft(ifftshift(ζf[i,:])))

    end

    return ζy

end


text:
abstract type Constraint end

mutable struct EpipolarConstraint <: Constraint
end

# Triangulation from Two Views Revisited: Hartley-Sturm vs. Optimal Correction
function satisfy(entity::FundamentalMatrix, constraint::EpipolarConstraint, 𝐅::AbstractArray, 𝒟::Tuple{AbstractArray, Vararg{AbstractArray}})
    ℳ, ℳʹ = 𝒟
    𝒪 = similar(ℳ)
    𝒪ʹ = similar(ℳʹ)

    N = length(ℳ)
    𝐏₂ = SMatrix{3,3,Float64,3^2}([1 0 0; 0 1 0; 0 0 0])

    I = 10
    for n = 1:N
        𝐦 = hom(ℳ[n])
        𝐦ʹ = hom(ℳʹ[n])
        𝐦ₕ = init_correction_view_1(𝐅, 𝐦, 𝐦ʹ, 𝐏₂)
        𝐦ₕʹ= init_correction_view_2(𝐅, 𝐦, 𝐦ʹ, 𝐏₂)
        for i = 1:I
            𝐦ₜ =  𝐦 - 𝐦ₕ
            𝐦ₜʹ =  𝐦ʹ - 𝐦ₕʹ
            𝐦ₕ = update_correction_view_1(𝐅, 𝐦, 𝐦ₕ, 𝐦ₜ, 𝐦ʹ, 𝐦ₕʹ, 𝐦ₜʹ, 𝐏₂)
            𝐦ₕʹ = update_correction_view_2(𝐅, 𝐦, 𝐦ₕ, 𝐦ₜ, 𝐦ʹ, 𝐦ₕʹ, 𝐦ₜʹ, 𝐏₂)
        end
        𝒪[n] = hom⁻¹(𝐦ₕ)
        𝒪ʹ[n] = hom⁻¹(𝐦ₕʹ)
    end
    𝒪 ,𝒪ʹ
end

function init_correction_view_1(𝐅::AbstractArray, 𝐦::AbstractVector, 𝐦ʹ::AbstractVector, 𝐏₂::AbstractArray)
    𝐦 -  dot(𝐦,𝐅*𝐦ʹ)*𝐏₂*𝐅*𝐦ʹ / ( dot(𝐅*𝐦ʹ,𝐏₂*𝐅*𝐦ʹ) + dot(𝐅'*𝐦, 𝐏₂*𝐅'*𝐦) )
end

function init_correction_view_2(𝐅::AbstractArray, 𝐦::AbstractVector, 𝐦ʹ::AbstractVector, 𝐏₂::AbstractArray)
    𝐦ʹ -  dot(𝐦,𝐅*𝐦ʹ)*𝐏₂*𝐅'*𝐦 / ( dot(𝐅*𝐦ʹ,𝐏₂*𝐅*𝐦ʹ) + dot(𝐅'*𝐦, 𝐏₂*𝐅'*𝐦) )
end

function update_correction_view_1(𝐅::AbstractArray, 𝐦::AbstractVector, 𝐦ₕ::AbstractVector, 𝐦ₜ::AbstractVector,  𝐦ʹ::AbstractVector, 𝐦ₕʹ::AbstractVector, 𝐦ₜʹ::AbstractVector,  𝐏₂::AbstractArray)
    𝐦 -  ( ( dot(𝐦ₕ,𝐅*𝐦ₕʹ) + dot(𝐅*𝐦ₕʹ, 𝐦ₜ) + dot(𝐅'*𝐦ₕ, 𝐦ₜʹ)  ) * 𝐏₂*𝐅*𝐦ₕʹ)   / (dot(𝐅*𝐦ₕʹ, 𝐏₂*𝐅*𝐦ₕʹ) + dot(𝐅'*𝐦ₕ, 𝐏₂*𝐅'*𝐦ₕ) )
end

function update_correction_view_2(𝐅::AbstractArray, 𝐦::AbstractVector, 𝐦ₕ::AbstractVector, 𝐦ₜ::AbstractVector,  𝐦ʹ::AbstractVector, 𝐦ₕʹ::AbstractVector, 𝐦ₜʹ::AbstractVector,  𝐏₂::AbstractArray)
    𝐦ʹ -  ( ( dot(𝐦ₕ,𝐅*𝐦ₕʹ) + dot(𝐅*𝐦ₕʹ, 𝐦ₜ) + dot(𝐅'*𝐦ₕ, 𝐦ₜʹ)  ) * 𝐏₂*𝐅'*𝐦ₕ)   / (dot(𝐅*𝐦ₕʹ, 𝐏₂*𝐅*𝐦ₕʹ) + dot(𝐅'*𝐦ₕ, 𝐏₂*𝐅'*𝐦ₕ) )
end


text:
using MicroMamba
using Test

@testset "MicroMamba" begin
    @test MicroMamba.available()
    @test isfile(MicroMamba.executable())
    MicroMamba.STATE.version = VersionNumber(0)
    MicroMamba.STATE.executable = ""
    @test MicroMamba.version() isa VersionNumber
    help = read(MicroMamba.cmd(`--help`), String)
    @test occursin("--help", help)
end


text:
function countUniqueValues_test()::Bool
  a=[1,1,1,2,2,1]*1.
  (uCounts,uValues)=MetoceanTools.countUniqueValues(a)
  if uCounts[1]==4
    return true
  else
    return false
  end
end

max_stars_repo_name: jadamc/MetoceanTools

max_issues_repo_name: jadamc/MetoceanTools

max_forks_repo_name: jadamc/MetoceanTools


text:
# Coordinate descent method, translated from the Python/Cython implementation
#  in scikit-learn and modified to comply with the interfaces of the NMF package

# Original files
# https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/decomposition/nmf.py
# https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/decomposition/cdnmf_fast.pyx

# Original implementation authors:
# Vlad Niculae
# Lars Buitinck
# Mathieu Blondel <mathieu@mblondel.org>
# Tom Dupre la Tour

# Original license: BSD 3 clause

# Julia translation: Vilim Štih

# Reference: Cichocki, Andrzej, and P. H. A. N. Anh-Huy. "Fast local algorithms for
#  large scale nonnegative matrix and tensor factorizations."
#  IEICE transactions on fundamentals of electronics, communications and
#  computer sciences 92.3: 708-721, 2009.

mutable struct CoordinateDescent{T}
    maxiter::Int           # maximum number of iterations (in main procedure)
    verbose::Bool          # whether to show procedural information
    tol::T                 # tolerance of changes on W and H upon convergence
    update_H::Bool         # whether to update H
    α::T                   # constant that multiplies the regularization terms
    l₁ratio::T             # select whether the regularization affects the components (H),
                           # the transformation (W), both or none of them
                           # (:components, :transformation, :both, :none)
    regularization::Symbol # l1 / l2 regularization mixing parameter (in [0; 1])
    shuffle::Bool          # # if true, randomize the order of coordinates in the CD solver

    function CoordinateDescent{T}(;maxiter::Integer=100,
                              verbose::Bool=false,
                              tol::Real=cbrt(eps(T)),
                              update_H::Bool=true,
                              α::Real=zero(T),
                              regularization=:both,
                              l₁ratio::Real=zero(T),
                              shuffle::Bool=false) where T
        new{T}(maxiter, verbose, tol, update_H, α, l₁ratio, regularization, shuffle)
    end
end

solve!(alg::CoordinateDescent{T}, X, W, H) where {T} =
    nmf_skeleton!(CoordinateDescentUpd{T}(alg.α, alg.l₁ratio, alg.regularization, alg.shuffle, alg.update_H),
                  X, W, H, alg.maxiter, alg.verbose, alg.tol)

struct CoordinateDescentUpd{T} <: NMFUpdater{T}
    l₁W::T
    l₂W::T
    l₁H::T
    l₂H::T
    shuffle::Bool
    update_H::Bool
    function CoordinateDescentUpd{T}(α::T, l₁ratio::T, regularization::Symbol, shuffle::Bool, update_H::Bool) where {T}
        αW = zero(T)
        αH = zero(T)

        if (regularization == :both) || (regularization == :components)
            αH = α
        end

        if (regularization == :both) || (regularization == :transformation)
            αW = α
        end

        new{T}(αW*l₁ratio,
               αW*(1-l₁ratio),
               αH*l₁ratio,
               αH*(1-l₁ratio),
               shuffle,
               update_H)
    end
end

mutable struct CoordinateDescentState{T}
    WH::Matrix{T}
    HHt::Matrix{T}
    XHt::Matrix{T}
    XtW::Matrix{T}
    violation::T
    violation_init::Union{Nothing, T}

    function CoordinateDescentState{T}(X, W, H, violation, violation_init) where T
        p, n, k = nmf_checksize(X, W, H)
        new{T}(W * H,
               Matrix{T}(undef, k, k),
               Matrix{T}(undef, p, k),
               Matrix{T}(undef, n, k),
               violation,
               violation_init)
    end
end

prepare_state(::CoordinateDescentUpd{T}, X, W, H) where T = CoordinateDescentState{T}(X, W, H, zero(T), nothing)

function evaluate_objv(::CoordinateDescentUpd{T}, s::CoordinateDescentState{T}, X, W, H) where T
    mul!(s.WH, W, H)
    convert(T, 0.5) * sqL2dist(X, s.WH)
end

"Updates W only"
function _update_coord_descent!(s::CoordinateDescentState{T}, X, W, H,
                                l1_reg, l2_reg, shuffle::Bool, W_flag::Bool) where T
    Ht = transpose(H)
    HHt = s.HHt
    mul!(HHt, H, Ht)
    if W_flag
        XHt = s.XHt
    else
        XHt = s.XtW
    end
    mul!(XHt, X, Ht)

    n_components = size(H, 1)
    n_samples = size(W, 1)

    if l2_reg > 0.
        HHt[diagind(HHt)] .+= l2_reg
    end
    if l1_reg > 0.
        XHt .-= l1_reg
    end
    if shuffle
        permutation = randperm(n_components)
    else
        permutation = 1:n_components
    end

    violation = zero(eltype(X))

    for t in permutation
        for i in 1:n_samples
             # gradient = GW[t, i] where GW = np.dot(W, HHt) - XHt
            grad = -XHt[i, t]

            for r in 1:n_components
                grad += HHt[t, r] * W[i, r]
            end

            # projected gradient
            pg = W[i, t] == 0 ? min(zero(grad), grad) : grad
            violation += abs(pg)

            # Hessian
            hess = HHt[t, t]
            if hess != 0
                W[i, t] = max(W[i, t] - grad / hess, zero(grad))
            end
        end
    end
    return violation
end

function update_wh!(upd::CoordinateDescentUpd{T}, s::CoordinateDescentState{T},
                    X::AbstractArray{T}, W::AbstractArray{T}, H::AbstractArray{T}) where T
    violation = zero(T)

    # update W
    violation += _update_coord_descent!(s, X, W, H, upd.l₁W, upd.l₂W, upd.shuffle, true)

    # update H
    if upd.update_H
        Wt = transpose(W)
        Ht = transpose(H)
        Xt = transpose(X)
        violation += _update_coord_descent!(s, Xt, Ht, Wt, upd.l₁H, upd.l₂H, upd.shuffle, false)
    end

    s.violation = violation
    if s.violation_init !== nothing
        s.violation_init = violation
    end
end


text:
using Printf
using Test: AbstractTestSet, Result, Pass, Fail, Error, Broken,
    get_testset_depth, get_testset

import Test: record, finish

const WIDTH = 70
const column = Ref(0)

mutable struct TodyTestSet <: AbstractTestSet
    desc::AbstractString
    passes::Vector
    fails::Vector
    errors::Vector
    starttime::Float64
    opts
end
TodyTestSet(desc::AbstractString; opts...) =
    TodyTestSet(desc, [], [], [], time_ns(), opts)

function log(ts::TodyTestSet, str::AbstractString; err=false, fail=false)
    w = textwidth(str)
    column[] = column[] + w
    column[] > WIDTH && (println(); column[] = w)
    (err || fail) ? printstyled(str, color=:red, bold=true) : print(str)
end

record(ts::TodyTestSet, res::Pass) = (log(ts, "."); push!(ts.passes, res))

record(ts::TodyTestSet, res::Fail) =
    (log(ts, "F", fail=true); push!(ts.fails, res))

record(ts::TodyTestSet, res::Error) =
    (log(ts, "E", err=true); push!(ts.errors, res))

record(ts::TodyTestSet, res::Broken) =
    (log(ts, "B"); push!(ts.passes, res))

function record(ts::TodyTestSet, kid::TodyTestSet)
    append!(ts.passes, kid.passes)
    append!(ts.fails, kid.fails)
    append!(ts.errors, kid.errors)
end

function finish(ts::TodyTestSet)
    sep() = println("-"^WIDTH)
    if get_testset_depth() == 0
        column[] = 0
        println()
        for err in ts.errors
            sep()
            println(err)
        end
        for fail in ts.fails
            sep()
            println(fail)
        end
        sep()
        n = length(ts.passes) + length(ts.fails) + length(ts.errors)
        s = (time_ns() - ts.starttime)/1e9
        @printf("Ran %i tests in %.3f s\n\n", n, s)

        nfails = length(ts.fails)
        nerrors = length(ts.errors)

        if nfails == nerrors == 0
            println("OK")
        else
            printstyled("FAILED", color=:red, bold=true)
            print(" (")

            problems = []
            nfails > 0 && push!(problems, "failures=$nfails")
            nerrors > 0 && push!(problems, "errors=$nerrors")
            print(join(problems, ", "))

            println(")")
        end
    else
        record(get_testset(), ts)
    end
end


text:
using Test

# Symmetric monoidal category
#############################

A, B = Ob(FreeSymmetricMonoidalCategoryAdditive, :A, :B)
f, g = Hom(:f, A, B), Hom(:g, B, A)

# Domains and codomains
@test dom(oplus(f,g)) == oplus(dom(f),dom(g))
@test codom(oplus(f,g)) == oplus(codom(f),codom(g))
@test dom(swap(A,B)) == oplus(A,B)
@test codom(swap(A,B)) == oplus(B,A)

# Associativity and unit
O = mzero(FreeSymmetricMonoidalCategoryAdditive.Ob)
@test oplus(A,O) == A
@test oplus(O,A) == A
@test oplus(oplus(A,B),A) == oplus(A,oplus(B,A))
@test oplus(oplus(f,g),f) == oplus(f,oplus(g,f))

# Extra syntax
@test oplus(A,B,A) == oplus(oplus(A,B),A)
@test oplus([A,B,A]) == oplus(oplus(A,B),A)
@test oplus(f,f,f) == oplus(oplus(f,f),f)
@test oplus([f,f,f]) == oplus(oplus(f,f),f)
@test oplus(FreeSymmetricMonoidalCategoryAdditive.Ob[]) == O
@test_throws MethodError oplus([])
@test A⊕B == oplus(A,B)
@test f⊕g == oplus(f,g)
@test A⊕B⊕A == oplus(A,B,A)

# Extra functions
@test collect(A) == [A]
@test collect(oplus(A,B)) == [A,B]
@test collect(O) == []
@test typeof(collect(O)) == Vector{FreeSymmetricMonoidalCategoryAdditive.Ob}
@test ndims(A) == 1
@test ndims(oplus(A,B)) == 2
@test ndims(O) == 0

# String format
@test string(oplus(f,g)) == "oplus(f,g)"
@test string(compose(oplus(f,f),oplus(g,g))) == "compose(oplus(f,f),oplus(g,g))"

# S-expressions
@test sexpr(O) == "(mzero)"
@test sexpr(oplus(A,B)) == "(oplus :A :B)"
@test sexpr(oplus(f,g)) == "(oplus :f :g)"
@test sexpr(compose(oplus(f,f),oplus(g,g))) == "(compose (oplus :f :f) (oplus :g :g))"

# Infix notation (Unicode)
@test unicode(O) == "O"
@test unicode(oplus(A,B)) == "A⊕B"
@test unicode(oplus(f,g)) == "f⊕g"
@test unicode(compose(oplus(f,f),oplus(g,g))) == "(f⊕f)⋅(g⊕g)"
@test unicode(oplus(compose(f,g),compose(g,f))) == "(f⋅g)⊕(g⋅f)"

# Infix notation (LaTeX)
@test latex(O) == "O"
@test latex(oplus(A,B)) == "A \\oplus B"
@test latex(oplus(f,g)) == "f \\oplus g"
@test latex(compose(oplus(f,f),oplus(g,g))) ==
  "\\left(f \\oplus f\\right) \\cdot \\left(g \\oplus g\\right)"
@test latex(oplus(compose(f,g),compose(g,f))) ==
  "\\left(f \\cdot g\\right) \\oplus \\left(g \\cdot f\\right)"
@test latex(swap(A,B)) == "\\sigma_{A,B}"

# Cocartesian category
######################

A, B = Ob(FreeCocartesianCategory, :A, :B)
f, g = Hom(:f, A, B), Hom(:g, B, A)

# Domains and codomains
@test dom(plus(A)) == oplus(A,A)
@test codom(plus(A)) == A
@test dom(zero(A)) == O
@test codom(zero(A)) == A

# Derived syntax
@test copair(f,f) == compose(oplus(f,f), plus(B))
@test coproj1(A,B) == oplus(id(A), zero(B))
@test coproj2(A,B) == oplus(zero(A), id(B))

# LaTeX notation
@test latex(plus(A)) == "\\nabla_{A}"
@test latex(zero(A)) == "0_{A}"


text:
println("\nRunLargeTest: Demonstrates arrays of components and timing.")

using Modia
using Modia.Electric

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@testset "RunLargeTests" begin

    @model LPfilter begin
        R = Resistor(R=100.0)
        C = Capacitor(C=2.5E-3, v=Float(start=0.0))
        V = ConstantVoltage(V=10.0)
        ground = Ground()
        @equations begin
            connect(V.n, ground.p)
            connect(V.p, R.p)
            connect(R.n, C.p)
            connect(C.n, V.n)
        end
    end

    nFilters = 100
    @model ManyFilters begin
        F = [LPfilter() for i in 1:nFilters]
    end

    @time simulate(ManyFilters, 2, logTiming=true, storeEliminated=false, removeSingularities=false)

end # testset
nothing


text:
using Pathways, Base.Iterators, Test, Statistics, Printf

const GROUNDTRUTH = joinpath(@__DIR__, "stringpa")
const DISPLAYOFF = "--DisplayOff"
const TMPFILE = joinpath(@__DIR__, "test.txt")
const PATHFILE = joinpath(@__DIR__, "test_pathway.txt")

const NUM_PATHS = r"Number of pathways: (\d+)"
const ASSEMBLY_INDEX = r"Assembly Index = (\d+)"

function stringpa(str::AbstractString)
    expected = -1

    open(TMPFILE, "w") do io
        println(io, str)
    end

    run(`$GROUNDTRUTH $TMPFILE $DISPLAYOFF`)

    expected = -1
    for line in readlines(PATHFILE)
        m = match(NUM_PATHS, line)
        if m !== nothing && parse(Int, m[1]) == 0
            expected = length(str) - 1;
            break
        end

        m = match(ASSEMBLY_INDEX, line)
        if m !== nothing
            expected = parse(Int, m[1])
        end
    end

    expected
end

outer(f, iter, n) = map(f ∘ collect, product(fill(iter, n)...))
outer(iter, n) = map(collect, product(fill(iter, n)...))

@time @testset "Assembly" begin
    include("assembly.jl")
end

@time @testset "Shortest Chain" begin
    include("thurber.jl")
end


text:
# Objects for keeping track of indices in models.
module VariationalModelIndices

VERSION < v"0.4.0-dev" && using Docile

export make_ud_index_matrix, linearize_matrix, linearized_matrix_size
export IndexDict, ModelIndices
export get_transformed_params, get_moment_params

@doc """
Return a matrix of indices that can be used to index into a linearization
of the upper diagonal of a matrix.
""" ->
function make_ud_index_matrix(k::Int64)
	ud_mat = Array(Int64, (k, k))
	for k1=1:k, k2=1:k
		ud_mat[k1, k2] =
			(k1 <= k2 ? (k1 + (k2 - 1) * k2 / 2) :
				        (k2 + (k1 - 1) * k1 / 2))
	end
	ud_mat
end

@doc """
Convert the upper diagonal of a matrix to a vector using the indieces of ud_ind.
""" ->
function linearize_matrix(mat::Matrix{Float64}, ud_ind::Matrix{Int64})
	k_ud = maximum(ud_ind)
	k_tot = size(ud_ind, 1)

	mat_lin = Array(Float64, k_ud)
	for k1=1:k_tot, k2=1:k1
		mat_lin[ud_ind[k1, k2]] = mat[k1, k2]
	end
	mat_lin
end

@doc """
Get a matrix size from the number of upper triangular elements.
""" ->
function linearized_matrix_size(k_ud::Int64)
	k_tot = 0.5 * (sqrt(1 + 8 * k_ud) - 1)
	k_tot_int = round(Int, k_tot)
	@assert round(Int, k_tot) == k_tot
	k_tot_int
end

# A dictionary from (symbol, extra information) to a vector of indices.
typealias IndexDict Dict{Any, Vector{Int64}}

@doc """
Index dictionaries that characterize a VB solution.

model_params: The parameters in the optimization problem
moment_params: The moment parameters of the variational solution
model_const: Constant parameters (e.g. prior parameters) in the model
moment_const: Constant parameters (e.g. prior parameters) in the LRVB matrix
trans_input: Inputs of transformed parameters (must be in model_params)
trans_output: Outputs of transformed parameteres (must be in moment_params)
""" ->
type ModelIndices
	# Indices into the objective (JuMP) model
	model_params::IndexDict

	# Indices into the vector of moment parameters.
	moment_params::IndexDict

	# Indices for constants (e.g. prior parameters) in the model
	# and moment matrix, respectively.
	model_const::IndexDict
	moment_const::IndexDict

	# Model parameters that are transformed into moment parameters.
	trans_input::IndexDict
	trans_output::IndexDict

	num_moment_indices::Int64
	num_const_indices::Int64
	num_transformed_indices::Int64
end

ModelIndices(model_params::IndexDict, moment_params::IndexDict,
             model_const::IndexDict, moment_const::IndexDict,
             trans_input::IndexDict, trans_output::IndexDict) = begin

  for key in keys(trans_input)
    @assert haskey(model_params, key)  "Missing key $(key) from model_params"
  end

  for key in keys(trans_output)
    @assert haskey(moment_params, key) "Missing key $(key) from moment_params"
  end

  num_moment_indices = maximum([ maximum(val) for val in values(moment_params) ])
  num_const_indices = maximum([ maximum(val) for val in values(moment_const) ])
  num_transformed_indices =
    sum(Int64[ length(indices) for (param, indices) in trans_output ])

  ModelIndices(model_params, moment_params,
	             model_const, moment_const, trans_input, trans_output,
               num_moment_indices, num_const_indices, num_transformed_indices)
end

function get_transformed_params(model_params::Vector{Float64},
                                indices::ModelIndices, transform::Function)
  # Get the transformed parameter vector.
  trans_input_params = zeros(Float64, indices.num_transformed_indices)
  for trans_param in keys(indices.trans_input)
  	trans_input_params[indices.trans_input[trans_param]] =
  		model_params[indices.model_params[trans_param]]
  end
  transform(trans_input_params)
end

function get_moment_params(model_params::Vector{Float64},
                           indices::ModelIndices, transform::Function)

	transformed_params = keys(indices.trans_input);
	untransformed_params = setdiff(keys(indices.model_params), transformed_params);

  trans_moment_params = get_transformed_params(model_params, indices, transform)

	moment_params = zeros(Float64, indices.num_moment_indices);
	for param0 in untransformed_params
		param0_indices = indices.model_params[param0]
		moment0_indices = indices.moment_params[param0]
		moment_params[moment0_indices] = model_params[param0_indices]
	end
	for param1 in keys(indices.trans_output)
		param1_indices = indices.trans_output[param1]
		moment1_indices = indices.moment_params[param1]
		moment_params[moment1_indices] = trans_moment_params[param1_indices]
	end

	moment_params
end

end # Module


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -32.35079088597315 <= q <= 33.04707617725232)

                     add_NL_constraint(m, :(softplus(0.7516814527569986 + 0.07211062597860263*softplus(0.800956367833014 + 0.1575422515944389*softplus(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.30181017969165413*softplus(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.5024447348308962*softplus(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + -0.7463877606567912*softplus(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.28099609060437203*softplus(-0.8241959750660866 + 0.48752067922981546*softplus(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + -0.8942226763869914*softplus(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + -0.1882956347370066*softplus(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.0613423792438379*softplus(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + 0.19401527858092393*softplus(-0.014178711859067938 + -0.6496362145686132*softplus(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5744456236018198*softplus(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.06990002879133872*softplus(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.11738131892224501*softplus(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.43991506890783594*softplus(-0.38559964579408224 + -0.6882459841770032*softplus(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5869261088566509*softplus(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.7286086900917508*softplus(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.9385542676650069*softplus(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4])))) + softplus(0.5442086655639207 + -0.859033717055464*softplus(0.800956367833014 + 0.1575422515944389*softplus(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.30181017969165413*softplus(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.5024447348308962*softplus(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + -0.7463877606567912*softplus(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.6994107397450926*softplus(-0.8241959750660866 + 0.48752067922981546*softplus(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + -0.8942226763869914*softplus(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + -0.1882956347370066*softplus(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.0613423792438379*softplus(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.21447765318835765*softplus(-0.014178711859067938 + -0.6496362145686132*softplus(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5744456236018198*softplus(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.06990002879133872*softplus(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.11738131892224501*softplus(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + 0.9247182813449877*softplus(-0.38559964579408224 + -0.6882459841770032*softplus(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5869261088566509*softplus(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.7286086900917508*softplus(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.9385542676650069*softplus(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4])))) + softplus(-0.33457802196328057 + 0.8511359109817849*softplus(0.800956367833014 + 0.1575422515944389*softplus(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.30181017969165413*softplus(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.5024447348308962*softplus(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + -0.7463877606567912*softplus(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + 0.9056512987027685*softplus(-0.8241959750660866 + 0.48752067922981546*softplus(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + -0.8942226763869914*softplus(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + -0.1882956347370066*softplus(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.0613423792438379*softplus(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.901983268277478*softplus(-0.014178711859067938 + -0.6496362145686132*softplus(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5744456236018198*softplus(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.06990002879133872*softplus(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.11738131892224501*softplus(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.02854306019224717*softplus(-0.38559964579408224 + -0.6882459841770032*softplus(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5869261088566509*softplus(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.7286086900917508*softplus(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.9385542676650069*softplus(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4])))) + softplus(0.85223750045835 + -0.7468516096711739*softplus(0.800956367833014 + 0.1575422515944389*softplus(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.30181017969165413*softplus(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.5024447348308962*softplus(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + -0.7463877606567912*softplus(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.1269236474608202*softplus(-0.8241959750660866 + 0.48752067922981546*softplus(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + -0.8942226763869914*softplus(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + -0.1882956347370066*softplus(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.0613423792438379*softplus(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.5407492362268407*softplus(-0.014178711859067938 + -0.6496362145686132*softplus(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5744456236018198*softplus(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.06990002879133872*softplus(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.11738131892224501*softplus(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.8670343803826994*softplus(-0.38559964579408224 + -0.6882459841770032*softplus(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5869261088566509*softplus(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.7286086900917508*softplus(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.9385542676650069*softplus(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Test
using PolynomialBasis
using Roots
using IntervalArithmetic
# using Revise
using ImplicitDomainQuadrature

IDQ = ImplicitDomainQuadrature

function allapprox(v1, v2; tol = 1e-14)
    @assert length(v1) == length(v2)
    flags = [isapprox(v1[i], v2[i], atol = tol) for i = 1:length(v1)]
    return all(flags)
end

f(x) = (x - 1.0) * (x - 2.0) * (x - 3.0)
g(x) = (x - 0.1) * (x - 0.2) * (x - 2.5)
@test IDQ.sign_conditions_satisfied([f, g], 0.5, [-1, -1])
@test IDQ.sign_conditions_satisfied([f, g], 1.5, [+1, -1])

quad1d = IDQ.ReferenceQuadratureRule(3)
@test_throws AssertionError IDQ.one_dimensional_quadrature(
    [f, g],
    [+1, -1, +1],
    0.0,
    3.5,
    quad1d,
)

quad = IDQ.one_dimensional_quadrature([f, g], [+1, -1], 0.0, 3.5, quad1d)
p, w = IDQ.transform(quad1d, 1.0, 2.0)
@test allapprox(quad.points, p)
@test allapprox(quad.weights, w)

quad =
    IDQ.one_dimensional_quadrature([f, g], [+1, +1], Interval(0.0, 3.5), quad1d)
p, w = IDQ.transform(quad1d, 3.0, 3.5)
@test allapprox(quad.points, p)
@test allapprox(quad.weights, w)

max_stars_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_issues_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_forks_repo_name: ArjunNarayanan/ImplicitDomainQuadrature


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
using BenchmarkTools

using SHA
using Blake3Hash

const SUITE     = BenchmarkGroup()
const DATA_SIZE = 2*8192

SUITE["blake3"] = BenchmarkGroup()

SUITE["blake3"]["RefImpl"] = @benchmarkable begin
    Blake3Hash.RefImpl.update!(h, d)
    Blake3Hash.RefImpl.digest(h)
end setup = begin
    h = Blake3Hash.RefImpl.Blake3Ctx()
    d = rand(UInt8, DATA_SIZE)
end

SUITE["blake3"]["SAImpl"] = @benchmarkable begin
    Blake3Hash.SAImpl.update!(h, d)
    Blake3Hash.SAImpl.digest(h)
end setup = begin
    h = Blake3Hash.SAImpl.Blake3Ctx()
    d = rand(UInt8, DATA_SIZE)
end

SUITE["blake3"]["Simd128"] = @benchmarkable begin
    Blake3Hash.Simd128.update!(h, d)
    Blake3Hash.Simd128.digest(h)
end setup = begin
    h = Blake3Hash.Simd128.Blake3Ctx()
    d = rand(UInt8, DATA_SIZE)
end

SUITE["ref"] = BenchmarkGroup()
SUITE["ref"]["sha2_256"] = @benchmarkable begin
    SHA.update!(c, d)
    SHA.digest!(c)
end setup = begin
    c = SHA2_256_CTX()
    d = rand(UInt8, DATA_SIZE)
end

SUITE["ref"]["sha3_256"] = @benchmarkable begin
    SHA.update!(c, d)
    SHA.digest!(c)
end setup = begin
    c = SHA3_256_CTX()
    d = rand(UInt8, DATA_SIZE)
end


text:
export
    CrossEntropyMethod

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

# TODO: add this to MLPlots??

# export
#     AnimationStrategy

# # add this to your MasterLearner to save Plots animations of your learning process
# type AnimationStrategy <: LearningStrategy
#     anim::Animation
#     f::Function
# end
# AnimationStrategy(f::Function) = AnimationStrategy(Animation(), f)
# hook(strat::AnimationStrategy, policy, i) = (strat.f(policy, i); frame(strat.anim))
# cleanup!(strat::AnimationStrategy, policy) = gif(strat.anim)

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

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

# TODO: this is a strategy and policy combined... they should really be split into a Learnable and a LearningStrategy
mutable struct CrossEntropyMethod <: LearningStrategy
    noise_func::Function    # additional deviation at each timestep
    maxsteps::Int           # max num steps in one episode
    cem_batchsize::Int
    cem_elitefrac::Float64
    stopping_norm::Float64
    t::Int
    n::Int
    μ::Vector{Float64}
    last_μ::Vector{Float64}
    σ::Vector{Float64}
    Z::Vector{Float64}      # extra variance
    CrossEntropyMethod(nf, iter, bs, ef, sn) = new(nf, iter, bs, ef, sn, 0)
end
function CrossEntropyMethod(; #f::Function;
           noise_func = t->0.0,
           maxsteps::Int = 100,
           cem_batchsize::Int = 20,
           cem_elitefrac::Float64 = 0.2,
           stopping_norm::Float64 = 1e-2)
    CrossEntropyMethod(noise_func, maxsteps, cem_batchsize, cem_elitefrac, stopping_norm)
end

function setup!(strat::CrossEntropyMethod, policy)
    n = length(params(policy))
    strat.n = n
    strat.μ = zeros(n)
    strat.last_μ = zeros(n)
    strat.σ = ones(n)
    strat.Z = zeros(n)
    return
end

# # the core loop, act/step in a simulation, update the env state, get a reward, update a policy, etc
function learn!(policy::AbstractPolicy, strat::CrossEntropyMethod, env::AbstractEnvironment)
    strat.last_μ = copy(strat.μ)
    strat.t += 1

    # sample thetas from a multivariate normal distribution
    N = MultivariateNormal(strat.μ, strat.σ)
    θs = [rand(N) for k=1:strat.cem_batchsize]

    # overwrite the parameters of the policy and run an episode for each θ
    Rs = map(θ -> begin
        params(policy)[:] = θ
        R = 0
        for (i,sars) in enumerate(Episode(env,policy))
            R += sars[3]
            i < strat.maxsteps || break
        end
        # R, T = episode!(env, policy; maxsteps = strat.maxsteps)
        R
    end, θs)

    # pick out the elite set
    n_elite = round(Int, strat.cem_batchsize * strat.cem_elitefrac)
    elite_indices = sortperm(Rs, rev=true)[1:n_elite]
    elite_θs = θs[elite_indices]
    info("Iteration $(strat.t). mean(R): $(mean(Rs)) max(R): $(maximum(Rs)) ‖μ‖²: $(norm(strat.μ)) ‖σ‖²: $(norm(strat.σ))")

    # update the policy from the empirical statistics of the elite set
    for j=1:length(strat.μ)
        θj = [θ[j] for θ in elite_θs]
        strat.μ[j] = mean(θj)
        strat.Z[j] = strat.noise_func(strat.t)
        strat.σ[j] = sqrt(var(θj) + strat.Z[j])
    end
    # @show strat.μ strat.σ strat.Z
end

# are we done iterating?  check for convergence, etc
function finished(strat::CrossEntropyMethod, policy::AbstractPolicy, i::Int)
    strat.t > 0 || return false
    normdiff = norm(strat.μ - strat.last_μ)
    if normdiff < strat.stopping_norm
        info("Converged after $(strat.t * strat.cem_batchsize) episodes.")
        return true
    end
    false
end


text:
function correctness_check!(coeffaccum, coeffbuffer, primetracker, meta,
                                ring, exps, coeffs, coeffs_zz, gens_temp_ff, gb_ff, ht)

    # first we check coefficients only
    if meta.heuristiccheck
        if !heuristic_correctness_check(coeffaccum.gb_coeffs_qq, primetracker.modulo)
            @info "Heuristic check failed."
            return false
        end
        @info "Heuristic check passed!"
    end

    if meta.randomizedcheck
        if !randomized_correctness_check!(coeffbuffer, coeffaccum, ring,
                                            coeffs_zz, gens_temp_ff, gb_ff, primetracker, ht)
            @info "Randomized check failed."
            return false
        end
        @info "Randomized check passed!"
    end

    if meta.guaranteedcheck
        return guaranteed_correctness_check!(ring, gb_ff.gens, coeffaccum.gb_coeffs_qq,
                                            exps, coeffs, gens_temp_ff, ht)
    end

    return true
end

@inline function threshold_in_heuristic(sznum, szden, szmod)
    # coefficient 1.10
    1.30*(sznum + sznum) >= szmod
end

# ln(num) + ln(den) < c ln p
function heuristic_correctness_check(gbcoeffs_qq, modulo)
    #=

    =#

    lnm = Base.GMP.MPZ.sizeinbase(modulo, 2)
    for i in 1:length(gbcoeffs_qq)
        @inbounds for j in 1:length(gbcoeffs_qq[i])
            n = numerator(gbcoeffs_qq[i][j])
            d = denominator(gbcoeffs_qq[i][j])
            # println(Base.GMP.MPZ.sizeinbase(c, 2), " ", lnm)
            if threshold_in_heuristic(Base.GMP.MPZ.sizeinbase(n, 2), Base.GMP.MPZ.sizeinbase(d, 2), lnm)
                return false
            end
        end
    end
    return true
end

# TODO: cleanup
function randomized_correctness_check!(
            coeffbuffer, coeffaccum,
            ring, coeffs_zz, gens_temp_ff, gb_ff, primetracker, ht)

    goodprime = nextgoodprime!(primetracker)

    reduce_modulo!(coeffbuffer, coeffs_zz, gens_temp_ff.coeffs, goodprime)
    gens_ff_copy = copy_basis_thorough(gens_temp_ff)
    cleanup_gens!(ring, gens_ff_copy, goodprime)

    # TODO: can encapsulate BigInt coefficients into coeffaccum
    gb_coeffs_zz = scale_denominators(coeffbuffer, coeffaccum.gb_coeffs_qq)
    gb_ff_copy = copy_basis_thorough(gb_ff)
    reduce_modulo!(coeffbuffer, gb_coeffs_zz, gb_ff_copy.coeffs, goodprime)
    cleanup_gens!(ring, gb_ff_copy, goodprime)

    @assert ring.ch == gb_ff_copy.ch == gens_ff_copy.ch

    # check that initial ideal contains in the computed groebner basis modulo goodprime
    normal_form_f4!(ring, gb_ff_copy, ht, gens_ff_copy)

    for i in 1:gens_ff_copy.ndone
        # meaning that it is not reduced
        if !isempty(gens_ff_copy.coeffs[i])
            return false
        end
    end

    # check that the basis is groebner basis modulo goodprime
    if !isgroebner_f4!(ring, gb_ff_copy, ht)
        return false
    end

    return true
end

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

function crt_issame_check!(coeffbuffer, coeffaccum, gbcoeffs_ff,
                                        modulo, prime)
    prevcoeffs = coeffaccum.prev_gb_coeffs_zz
    currcoeffs = coeffaccum.gb_coeffs_zz
    prevmodulo = div(modulo, prime)
    println("In same check")
    println(prevcoeffs)
    println(prevmodulo)
    println(currcoeffs)
    println(modulo)
    for i in 1:length(prevcoeffs)
        @inbounds for j in 1:length(prevcoeffs[i])
            if currcoeffs[i][j] != prevcoeffs[i][j]
                if (modulo - currcoeffs[i][j]) != (prevmodulo - prevcoeffs[i][j])
                    # if (currcoeffs[i][j] % prevmodulo) != prevcoeffs[i][j]
                        return false
                    # end
                end
            end
        end
    end
    return true
end

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

function guaranteed_correctness_check!(ring, gbexps, gb_coeffs_qq,
                                        exps, coeffs, gens_tmp_ff, ht)

    @warn "Setting parameter certify=true in groebner is not recommended."

    #=
    println("GB\n", gbexps, "\n", gb_coeffs_qq)
    println("GENS\n", gens_tmp_ff, "\n", coeffs)

    println("HT $(gens_tmp_ff.ntotal)")
    println(ht.exponents[1:3])
    =#

    # tmp_exps = [copy(init_gens.gens[i]) for i in 1:init_gens.ntotal]
    # tmp_coeffs = [
    #                copy(coeffs[i])
    #                for i in 1:init_gens.ntotal]

    gens_qq, _ = initialize_structures(ring, gens_tmp_ff.gens[1:gens_tmp_ff.ntotal], coeffs, ht)
    gb_qq, _   = initialize_structures(ring, gbexps, gb_coeffs_qq, ht)

    normalize_basis!(gb_qq)
    normalize_basis!(gens_qq)

    gens_qq_copy = copy_basis_thorough(gens_qq)

    # check that initial ideal contains in the computed groebner basis modulo goodprime
    normal_form_f4!(ring, gb_qq, ht, gens_qq_copy)

    # TODO

    # @error "normal form"
    # println(initial_ff)
    for i in 1:gens_qq_copy.ndone
        # meaning that it is not reduced
        if !isempty(gens_qq_copy.coeffs[i])
            return false
        end
    end

    # check that the basis is groebner basis modulo goodprime
    if !isgroebner_f4!(ring, gb_qq, ht)
        return false
    end

    return true
end

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

function linear_algebra_error(reason)
    throw(AssertionError("Probabilistic part of the algorithm failed, sorry. Please, submit a github issue. You can try to rerun with another random seed.\nReason estimated: $reason"))
end

function shape_control(gb_coeffs1, gb_coeffs2)
    if length(gb_coeffs1) != length(gb_coeffs2)
        linear_algebra_error("Unlucky reduction to zero in probabilistic linear algebra.")
    end

    for p in 1:length(gb_coeffs1)
        if length(gb_coeffs1[p]) != length(gb_coeffs2[p])
            linear_algebra_error("Unlucky cancellation of basis coefficients.")
        end
    end

    nothing
end

function rational_correctness_bound(modulo::BigInt)
    setprecision(2*Base.GMP.MPZ.sizeinbase(modulo, 2)) do
        ceil(BigInt, BigFloat(modulo)^(1/10))
    end
end

function coeff_control(primetracker, coeffaccum)
    modulo = primetracker.modulo
    bound = rational_correctness_bound(modulo)
    for c in coeffaccum.gb_coeffs_qq
        for cc in c
            if cc > bound
                return true
            end
        end
    end

    linear_algebra_error("Unlucky reduction to zero in probabilistic linear algebra.")

    return false
end


text:
### testing reflect method
using Test
using CUDA
using Random
# using Revise
using Flux
using HighDimPDE

@testset "CPU reflect methods" begin
    # testing reflection on batchsize 0
    d = 10
    X0 = fill(0.0,d)
    X1 = X0 + randn(d)
    X11 = HighDimPDE._reflect(X0, X1, fill(-1.,d), fill(1.,d))
    @test all((X11 .< 1) .& (X11 .> -1))

    @testset "testing NaN reflect" begin
        # testing reflection on batchsize
        y0 = X0[:]
        y1 = X0[:]
        dWall = zeros(Float32, d, 1000 )
        randn!(dWall) .* 10
        for i in 1:1000
            dW = @view dWall[:,i]
            y0 .= y1
            y1 .= y0 .+  dW
            y1 .= HighDimPDE._reflect(y0, y1, fill(-1.,d), fill(1.,d))
        end
        @test count(isnan.(y1)) == 0
    end
end

# testing reflection on batchsize
if CUDA.functional()
    CUDA.allowscalar(false)

    @testset "GPU reflect methods" begin
        d = 10
        X0 = fill(0.0f0,d)
        X1 = X0 + randn(Float32,d)
        X11 = HighDimPDE._reflect(X0, X1, fill(-1f0,d), fill(1f0,d))
        @testset "testing equivalence cpu gpu" begin
            batch_size = 1000
            a = repeat(X0[:],1,batch_size)
            b = a + 2 .* randn(size(a))
            a = hcat(X0,a) |> gpu
            b = hcat(X1,b) |> gpu
            s = -1f0; e = 1f0;
            b_gpu = HighDimPDE._reflect(a,b,s,e)
            b_cpu = b_gpu |> cpu
            @test prod(b_cpu .< 1) && prod(b_cpu .> -1)
            @test prod((b_cpu[:,1] .≈ X11[:,1]))
        end

        @testset "testing NaN reflect" begin
            # testing reflection on batchsize
            batch_size = 10000
            y0 = repeat(X0[:],1,batch_size) |> gpu
            y1 = repeat(X0[:],1,batch_size) |> gpu
            dWall = zeros(Float32, d, batch_size, 100 ) |> gpu
            randn!(dWall)
            for i in 1:100
                dW = @view dWall[:,:,i]
                y0 .= y1
                y1 .= y0 .+  dW
                y1 .= HighDimPDE._reflect(y0,y1,-1f0,1f0)
            end
            @test count(isnan.(y1)) == 0
        end
    end
end


text:
#="""
```
Setting{T}
```

The `Setting` type is an interface for computational settings that affect how the code runs
without affecting the mathematical definition of the model. It also provides support for
non-conflicting file names for output of 2 models that differ only in the values of their
computational settings.

### Fields
- `key::Symbol`: Name of setting
- `value::T`: Value of setting
- `print::Bool`: Indicates whether to append this setting's code and value to output file
  names. If true, output file names will include a suffix of the form `_code1=val1_code2=val2`
  etc. where codes are listed in alphabetical order.
- `code::String`: string of <=4 characters to print to output file suffixes when
  `print=true`.
- `description::String`: Short description of what the setting is used for.
"""
mutable struct Setting{T}
    key::Symbol                  # name of setting
    value::T                     # whatever the setting is
    print::Bool                  # whether or not to add this setting to the print
    code::String         # what gets printed to the print
    description::String  # description of what the setting is for
end

# for printing codes to filename string
Base.convert(::Type{T}, s::Setting{U}) where {T<:Number, U<:Number} = convert(T, s.value)
Base.convert(::Type{String}, s::Setting{String}) = convert(String, s.value)

Base.promote_rule(::Type{Setting{T}}, ::Type{U}) where {T<:Number,U<:Number} = promote_rule(T,U)
Base.promote_rule(::Type{Setting{String}}, ::Type{String}) = String # promote_rule(String, String)
Base.promote_rule(::Type{Setting{Bool}}, ::Type{Bool})  = Bool #promote_rule(Bool, Bool)

Base.string(s::Setting{String}) = string(s.value)

to_filestring(s::Setting) = string(s.code) * "=" * string(s.value)

# key, value constructor
Setting(key, value) = Setting(key, value, false, "", "")
Setting(key, value, description) = Setting(key, value, false, "", description)

function Base.show(io::IO, s::Setting)
    @printf io "%s: %s" s.key s.value
end

"""
```
(<=)(m::AbstractDSGEModel, s::Setting)
```
Syntax for adding a setting to a model/overwriting a setting via `m <= Setting(...)`
"""
function (<=)(m::AbstractDSGEModel, s::Setting)
    if !m.testing
        setting_field_name = :settings
    else
        setting_field_name = :test_settings
    end

    if !haskey(getfield(m, setting_field_name), s.key) ||
        isa(s.value, Function)
        getfield(m, setting_field_name)[s.key] = s
    else
        update!(getfield(m, setting_field_name)[s.key], s)
    end
end

"""
```
update!(a::Setting, b::Setting)
```

Update `a` with the fields of `b` if:

- The `key` field is updated if `a.key == b.key`
- The `print` boolean and `code` string are overwritten if `a.print` is false and
  `b.print` is true, or `a.print` is true, `b.print` is false, and
  b.code is non-empty.
- The `description` field is updated if `b.description` is nonempty
"""
function update!(a::Setting, b::Setting)
    # Make sure Setting a can appropriately be updated by b.
    if a.key ≠ b.key
        return
    end

    a.value = b.value

    # b overrides the print boolean of a if:
    # - a.print is false and b.print is true.
    # - a.print is true, b.print is false, and b.code is non-empty.
    # We must be able to tell if b was created via a constructor like `Setting(:key,
    # value)`, in which case the print, code, and description values are set to defaults. We
    # do not overwrite if we can't determine whether or not those fields are just the
    # defaults.
    if (a.print == false && b.print == true) ||
        (a.print == true && b.print == false && !isempty(b.code))
        a.print = b.print
    end
    if !isempty(b.code) && b.code ≠ a.code
        a.code = b.code
    end
    if !isempty(b.description) && b.description ≠ a.description
        a.description = b.description
    end
end

"""
```
get_setting(m::AbstractDSGEModel, setting::Symbol)
```

Returns the value of the setting
"""
function get_setting(m::AbstractDSGEModel, s::Symbol)
    if m.testing && in(s, keys(m.test_settings))
        return m.test_settings[s].value
    else
        return m.settings[s].value
    end
end
=#