text:
#
#      S1 - The manifold of the 1-dimensional sphere represented by angles
#
# Manopt.jl, R. Bergmann, 2018-06-26
import Base: exp, log, show

export Circle, S1Point, S1TVector
export distance, dot, embed, exp, log, manifoldDimension, norm, opposite
export parallelTransport
export parallelTransport, randomMPoint, randomTVector, typeofMPoint, typeofTVector
export zeroTVector
export show, getValue
export validateMPoint, validateTVector

export symRem
# Types
# ---
@doc doc"""
    Circle <: Manifold

The one-dimensional manifold $\mathbb S^1$ represented by angles.
Note that one can also use the $n$-dimensional sphere with $n=1$ to obtain the
same manifold represented by unit vectors in $\mathbb R^2$.

# Abbreviation

`S1`

# Constructor

    Circle()

construct a circle
"""
struct Circle <: Manifold
  name::String
  dimension::Int
  abbreviation::String
  Circle() = new("1-Sphere as angles",1,"S1")
end
@doc doc"""
    S1Point <: MPoint

a point $x\in\mathbb S^1$ represented by an angle `getValue(x)`$\in[-\pi,\pi)$,
usually referred to as “cyclic“ or “phase” data.
"""
struct S1Point <: MPoint
  value::Float64
  S1Point(value::Float64) = new(value)
end
getValue(x::S1Point) = x.value

@doc doc"""
    S1TVector <: TVector

a tangent vector $\xi\in\mathbb S^1$ represented by a real valiue
`getValue(ξ)`$\in\mathbb R$.
"""
struct S1TVector <: TVector
  value::Float64
  S1TVector(value::Float64) = new(value)
end
getValue(ξ::S1TVector) = ξ.value
# Traits
# ---
#(a) S1 is a matrix manifold
@traitimpl IsMatrixM{Circle}
@traitimpl IsMatrixP{S1Point}
@traitimpl IsMatrixTV{S1TVector}

# Functions
# ---

@doc doc"""
    distance(M,x,y)

return the distance two cyclic data items, which is given by
$\lvert (x-y)_{2\pi} \rvert $,
where $(\cdot)_{2\pi}$ is the symmetric remainder modulo $2\pi$,
see [`symRem`](@ref).
"""
distance(M::Circle, x::S1Point,y::S1Point) = abs( symRem(getValue(y) - getValue(x)) )
@doc doc"""
    dot(M,x,ξ,ν)

compute the inner product of two [`S1TVector`](@ref)s in the tangent space $T_x\mathbb S^1$
of the [`S1Point`](@ref) `x`. Since the values are angles, we
obtain $\langle \xi,\nu\rangle_x = \xi\nu$.
"""
dot(M::Circle, x::S1Point, ξ::S1TVector, ν::S1TVector) = getValue(ξ)*getValue(ν)

@doc doc"""
    embed(M,x)

embed the [`Circle`](@ref)` `[`Manifold`](@ref), i.e. turn the [`S1Point`](@ref)
into an [`SnPoint`](@ref) on the manifold [`Sphere`](@ref)`(1)` embedded in
$\mathbb R^2$.
"""
embed(M::Circle,x::S1Point) = SnPoint([cos(getValue(x)),sin(getValue(x))])

@doc doc"""
    exp(M,x,ξ,[t=1.0])

compute the exponential map on the [`Circle`](@ref) $\mathbb S^1$ with
respect to the [`S1Point`](@ref) `x` and the [`S1TVector`](@ref) `ξ`, which can
be shortened with `t` to `tξ`. The formula reads
```math
y = (x+\xi)_{2\pi}
```
where $(\cdot)_{2\pi}$ is the symmetric remainder modulo $2\pi$,
see [`symRem`](@ref).
"""
exp(M::Circle, x::S1Point,ξ::S1TVector,t::Float64=1.0) = S1Point( symRem(getValue(x) + t*getValue(ξ)) )

@doc doc"""
    log(M,x,y)

compute the logarithmic map on the [`Circle`](@ref) $\mathbb S^1$,
i.e., the [`S1TVector`](@ref) `ξ` whose corresponding
[`geodesic`](@ref) starting from [`S1Point`](@ref) `x` reaches the
[`S1Point`](@ref)` y` after time 1. The formula reads

```math
\xi = (y-x)_{2\pi},
```

where $(\cdot)_{2\pi}$ is the symmetric remainder modulo $2\pi$,
see [`symRem`](@ref).
"""
log(M::Circle, x::S1Point,y::S1Point)::S1TVector = S1TVector(symRem( getValue(y) - getValue(x) ))

"""
    manifoldDimension(x)

return the dimension of the manifold the [`S1Point`](@ref) `x` belongs to, i.e.
of the [`Circle`](@ref), which is 1.
"""
manifoldDimension(x::S1Point) = 1

"""
    manifoldDimension(M)

return the dimension of the [`Circle`](@ref) manifold, i.e., 1.
"""
manifoldDimension(M::Circle) = 1

@doc doc"""
    norm(M,x,ξ)

compute the norm of the [`S1TVector`](@ref) `ξ` in the tangent space
$T_x\mathcal M$ at [`S1Point`](@ref) `x` of the
[`Circle`](@ref) $\mathbb S^1$, which is just its absolute value $\lvert\xi\rvert$.
"""
norm(M::Circle, x::S1Point, ξ::S1TVector)::Float64 = abs( getValue(ξ) )

@doc doc"""
    opposite(M,x)

return the antipodal [`S1Point`](@ref) of `x` on the [`Circle`](@ref) `M`,
i.e. $y = (x+\pi)_{2\pi}$.
"""
opposite(M::Circle, x::S1Point) = S1Point( symRem(getValue(x)+π) )
@doc doc"""
    parallelTransport(M,x,y,ξ)

compute the parallel transport of the [`S1TVector`](@ref) `ξ` from the tangent
space $T_x\mathbb S^1$ at the [`S1Point`](@ref) `x` to $T_y\mathbb S^1$
at the [`S1Point`](@ref)` y`.
Since the [`Sphere`](@ref) `M` is represented in angles this is the identity.
"""
parallelTransport(M::Circle, x::S1Point, y::S1Point, ξ::S1TVector) = ξ
@doc doc"""
    randomMPoint(M,:Uniform)

return a random [`S1Point`](@ref) on the [`Circle`](@ref) $\mathbb S^1$ by
picking a random element from $[-\pi,\pi)$ uniformly.
"""
randomMPoint(M::Circle, ::Val{:Uniform}) = S1Point((rand()-0.5)*2*π)
randomMPoint(M::Circle) = randomMPoint(M,Val(:Uniform)) # introduce different default

@doc doc"""
    randomTVector(M,x [,Gaussian,σ=1.0])

returns a random tangent vector from the tangent space of the [`S1Point`](@ref)
 `x` on the [`Circle`](@ref) $\mathbb S^1$ by using a normal distribution with
mean 0 and standard deviation 1.
"""
randomTVector(M::Circle, x::S1Point, ::Val{:Gaussian}, σ::Real=1.0) = S1TVector(σ*randn())

@doc doc"""
    (Ξ,κ) = tangentONB(M,x,ξ)

compute an ONB within the tangent space $T_x\mathcal M$ such that $\xi$ is the
first vector and compute the eigenvalues of the curvature tensor
$R(\Xi,\dot g)\dot g$, where $g=g_{x,\xi}$ is the geodesic with $g(0)=x$,
$\dot g(0) = \xi$, i.e. $\kappa_1$ corresponding to $\Xi_1=\xi$ is zero.

# See also
 [`jacobiField`](@ref), [`adjointJacobiField`](@ref).
"""
tangentONB(M::Circle,x::S1Point,ξ::S1TVector) = [S1TVector(sign(getValue(ξ))==0 ? 1.0 : sign(getValue(ξ)))], [0.]

@doc doc"""
    (Ξ,κ) = tangentONB(M,x,y)

compute an ONB within the tangent space $T_x\mathcal M$ such that $\xi=\log_xy$ is the
first vector and compute the eigenvalues of the curvature tensor
$R(\Xi,\dot g)\dot g$, where $g=g_{x,\xi}$ is the geodesic with $g(0)=x$,
$\dot g(0) = \xi$, i.e. $\kappa_1$ corresponding to $\Xi_1=\xi$ is zero.

# See also
 [`jacobiField`](@ref), [`adjointJacobiField`](@ref).
"""
tangentONB(M::Circle, x::S1Point, y::S1Point) = tangentONB(M,x,log(M,x,y))

typeofTVector(::Type{S1Point}) = S1TVector
typeofMPoint(::Type{S1TVector}) = S1Point

"""
    typicalDistance(M)

returns the typical distance on the [`Circle`](@ref) `M`: π.
"""
typicalDistance(M::Circle) = π/2;

@doc doc"""
    ξ = zeroTVector(M,x)

returns a zero vector in the tangent space $T_x\mathcal M$ of the
[`S1Point`](@ref) $x\in\mathbb S^1$ on the [`Circle`](@ref)` S1`.
"""
zeroTVector(M::Circle, x::S1Point) = S1TVector(  zero( getValue(x) )  );
# Display
# ---
show(io::IO, M::Circle) = print(io, "The manifold S1 consisting of angles");
show(io::IO, x::S1Point) = print(io, "S1($( getValue(x) ))");
show(io::IO, ξ::S1TVector) = print(io, "S1T($( getValue(ξ) ))");
# little Helpers
# ---
@doc doc"""
    symRem(x,[T=π])

symmetric remainder of `x` with respect to the interall 2*`T`, i.e.
`(x+T)%2T`, where the default for `T` is $\pi$
"""
function symRem(x::Float64, T::Float64=Float64(π))::Float64
  return rem(x, 2*T,RoundNearest)
end
@doc doc"""
    validateMPoint(M,x)

validate, that a [`S1Point`](@ref) `x` has a valid value for a point on the
[`Circle`](@ref) `M`$=\mathbb S^1$, i.e. is within $[-\pi,\pi)$.
"""
function validateMPoint(M::Circle, x::S1Point)
    if (getValue(x) < - π) || (getValue(x) >= π) # out of range
        throw( ErrorException(
            "The Point $x is out of range for the Circle represented by anfgles in radians [-π,π),"
        ))
    end
    return true
end
@doc doc"""
    validateTVector(M,x,ξ)

validate, that the [`S1TVector`](@ref) `ξ` is a valid tangent vector in the
tangent space of the [`S1Point`](@ref) `x` ont the [`Circle`](@ref) `M`$=\mathbb S^1$,
though this is always the case since all real values are valid.
"""
function validateTVector(M::Circle,x::S1Point,ξ::S1TVector)
    return true
end


text:
using FFTW

test_fftwavenums() = FourierFlows.fftwavenums(6; L=2π) == [0, 1, 2, 3, -2, -1]

function test_domainaverage(n; xdir=true)
  g = TwoDGrid(n, 2π)
  if xdir; c = cos.(g.X).^2
  else;    c = cos.(g.Y).^2
  end
  0.5 ≈ FourierFlows.domainaverage(c, g)
end

function test_rms(n)
  g = TwoDGrid(n, 2π)
  q = cos.(g.X)
  isapprox(FourierFlows.rms(q), sqrt(1/2))
end

# This test could use some further work.
function test_radialspectrum(n, ahkl, ahρ; debug=false, atol=0.1)
  g = TwoDGrid(n, 2π)
  ah = ahkl.(g.K, g.L)
  ah[1, 1] = 0.0

  ρ, ahρ_estimate = FourierFlows.radialspectrum(ah, g; refinement=16)

  if debug
    println(sum(ahρ.(ρ) - ahρ_estimate))
    return ah, ahρ_estimate, ahρ.(ρ)
  else
    normalizeddiff = sum(abs.(ahρ.(ρ) - ahρ_estimate)) / length(ρ)
    return isapprox(normalizeddiff, 0.0, atol=atol)
  end
end

function integralsquare(func, grid)
    sum(abs2.(func))*grid.dx*grid.dy
end

function test_parsevalsum(func, grid; realvalued=true)
  # Compute integral in physical space
  integral = integralsquare(func, grid)
  # Compute integral in wavenumber space using Parseval's theorem
  if realvalued==true;      funch = rfft(func)
  elseif realvalued==false; funch = fft(func)
  end
  parsevalsum = FourierFlows.parsevalsum(abs2.(funch), grid)
  isapprox(integral, parsevalsum; atol=1.0e-14)
end

function test_parsevalsum2(func, grid; realvalued=true)
  # Compute integral in physical space
  integral = integralsquare(func, grid)
  # Compute integral in wavenumber space using Parseval's theorem
  if realvalued==true;      funch = rfft(func)
  elseif realvalued==false; funch = fft(func)
  end
  parsevalsum2 = FourierFlows.parsevalsum2(funch, grid)
  isapprox(integral, parsevalsum2; atol=1.0e-14)
end

"""
Compute the J(a, b) and compare with analytic_answer.
"""
function test_jacobian(a, b, analytic_answer, grid)
    # it's important to use atol for this test since when analytic_answer=0
    # the rtol is not conclusive (e.g., isapprox(1e-5, 0, rtol=1e-10) is false)
    isapprox(FourierFlows.jacobian(a, b, grid), analytic_answer;
             atol=g.nx*g.ny*1e-14)
end

"""
Test the createarrays macro.
"""
function test_createarrays(T=Float64, dims=(13, 45))
  a1, b1 = zeros(T, dims), zeros(T, dims)
  FourierFlows.@createarrays T dims a2 b2
  a1 == a2 && b1 == b2
end

"""
Test the peakedisotropicspectrum function.
"""
function testpeakedisotropicspectrum()
  n, L = 128, 2π
  gr = TwoDGrid(n, L)
  k0, E0 = 6, 0.5
  qi = FourierFlows.peakedisotropicspectrum(gr, k0, E0; allones=true)
  ρ, qhρ = FourierFlows.radialspectrum(rfft(qi).*gr.invKKrsq, gr)

  ρtest = ρ[ (ρ.>15.0) .& (ρ.<=17.5)]
  qhρtest = qhρ[ (ρ.>15.0) .& (ρ.<=17.5)]

  isapprox(abs.(qhρtest)/abs(qhρtest[1]), (ρtest/ρtest[1]).^(-2), rtol=5e-3)
end

abstract type TestVars <: AbstractVars end

physicalvars = [:a, :b]
transformvars = [ Symbol(var, :h) for var in physicalvars ]
 forcedvar = [:Fh]

varspecs = cat(FourierFlows.getfieldspecs(physicalvars, :(Array{T,2})),
  FourierFlows.getfieldspecs(transformvars, :(Array{Complex{T},2})), dims=1)

eval(FourierFlows.structvarsexpr(:VarsFields, physicalvars, transformvars))
eval(FourierFlows.structvarsexpr(:VarsFieldsParent, physicalvars, transformvars; parent=:TestVars))
eval(FourierFlows.structvarsexpr(:VarsSpecs, varspecs))
eval(FourierFlows.structvarsexpr(:VarsSpecsParent, varspecs; parent=:TestVars))

function VarsFields(g)
  @createarrays typeof(g.Lx) (g.nx, g.ny) a b
  @createarrays Complex{typeof(g.Lx)} (g.nkr, g.nl) ah bh
  VarsFields(a, b, ah, bh)
end

function VarsFieldsParent(g)
  @createarrays typeof(g.Lx) (g.nx, g.ny) a b
  @createarrays Complex{typeof(g.Lx)} (g.nkr, g.nl) ah bh
  VarsFieldsParent(a, b, ah, bh)
end

function VarsSpecs(g)
  @createarrays typeof(g.Lx) (g.nx, g.ny) a b
  @createarrays Complex{typeof(g.Lx)} (g.nkr, g.nl) ah bh
  VarsSpecs(a, b, ah, bh)
end

function VarsSpecsParent(g)
  @createarrays typeof(g.Lx) (g.nx, g.ny) a b
  @createarrays Complex{typeof(g.Lx)} (g.nkr, g.nl) ah bh
  VarsSpecsParent(a, b, ah, bh)
end

function test_structvarsexprFields(g)
  v1 = VarsFields(g)
  (
    typeof(v1.a)==Array{Float64,2} &&
    typeof(v1.ah)==Array{Complex{Float64},2} &&
    size(v1.a)==size(v1.b) &&
    size(v1.ah)==size(v1.bh) &&
    size(v1.a)==(g.nx, g.ny) &&
    size(v1.ah) == (g.nkr, g.nl)
  )
  v2 = VarsFieldsParent(g)
  (
    typeof(v2.a)==Array{Float64,2} &&
    typeof(v2.ah)==Array{Complex{Float64},2} &&
    size(v2.a)==size(v2.b) &&
    size(v2.ah)==size(v2.bh) &&
    size(v2.a)==(g.nx, g.ny) &&
    size(v2.ah) == (g.nkr, g.nl)
  )
end

function test_structvarsexprSpecs(g)
  v1 = VarsSpecs(g)
  (
    typeof(v1.a)==Array{Float64,2} &&
    typeof(v1.ah)==Array{Complex{Float64},2} &&
    size(v1.a)==size(v1.b) &&
    size(v1.ah)==size(v1.bh) &&
    size(v1.a)==(g.nx, g.ny) &&
    size(v1.ah) == (g.nkr, g.nl)
  )
  v2 = VarsSpecsParent(g)
  (
    typeof(v2.a)==Array{Float64,2} &&
    typeof(v2.ah)==Array{Complex{Float64},2} &&
    size(v2.a)==size(v2.b) &&
    size(v2.ah)==size(v2.bh) &&
    size(v2.a)==(g.nx, g.ny) &&
    size(v2.ah) == (g.nkr, g.nl)
  )
end

# Run tests -------------------------------------------------------------------

# Test on a rectangular grid
nx, ny = 64, 128   # number of points
Lx, Ly = 2π, 3π    # Domain width
g = TwoDGrid(nx, Lx, ny, Ly)
x, y = g.X, g.Y
k0, l0 = 2π/Lx, 2π/Ly

# Real and complex-valued functions
σ = 0.5
f1 = exp.(-(x.^2 + y.^2)/(2σ^2))
f2 = exp.( im*(2k0*x + 3l0*y.^2) ).*(
      exp.(-(x.^2 + y.^2)/(2σ^2)) + 2im*exp.(-(x.^2 + y.^2)/(5σ^2)) )

# Sine/Exp waves
k1, l1 = 2*k0, 6*l0
k2, l2 = 3*k0, -3*l0

sinkl1 = sin.(k1*x + l1*y)
sinkl2 = sin.(k2*x + l2*y)
expkl1 = exp.(im*(k1*x + l1*y))
expkl2 = exp.(im*(k2*x + l2*y))

# Analytical expression for the Jacobian of sin1 and sin2 and of exp1 and exp2
Jsinkl1sinkl2 = (k1*l2-k2*l1)*cos.(k1*x + l1*y).*cos.(k2*x + l2*y)
Jexpkl1expkl2 = (k2*l1-k1*l2)*exp.(im*((k1+k2)*x + (l1+l2)*y))

@test test_parsevalsum(f1, g; realvalued=true)   # Real valued f with rfft
@test test_parsevalsum(f1, g; realvalued=false)  # Real valued f with fft
@test test_parsevalsum(f2, g; realvalued=false)  # Complex valued f with fft
@test test_parsevalsum2(f1, g; realvalued=true)  # Real valued f with rfft
@test test_parsevalsum2(f1, g; realvalued=false) # Real valued f with fft
@test test_parsevalsum2(f2, g; realvalued=false) # Complex valued f with fft

@test test_jacobian(sinkl1, sinkl1, 0*sinkl1, g)  # Test J(a, a) = 0
@test test_jacobian(sinkl1, sinkl2, Jsinkl1sinkl2, g) # Test J(sin1, sin2) = Jsin1sin2
@test test_jacobian(expkl1, expkl2, Jexpkl1expkl2, g) # Test J(exp1, exp2) = Jexp1exps2

@test test_createarrays()
@test test_fftwavenums()
@test test_rms(32)
@test test_domainaverage(32; xdir=true)
@test test_domainaverage(32; xdir=false)
@test test_structvarsexprFields(g)
@test test_structvarsexprSpecs(g)

# Radial spectrum tests. Note that ahρ = ∫ ah ρ dθ.
n = 128; δ = n/10                 # Parameters
ahkl(k, l) = exp(-(k^2+l^2)/2δ^2) #  a = exp(-ρ²/2δ²)
    ahρ(ρ) = 2π*ρ*exp(-ρ^2/2δ^2)  # aᵣ = 2π ρ exp(-ρ²/2δ²)
@test test_radialspectrum(n, ahkl, ahρ)

ahkl(k, l) = exp(-(k^2+l^2)/2δ^2) * k^2/(k^2+l^2) #  a = exp(-ρ²/2δ²)*cos(θ)²
    ahρ(ρ) = π*ρ*exp(-ρ^2/2δ^2)                   # aᵣ = π ρ exp(-ρ²/2δ²)
@test test_radialspectrum(n, ahkl, ahρ)

@test testpeakedisotropicspectrum()

ext: jl


text:
module Juqst

export 	Tableau,
		cnot,
		hadamard,
		X,
		Z,
		phase,
		measure,
		kets,
		reinitialise,
		initialise,
		initialize

using LinearAlgebra

mutable struct Tableau
	state::Array{Int32,2}
	qubits::Int32
	showRaw::Bool
	commands::Array{String,1}
	executeCommands::Array{Expr,1}
	trackCommands::Bool
end

include("Initial.jl")

export 	cliffordToTableau,
		tableauToClifford,
		decomposeState,
		drawCircuit,
		qiskitCircuit,
		getNumberOfCliffords,
		getNumberOfSymplecticCliffords,
		getNumberOfBitStringsCliffords,
		generateRawCliffords,
		makeFromCommand

include("Symplectic.jl")
include("open-systems.jl")

export mat,
       liou,
       choi_liou_involution,
       swap_involution,
       choi2liou,
       choiX2liou,
       choi2kraus,
       kraus2choi,
       kraus2liou,
       liou2choi,
       liou2choiX,
       liou2kraus,
       liou2pauliliou,
       pauliliou2liou,
       depol,
       istp,
       iscp,
       ischannel,
       isunital,
       nearestu,
       nicePrint,
       writemime,
	   makeSuper,
	   choi2chi

include("rchannels.jl")

export 	fidelity,
		unitarity,
		unitarityPercent,
		genChannel,
		randomFidelityNoise,
		randomPrepNoise,
		randomMeasureNoise,
		genChannelMap

include("marginal.jl")

export readInCSVFile,transformToFidelity,fitTheFidelities,convertAndProject,gibbsRandomField,marginalise,projectSimplex,getIndices,reconstructedJS
export reconstruct,covarianceMatrix
export getGrainedP,mutualInformation,relativeEntropy,conditionalMutualInfo,covarianceMatrix,JSD
export correlationMatrix,marginaliseFromRawData

end # module


text:
function bus_gen_values(data, solution, value_key)
    bus_pg = Dict(i => 0.0 for (i, bus) in data["bus"])
    for (i, gen) in data["gen"]
        bus_pg["$(gen["gen_bus"])"] += solution["gen"][i][value_key]
    end
    return bus_pg
end

function parse_pm_pmi(file)
    pm_data = PowerModels.parse_file(file)
    pmi_data = PMI.get_pm_data(System(file))
    #TODO: remove the following for loop once PM#v0.20.0 is tagged
    for id in keys(pmi_data["bus"])
        pm_data["bus"][id]["vm"] = pmi_data["bus"][id]["vm"]
    end
    return pm_data, pmi_data
end

@testset "Test Power Flow" begin
    file = joinpath(PM_DATA_DIR, "matpower", "case5.m")
    for (model, solver) in MODEL_SOLVER_MAP
        @testset "Testing $model with $(solver.optimizer_constructor.name.module)" begin
            pm_result = run_pf(file, model, solver)
            pmi_result = run_pf(System(file), model, solver)

            @test pm_result["termination_status"] == pmi_result["termination_status"]
            @test isapprox(pm_result["objective"], pmi_result["objective"]; atol = 1e-2)
        end
    end
end

@testset "Test Optimal PowerFlow" begin
    file = joinpath(PM_DATA_DIR, "matpower", "case5.m")
    data = PowerModels.parse_file(file)
    pmi_data = PMI.get_pm_data(System(file))

    for (model, solver) in MODEL_SOLVER_MAP
        pm_result = run_opf(data, model, solver)
        pmi_result = run_opf(pmi_data, model, solver)

        @test pm_result["termination_status"] == pmi_result["termination_status"]
        @test isapprox(pm_result["objective"], pmi_result["objective"]; atol = 1.8e1)
    end
end

@testset "Test DC PowerFlow" begin
    file = joinpath(PM_DATA_DIR, "matpower", "case5.m")
    data, pmi_data = parse_pm_pmi(file)
    pm_native = compute_dc_pf(data)
    pm_opt = run_dc_pf(data, ipopt_solver)
    ps_native = compute_dc_pf(System(file))
    ps_opt = run_dc_pf(System(file), ipopt_solver)

    @test ps_opt["termination_status"] == pm_opt["termination_status"] == LOCALLY_SOLVED

    for (i, bus) in data["bus"]
        @test isapprox(
            pm_opt["solution"]["bus"][i]["va"],
            ps_opt["solution"]["bus"][i]["va"];
            atol = 1e-3,
        )
        @test isapprox(
            pm_native["solution"]["bus"][i]["va"],
            ps_native["solution"]["bus"][i]["va"];
            atol = 1e-3,
        )
    end
end

@testset "5-bus case" begin
    file = joinpath(PM_DATA_DIR, "matpower", "case5.m")
    data, pmi_data = parse_pm_pmi(file)

    pm_native = compute_ac_pf(data)
    pm_opt = run_ac_pf(data, ipopt_solver)

    ps_native = compute_ac_pf(pmi_data)
    ps_opt = run_ac_pf(pmi_data, ipopt_solver)

    @test ps_opt["termination_status"] == pm_opt["termination_status"] == LOCALLY_SOLVED

    bus_pg_pm_opt = bus_gen_values(data, pm_opt["solution"], "pg")
    bus_qg_pm_opt = bus_gen_values(data, pm_opt["solution"], "qg")

    bus_pg_pm_native = bus_gen_values(data, pm_native["solution"], "pg")
    bus_qg_pm_native = bus_gen_values(data, pm_native["solution"], "qg")

    bus_pg_ps_opt = bus_gen_values(pmi_data, ps_opt["solution"], "pg")
    bus_qg_ps_opt = bus_gen_values(pmi_data, ps_opt["solution"], "qg")

    bus_pg_ps_native = bus_gen_values(pmi_data, ps_native["solution"], "pg")
    bus_qg_ps_native = bus_gen_values(pmi_data, ps_native["solution"], "qg")

    for (i, bus) in data["bus"]
        @test isapprox(
            pm_opt["solution"]["bus"][i]["va"],
            ps_opt["solution"]["bus"][i]["va"];
            atol = 1e-4,
        )
        @test isapprox(
            pm_opt["solution"]["bus"][i]["vm"],
            ps_opt["solution"]["bus"][i]["vm"];
            atol = 1e-3,
        )
        @test isapprox(
            pm_native["solution"]["bus"][i]["va"],
            ps_native["solution"]["bus"][i]["va"];
            atol = 1e-7,
        )
        @test isapprox(
            pm_native["solution"]["bus"][i]["vm"],
            ps_native["solution"]["bus"][i]["vm"];
            atol = 1e-7,
        )

        @test isapprox(bus_pg_pm_opt[i], bus_pg_ps_opt[i]; atol = 1e-4)
        if !Sys.iswindows() #this is actually pretty different on Windows, not sure why
            @test isapprox(bus_qg_pm_opt[i], bus_qg_ps_opt[i]; atol = 1e-2)
        end
        @test isapprox(bus_pg_pm_native[i], bus_pg_ps_native[i]; atol = 1e-6)
        @test isapprox(bus_qg_pm_native[i], bus_qg_ps_native[i]; atol = 1e-2)
    end
end

max_stars_repo_name: NREL-SIIP/PowerModelsInterface

max_issues_repo_name: NREL-SIIP/PowerModelsInterface

max_forks_repo_name: NREL-SIIP/PowerModelsInterface


text:
"""
    slidingwindow(L::LT, f::FT, d::IT) where {LT <: SimpleSDMLayer, FT <: Function, IT <: Number}

This function will replace the value at any cell by applying the function `f` to
the array of cells values that are within a distance `d` (in kilometers) from
the focal cell. This is, for example, useful to use an average to smooth out the
layers. The distance is estimated using the haversine distance, assuming that
the radius of the Earth is 6371.0 km. This means that the size of the window
will vary a little bit across latitudes, but this is far better than using a
number of cells, which would have dramatic consequences near the poles.

It *always* returns a `SimpleSDMResponse`, and the cells containing `nothing`
will also not contain a value in the output. This is *different* from the
behavior of `coarsen`, which tends to expand the area of the layer in which we
have data.

This function is currently relatively slow. Performance improvements will arrive
at some point.
"""
function slidingwindow(layer::LT, f::FT, d::IT; threaded::Bool=Threads.nthreads()>1) where {LT <: SimpleSDMLayer, FT <: Function, IT <: Number}
    # We infer the return type from a call to the function on the first three elements
    return_type = typeof(f(collect(layer)[1:min(3, length(layer))]))

    # New layer using typed similar
    N = similar(layer, return_type)

    # Vector of all positions with a value
    filled_positions = keys(layer)

    # We then filter in the occupied positions
    if threaded
        Threads.@threads for pos in filled_positions
            N[pos] = f(_sliding_values(layer, pos, d))
        end
    else
        for pos in filled_positions
           N[pos] = f(_sliding_values(layer, pos, d))
        end
    end

    # And we return the object
    return N
end

function _sliding_values(layer, pt, d; R=6371.0)
    # Bounding box (approx.) for the sliding window of length d at the given point
    max_lat = min(layer.top, pt[2]+(180.0*d)/(π*R))
    min_lat = max(layer.bottom, pt[2]-(180.0*d)/(π*R))
    max_lon = min(layer.right, pt[1]+(360.0*d)/(π*R))
    min_lon = max(layer.left, pt[1]-(360.0*d)/(π*R))

    # Extracted layer for the sliding window
    _tmp = clip(layer; left=min_lon, right=max_lon, top=max_lat, bottom=min_lat)

    # Filter the correct positions
    filled_positions = keys(_tmp)
    neighbors = filter(p -> Distances.haversine((pt[1], pt[2]), (p[1], p[2]), R) < d, filled_positions)

    # Return
    return _tmp[neighbors]
end


text:
@doc raw"""
    PeriodicKernel([α=1 [,p=π]])

The periodic kernel is given by:

```math
\kappa(\mathbf{x},\mathbf{y}) =
\exp\left(-\alpha \sum_{i=1}^n \sin(p(x_i - y_i))^2\right)
\qquad p >0, \; \alpha > 0
```

where ``\mathbf{x}`` and ``\mathbf{y}`` are ``n`` dimensional vectors. The parameters ``p``
and ``\alpha`` are scaling parameters for the periodicity and the magnitude, respectively.
This kernel is useful when data has periodicity to it.
"""
struct PeriodicKernel{T<:AbstractFloat} <: MercerKernel{T}
    alpha::HyperParameter{T}
    PeriodicKernel{T}(α::Real) where {T<:AbstractFloat} = new{T}(
        HyperParameter(convert(T,α), interval(OpenBound(zero(T)), nothing))
    )
end
PeriodicKernel(α::T1 = 1.0) where {T1<:Real} = PeriodicKernel{floattype(T1)}(α)

@inline basefunction(::PeriodicKernel) = SineSquared()
@inline kappa(κ::PeriodicKernel{T}, z::T) where {T} = squaredexponentialkernel(z, getvalue(κ.alpha))


text:
function update!(dQ, Q, Qstage, rka1, rka2, rkb, dt, slow_δ, slow_dQ,
                 slow_scaling)
  @inbounds @loop for i = (1:length(Q);
                           (blockIdx().x - 1) * blockDim().x + threadIdx().x)
    if slow_δ !== nothing
      dQ[i] += slow_δ * slow_dQ[i]
    end
    Qstage[i] = rka1 * Q[i] + rka2 * Qstage[i] + dt * rkb * dQ[i]
    if slow_scaling !== nothing
      slow_dQ[i] *= slow_scaling
    end
  end
end

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
@testset "897.increasing-order-search-tree.jl" begin
    @test increasing_bst(TreeNode(5,
                                  TreeNode(3, TreeNode(2, TreeNode(1), nothing),
                                           TreeNode(4)),
                                  TreeNode(6, nothing,
                                           TreeNode(8, TreeNode(7), TreeNode(9))))) ==
          TreeNode(1, nothing,
                   TreeNode(2, nothing,
                            TreeNode(3, nothing,
                                     TreeNode(4, nothing,
                                              TreeNode(5, nothing,
                                                       TreeNode(6, nothing,
                                                                TreeNode(7, nothing,
                                                                         TreeNode(8,
                                                                                nothing,
                                                                                TreeNode(9)))))))))
end


text:
"""Given a matrix of energies `E` (in Hartrees), where each column
corresponds to a different angular symmetry, a vector of matrices `V`
containing the basis vectors for each angular symmetry, generate a
field-free Hamiltonian on diagonal form and transform all operators
`O` to the field-free basis. Optionally, `cutoff` (in eV) limits the
spectrum. Also returned is the index of the lowest, i.e., groundstate
energy."""
function hamiltonian{T<:AbstractFloat,
                     M<:AbstractMatrix}(E::Vector{Vector{T}},
                                        V::Vector{M},
                                        O::Vector{Operator},
                                        cutoff::Real = 0,
                                        verbose::Bool = false)
    n = length(E)
    sel = [(1:length(E[j])) for j = 1:n]

    # Find range of elements with energies below the cutoff energy.
    if cutoff != 0
        cutoff /= 27.211
        for j = 1:n
            i = findlast(e -> e <= cutoff, E[j])
            verbose && i == sel[j][end] && warn("Largest energy $(E[j][end]) smaller than cutoff $(cutoff).")
            i == 0 && error("No energy smaller than cutoff in partial wave $j.")
            sel[j] = 1:i
        end
    end

    # Find block coordinates
    coords = Vector{UnitRange{Int64}}(n)
    coords[1] = sel[1]
    for j = 2:n
        coords[j] = sel[j] + coords[j-1][end]
    end

    # Reduce the space
    E = vcat([E[j][sel[j]] for j = 1:n]...)
    V = Matrix{T}[V[j][:,sel[j]] for j in 1:n]

    # Field-free Hamiltonian in its eigenbasis
    H = sparse(Diagonal(E))

    # Transform all operators O into the eigenbasis
    O = map(O) do o
        ot = spzeros(eltype(first(o.o).m), size(H)...)
        for oo in o.o
            o12 = sparse(V[oo.k1]'*oo.m*V[oo.k2])
            ot[coords[oo.k1],coords[oo.k2]] += o12
            ot[coords[oo.k2],coords[oo.k1]] += o12'
        end
        ot
    end

    # Find groundstate
    gst = indmin(E)

    H,O,gst
end

function hamiltonian{T<:AbstractFloat,
                     M<:AbstractMatrix}(E::Vector{Vector{T}},
                                        V::Vector{M},
                                        O::Operator,
                                        cutoff::Real = 0,
                                        verbose::Bool = false)
    H,O,gst = hamiltonian(E, V, [O], cutoff, verbose)
    H,O[1],gst
end

export hamiltonian


text:
@testset "hyperkron_tests" begin
    @testset "num_multiset_permutations" begin
        @test num_multiset_permutations(Integer[]) == 1
        @test num_multiset_permutations([0]) == 1
        @test num_multiset_permutations([0,0,0]) == 1
        @test num_multiset_permutations([0,0,1]) == 3
        @test num_multiset_permutations([1,2,3,3]) == 12
        @test num_multiset_permutations([1,2,3,4]) == 24
        @test num_multiset_permutations([1,1,5,6,7,7,8]) == 1260
    end

    @testset "_next_distinct_character" begin
        @test _next_distinct_character(Integer[], 0, 0) == 1
        @test _next_distinct_character([0], 1, 1) == 2
        @test _next_distinct_character([0], 4, 1) == 2
        @test _next_distinct_character([0,0,0], 1, 3) == 4
        @test _next_distinct_character([0,0,0], 6, 3) == 4
        @test _next_distinct_character([0,0,1], 1, 3) == 3
        @test _next_distinct_character([0,0,1], 2, 3) == 3
        @test _next_distinct_character([0,0,1], 3, 3) == 4
        @test _next_distinct_character([1,1,1,3,5], 2, 5) == 4
        @test _next_distinct_character([4,5,7,7,7], 3, 5) == 6
        @test _next_distinct_character([4,5,7,7,7,8], 1, 6) == 2
    end

    @testset "unrank5!" begin
        @test unrank5!(Integer[], 1, 1, 0) == Integer[]
        @test unrank5!([0], 1, 1, 1) == [0]
        @test unrank5!([0], 2, 1, 1) == [0]
        @test unrank5!([0], 0, 1, 1) == [0]
        @test unrank5!([0,0,1], 1, 3, 3) == [0,0,1]
        @test unrank5!([0,0,1], 2, 3, 3) == [0,1,0]
        @test unrank5!([0,1,2,2], 7, 12, 4) == [2,0,1,2]
    end

    @testset "_morton_decode3" begin
        @test _morton_decode3(Integer[], 2, 8) == (1,1,1)
        @test _morton_decode3([0], 2, (2,2,2)) == (1,1,1)
        @test _morton_decode3([1], 2, (2,2,2)) == (1,1,2)
        @test _morton_decode3([7], 2, (2,2,2)) == (2,2,2)
        @test _morton_decode3([0,0], 2, (4,4,4)) == (1,1,1)
        @test _morton_decode3([3,4], 2, (4,4,4)) == (4,1,2)
        @test _morton_decode3([7,1,2], 2, (8,8,8)) == (4,5,7)
    end

    @testset "hyperkron_graph" begin
        G, _ = hyperkron_graph(kron_params(0.4, 0.8, 0.2, 0.5), 3)
        @test all([G.colptr[i] != G.rowval[i] for i = 1:size(G,1)]) # no self-loops
        @test maximum(G.nzval) == 1
        @test minimum(G.nzval) == 1
        @test typeof(G) == SparseMatrixCSC{Float64, Int64}
        @test G.m == G.n
    end
end


text:
## mapreduce implementation that skips NA

function skipna_init(f, op, na::BitArray, data::Array, ifirst::Int, ilast::Int)
    # Get first non-NA element
    ifirst = Base.findnextnot(na, ifirst)
    @inbounds d1 = data[ifirst]

    # Get next non-NA element
    ifirst = Base.findnextnot(na, ifirst+1)
    @inbounds d2 = data[ifirst]

    # Reduce first two elements
    (op(f(d1), f(d2)), ifirst)
end

function mapreduce_seq_impl_skipna(f, op, T, A::DataArray, ifirst::Int, ilast::Int)
    data = A.data
    na = A.na
    chunks = na.chunks

    v, i = skipna_init(f, op, na, data, ifirst, ilast)

    while i < ilast
        i += 1
        @inbounds na = Base.unsafe_bitgetindex(chunks, i)
        na && continue
        @inbounds d = data[i]
        v = op(v, f(d))
    end
    v
end

# Pairwise map-reduce
function mapreduce_pairwise_impl_skipna{T}(f, op, A::DataArray{T}, bytefirst::Int, bytelast::Int, n_notna::Int, blksize::Int)
    if n_notna <= blksize
        ifirst = 64*(bytefirst-1)+1
        ilast = min(64*bytelast, length(A))
        # Fall back to Base implementation if no NAs in block
        return ilast - ifirst + 1 == n_notna ? Base.mapreduce_seq_impl(f, op, A.data, ifirst, ilast) :
                                               mapreduce_seq_impl_skipna(f, op, T, A, ifirst, ilast)
    end

    # Find byte in the middle of range
    # The block size is restricted so that there will always be at
    # least two non-NA elements in the returned range
    chunks = A.na.chunks
    nmid = 0
    imid = bytefirst-1
    while nmid < (n_notna >> 1)
        imid += 1
        @inbounds nmid += count_zeros(chunks[imid])
    end

    v1 = mapreduce_pairwise_impl_skipna(f, op, A, bytefirst, imid, nmid, blksize)
    v2 = mapreduce_pairwise_impl_skipna(f, op, A, imid+1, bytelast, n_notna-nmid, blksize)
    op(v1, v2)
end

if isdefined(Base, :pairwise_blocksize)
    sum_pairwise_blocksize(f) = Base.pairwise_blocksize(f, +)
else
    const sum_pairwise_blocksize = Base.sum_pairwise_blocksize
end

mapreduce_impl_skipna{T}(f, op, A::DataArray{T}) =
    mapreduce_seq_impl_skipna(f, op, T, A, 1, length(A.data))
mapreduce_impl_skipna(f, op::typeof(@functorize(+)), A::DataArray) =
    mapreduce_pairwise_impl_skipna(f, op, A, 1, length(A.na.chunks),
                                   length(A.na)-countnz(A.na),
                                   max(128, sum_pairwise_blocksize(f)))

## general mapreduce interface

function _mapreduce_skipna{T}(f, op, A::DataArray{T})
    n = length(A)
    na = A.na

    nna = countnz(na)
    nna == n && return Base.mr_empty(f, op, T)
    nna == n-1 && return Base.r_promote(op, f(A.data[Base.findnextnot(na, 1)]))
    nna == 0 && return Base.mapreduce_impl(f, op, A.data, 1, n)

    mapreduce_impl_skipna(f, op, A)
end

# This is only safe when we can guarantee that if a function is passed
# NA, it returns NA. Otherwise we will fall back to the implementation
# in Base, which is slow because it's type-unstable, but guarantees the
# correct semantics
typealias SafeMapFuns @compat Union{typeof(@functorize(identity)), typeof(@functorize(abs)), typeof(@functorize(abs2)),
                            typeof(@functorize(exp)), typeof(@functorize(log)), typeof(@functorize(centralizedabs2fun))}
typealias SafeReduceFuns @compat Union{typeof(@functorize(+)), typeof(@functorize(*)), typeof(@functorize(max)), typeof(@functorize(min))}
function Base._mapreduce(f::SafeMapFuns, op::SafeReduceFuns, A::DataArray)
    any(A.na) && return NA
    Base._mapreduce(f, op, A.data)
end

function Base.mapreduce(f, op::Function, A::DataArray; skipna::Bool=false)
    is(op, +) ? (skipna ? _mapreduce_skipna(f, @functorize(+), A) : Base._mapreduce(f, @functorize(+), A)) :
    is(op, *) ? (skipna ? _mapreduce_skipna(f, @functorize(*), A) : Base._mapreduce(f, @functorize(*), A)) :
    is(op, &) ? (skipna ? _mapreduce_skipna(f, @functorize(&), A) : Base._mapreduce(f, @functorize(&), A)) :
    is(op, |) ? (skipna ? _mapreduce_skipna(f, @functorize(|), A) : Base._mapreduce(f, @functorize(|), A)) :
    skipna ? _mapreduce_skipna(f, op, A) : Base._mapreduce(f, op, A)
end

# To silence deprecations, but could be more efficient
Base.mapreduce(f, op::(@compat Union{typeof(@functorize(|)), typeof(@functorize(&))}), A::DataArray; skipna::Bool=false) =
    skipna ? _mapreduce_skipna(f, op, A) : Base._mapreduce(f, op, A)

Base.mapreduce(f, op, A::DataArray; skipna::Bool=false) =
    skipna ? _mapreduce_skipna(f, op, A) : Base._mapreduce(f, op, A)

Base.reduce(op, A::DataArray; skipna::Bool=false) =
    mapreduce(@functorize(identity), op, A; skipna=skipna)

## usual reductions

for (fn, op) in ((:(Base.sum), @functorize(+)),
                 (:(Base.prod), @functorize(*)),
                 (:(Base.minimum), @functorize(min)),
                 (:(Base.maximum), @functorize(max)))
    @eval begin
        $fn(f::(@compat Union{Function,$(supertype(typeof(@functorize(abs))))}), a::DataArray; skipna::Bool=false) =
            mapreduce(f, $op, a; skipna=skipna)
        $fn(a::DataArray; skipna::Bool=false) =
            mapreduce(@functorize(identity), $op, a; skipna=skipna)
    end
end

for (fn, f, op) in ((:(Base.sumabs), @functorize(abs), @functorize(+)),
                    (:(Base.sumabs2), @functorize(abs2), @functorize(+)))
    @eval $fn(a::DataArray; skipna::Bool=false) = mapreduce($f, $op, a; skipna=skipna)
end

## mean

Base.mean(a::DataArray; skipna::Bool=false) =
    sum(a; skipna=skipna) / (length(a.na)-(skipna ? countnz(a.na) : 0))

## variance

function Base.varm{T}(A::DataArray{T}, m::Number; corrected::Bool=true, skipna::Bool=false)
    if skipna
        n = length(A)
        na = A.na

        nna = countnz(na)
        nna == n && return convert(Base.momenttype(T), NaN)
        nna == n-1 && return convert(Base.momenttype(T),
                                     abs2(A.data[Base.findnextnot(na, 1)] - m)/(1 - @compat(Int(corrected))))

        /(nna == 0 ? Base.centralize_sumabs2(A.data, m, 1, n) :
                     mapreduce_impl_skipna(@functorize(centralizedabs2fun)(m), @functorize(+), A),
          n - nna - @compat(Int(corrected)))
    else
        any(A.na) && return NA
        Base.varm(A.data, m; corrected=corrected)
    end
end
Base.varm{T}(A::DataArray{T}, m::NAtype; corrected::Bool=true, skipna::Bool=false) = NA

function Base.var(A::DataArray; corrected::Bool=true, mean=nothing, skipna::Bool=false)
    mean == 0 ? Base.varm(A, 0; corrected=corrected, skipna=skipna) :
    mean == nothing ? varm(A, Base.mean(A; skipna=skipna); corrected=corrected, skipna=skipna) :
    isa(mean, (@compat Union{Number, NAtype})) ?
        varm(A, mean; corrected=corrected, skipna=skipna) :
        throw(ErrorException("Invalid value of mean."))
end

Base.stdm(A::DataArray, m::Number; corrected::Bool=true, skipna::Bool=false) =
    sqrt(varm(A, m; corrected=corrected, skipna=skipna))

Base.std(A::DataArray; corrected::Bool=true, mean=nothing, skipna::Bool=false) =
    sqrt(var(A; corrected=corrected, mean=mean, skipna=skipna))

## weighted mean

function Base.mean(a::DataArray, w::WeightVec; skipna::Bool=false)
    if skipna
        v = a .* w.values
        sum(v; skipna=true) / sum(DataArray(w.values, v.na); skipna=true)
    else
        anyna(a) ? NA : mean(a.data, w)
    end
end

function Base.mean{W,V<:DataArray}(a::DataArray, w::WeightVec{W,V}; skipna::Bool=false)
    if skipna
        v = a .* w.values
        sum(v; skipna=true) / sum(DataArray(w.values.data, v.na); skipna=true)
    else
        anyna(a) || anyna(w.values) ? NA : wsum(a.data, w.values.data) / w.sum
    end
end


text:
"""
    $(SIGNATURES)

Discretize given percentiles (unweighted). The smallest `pctV[1]` of values in `inV` end up in bin 1, etc.
"""
function discretize_given_percentiles(inV :: AbstractVector{F1},
    pctV :: AbstractVector{F2}) where {F1, F2}

    edgeV = bin_edges_from_percentiles(inV, pctV);
    return discretize(inV, edgeV)
end

"""
    $(SIGNATURES)

Discretize given percentiles (weighted). The smallest `pctV[1]` of values in `inV` end up in bin 1, etc.
"""
function discretize_given_percentiles(inV :: AbstractVector{F1},
    wtV :: AbstractVector{F2},
    pctV :: AbstractVector{F3}) where {F1, F2, F3}

   edgeV = bin_edges_from_percentiles(inV, wtV, pctV);
   return discretize(inV, edgeV)
end

## Discretize given bounds
# Equivalent to Discretizers LinearDiscretizer
function discretize(inV :: AbstractVector{F1},
    edgeV :: AbstractVector{F2}) where {F1, F2}

    n = length(edgeV);
    outV = zeros(Int64, size(inV));
    for i1 = n : -1 : 2
        outV[inV .<= edgeV[i1]] .= i1 - 1;
    end
    outV[inV .<= edgeV[1]] .= 0;
    return outV
end

"""
Bin edges from percentiles (unweighted)
Lowest bin includes minimum of inV
"""
function bin_edges_from_percentiles(inV :: AbstractVector{F1},
    pctV :: AbstractVector{F2}) where {F1, F2}

    edgeV = quantile(inV, [zero(F2); pctV]);
    # Ensure that lowest point is inside edges
    edgeV[1] -= F1(1e-8);
    return edgeV
end

## Bin edges from percentiles (weighted)
function bin_edges_from_percentiles(inV :: AbstractVector{F1},
    wtV :: AbstractVector{F2}, pctV :: AbstractVector{F3}) where {F1, F2, F3}

    @assert check_weights(wtV)  "Invalid weights"
    # Quantile is also a built in function
    edgeV = [minimum(inV) - F1(1e-8); quantile(inV, Weights(wtV), pctV)];
    edgeV[end] += F1(1e-8);
    return edgeV
end

"""
	$(SIGNATURES)

Count how many cases occur in each bin. Unweighted.
"""
function count_bins(inV :: AbstractVector{F1},
    edgeV :: AbstractVector{F2}) where {F1, F2}

    nBins = length(edgeV) - 1;
    cntV = Vector{Int}(undef, nBins);
    for iBin = 1 : nBins
        cntV[iBin] = sum(x -> inbin(x, edgeV, iBin), inV);
    end
    return cntV
end

"""
	$(SIGNATURES)

Count how many cases occur in each bin. Weighted.
"""
function count_bins(inV :: AbstractVector{F1},
    wtV :: AbstractVector{F2},
    edgeV :: AbstractVector{F3}) where {F1, F2, F3}

    @assert check_weights(wtV)  "Invalid weights"
    nBins = length(edgeV) - 1;
    massV = Vector{Float64}(undef, nBins);
    for iBin = 1 : nBins
        massV[iBin] = mapreduce(
            (x, wt) -> (inbin(x, edgeV, iBin) ? wt : zero(F2)),
            +, inV, wtV);
    end
    return massV
end

inbin(x, edgeV, iBin) = (x > edgeV[iBin])  &&  (x <= edgeV[iBin+1]);

function check_weights(wtV :: AbstractVector{F1}) where F1
    isValid = true;
    if any_less(wtV, zero(F1))
        @warn "Negative weights"
        isValid = false;
    end
    return isValid
end

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


text:
using JultraDark

@testset "Initialise Grids" begin
    grids = Grids(zeros(Complex{Float64}, 16, 16, 16), 1)
    @test typeof(grids) == Grids

    grids = Grids(1.0, 4)
    @test typeof(grids) == Grids

end

@testset "Initialise PencilGrids" begin
    grids = PencilGrids(zeros(Complex{Float64}, 16, 16, 16), 1)
    @test typeof(grids) == PencilGrids

    grids = PencilGrids(1.0, 4)
    @test typeof(grids) == PencilGrids

end


text:
@valueequals mutable struct ObjectTiles
    data::Array{Int, 2}

    ObjectTiles() = new(fill(-1, (40, 23)))

    ObjectTiles(d::Array{Int, 2}) = new(d)

    function ObjectTiles(s::String)
        s = chomp(replace(s, "\r\n" => "\n"))
        lines = split(s, "\n")
        lineLengths = length.(split.(lines, ','))

        cols = maximum(lineLengths)
        rows = length(lines)

        res = fill(-1, rows, cols)

        for (i, line) in enumerate(lines)
            if isempty(line)
                continue
            end

            res[i, 1:lineLengths[i]] = parse.(Int, split(line, ','))
        end

        return new(res)
    end
end

function Base.string(t::ObjectTiles)
    return join(join.([t.data[i, :] for i in 1:size(t.data, 1)], ','), "\n")
end

max_stars_repo_name: CelestialCartographers/Maple

max_issues_repo_name: CelestialCartographers/Maple

max_forks_repo_name: CelestialCartographers/Maple


text:
using MOCNeutronTransport
using BenchmarkTools
using Plots
pyplot()

import MOCNeutronTransport.area

struct UM_Tuple{P, F, T}
    points::NTuple{P, Point_2D{T}}
    faces::NTuple{F, Tuple{Vararg{Int64}}}
end
struct UM_Vector{P, T}
    points::NTuple{P, Point_2D{T}}
    faces::Vector{Vector{Int64}}
end

function area(mesh::UM_Vector)
    mesh_area = 0.0
    for face in mesh.faces
        if face[1] == 9
            mesh_area += area(
                        Quadrilateral_2D(
                                map(x -> mesh.points[x], Tuple(face[2:5]))
                            )
                        )
        end
    end
    return mesh_area
end

function area(mesh::UM_Tuple)
    mesh_area = 0.0
    for face in mesh.faces
        if face[1] == 9
            mesh_area += area(
                        Quadrilateral_2D(
                                map(x -> mesh.points[x], face[2:5])
                            )
                        )
        end
    end
    return mesh_area
end

tuple_times = Float64[]
vector_times = Float64[]
Nset = Int64.([1E1, 1E2, 1E3, 1E4, 1E5, 1E6])
Nctr = 1
for N = Nset
    # Generate the points
    points = Tuple(reduce(vcat, [[Point_2D(1.0*i, 0.0), Point_2D(1.0*i, 1.0)] for i = 0:N]))
    # Generate the faces for quadrilaterals
    quad_faces_vec = [[9, i, i+2, i+3, i+1] for i = 1:2:2N]
    quad_faces_tuple = Tuple([ Tuple(v) for v in quad_faces_vec ])

    faces_tuple = quad_faces_tuple
    faces_vec = quad_faces_vec

    tuple_mesh = UM_Tuple(points, faces_tuple)
    vector_mesh = UM_Vector(points, faces_vec)

    # Test speed of converting to quadrilaterals and getting the area
    println("N = 1E$Nctr")
    time = @belapsed area($vector_mesh)
    push!(vector_times, time/N)
    println("  Vector: $time")
    time = @belapsed area($tuple_mesh)
    push!(tuple_times, time/N)
    println("  Tuple: $time")
    global Nctr += 1
end
println("tuple_times ./ vector_times = ", tuple_times ./ vector_times)
plot(Nset,
     vector_times,
     label = "Vector",
     title = "Number of mesh cells vs area time per cell",
     xaxis=:log,
     yaxis=:log
    )
plot!(Nset, tuple_times, label = "Tuple")
xlabel!("Mesh cells")
ylabel!("Area computation time per cell")

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
# Autogenerated wrapper script for SBML_jll for x86_64-apple-darwin
export libsbml

using XML2_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("SBML")
JLLWrappers.@declare_library_product(libsbml, "@rpath/libsbml.5.dylib")
function __init__()
    JLLWrappers.@generate_init_header(XML2_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libsbml,
        "lib/libsbml.5.19.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# const VTKtypemap = Dict{DataType, Int}(P1=>1, L2=>3, T3=>5,
#     Q4=>9, T4=>10, H8=>12, Q8=>23,
#     L3=>21, T6=>22,
#     T10=>24, H20=>25)
const _VTK_TYPE_MAP = Dict{AbsShapeDesc, Int}(P1=>1, L2=>3, T3=>5, Q4=>9, T4=>10, H8=>12, T6=>22, Q8=>23)

"""
    vtkwrite(filename, connectivity)

Write VTK file with the mesh.

- `connectivity` = incidence relation of the type `d -> 0`. It must have a "geom"
  attribute for access to the locations of all the vertices that the
  connectivity incidence relation references.
- `data` = array of named tuples. Names of attributes of either the left or the
  right shape collection of the `connectivity` incidence relation. Property names: `:name` required, name of the attribute. `:allxyz` optional, pad the attribute to be a three-dimensional quantity (in the global Cartesian coordinate system).
"""
function vtkwrite(filename, connectivity, data = [])
    locs = attribute(connectivity.right, "geom")
    # Fill in a matrix of point coordinates
    points = fill(zero(eltype(locs[1])), length(locs[1]), length(locs))
    for i in 1:length(locs)
        points[:, i] .= locs[i]
    end #
    # Figure out the cell type
    celltype = WriteVTK.VTKCellTypes.VTKCellType(_VTK_TYPE_MAP[shapedesc(connectivity.left)])
    # Prepare an array of the cells
    cells = [MeshCell(celltype, [j for j in connectivity[i]]) for i in 1:nrelations(connectivity)]
    vtkfile = vtk_grid(filename, points, cells, compress=3)
    for nt in data
        an = nt.name # name of the shape collection
        pn = propertynames(nt)
        allxyz = (:allxyz in pn ? nt.allxyz : false)
        if an in keys(connectivity.right.attributes)
            a = attribute(connectivity.right, an)
            nc = length(a[1])
            ncz = ((nc < 3) && allxyz ? ncz = 3 : nc)
            pdata = fill(0.0, ncz, length(a))
            for j in 1:nc
                for i in 1:length(a)
                    pdata[j, i] = a[i][j]
                end
            end
            vtkfile[an] = pdata
        elseif an in keys(connectivity.left.attributes)
            a = attribute(connectivity.left, an)
            nc = length(a[1])
            cdata = fill(0.0, nc, length(a))
            for j in 1:nc
                for i in 1:length(a)
                    cdata[j, i] = a[i][j]
                end
            end
            vtkfile[an] = cdata
        end
    end
    return vtk_save(vtkfile)
end

"""
    export_MESH(meshfile, mesh)

Import vertices and shapes in the MESH format.

# Output
Data dictionary, with keys
- "`vertices`" (vertices),
- "`shapes`" (array of shape collections).
"""
function export_MESH(meshfile, connectivity)
    meshfilebase, ext = splitext(meshfile)
    # Name of the shape descriptor
    t = shapedesc(connectivity.left).name
    datinfo = [meshfilebase * "-xyz.dat", t, meshfilebase * "-conn.dat"]
    locs = attribute(connectivity.right, "geom")
    open(datinfo[1], "w") do file
        X = [locs[idx] for idx in 1:length(locs)]
        writedlm(file, X, ' ')
    end
    open(datinfo[3], "w") do file
        c = [connectivity[idx] for idx in 1:nrelations(connectivity)]
        writedlm(file, c, ' ')
    end
    open(meshfilebase * ".mesh", "w") do file
        writedlm(file, datinfo)
    end

    return true
end


text:
# Now that the initialization procedure is complete, you can reboot the
# reactor.
#
# Starting with all cubes off, run all of the reboot steps for all cubes in the
# reactor.
#
# Consider the following reboot steps:
#
# on x=-5..47,y=-31..22,z=-19..33
# on x=-44..5,y=-27..21,z=-14..35
# on x=-49..-1,y=-11..42,z=-10..38
# on x=-20..34,y=-40..6,z=-44..1
# off x=26..39,y=40..50,z=-2..11
# on x=-41..5,y=-41..6,z=-36..8
# off x=-43..-33,y=-45..-28,z=7..25
# on x=-33..15,y=-32..19,z=-34..11
# off x=35..47,y=-46..-34,z=-11..5
# on x=-14..36,y=-6..44,z=-16..29
# on x=-57795..-6158,y=29564..72030,z=20435..90618
# on x=36731..105352,y=-21140..28532,z=16094..90401
# on x=30999..107136,y=-53464..15513,z=8553..71215
# on x=13528..83982,y=-99403..-27377,z=-24141..23996
# on x=-72682..-12347,y=18159..111354,z=7391..80950
# on x=-1060..80757,y=-65301..-20884,z=-103788..-16709
# on x=-83015..-9461,y=-72160..-8347,z=-81239..-26856
# on x=-52752..22273,y=-49450..9096,z=54442..119054
# on x=-29982..40483,y=-108474..-28371,z=-24328..38471
# on x=-4958..62750,y=40422..118853,z=-7672..65583
# on x=55694..108686,y=-43367..46958,z=-26781..48729
# on x=-98497..-18186,y=-63569..3412,z=1232..88485
# on x=-726..56291,y=-62629..13224,z=18033..85226
# on x=-110886..-34664,y=-81338..-8658,z=8914..63723
# on x=-55829..24974,y=-16897..54165,z=-121762..-28058
# on x=-65152..-11147,y=22489..91432,z=-58782..1780
# on x=-120100..-32970,y=-46592..27473,z=-11695..61039
# on x=-18631..37533,y=-124565..-50804,z=-35667..28308
# on x=-57817..18248,y=49321..117703,z=5745..55881
# on x=14781..98692,y=-1341..70827,z=15753..70151
# on x=-34419..55919,y=-19626..40991,z=39015..114138
# on x=-60785..11593,y=-56135..2999,z=-95368..-26915
# on x=-32178..58085,y=17647..101866,z=-91405..-8878
# on x=-53655..12091,y=50097..105568,z=-75335..-4862
# on x=-111166..-40997,y=-71714..2688,z=5609..50954
# on x=-16602..70118,y=-98693..-44401,z=5197..76897
# on x=16383..101554,y=4615..83635,z=-44907..18747
# off x=-95822..-15171,y=-19987..48940,z=10804..104439
# on x=-89813..-14614,y=16069..88491,z=-3297..45228
# on x=41075..99376,y=-20427..49978,z=-52012..13762
# on x=-21330..50085,y=-17944..62733,z=-112280..-30197
# on x=-16478..35915,y=36008..118594,z=-7885..47086
# off x=-98156..-27851,y=-49952..43171,z=-99005..-8456
# off x=2032..69770,y=-71013..4824,z=7471..94418
# on x=43670..120875,y=-42068..12382,z=-24787..38892
# off x=37514..111226,y=-45862..25743,z=-16714..54663
# off x=25699..97951,y=-30668..59918,z=-15349..69697
# off x=-44271..17935,y=-9516..60759,z=49131..112598
# on x=-61695..-5813,y=40978..94975,z=8655..80240
# off x=-101086..-9439,y=-7088..67543,z=33935..83858
# off x=18020..114017,y=-48931..32606,z=21474..89843
# off x=-77139..10506,y=-89994..-18797,z=-80..59318
# off x=8476..79288,y=-75520..11602,z=-96624..-24783
# on x=-47488..-1262,y=24338..100707,z=16292..72967
# off x=-84341..13987,y=2429..92914,z=-90671..-1318
# off x=-37810..49457,y=-71013..-7894,z=-105357..-13188
# off x=-27365..46395,y=31009..98017,z=15428..76570
# off x=-70369..-16548,y=22648..78696,z=-1892..86821
# on x=-53470..21291,y=-120233..-33476,z=-44150..38147
# off x=-93533..-4276,y=-16170..68771,z=-104985..-24507
#
# After running the above reboot steps, 2758514936282235 cubes are on. (Just
# for fun, 474140 of those are also in the initialization procedure region.)
#
# Starting again with all cubes off, execute all reboot steps. Afterward,
# considering all cubes, how many cubes are on?

Cube = NTuple{3,UnitRange{Int}}

"""
    splitrange(a, b)

Given two ranges, split the first on the second. If the second is entirely
contained within the first, this will return two new ranges. Otherwise, if the
second is at the beginning or end of the first, this will only return one
range.

It is assumed that the first range completely contains the second. If the
ranges are equal, this function returns an empty vector.
"""
function splitrange(a::UnitRange{Int}, b::UnitRange{Int})
  if a.start < b.start
    if a.stop > b.stop
      return [a.start:(b.start-1), (b.stop+1):a.stop]
    else
      return [a.start:(b.start-1)]
    end
  elseif a.stop > b.stop
    return [(b.stop+1):a.stop]
  else
    return []
  end
end

"""
    cubediff(a, b)

Given two cubes, computes the cube(s) that remain after subtracting the second
cube from the first. Returns nothing if the second cube completely contains the
first. Otherwise, returns a vector of cubes.
"""
function cubediff(a::Cube, b::Cube)
  intersects = intersect.(a, b)
  any(isequal(0), length.(intersects)) && return [a]
  length.(a) == length.(intersects) && return nothing

  results = Cube[]
  (xpieces, ypieces, zpieces) = splitrange.(a, intersects)
  for zpiece in zpieces
    push!(results, (a[1], a[2], zpiece))
  end
  for ypiece in ypieces
    push!(results, (a[1], ypiece, intersects[3]))
  end
  for xpiece in xpieces
    push!(results, (xpiece, intersects[2], intersects[3]))
  end
  return results
end

# do eet
on = Cube[]
foreach(eachline("input.txt")) do l
  # parse line
  (instr,coords...) = match(
    r"(on|off) x=(-?\d+)..(-?\d+),y=(-?\d+)..(-?\d+),z=(-?\d+)..(-?\d+)",
    l
  )
  (x1,x2,y1,y2,z1,z2) = parse.(Int, coords)

  if instr == "on"
    # given this new cube to turn on, subtract all regions that are already on
    # and add the resulting cubes to the list.
    newon = Cube[(x1:x2, y1:y2, z1:z2)]
    for alreadyon in on
      nextnewon = Cube[]
      for teston in newon
        cubes = cubediff(teston, alreadyon)
        !isnothing(cubes) && append!(nextnewon, cubes)
      end
      newon = nextnewon
      length(newon) == 0 && break
    end
    append!(on, newon)
  else
    # given this new cube to turn off, subtract it from all regions that are
    # already on, replacing the list of on cubes.
    newon = Cube[]
    off = (x1:x2, y1:y2, z1:z2)
    for alreadyon in on
      cubes = cubediff(alreadyon, off)
      !isnothing(cubes) && append!(newon, cubes)
    end
    global on = newon
  end
end

println("Result: ", mapreduce(o -> prod(length.(o)), +, on))


text:
using PkgTemplateXKDR
using Documenter

DocMeta.setdocmeta!(PkgTemplateXKDR, :DocTestSetup, :(using PkgTemplateXKDR); recursive=true)

makedocs(;
    modules=[PkgTemplateXKDR],
    authors="xKDR Forum",
    repo="https://github.com/xKDR/PkgTemplateXKDR.jl/blob/{commit}{path}#{line}",
    sitename="$PkgTemplateXKDR.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://xKDR.github.io/PkgTemplateXKDR.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/xKDR/PkgTemplateXKDR.jl",
    target = "build",
    devbranch="main"
)


text:
function clustercoeffs_test()
    A = load_matrix_network("clique-10")
    cc = clustercoeffs(MatrixNetwork(A))

    v = ones(Int64,10)
    if v != cc
        error("clustercoeffs failed")
    end

    return true
end


text:
export dbc_icons

dbc_icons = (
    BOOTSTRAP = "https://cdn.jsdelivr.net/npm/bootstrap-icons@1.8.1/font/bootstrap-icons.css",
    FONT_AWESOME = "https://use.fontawesome.com/releases/v6.1.1/css/all.css",
)


text:
using ECOS
const ECOS_PRIMAL_FACTORY = MOI.OptimizerWithAttributes(ECOS.Optimizer, MOI.Silent() => true)
const ECOS_DUAL_FACTORY = dual_optimizer(ECOS_PRIMAL_FACTORY)
const ECOS_PRIMAL_OPT = MOI.instantiate(ECOS_PRIMAL_FACTORY)
const ECOS_DUAL_OPT = MOI.instantiate(ECOS_DUAL_FACTORY)

# push!(dual_linear_optimizer, DualOptimizer(ECOS.Optimizer(verbose = 0)))

@testset "ECOS conic Problems" begin
    @testset "ECOS Conic linear problems" begin
        list_of_soc_problems = [
            conic_linear1_test,
            conic_linear2_test,
            conic_linear3_test,
            conic_linear4_test
        ]
        test_strong_duality(list_of_soc_problems, ECOS_PRIMAL_FACTORY)
    end
    @testset "ECOS SOC problems" begin
        list_of_soc_problems = [
            soc1_test,
            soc2_test,
            soc3_test,
            soc4_test,
            soc5_test,
            soc6_test
        ]
        test_strong_duality(list_of_soc_problems, ECOS_PRIMAL_FACTORY)
    end
    @testset "ECOS RotatedSOC problems" begin
        list_of_rsoc_problems = [
            rsoc1_test,
            rsoc2_test,
            rsoc3_test,
            rsoc4_test
        ]
        test_strong_duality(list_of_rsoc_problems, ECOS_PRIMAL_FACTORY)
    end
end


text:
"""
    struct DiagnosticOptions <: Any
Hold a pointer to a `clang::DiagnosticOptions` object.
"""
struct DiagnosticOptions
    ptr::CXDiagnosticOptions
end

Base.unsafe_convert(::Type{CXDiagnosticOptions}, x::DiagnosticOptions) = x.ptr
Base.cconvert(::Type{CXDiagnosticOptions}, x::DiagnosticOptions) = x


text:
#
# interp/sparse.jl --
#
# Implement sparse linear interpolator.
#
#------------------------------------------------------------------------------
#
# This file is part of the LinearInterpolators package licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2015-2016, Éric Thiébaut, Jonathan Léger & Matthew Ozon.
# Copyright (C) 2016-2021, Éric Thiébaut.
#

# All code is in a module to "hide" private methods.
module SparseInterpolators

export
    SparseInterpolator,
    SparseUnidimensionalInterpolator

using InterpolationKernels

using LazyAlgebra
using LazyAlgebra.Foundations
import LazyAlgebra: apply, apply!, vcreate, output_size, input_size

import Base: axes, eltype, size
import SparseArrays: sparse

using ..LinearInterpolators
using ..LinearInterpolators: limits, getcoefs
import ..LinearInterpolators.Meta
import ..LinearInterpolators: coefficients, columns, rows,
    fit, regularize, regularize!

abstract type AbstractSparseInterpolator{T<:AbstractFloat} <: LinearMapping end

eltype(A::AbstractSparseInterpolator) = eltype(typeof(A))
eltype(::Type{<:AbstractSparseInterpolator{T}}) where {T} = T

struct SparseInterpolator{T<:AbstractFloat,S,N} <: AbstractSparseInterpolator{T}
    C::Vector{T}
    J::Vector{Int}
    nrows::Int
    ncols::Int
    dims::Dims{N} # dimensions of result
    function SparseInterpolator{T,S,N}(C::Vector{T},
                                       J::Vector{Int},
                                       dims::Dims{N},
                                       ncols::Int) where {T,S,N}
        @assert S ≥ 1
        @assert minimum(dims) ≥ 1
        nrows = prod(dims)
        nvals = S*nrows       # number of non-zero coefficients
        @assert length(C) == nvals
        @assert length(J) == nvals
        new{T,S,N}(C, J, nrows, ncols, dims)
    end
end

# Interpolator can be used as a function.
(A::SparseInterpolator)(x::AbstractVector) = apply(A, x)

output_size(A::SparseInterpolator) = A.dims
input_size(A::SparseInterpolator) = (A.ncols,)
width(A::SparseInterpolator{T,S,N}) where {T,S,N} = S
coefficients(A::SparseInterpolator) = A.C
columns(A::SparseInterpolator) = A.J
function rows(A::SparseInterpolator{T,S,N}) where {T,S,N}
    nrows = A.nrows
    nvals = S*nrows       # number of non-zero coefficients
    @assert length(A.C) == nvals
    @assert length(A.J) == nvals
    I = Array{Int}(undef, nvals)
    k0 = 0
    for i in 1:nrows
        for s in 1:S
            k = k0 + s
            @inbounds I[k] = i
        end
        k0 += S
    end
    return I
end

# Convert to a sparse matrix.
sparse(A::SparseInterpolator) =
    sparse(rows(A), columns(A), coefficients(A), A.nrows, A.ncols)

"""
    A = SparseInterpolator{T=eltype(ker)}(ker, pos, grd)

yields a sparse linear interpolator suitable for interpolating with kernel
`ker` at positions `pos` a function sampled on the grid `grd`.  Optional
parameter `T` is the floating-point type of the coefficients of the operator
`A`.  Call `eltype(A)` to query the type of the coefficients of the sparse
interpolator `A`.

Then `y = apply(A, x)` or `y = A(x)` or `y = A*x` yield the result of
interpolation array `x`.  The shape of `y` is the same as that of `pos`.
Formally, this amounts to computing:

    y[i] = sum_j ker((pos[i] - grd[j])/step(grd))*x[j]

with `step(grd)` the (constant) step size between the nodes of the grid `grd`
and `grd[j]` the `j`-th position of the grid.

"""
SparseInterpolator(ker::Kernel{T}, args...) where {T<:AbstractFloat} =
    SparseInterpolator{T}(ker, args...)

@deprecate SparseInterpolator(T::Type{<:AbstractFloat}, ker::Kernel, args...) SparseInterpolator{T}(ker, args...)

SparseInterpolator{T}(ker::Kernel, args...) where {T<:AbstractFloat} =
    SparseInterpolator{T}(T(ker), args...)

function SparseInterpolator{T}(ker::Kernel{T},
                               pos::AbstractArray{<:Real},
                               grd::AbstractRange) where {T<:AbstractFloat}
    SparseInterpolator{T}(ker, fractional_index(T, pos, grd),
                          CartesianIndices(axes(pos)), length(grd))
end

function SparseInterpolator{T}(ker::Kernel{T},
                               pos::AbstractArray{<:Real},
                               len::Integer) where {T<:AbstractFloat}
    SparseInterpolator{T}(ker, fractional_index(T, pos),
                          CartesianIndices(axes(pos)), len)
end

function SparseInterpolator{T}(ker::Kernel{T,S},
                               f::Function,
                               R::CartesianIndices{N},
                               ncols::Integer) where {T<:AbstractFloat,S,N}
    C, J = _sparsecoefs(R, Int(ncols), ker, f)
    return SparseInterpolator{T,S,N}(C, J, size(R), ncols)
end

@generated function _sparsecoefs(R::CartesianIndices{N},
                                 ncols::Int,
                                 ker::Kernel{T,S},
                                 f::Function) where {T,S,N}

    J_ = [Symbol(:j_,s) for s in 1:S]
    C_ = [Symbol(:c_,s) for s in 1:S]
    code = (Meta.generate_getcoefs(J_, C_, :ker, :lim, :x),
            [:( J[k+$s] = $(J_[s]) ) for s in 1:S]...,
            [:( C[k+$s] = $(C_[s]) ) for s in 1:S]...)

    quote
        lim = limits(ker, ncols)
        nvals = S*length(R)
        J = Array{Int}(undef, nvals)
        C = Array{T}(undef, nvals)
        k = 0
        @inbounds for i in R
            x = convert(T, f(i))
            $(code...)
            k += S
        end
        return C, J
    end
end

function _check(A::SparseInterpolator{T,S,N},
                out::AbstractArray{T,N},
                inp::AbstractVector{T}) where {T,S,N}
    nvals = S*A.nrows # number of non-zero coefficients
    J, ncols = A.J, A.ncols
    length(A.C) == nvals ||
        error("corrupted sparse interpolator (bad number of coefficients)")
    length(J) == nvals ||
        error("corrupted sparse interpolator (bad number of indices)")
    length(inp) == ncols ||
        error("bad vector length (expecting $(A.ncols), got $(length(inp)))")
    size(out) == A.dims ||
        error("bad output array size (expecting $(A.dims), got $(size(out)))")
    length(out) == A.nrows ||
        error("corrupted sparse interpolator (bad number of \"rows\")")
    @inbounds for k in 1:nvals
        1 ≤ J[k] ≤ ncols ||
            error("corrupted sparse interpolator (out of bound indices)")
    end
end

function vcreate(::Type{Direct},
                 A::SparseInterpolator{T,S,N},
                 x::AbstractVector{T},
                 scratch::Bool=false) where {T,S,N}
    return Array{T}(undef, output_size(A))
end

function vcreate(::Type{Adjoint},
                 A::SparseInterpolator{T,S,N},
                 x::AbstractArray{T,N},
                 scratch::Bool=false) where {T,S,N}
    return Array{T}(undef, input_size(A))
end

function apply!(α::Real,
                ::Type{Direct},
                A::SparseInterpolator{Ta,S,N},
                x::AbstractVector{Tx},
                scratch::Bool,
                β::Real,
                y::AbstractArray{Ty,N}) where {Ta,Tx<:Real,
                                               Ty<:AbstractFloat,S,N}
    _check(A, y, x)
    if α == 0
        vscale!(y, β)
    else
        T = float(promote_type(Ta, Tx))
        alpha = convert(T, α)
        nrows, ncols = A.nrows, A.ncols
        C, J = coefficients(A), columns(A)
        k0 = 0
        if β == 0
            @inbounds for i in 1:nrows
                sum = zero(T)
                @simd for s in 1:S
                    k = k0 + s
                    j = J[k]
                    sum += C[k]*x[j]
                end
                y[i] = alpha*sum
                k0 += S
            end
        else
            beta = convert(Ty, β)
            @inbounds for i in 1:nrows
                sum = zero(T)
                @simd for s in 1:S
                    k = k0 + s
                    j = J[k]
                    sum += C[k]*x[j]
                end
                y[i] = alpha*sum + beta*y[i]
                k0 += S
            end
        end
    end
    return y
end

function apply!(α::Real,
                ::Type{Adjoint},
                A::SparseInterpolator{Ta,S,N},
                x::AbstractArray{Tx,N},
                scratch::Bool,
                β::Real,
                y::AbstractVector{Ty}) where {Ta,Tx<:Real,
                                              Ty<:AbstractFloat,S,N}
    _check(A, x, y)
    vscale!(y, β)
    if α != 0
        T = float(promote_type(Ta, Tx))
        alpha = convert(T, α)
        nrows, ncols = A.nrows, A.ncols
        C, J = coefficients(A), columns(A)
        k0 = 0
        @inbounds for i in 1:nrows
            c = alpha*x[i]
            if c != 0
                @simd for s in 1:S
                    k = k0 + s
                    j = J[k]
                    y[j] += C[k]*c
                end
            end
            k0 += S
        end
    end
    return y
end

"""

`AtWA(A,w)` yields the matrix `A'*W*A` from a sparse linear operator `A` and
weights `W = diag(w)`.

"""
function AtWA(A::SparseInterpolator{T,S,N},
              w::AbstractArray{T,N}) where {T,S,N}
    ncols = A.ncols
    AtWA!(Array{T}(undef, ncols, ncols), A, w)
end

"""

`AtA(A)` yields the matrix `A'*A` from a sparse linear operator `A`.

"""
function AtA(A::SparseInterpolator{T,S,N}) where {T,S,N}
    ncols = A.ncols
    AtA!(Array{T}(undef, ncols, ncols), A)
end

# Build the `A'*A` matrix from a sparse linear operator `A`.
function AtA!(dst::AbstractArray{T,2},
              A::SparseInterpolator{T,S,N}) where {T,S,N}
    nrows, ncols = A.nrows, A.ncols
    @assert size(dst) == (ncols, ncols)
    fill!(dst, zero(T))
    C, J = coefficients(A), columns(A)
    k0 = 0
    @assert length(J) == length(C)
    @inbounds for i in 1:nrows
        for s in 1:S
            k = k0 + s
            1 ≤ J[k] ≤ ncols || error("corrupted interpolator table")
        end
        for s1 in 1:S
            k1 = k0 + s1
            j1, c1 = J[k1], C[k1]
            @simd for s2 in 1:S
                k2 = k0 + s2
                j2, c2 = J[k2], C[k2]
                dst[j1,j2] += c1*c2
            end
        end
        k0 += S
    end
    return dst
end

# Build the `A'*W*A` matrix from a sparse linear operator `A` and weights `W`.
function AtWA!(dst::AbstractArray{T,2}, A::SparseInterpolator{T,S,N},
               wgt::AbstractArray{T,N}) where {T,S,N}
    nrows, ncols = A.nrows, A.ncols
    @assert size(dst) == (ncols, ncols)
    @assert size(wgt) == output_size(A)
    fill!(dst, zero(T))
    C, J = coefficients(A), columns(A)
    k0 = 0
    @assert length(J) == length(C)
    @inbounds for i in 1:nrows
        for s in 1:S
            k = k0 + s
            1 ≤ J[k] ≤ ncols || error("corrupted interpolator table")
        end
        w = wgt[i]
        for s1 in 1:S
            k1 = k0 + s1
            j1 = J[k1]
            wc1 = w*C[k1]
            @simd for s2 in 1:S
                k2 = k0 + s2
                j2 = J[k2]
                dst[j1,j2] += C[k2]*wc1
            end
        end
        k0 += S
    end
    return dst
end

# Default regularization levels.
const RGL_EPS = 1e-9
const RGL_MU = 0.0

"""
    fit(A, y [, w]; epsilon=1e-9, mu=0.0) -> x

performs a linear fit of `y` by the model `A*x` with `A` a linear interpolator.
The returned value `x` minimizes:

    sum(w.*(A*x - y).^2)

where `w` are given weights.  If `w` is not specified, all weights are assumed
to be equal to one; otherwise `w` must be an array of nonnegative values and of
same size as `y`.

Keywords `epsilon` and `mu` may be specified to regularize the solution and
minimize:

    sum(w.*(A*x - y).^2) + rho*(epsilon*norm(x)^2 + mu*norm(D*x)^2)

where `D` is a finite difference operator, `rho` is the maximum diagonal
element of `A'*diag(w)*A` and `norm` is the Euclidean norm.

"""
function fit(A::SparseInterpolator{T,S,N},
             y::AbstractArray{T,N},
             w::AbstractArray{T,N};
             epsilon::Real = RGL_EPS,
             mu::Real = RGL_MU) where {T,S,N}
    @assert size(y) == output_size(A)
    @assert size(w) == size(y)

    # Compute RHS vector A'*W*y with W = diag(w).
    rhs = A'*(w.*y)

    # Compute LHS matrix A'*W*A with W = diag(w).
    lhs = AtWA(A, w)

    # Regularize a bit.
    regularize!(lhs, epsilon, mu)

    # Solve the linear equations.
    cholfact!(lhs,:U,Val{true})\rhs
end

function fit(A::SparseInterpolator{T,S,N},
             y::AbstractArray{T,N};
             epsilon::Real = RGL_EPS,
             mu::Real = RGL_MU) where {T,S,N}
    @assert size(y) == output_size(A)
    @assert size(w) == size(y)

    # Compute RHS vector A'*y.
    rhs = A'*y

    # Compute LHS matrix A'*W*A with W = diag(w).
    lhs = AtA(A)

    # Regularize a bit.
    regularize!(lhs, epsilon, mu)

    # Solve the linear equations.
    cholfact!(lhs,:U,Val{true})\rhs
end

"""
    regularize(A, ϵ, μ) -> R

regularizes the symmetric matrix `A` to produce the matrix:

    R = A + ρ*(ϵ*I + μ*D'*D)

where `I` is the identity, `D` is a finite difference operator and `ρ` is the
maximum diagonal element of `A`.

"""
regularize(A::AbstractArray{T,2}, args...) where {T<:AbstractFloat} =
    regularize!(copyto!(Array{T}(undef, size(A)), A), args...)

"""
    regularize!(A, ϵ, μ) -> A

stores the regularized matrix in `A` (and returns it).  This is the in-place
version of [`LinearInterpolators.SparseInterpolators.regularize`].

"""
function regularize!(A::AbstractArray{T,2},
                     eps::Real = RGL_EPS,
                     mu::Real = RGL_MU) where {T<:AbstractFloat}
    regularize!(A, T(eps), T(mu))
end

function regularize!(A::AbstractArray{T,2},
                     eps::T, mu::T) where {T<:AbstractFloat}
    local rho::T
    @assert eps ≥ zero(T)
    @assert mu ≥ zero(T)
    @assert size(A,1) == size(A,2)
    n = size(A,1)
    if eps > zero(T) || mu > zero(T)
        rho = A[1,1]
        for j in 2:n
            d = A[j,j]
            rho = max(rho, d)
        end
        rho > zero(T) || error("we have a problem!")
    end
    if eps > zero(T)
        q = eps*rho
        for j in 1:n
            A[j,j] += q
        end
    end
    if mu > zero(T)
        q = rho*mu
        if n ≥ 2
            r = q + q
            A[1,1] += q
            A[2,1] -= q
            for i in 2:n-1
                A[i-1,i] -= q
                A[i,  i] += r
                A[i+1,i] -= q
            end
            A[n-1,n] -= q
            A[n,  n] += q
        elseif n == 1
            A[1,1] += q
        end
    end
    return A
end

# Yields a function that takes an index and returns the corresponding
# interpolation position as fractional index into the source array.
function fractional_index(T::Type{<:AbstractFloat},
                          pos::AbstractArray{<:Real},
                          grd::AbstractRange)
    # Use the central position of the grid to minimize rounding errors.
    c = (convert(T, first(grd)) + convert(T, last(grd)))/2
    q = 1/convert(T, step(grd))
    r = convert(T, 1 + length(grd))/2
    return i -> q*(convert(T, pos[i]) - c) + r
end

function fractional_index(T::Type{<:AbstractFloat},
                          pos::AbstractArray{<:Real})
    return i -> T(pos[i])
end

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

"""
    SparseUnidimensionalInterpolator{T<:AbstractFloat,S,D} <: AbstractSparseInterpolator{T}
* `T` is the floating-point type of the coefficients,
* `S` is the size of the kernel
  (number of nodes to combine for a single interpolator)
* `D` is the dimension of interpolation.
"""
struct SparseUnidimensionalInterpolator{T<:AbstractFloat,S,D} <: AbstractSparseInterpolator{T}
    nrows::Int     # number of rows
    ncols::Int     # number of columns
    C::Vector{T}   # coefficients along the dimension of interpolation
    J::Vector{Int} # columns indices along the dimension of interpolation
end

(A::SparseUnidimensionalInterpolator)(x) = apply(A, x)

interp_dim(::SparseUnidimensionalInterpolator{T,S,D}) where {T,S,D} = D

coefficients(A::SparseUnidimensionalInterpolator) = A.C
columns(A::SparseUnidimensionalInterpolator) = A.J
size(A::SparseUnidimensionalInterpolator) = (A.nrows, A.ncols)
size(A::SparseUnidimensionalInterpolator, i::Integer) =
    (i == 1 ? A.nrows :
     i == 2 ? A.ncols : error("out of bounds dimension"))

"""
    SparseUnidimensionalInterpolator{T=eltype(ker)}(ker, d, pos, grd)

yields a linear mapping which interpolates the `d`-th dimension of an array
with kernel `ker` at positions `pos` along the dimension of interpolation `d`
and assuming the input array has grid coordinates `grd` along the the `d`-th
dimension of interpolation.  Argument `pos` is a vector of positions, argument
`grd` may be a range or the length of the dimension of interpolation.  Optional
parameter `T` is the floating-point type of the coefficients of the operator.

This kind of interpolator is suitable for separable multi-dimensional
interpolation with precomputed interpolation coefficients.  Having precomputed
coefficients is mostly interesting when the operator is to be applied multiple
times (for instance in iterative methods).  Otherwise, separable operators
which compute the coefficients *on the fly* may be preferable.

A combination of instances of `SparseUnidimensionalInterpolator` can be built
to achieve sperable multi-dimensional interpolation.  For example:

    using LinearInterpolators
    ker = CatmullRomSpline()
    n1, n2 = 70, 50
    x1 = linspace(1, 70, 201)
    x2 = linspace(1, 50, 201)
    A1 = SparseUnidimensionalInterpolator(ker, 1, x1, 1:n1)
    A2 = SparseUnidimensionalInterpolator(ker, 2, x2, 1:n2)
    A = A1*A2

"""
SparseUnidimensionalInterpolator(ker::Kernel{T}, args...) where {T<:AbstractFloat} =
    SparseUnidimensionalInterpolator{T}(ker, args...)

@deprecate SparseUnidimensionalInterpolator(T::Type{<:AbstractFloat}, ker::Kernel, args...) SparseUnidimensionalInterpolator{T}(ker, args...)

SparseUnidimensionalInterpolator{T}(ker::Kernel, args...) where {T<:AbstractFloat} =
    SparseUnidimensionalInterpolator{T}(T(ker), args...)

function SparseUnidimensionalInterpolator{T}(ker::Kernel{T},
                                             d::Integer,
                                             pos::AbstractVector{<:Real},
                                             len::Integer) where {T<:AbstractFloat}
    len ≥ 1 || throw(ArgumentError("invalid dimension length"))
    return SparseUnidimensionalInterpolator{T}(ker, d, pos, 1:Int(len))
end

# FIXME: not type-stable
function SparseUnidimensionalInterpolator{T}(ker::Kernel{T,S},
                                             d::Integer,
                                             pos::AbstractVector{<:Real},
                                             grd::AbstractRange
                                             ) where {T<:AbstractFloat,S}
     SparseUnidimensionalInterpolator{T,S,Int(d)}(ker, pos, grd)
end

function SparseUnidimensionalInterpolator{T,S,D}(ker::Kernel{T,S},
                                                 pos::AbstractVector{<:Real},
                                                 grd::AbstractRange
                                                 ) where {D,T<:AbstractFloat,S}
    isa(D, Int) || throw(ArgumentError("invalid type for dimension of interpolation"))
    D ≥ 1 || throw(ArgumentError("invalid dimension of interpolation"))
    nrows = length(pos)
    ncols = length(grd)
    C, J = _sparsecoefs(CartesianIndices((nrows,)), ncols, ker,
                        fractional_index(T, pos, grd))
    return SparseUnidimensionalInterpolator{T,S,D}(nrows, ncols, C, J)
end

function vcreate(::Type{Direct},
                 A::SparseUnidimensionalInterpolator,
                 x::AbstractArray,
                 scratch::Bool=false)
    nrows, ncols = size(A)
    return _vcreate(nrows, ncols, A, x)
end

function vcreate(::Type{Adjoint},
                 A::SparseUnidimensionalInterpolator,
                 x::AbstractArray,
                 scratch::Bool=false)
    nrows, ncols = size(A)
    return _vcreate(ncols, nrows, A, x)
end

function _vcreate(ny::Int, nx::Int,
                  A::SparseUnidimensionalInterpolator{Ta,S,D},
                  x::AbstractArray{Tx,N}) where {Ta,Tx<:Real,S,D,N}
    xdims = size(x)
    1 ≤ D ≤ N ||
        throw(DimensionMismatch("out of range dimension of interpolation"))
    xdims[D] == nx ||
        throw(DimensionMismatch("dimension $D of `x` must be $nx"))
    Ty = float(promote_type(Ta, Tx))
    ydims = [(d == D ? ny : xdims[d]) for d in 1:N]
    return Array{Ty,N}(undef, ydims...)
end

function apply!(α::Real, ::Type{Direct},
                A::SparseUnidimensionalInterpolator{Ta,S,D},
                x::AbstractArray{Tx,N},
                scratch::Bool,
                β::Real,
                y::AbstractArray{Ty,N}) where {Ta<:AbstractFloat,
                                               Tx<:Real,
                                               Ty<:AbstractFloat,S,D,N}
    # Check arguments.
    _check(A, N)
    xdims = size(x)
    ydims = size(y)
    nrows, ncols = size(A)
    xdims[D] == ncols ||
        throw(DimensionMismatch("dimension $D of `x` must be $ncols"))
    ydims[D] == nrows ||
        throw(DimensionMismatch("dimension $D of `y` must be $nrows"))
    for k in 1:N
        k == D || xdims[k] == ydims[k] ||
            throw(DimensionMismatch("`x` and `y` have incompatible dimensions"))
    end

    # Apply operator.
    if α == 0
        vscale!(y, β)
    else
        C = coefficients(A)
        J = columns(A)
        I_pre = CartesianIndices(xdims[1:D-1])
        I_post = CartesianIndices(xdims[D+1:N])
        T = promote_type(Ta,Tx)
        alpha = convert(T, α)
        if β == 0
            _apply_direct!(T, Val{S}, C, J, alpha, x, y,
                           I_pre, nrows, I_post)
        else
            beta = convert(Ty, β)
            _apply_direct!(T, Val{S}, C, J, alpha, x, beta, y,
                           I_pre, nrows, I_post)
        end
    end
    return y
end

function apply!(α::Real, ::Type{Adjoint},
                A::SparseUnidimensionalInterpolator{Ta,S,D},
                x::AbstractArray{Tx,N},
                scratch::Bool,
                β::Real,
                y::AbstractArray{Ty,N}) where {Ta<:AbstractFloat,
                                               Tx<:Real,
                                               Ty<:AbstractFloat,S,D,N}
    # Check arguments.
    _check(A, N)
    xdims = size(x)
    ydims = size(y)
    nrows, ncols = size(A)
    xdims[D] == nrows ||
        throw(DimensionMismatch("dimension $D of `x` must be $nrows"))
    ydims[D] == ncols ||
        throw(DimensionMismatch("dimension $D of `y` must be $ncols"))
    for k in 1:N
        k == D || xdims[k] == ydims[k] ||
            throw(DimensionMismatch("`x` and `y` have incompatible dimensions"))
    end

    # Apply adjoint operator.
    vscale!(y, β)
    if α != 0
        T = promote_type(Ta,Tx)
        _apply_adjoint!(Val{S}, coefficients(A), columns(A),
                        convert(T, α), x, y,
                        CartesianIndices(xdims[1:D-1]), nrows,
                        CartesianIndices(xdims[D+1:N]))
    end
    return y
end

# The 3 following private methods are needed to achieve type invariance and win
# a factor ~1000 in speed!  Also note the way the innermost loop is written
# with a constant range and an offset k0 which is updated; this is critical for
# saving a factor 2-3 in speed.
#
# The current version takes ~ 4ms (7 iterations of linear conjugate gradients)
# to fit a 77×77 array of weights interpolated by Catmull-Rom splines to
# approximate a 256×256 image.

function _apply_direct!(::Type{T},
                        ::Type{Val{S}},
                        C::Vector{<:AbstractFloat},
                        J::Vector{Int},
                        α::AbstractFloat,
                        x::AbstractArray{<:Real,N},
                        y::AbstractArray{<:AbstractFloat,N},
                        I_pre::CartesianIndices{N_pre},
                        len::Int,
                        I_post::CartesianIndices{N_post}
                        ) where {T<:AbstractFloat,S,N,N_post,N_pre}
    @assert N == N_post + N_pre + 1
    @inbounds for i_post in I_post
        for i_pre in I_pre
            k0 = 0
            for i in 1:len
                sum = zero(T)
                @simd for s in 1:S
                    k = k0 + s
                    sum += C[k]*x[i_pre,J[k],i_post]
                end
                y[i_pre,i,i_post] = α*sum
                k0 += S
            end
        end
    end
end

function _apply_direct!(::Type{T},
                        ::Type{Val{S}},
                        C::Vector{<:AbstractFloat},
                        J::Vector{Int},
                        α::AbstractFloat,
                        x::AbstractArray{<:Real,N},
                        β::AbstractFloat,
                        y::AbstractArray{<:AbstractFloat,N},
                        I_pre::CartesianIndices{N_pre},
                        len::Int,
                        I_post::CartesianIndices{N_post}
                        ) where {T<:AbstractFloat,S,N,N_post,N_pre}
    @assert N == N_post + N_pre + 1
    @inbounds for i_post in I_post
        for i_pre in I_pre
            k0 = 0
            for i in 1:len
                sum = zero(T)
                @simd for s in 1:S
                    k = k0 + s
                    sum += C[k]*x[i_pre,J[k],i_post]
                end
                y[i_pre,i,i_post] = α*sum + β*y[i_pre,i,i_post]
                k0 += S
            end
        end
    end
end

function _apply_adjoint!(::Type{Val{S}},
                         C::Vector{<:AbstractFloat},
                         J::Vector{Int},
                         α::AbstractFloat,
                         x::AbstractArray{<:Real,N},
                         y::AbstractArray{<:AbstractFloat,N},
                         I_pre::CartesianIndices{N_pre},
                         len::Int,
                         I_post::CartesianIndices{N_post}
                         ) where {S,N,N_post,N_pre}
    @assert N == N_post + N_pre + 1
    @inbounds for i_post in I_post
        for i_pre in I_pre
            k0 = 0
            for i in 1:len
                c = α*x[i_pre,i,i_post]
                @simd for s in 1:S
                    k = k0 + s
                    y[i_pre,J[k],i_post] += C[k]*c
                end
                k0 += S
            end
        end
    end
end

function _check(A::SparseUnidimensionalInterpolator{T,S,D},
                N::Int) where {T<:AbstractFloat,S,D}
    1 ≤ D ≤ N ||
        throw(DimensionMismatch("out of range dimension of interpolation"))
    nrows, ncols = size(A)
    nvals = S*nrows
    C = coefficients(A)
    J = columns(A)
    length(C) == nvals ||
        throw(DimensionMismatch("array of coefficients must have $nvals elements (has $(length(C)))"))
    length(J) == nvals ||
        throw(DimensionMismatch("array of indices must have $nvals elements (has $(length(C)))"))
    for k in eachindex(J)
        1 ≤ J[k] ≤ ncols || throw(ErrorException("out of bounds indice(s)"))
    end
end

end # module


text:
function _power_flow_solution!(
    initial_guess::Vector{Float64},
    sys::PSY.System,
    inputs::SimulationInputs,
)
    res = PSY.solve_powerflow!(sys)
    if !res
        @error("PowerFlow failed to solve")
        return BUILD_FAILED
    end
    bus_size = length(PSY.get_bus_numbers(sys))
    @debug "Updating Voltage guess"
    for bus in PSY.get_components(PSY.Bus, sys)
        @debug PSY.get_name(bus)
        @debug "V_r" PSY.get_magnitude(bus) * cos(PSY.get_angle(bus))
        @debug "V_i" PSY.get_magnitude(bus) * sin(PSY.get_angle(bus))
        bus_n = PSY.get_number(bus)
        bus_ix = get_lookup(inputs)[bus_n]
        initial_guess[bus_ix] = PSY.get_magnitude(bus) * cos(PSY.get_angle(bus))
        initial_guess[bus_ix + bus_size] = PSY.get_magnitude(bus) * sin(PSY.get_angle(bus))
    end
    return BUILD_INCOMPLETE
end

function _initialize_static_injection!(inputs::SimulationInputs)
    @debug "Updating Source internal voltages"
    static_injection_devices = get_static_injections_data(inputs)
    if !isempty(static_injection_devices)
        try
            for s in static_injection_devices
                initialize_static_device!(s)
            end
        catch e
            bt = catch_backtrace()
            @error "Static Injection Failed to Initialize" exception = e, bt
            return BUILD_FAILED
        end
    end
    return BUILD_INCOMPLETE
end

function _initialize_dynamic_injection!(
    initial_guess::Vector{Float64},
    inputs::SimulationInputs,
)
    @debug "Updating Dynamic Injection Component Initial Guess"
    injection_start = get_injection_pointer(inputs)
    try
        for d in get_injectors_data(inputs)
            dynamic_device = PSY.get_dynamic_injector(d)
            @debug PSY.get_name(d) typeof(d) typeof(dynamic_device)
            n_states = PSY.get_n_states(dynamic_device)
            ix_range = range(injection_start, length = n_states)
            injection_start = injection_start + n_states
            x0_device = initialize_dynamic_device!(dynamic_device, d)
            @assert length(x0_device) == n_states
            initial_guess[ix_range] = x0_device
        end
    catch e
        bt = catch_backtrace()
        @error "Dynamic Injection Failed to Initialize" exception = e, bt
        return BUILD_FAILED
    end
    return BUILD_INCOMPLETE
end

function _initialize_dynamic_branches!(
    initial_guess::Vector{Float64},
    inputs::SimulationInputs,
)
    @debug "Updating Component Initial Guess"
    branches_start = get_branches_pointer(inputs)
    try
        for br in get_dynamic_branches(inputs)
            @debug PSY.get_name(br) typeof(br)
            n_states = PSY.get_n_states(br)
            ix_range = range(branches_start, length = n_states)
            branches_start = branches_start + n_states
            x0_branch = initialize_dynamic_device!(br)
            @assert length(x0_branch) == n_states
            initial_guess[ix_range] = x0_branch
        end
    catch e
        bt = catch_backtrace()
        @error "Dynamic Branches Failed to Initialize" exception = e, bt
        return BUILD_FAILED
    end
    return BUILD_INCOMPLETE
end

function _check_valid_values(initial_guess::Vector{Float64}, inputs::SimulationInputs)
    if any(!isfinite, initial_guess)
        i = findall(!isfinite, initial_guess)
        invalid_initial_guess = String[]
        for (device, states) in get_global_index(inputs)
            for state in states
                if state.second ∈ i
                    push!(invalid_initial_guess, "$device - $(p.first)")
                end
            end
        end
        @error("Invalid initial guess values $invalid_initial_guess")
        return BUILD_FAILED
    end
    return BUILD_INCOMPLETE
end

struct NLsolveWrapper
    zero::Vector{Float64}
    converged::Bool
    failed::Bool
end

NLsolveWrapper() = NLsolveWrapper(Vector{Float64}(), false, true)
converged(sol::NLsolveWrapper) = sol.converged
failed(sol::NLsolveWrapper) = sol.failed

function _nlsolve_call(initial_guess, inputs::SimulationInputs, tolerance::Float64)
    dx0 = zeros(length(initial_guess)) #Define a vector of zeros for the derivative
    inif! = (out, x) -> system_implicit!(
        out,    #output of the function
        dx0,    #derivatives equal to zero
        x,      #states
        inputs,    #Parameters
        -99.0,    #time val not relevant
    )
    return try
        sys_solve = NLsolve.nlsolve(
            inif!,
            initial_guess,
            xtol = tolerance,
            ftol = tolerance,
            method = :trust_region,
            autodiff = :forward,
        ) #Solve using initial guess x0
        NLsolveWrapper(sys_solve.zero, NLsolve.converged(sys_solve), false)
    catch e
        bt = catch_backtrace()
        @error "NLsolve failed to solve" exception = e, bt
        NLsolveWrapper()
    end
end

function _convergence_check(sys_solve::NLsolveWrapper, tol::Float64)
    if converged(sys_solve)
        @info("Initialization succeeded with a tolerance of $(tol). Saving solution")
    else
        @warn(
            "Initialization convergence failed, initial conditions do not meet conditions for an stable equilibrium.\nTrying to solve again reducing numeric tolerance"
        )
    end
    return converged(sys_solve)
end

function _refine_initial_condition!(
    initial_guess::Vector{Float64},
    inputs::SimulationInputs,
)
    @debug "Start NLSolve System Run"
    @debug initial_guess
    converged = false
    for tol in [STRICT_NL_SOLVE_TOLERANCE, RELAXED_NL_SOLVE_TOLERANCE]
        sys_solve = _nlsolve_call(initial_guess, inputs, tol)
        failed(sys_solve) && return BUILD_FAILED
        converged = _convergence_check(sys_solve, tol)
        @debug "Write result to initial guess vector under condition converged = $(converged)"
        initial_guess .= sys_solve.zero
        if converged
            break
        end
        @warn(
            "Initialization never converged to desired tolerances. Initial conditions do not meet conditions for an stable equilibrium. Simulation migth diverge"
        )
    end
    return converged ? SIMULATION_INITIALIZED : _check_valid_values(initial_guess, inputs)
end

# Default implementation for both models. This implementation is to future proof if there is
# a divergence between the required build methods
function _calculate_initial_conditions!(sim::Simulation)
    inputs = get_simulation_inputs(sim)
    simulation_pre_step!(inputs, Real)
    @debug "Start state intialization routine"
    while sim.status == BUILD_INCOMPLETE
        sim.status = _power_flow_solution!(sim.x0_init, get_system(sim), inputs)
        sim.status = _initialize_static_injection!(inputs)
        sim.status = _initialize_dynamic_injection!(sim.x0_init, inputs)
        if has_dyn_lines(inputs)
            sim.status = _initialize_dynamic_branches!(sim.x0_init, inputs)
        else
            @debug "No Dynamic Branches in the system"
        end
        sim.status = _check_valid_values(sim.x0_init, inputs)
        sim.status = _refine_initial_condition!(sim.x0_init, inputs)
    end
    return sim.status != BUILD_FAILED
end

function calculate_initial_conditions!(sim::Simulation{ImplicitModel})
    return _calculate_initial_conditions!(sim)
end

function calculate_initial_conditions!(sim::Simulation{MassMatrixModel})
    return _calculate_initial_conditions!(sim)
end

"""
Returns a Dictionary with the resulting initial conditions of the simulation
"""
function get_initial_conditions(sim::Simulation)
    system = get_system(sim)
    simulation_inputs = get_simulation_inputs(sim)
    bus_size = get_bus_count(simulation_inputs)
    V_R = Dict{Int, Float64}()
    V_I = Dict{Int, Float64}()
    Vm = Dict{Int, Float64}()
    θ = Dict{Int, Float64}()
    for bus in PSY.get_components(PSY.Bus, system)
        bus_n = PSY.get_number(bus)
        bus_ix = get_lookup(simulation_inputs)[bus_n]
        V_R[bus_n] = sim.x0_init[bus_ix]
        V_I[bus_n] = sim.x0_init[bus_ix + bus_size]
        Vm[bus_n] = sqrt(sim.x0_init[bus_ix]^2 + sim.x0_init[bus_ix + bus_size]^2)
        θ[bus_n] = angle(sim.x0_init[bus_ix] + sim.x0_init[bus_ix + bus_size] * 1im)
    end
    results = Dict{String, Any}("V_R" => V_R, "V_I" => V_I, "Vm" => Vm, "θ" => θ)
    for device in PSY.get_components(PSY.DynamicInjection, system)
        states = PSY.get_states(device)
        name = PSY.get_name(device)
        global_index = get_global_index(simulation_inputs)[name]
        x0_device = Dict{Symbol, Float64}()
        for s in states
            x0_device[s] = sim.x0_init[global_index[s]]
        end
        results[name] = x0_device
    end
    dyn_branches = PSY.get_components(PSY.DynamicBranch, system)
    if !isempty(dyn_branches)
        for br in dyn_branches
            states = PSY.get_states(br)
            name = PSY.get_name(br)
            global_index = get_global_index(simulation_inputs)[name]
            x0_br = Dict{Symbol, Float64}()
            for s in states
                x0_br[s] = sim.x0_init[global_index[s]]
            end
            printed_name = "Line " * name
            results[printed_name] = x0_br
        end
    end
    return results
end


text:
#####
##### `sum`
#####

function frule((_, ẋ), ::typeof(sum), x; dims=:)
    return sum(x; dims=dims), sum(ẋ; dims=dims)
end

function rrule(::typeof(sum), x::AbstractArray{T}; dims=:) where {T<:Number}
    y = sum(x; dims=dims)
    function sum_pullback(ȳ)
        # broadcasting the two works out the size no-matter `dims`
        x̄ = broadcast(x, ȳ) do xi, ȳi
            ȳi
        end
        return (NO_FIELDS, x̄)
    end
    return y, sum_pullback
end

function frule(
    (_, _, ẋ),
    ::typeof(sum),
    ::typeof(abs2),
    x::AbstractArray{T};
    dims=:,
) where {T<:Union{Real,Complex}}
    y = sum(abs2, x; dims=dims)
    ∂y = if dims isa Colon
        2 * real(dot(x, ẋ))
    elseif VERSION ≥ v"1.2" # multi-iterator mapreduce introduced in v1.2
        2 * mapreduce(_realconjtimes, +, x, ẋ; dims=dims)
    else
        2 * sum(_realconjtimes.(x, ẋ); dims=dims)
    end
    return y, ∂y
end

function rrule(
    ::typeof(sum),
    ::typeof(abs2),
    x::AbstractArray{T};
    dims=:,
) where {T<:Union{Real,Complex}}
    y = sum(abs2, x; dims=dims)
    function sum_abs2_pullback(ȳ)
        return (NO_FIELDS, DoesNotExist(), 2 .* real.(ȳ) .* x)
    end
    return y, sum_abs2_pullback
end


text:
module PoissonLagrangeMultiplierTests

using Test
using Gridap
using Gridap.Geometry

function is_left(coords)
  n = length(coords)
  x = (1/n)*sum(coords)
  ( x[1] <= 0.5 )
end

domain = (0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)

Ω = Interior(model)
Γ = Boundary(model)
Λ = Skeleton(Γ)

Γface_to_coords = get_cell_coordinates(Γ)
Γface_mask = lazy_map(is_left,Γface_to_coords)
Γlface_Γface = findall(Γface_mask)
Γrface_Γface = findall(!,Γface_mask)
Γl = Triangulation(Γ,Γlface_Γface)
Γr = Triangulation(Γ,Γrface_Γface)

get_glue(Ω,Val(0))
get_glue(Γ,Val(0))
get_glue(Λ,Val(0))
get_glue(Γl,Val(0))
get_glue(Γr,Val(0))

order = 2
reffe_u = ReferenceFE(lagrangian,Float64,order)
reffe_λ = ReferenceFE(lagrangian,Float64,order-1)
V = TestFESpace(Ω,reffe_u,conformity=:H1)
S = TestFESpace(Γ,reffe_λ,conformity=:L2)
U = TrialFESpace(V)
L = TrialFESpace(S)
Y = MultiFieldFESpace([V,S])
X = MultiFieldFESpace([U,L])

#uh, λh = FEFunction(Y,rand(num_free_dofs(Y)))
#writevtk(Ω,"t_Ω",nsubcells=10,cellfields=["uh"=>uh])
#writevtk(Γ,"t_Γ",nsubcells=10,cellfields=["uh"=>uh,"λh"=>λh])
#writevtk(Λ,"t_Λ",cellfields=["uh"=>mean(uh),"λh"=>mean(λh)])
#writevtk(Γr,"t_Γr",nsubcells=10,cellfields=["uh"=>uh,"λh"=>λh])
#writevtk(Γl,"t_Γl",nsubcells=10,cellfields=["uh"=>uh,"λh"=>λh])

degree = 2*order
dΩ = Measure(Ω,degree)
dΓl = Measure(Γl,degree)
dΓr = Measure(Γr,degree)
dΛ = Measure(Λ,degree)
n = get_normal_vector(Γ)
nΛ = get_normal_vector(Λ)

uₑ(x) = x[1]^2 + x[2]^2
f(x) = -Δ(uₑ)(x)

# Weak form. Additional non needed terms are added for testing purposes
a((u,λ),(v,η)) = ∫( ∇(u)⋅∇(v) )dΩ +
               ∫( (u+λ)*(η+v) - u*v - η*λ + 0.0*(n⋅∇(u)-λ)*(n⋅∇(v)-η) )dΓl +
               ∫( (u+λ)*(η+v) - u*v - η*λ + 0.0*(n⋅∇(u)-λ)*(n⋅∇(v)-η) )dΓr +
               ∫( 0.0*jump(n⋅∇(u))*mean(η) )dΛ +
               ∫( 0.0*jump(nΛ⋅∇(u))*mean(η) )dΛ
l((v,η)) = ∫( f*v )dΩ + ∫( uₑ*η )dΓl + ∫( uₑ*η )dΓr

op = AffineFEOperator(a,l,X,Y)
uh,λh = solve(op)

e = uh-uₑ
l2(v) = √(∑(∫(v*v)dΩ))

tol = 1.0e-10
@test l2(e) < tol

end


text:
abstract type Transmission end

struct ExogenousTransmission <: Transmission
  individual::Int64
end

struct EndogenousTransmission <: Transmission
  individual::Int64
  source::Int64
end

struct NoTransmission <: Transmission
end

function Base.show(io::IO, x::ExogenousTransmission)
  return println(io, "Exogenous infection transmission to individual $(x.individual)")
end

function Base.show(io::IO, x::EndogenousTransmission)
  return println(io, "Endogenous infection transmission to individual $(x.individual) by individual $(x.source)")
end

function Base.show(io::IO, x::NoTransmission)
  return println(io, "No infection transmission")
end


text:
#! /usr/bin/julia

# Rosetta Code, User input/Text

print("String? ")
y = chomp(readline())
println("Your input was \"", y, "\".\n")
print("Integer? ")
y = chomp(readline())
try
    y = parseint(y)
    println("Your input was \"", y, "\".\n")
catch
    println("Sorry, but \"", y, "\" does not compute as an integer.")
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code