text:
struct StateSpaceNoise{T, MT<:AbstractMatrix{T}} <: ControlSystems.AbstractStateSpace
	A::MT
	B::MT
	K::MT
	Ts::Float64
	nx::Int
	nu::Int
	ny::Int
	function StateSpaceNoise(A::MT, B::MT, K::MT, Ts::Float64) where MT
		nx,nu,ny = ControlSystems.state_space_validation(A,B,K',zeros(size(K',1), size(B,2)),Ts)
		new{eltype(A), typeof(A)}(A, B, K, Ts, nx, nu, ny)
	end
end

ControlSystems.isstable(s::StateSpaceNoise) = all(abs(e) <= 1 for e in eigvals(s.A-s.K*s.C))

# Funtions for number of intputs, outputs and states
# ControlSystems.ninputs(sys::StateSpaceNoise) = sys.nu
# ControlSystems.noutputs(sys::StateSpaceNoise) = sys.ny
# ControlSystems.nstates(sys::StateSpaceNoise) = sys.nx
#
# Base.ndims(::StateSpaceNoise) = 2 # NOTE: Also for SISO systems?
# Base.size(sys::StateSpaceNoise) = (noutputs(sys), ninputs(sys)) # NOTE: or just size(get_D(sys))
# Base.size(sys::StateSpaceNoise, d) = d <= 2 ? size(sys)[d] : 1
Base.eltype(::Type{S}) where {S<:StateSpaceNoise} = S
Base.convert(::Type{StateSpace}, sys::StateSpaceNoise) = ss(sys.A, sys.B, sys.C, 0, sys.Ts)
ControlSystems.ss(sys::StateSpaceNoise) = convert(StateSpace,sys)
ControlSystems.tf(sys::StateSpaceNoise) = tf(ss(sys))

function Base.getproperty(sys::StateSpaceNoise, p::Symbol)
	if p == :C
		return [I zeros(sys.ny,sys.nx-sys.ny)]
	elseif p == :D
		return zeros(sys.ny, sys.nu)
	end
	return getfield(sys,p)
end

ControlSystems.innovation_form(sys::StateSpaceNoise) = ss(sys.A, sys.K, sys.C, Matrix(Eye(sys.ny)), sys.Ts) # innovation model

function Base.getindex(sys::StateSpaceNoise, inds...)
	if size(inds, 1) != 2
		error("Must specify 2 indices to index statespace model")
	end
	rows, cols = ControlSystems.index2range(inds...) # FIXME: ControlSystems.index2range(inds...)
	return StateSpaceNoise(copy(sys.A), sys.B[:, cols], sys.K[:, rows], sys.Ts)
end

struct SysFilter{T<:Union{StateSpaceNoise, StateSpace}, FT}
	sys::T
	state::Vector{FT}
	yh::Vector{FT}
end
SysFilter(sys::LTISystem,x0=zeros(sys.nx)) = SysFilter(sys,x0,zeros(eltype(x0), sys.ny))

(s::SysFilter)(y, u) = sysfilter!(s.state, s.sys, y, u)
(s::SysFilter)(u) = sysfilter!(s.state, s.sys, u)
sysfilter!(s::SysFilter, y, u) = sysfilter!(s.state, s.sys, y, u)
sysfilter!(s::SysFilter, u) = sysfilter!(s.state, s.sys, u)

function sysfilter!(state::AbstractVector, sys::StateSpaceNoise, y, u)
	@unpack A,B,K,ny = sys
	yh     = state[1:ny] #vec(sys.C*state)
	e      = y - yh
	state .= vec(A*state + B*u + K*e)
	yh
end

function sysfilter!(state::AbstractVector, sys::StateSpaceNoise, u)
	@unpack A,B,K,ny = sys
	yh     = state[1:ny] #vec(C*state)
	state .= vec(A*state + B*u)
	yh
end

function sysfilter!(state::AbstractVector, sys::StateSpace, y, u)
	@unpack A,B,C,D = sys
	yh     = vec(C*state + D*u)
	state .= vec(A*state + B*u)
	yh
end

function sysfilter!(state::AbstractVector, sys::StateSpace, u)
	@unpack A,B,C,D = sys
	yh     = vec(C*state + D*u)
	state .= vec(A*state + B*u)
	yh
end

struct OberservationIterator{T}
	y::T
	u::T
end

observations(y,u) = OberservationIterator(y,u)

function Base.iterate(i::OberservationIterator{<:AbstractMatrix}, state=1)
	state > length(i) && return nothing
	((i.y[:,state],i.u[:,state]),state+1)
end
Base.length(i::OberservationIterator{<:AbstractMatrix}) = size(i.y, 2)

function Base.iterate(i::OberservationIterator{<:AbstractVector{<:Union{AbstractVector, Number}}}, state=1)
	state > length(i) && return nothing
	((i.y[state],i.u[state]),state+1)
end
Base.length(i::OberservationIterator{<:AbstractVector{<:Union{AbstractVector, Number}}}) = length(i.y)

struct PredictionErrorIterator{T,MT}
	model::MT
	oi::OberservationIterator{T}
end

prediction_errors(model, y, u) = PredictionErrorIterator(model, observations(y,u))

function Base.iterate(i::PredictionErrorIterator, state=1)
	state >= length(i) && return nothing
	(y,u), state1 = iterate(i.oi, state)
	yh = i.model(y,u)
	(y-yh,state1)
end
Base.length(i::PredictionErrorIterator) = length(i.oi)

mutable struct SimulationErrorIterator{T,MT,YT}
	model::MT
	oi::OberservationIterator{T}
	yh::YT # Stores the last prediction
end

simulation_errors(model, y, u) = SimulationErrorIterator(model, observations(y,u), zeros(eltype(model.sys.A), obslength(y)))

function Base.iterate(i::SimulationErrorIterator, state=1)
	state >= length(i) && return nothing
	(y,u), state1 = iterate(i.oi, state)
	i.yh = i.model(u)
	(y-i.yh,state1)
end
Base.length(i::SimulationErrorIterator) = length(i.oi)


text:
import MathOptInterface
import MathOptInterface: Utilities

const MOI = MathOptInterface
const MOIU = MOI.Utilities

# indices
const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex

# function aliases
const SV = MOI.SingleVariable
const SAF = MOI.ScalarAffineFunction{Float64}
const SQF = MOI.ScalarQuadraticFunction{Float64}

# set aliases
const Bounds{T} = Union{
    MOI.EqualTo{T},
    MOI.GreaterThan{T},
    MOI.LessThan{T},
    MOI.Interval{T}
}

mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::BaronModel
    nlp_block_data::Union{Nothing, MOI.NLPBlockData}
end

Optimizer(; options...) = Optimizer(BaronModel(; options...), nothing)

MOI.get(model::Optimizer, ::MOI.SolverName) = "BARON"

# empty
function MOI.is_empty(model::Optimizer)
    BARON.is_empty(model.inner) && model.nlp_block_data === nothing
end

function MOI.empty!(model::Optimizer)
    model.inner = BaronModel(; model.inner.options...)
    model.nlp_block_data = nothing
end

# copy
MOIU.supports_default_copy_to(model::Optimizer, copy_names::Bool) = !copy_names
function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; kws...)
    return MOIU.automatic_copy_to(dest, src; kws...)
end

# optimize
function MOI.optimize!(model::Optimizer)
    write_bar_file(model.inner)
    run(`$baron_exec $(model.inner.problem_file_name)`)
    read_results(model.inner)
end

# TimeLimitSec
MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, val::Real)
    model.inner.options[:MaxTime] = Float64(val)
    return
end

# BARON's default time limit is 1000 seconds.
function MOI.get(model::Optimizer, ::MOI.TimeLimitSec)
    return get(model.inner.options, :MaxTime, 1000.0)
end

include(joinpath("moi", "util.jl"))
include(joinpath("moi", "variables.jl"))
include(joinpath("moi", "constraints.jl"))
include(joinpath("moi", "objective.jl"))
include(joinpath("moi", "nlp.jl"))
include(joinpath("moi", "results.jl"))

ext: jl


text:
# ---------------------------------------------------------------------------- #
#
#   quadratureTriangle.jl
#
#   Get Legendre quadrature points for Triangle
#
#   λυτέος
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

function quadratureTriangle( P::PG1 )

  nq = 1
  weight = Array{Float64}(nq)
  xc     = Array{Float64}(nq,1)
  yc     = Array{Float64}(nq,1)

  weight[1] = 1
  xc[1] = 1./3.
  yc[1] = 1./3.

  pts = hcat( xc, yc )

  weight = weight/2

  return pts, weight

end

function quadratureTriangle( P::PG2 )

  nq = 3

  weight = Array{Float64}(nq)
  xc     = Array{Float64}(nq,1)
  yc     = Array{Float64}(nq,1)

  weight[1] = 1./3.
  xc[1]     = 1./6.
  yc[1]     = 2./3.

  weight[2] = 1./3.
  xc[2]     = 2./3.
  yc[2]     = 1./6.

  weight[3] = 1./3.
  xc[3]     = 1./6.
  yc[3]     = 1./6.

  pts = hcat( xc, yc )

  weight = weight/2

  return pts, weight

end

function quadratureTriangle( P::PG3 )

  nq = 6
  weight = Array{Float64}(nq)
  xc     = Array{Float64}(nq,1)
  yc     = Array{Float64}(nq,1)

  w = 2./12.
  a = 0.65902762237409221517838077125540
  b = 0.23193336855303057249678456117469
  c = 0.10903900907287721232483466756991

  weight[1] = w
  xc[1]     = a
  yc[1]     = b

  weight[2] = w
  xc[2]     = b
  yc[2]     = c

  weight[3] = w
  xc[3]     = c
  yc[3]     = a

  weight[4] = w
  xc[4]     = a
  yc[4]     = c

  weight[5] = w
  xc[5]     = b
  yc[5]     = a

  weight[6] = w
  xc[6]     = c
  yc[6]     = b

  pts = hcat( xc, yc )

  weight = weight/2

  return pts, weight

end

function quadratureTriangle( P::PG4 )

  nq = 6
  weight = Array{Float64}(nq)
  xc     = Array{Float64}(nq,1)
  yc     = Array{Float64}(nq,1)

  x = 4./9. - sqrt(10)/18. + sqrt(950 - 220*sqrt(10))/90.
  w = 2.*(1./12. + 3*sqrt(95 - 22*sqrt(10))/496. - sqrt(950 - 220*sqrt(10))/7440.)
  xc[1]     = x
  yc[1]     = x
  weight[1] = w

  xc[2]     = x
  yc[2]     = 1 - 2*x
  weight[2] = w

  xc[3]     = 1 - 2*x
  yc[3]     = x
  weight[3] = w

  x         = 4./9. - sqrt(10)/18. - sqrt(950 - 220*sqrt(10))/90.
  w         = 2.*(1./12. - 3*sqrt(95 - 22*sqrt(10))/496. + sqrt(950 - 220*sqrt(10))/7440.)
  xc[4]     = x
  yc[4]     = x
  weight[4] = w

  xc[5]     = x
  yc[5]     = 1 - 2*x
  weight[5] = w

  xc[6]     = 1 - 2*x
  yc[6]     = x
  weight[6] = w

  pts = hcat( xc, yc )

  weight = weight/2

  return pts, weight

end

function quadratureTriangle( P::PG5 )

  nq = 7
  weight = Array{Float64}(nq)
  xc     = Array{Float64}(nq,1)
  yc     = Array{Float64}(nq,1)

  x         = 1./3.
  w         = 2.*9./80.
  xc[1]     = x
  yc[1]     = x
  weight[1] = w

  x         = 2./7. + sqrt(15)/21.
  w         = 2.*(31./480. + sqrt(15)/2400.)
  xc[2]     = x
  yc[2]     = x
  weight[2] = w

  xc[3]     = x
  yc[3]     = 1 - 2*x
  weight[3] = w

  xc[4]     = 1 - 2*x
  yc[4]     = x
  weight[4] = w

  x         = 2./7. - sqrt(15)/21.
  w         = 2.*(31./480. - sqrt(15)/2400.)
  xc[5]     = x
  yc[5]     = x
  weight[5] = w

  xc[6]     = x
  yc[6]     = 1 - 2*x
  weight[6] = w

  xc[7]     = 1 - 2*x
  yc[7]     = x
  weight[7] = w

  pts = hcat( xc, yc )

  weight = weight/2

  return pts, weight

end

function quadratureTriangle( P::PG8 )

  nq = 16
  weight = Array{Float64}(nq)
  xc     = Array{Float64}(nq,1)
  yc     = Array{Float64}(nq,1)

  x         = 1./3.
  w         = 2.*0.072157803838893584125545555244532
  xc[1]     = x
  yc[1]     = x
  weight[1] = w

  x         = 0.17056930775176020662229350149146
  w         = 2.*0.051608685267359125140895775146065
  xc[2]     = x
  yc[2]     = x
  weight[2] = w

  xc[3]     = x
  yc[3]     = 1 - 2*x
  weight[3] = w

  xc[4]     = 1 - 2*x
  yc[4]     = x
  weight[4] = w

  x         = 0.45929258829272315602881551449417
  w         = 2.*0.047545817133642312396948052194292
  xc[5]     = x
  yc[5]     = x
  weight[5] = w

  xc[6]     = x
  yc[6]     = 1 - 2*x
  weight[6] = w

  xc[7]     = 1 - 2*x
  yc[7]     = x
  weight[7] = w

  x         = 0.050547228317030975458423550596599
  w         = 2.*0.016229248811599040155462964170890
  xc[8]     = x
  yc[8]     = x
  weight[8] = w

  xc[9]     = x
  yc[9]     = 1 - 2*x
  weight[9] = w

  xc[10]     = 1 - 2*x
  yc[10]     = x
  weight[10] = w

  x          = 0.72849239295540428124100037917606
  y          = 0.26311282963463811342178578628464
  w          = 2.*0.0136151570872174971324223450369545
  xc[11]     = x
  yc[11]     = y
  weight[11] = w

  xc[12]     = y
  yc[12]     = x
  weight[12] = w

  xc[13]     = x
  yc[13]     = 1 - x - y
  weight[13] = w

  xc[14]     = 1 - x - y
  yc[14]     = x
  weight[14] = w

  xc[15]     = y
  yc[15]     = 1 - x - y
  weight[15] = w

  xc[16]     = 1 - x - y
  yc[16]     = y
  weight[16] = w

  pts = hcat( xc, yc )

  weight = weight/2

  return pts, weight

end


text:
# Test variational information

using Test
using Clustering

@testset "varinfo() (variational information)" begin

a1 = [1, 1, 1, 2, 2, 2, 3, 3, 3, 3]
a2 = [1, 1, 1, 1, 2, 2, 2, 2, 2, 2]

@test varinfo(a1, a1) ≈ 0.0 atol=1.0e-12
@test varinfo(a2, a2) ≈ 0.0 atol=1.0e-12

v = varinfo(a1, a2)
v_ = varinfo(a2, a1)
@test 0.0 < v < log(3)
@test v ≈ v_

a1 = [1, 2, 3, 4, 5]
a2 = [1, 1, 1, 1, 1]

@test varinfo(a1, a2) ≈ log(5)
@test varinfo(a2, a1) ≈ log(5)

a1 = [1, 1, 1, 2, 2, 2]
a2 = [2, 2, 2, 1, 1, 1]
@test varinfo(a1, a2) ≈ 0.0 atol=1.0e-12

end


text:
### Rollout policies for use with MCTS

import POMCP

## Open Loop

# # maybe I'll finish this later, but it doesn't seem that necessary for now
# type NonAdaptive <: Solver end
#
# function solve(s::NonAdaptive, p::OPCSP, pol::StaticPath=StaticPath(Array(Int,0)))
#
# end
#
# type StaticPath <: Policy
#     d::Dict{Int, Int} # maps current node to next
# end
# function StaticPath(Array
#
# function action(p::StaticPath, b::OPCSPDistribution, act::OPCSPAction=OPCSPAction(0))
#
# end

## Open loop value estimate
# set the value_estimate_method of the POMCPSolver to use this
function POMCP.estimate_value(pomcp::POMCP.POMCPPolicy, op::OPCSP, s::OPCSPState, h::POMCP.BeliefNode)
    @assert s.remaining == h.B.remaining
    openset = cat(1, s.i, collect(s.open))
    if pomcp.solver.value_estimate_method == :OPCSPMostLikely
        mean_future = SimpleOP(op.r[openset] + d_belief.mean[openset],
                          [],
                          s.remaining,
                          1,
                          findfirst(openset, op.stop),
                          op.distances[openset, openset])
        soln = solve_op(mean_future)
        path = build_path(mean_future, soln)
        return reward(mean_future, path)
    elseif pomcp.solver.value_estimate_method == :OPCSPCheat
        cheat_future = SimpleOP(op.r[openset] + s.mu[openset],
                                [],
                                s.remaining,
                                1,
                                findfirst(openset, op.stop),
                                op.distances[openset, openset])
        soln = solve_op(cheat_future)
        path = build_path(cheat_future, soln)
        return reward(cheat_future, path)
    else
        return POMCP.rollout(pomcp, s, h)
    end
end


text:
using Diban
using DataFrames

tbl = DataFrame(
    int32 = Int32[-1, 0, 1],
    int64 = Int64[-10, 0, 10],
    float32 = Float32[-0.5, 0, 0.5],
    float64 = Float64[-0.5, 0, 0.5],
    bool = [true, false, true],
    string = ["abc", "def", "ghi"],
    int32m = Union{Missing, Int32}[-1, missing, 1],
    int64m = Union{Missing, Int64}[-10, missing, 10],
    float32m = Union{Missing, Float32}[-0.5, missing, 0.5],
    float64m = Union{Missing, Float64}[-0.5, missing, 0.5],
    boolm = Union{Missing, Bool}[true, missing, false],
    stringm = Union{Missing, String}["abc", missing, "ghi"],
)

path = "c:/scratch/tmp.parquet"
write_parquet(path, tbl)
a = read_parquet(path)

Diban.metadata(path)

# using Diban: COMPRESSION_CODEC_CODE, write_col_dict, write_col_chunk,
#     create_schema_parent_node, create_col_schema, COL_TYPE_CODE, write_thrift, PAR2

using Thrift
using Tables
using Snappy, CodecZstd
using LittleEndianBase128

using LittleEndianBase128: encode

read_parquet("c:/scratch/abc.parquet")

Diban.metadata("c:/scratch/abc.parquet")

a.refs

a = categorical(["a",missing,"b"])

using DataAPI

DataAPI.refarray(a)
DataAPI.levels(a)
DataAPI.refpool(a)

DataAPI.defaultarray(eltype(a), 1)

io = open("c:/scratch/tmp.parquet")

sz = filesize("c:/scratch/tmp.parquet")

seek(io, sz - 4)

read(io, 4) |> String

seek(io, sz - 8)
metasize = read(io, UInt32)  |> Int

data_size = sz - 12 - metasize

seek(io, 4+data_size)

read_thrift(io, PAR2.FileMetaData)

close(io)

Diban.metadata(path)
ParFile(path)

write_parquet(path, tbl; compression_codec = "lz4")
write_parquet(path, tbl; compression_codec = "zstd")
write_parquet(path, tbl; compression_codec = "gzip")
write_parquet(path, tbl; compression_codec = "uncompressed")

Diban.metadata("c:/scratch/abc.parquet")

io=open("c:/scratch\\abc.parquet")
seek(io, 4)

ph = read_thrift(io, PAR2.PageHeader)

data = read(io, ph.compressed_page_size)
udata = Snappy.uncompress(data)

string(data[end], base=2)

close(io)

a = BitArray(rand(Bool, 65))

io=IOBuffer()
write(io, a)

seek(io, 0)

ok = read(io)
close(io)

string(ok[1], base=2, pad=8)
string(ok[2], base=2, pad=8)
string(ok[3], base=2, pad=8)

colname = "abc"
colvals = tbl[!, colname]

#close(fileio)

Diban.metadata(path)

Diban.metadata("c:/scratch/abc.parquet")

read_parquet(path)

### example
path = "c:/scratch/abc.parquet"
Diban.metadata(path)

ff = open(path)

String(read(ff, 4))
fsz = filesize(path)
seek(ff, fsz-8)

metadata_size = Int(read(ff, UInt32))

data_size = fsz - 12 - metadata_size
seek(ff, 4 + data_size)

example_schema = read_thrift(ff, PAR2.FileMetaData)
close(ff)

close(fileio)

ff = open("tmp.parquet")
String(read(ff, 4))
dh = read_thrift(ff, PAR2.PageHeader)
tmpio = read(ff, dh.compressed_page_size) |> Snappy.uncompress |> IOBuffer
read(tmpio, Int)
read(tmpio, Int)
read(tmpio, Int)
dh = read_thrift(ff, PAR2.PageHeader)

close(ff)

# write encoded data


text:
## Specialized type of nonlinear mixed-effects model with a single
## grouping factor and random effects for each parameter in the
## nonlinear model

type SimpleNLMM{T<:FP} <: NLMM{T}
    m::NLregModF{T}
    inds::Vector               # grouping factor indices
    nrep::Vector               # run-length encoding of inds
    lambda::AbstractMatrix{T}  # template block for Lambda
    L::Vector{Matrix{T}}       # should this be a 3-D array?
    beta::Vector{T}            # fixed-effects parameter vector
    u::Matrix{T}               # spherical random-effects values
    delu::Matrix{T}            # increment in the PNLS algorithm
    b::Matrix{T}               # random effects on original scale
    phi::Matrix{T}             # phi := beta + lambda * (u + fac*delu)
    nAGQ::Int                  # adaptive Gauss-Hermite quadrature points
    mxpnls::Int                # maximum number of PNLS iterations
    minfac::T                  # minimum step factor in PNLS
    tolsqr::T                  # squared tolerance for orthogonality conv.crit.
end
function SimpleNLMM{T<:FP}(m::NLregModF{T},inds::Vector,
                           lambda::AbstractMatrix{T},beta::Vector{T})
    n = size(m)[end]; p = npars(m); np = n*p; ui, nrep = rle(inds); ni = length(ui)
    length(inds) == n || error("length(inds) = $(length(inds)), should be $n")
    isperm(ui) ||
        error("similar inds must be adjacent and unique(inds) a permutation")
    size(lambda) == (p,p) ||
        error("size(lambda) = $(size(lambda)) should be $((p,p))")
    L = Matrix{T}[eye(T,p) for i in 1:ni]
    u = zeros(T,(p,ni))
    SimpleNLMM(m,inds,nrep,lambda,L,beta,u,copy(u),similar(u),similar(u),
               1,300,convert(T,0.5^9),convert(T,1e-8))
end
function SimpleNLMM(nl::Union(PLinearLS,NonlinearLS),inds::Vector,ltype::DataType)
    lambda = ltype(eye(length(coef(nl))))
    if ltype == Triangular
        lambda.uplo = 'L'
        lambda.unitdiag = 'N'
    end
    SimpleNLMM(deepcopy(nl.m),inds,lambda,coef(nl))
end

StatsBase.coef(nm::SimpleNLMM) = copy(nm.beta)

nlregmod(nm::SimpleNLMM) = nm.m

## penalized residual sum of squares
function prss!{T<:FP}(nm::SimpleNLMM{T},fac::T)
    b = nm.b                  # random effects on original scale
    copy!(b,nm.u)             # initialize with u
    fma!(b,nm.delu,fac)       # add fac*delu
    ssu = sumsq(b)            # record squared length of u + fac(delu)
    A_mul_B!(nm.lambda,b)     # convert to b scale
    updtmu!(nm.m,broadcast!(+,nm.phi,b,nm.beta),nm.inds) + ssu # rss + penalty
end
prss!{T<:FP}(nm::SimpleNLMM{T}) = prss!(nm,zero(T))

pwrss(nm::SimpleNLMM) = sumsq(nm.m.resid) + sumsq(vec(nm.u))

## logdet of L, the Cholesky factor
function ldL2{T<:FP}(m::Matrix{T})
    dd = zero(T)
    for i in 1:size(m,1) dd += log(m[i,i]) end
    dd + dd
end
ldL2{T<:FP}(nm::SimpleNLMM{T}) = sum(ldL2,nm.L)

## Update the L array using the current values of tgrad
function updtL!{T<:FP}(nm::SimpleNLMM{T})
    m = nm.m; u = nm.u; delu = copy!(nm.delu,nm.u); L = nm.L; ee = eye(T,size(L[1],1))
    mm = nm.lambda' * m.tgrad
    rr = residuals(m); nn = nm.nrep; offset = 0
    for i in 1:length(L)
        ni = nn[i]; ii = offset+(1:nn[i]); offset += ni;
        g = sub(mm,:,ii); du = sub(delu,:,i)
        BLAS.gemv!('N',1.,g,sub(rr,ii),-1.,du)
        ## Don't need to check for singularity in potrf! result b/c adding I
        ## TODO: use two applications of trsv and evaluate the numerator of the conv. crit.
        LAPACK.potrs!('L',LAPACK.potrf!('L',BLAS.syrk!('L','N',1.,g,1.,copy!(L[i],ee)))[1],du)
    end
    prss!(nm,one(T))                    # evaluate prss with full increment
end

## increment the spherical random effects as nm.u += nm.delu * fac and zero out nm.delu
function incr!{T<:FP}(nm::SimpleNLMM{T}, fac::T)
    fma!(nm.u,nm.delu,fac)
    fill!(nm.delu,zero(T))
end

## evaluate the approximate deviance using adaptive Gauss-Hermite quadrature
function StatsBase.deviance{T<:FP}(nm::SimpleNLMM{T})
    nm.nAGQ == 1 || error("General adaptive Gauss-Hermite quadrature not yet written")
    n = nobs(nm)
    ldL2(nm) + n * (one(T) + log(2pi*prss!(nm)/n))
end

## Determine the conditional mode of the random effects using penalized nonlinear least squares
function pnls!{T<:FP}(nm::SimpleNLMM{T},u0::Matrix{T};verbose=false)
    copy!(nm.u,u0); fill!(nm.delu,zero(T)) # start from a fixed position
    oldprss = prss!(nm,0.); conv = one(T);
    for i in 1:nm.mxpnls
        newprss = updtL!(nm); f = one(T)
        verbose && @printf(" %12f, %12f\n", oldprss, newprss)
        while oldprss < newprss && f >= nm.minfac
            f /= convert(T,2)
            newprss = prss!(nm,f)
            verbose && @printf("  f = %8e: %12f\n", f, newprss)
        end
        f < nm.minfac && error("Failure to reduce prss at u = ",nm.u)
        incr!(nm,f)
        (conv = (oldprss - newprss)/newprss) < nm.tolsqr && break
        oldprss = newprss
    end
    conv >= nm.tolsqr && error("Failure to converge in ", nm.mxpnls, " iterations")
    deviance(nm)
end
pnls!{T<:FP}(nm::SimpleNLMM{T};verbose=true) = pnls!(nm, zeros(T,size(nm.u));verbose=verbose)

function Base.show(io::IO, m::SimpleNLMM)
    mm = nlregmod(m)
    mstr = m.nAGQ == 1 ? "Laplace approximation" : "adaptive Gauss-Hermite quadrature ($(m.nAGQ))"
    println(io, "Simple, nonlinear mixed-effects model fit by ", mstr)
    oo = deviance(m)
    @printf(io, " logLik: %f, deviance: %f", -oo/2., oo)
    println(io); println(io)

    @printf(io, " Variance components:\n                Variance    Std.Dev.\n")
    stdm = std(m); fnms = vcat(pnames(mm),"Residual")
    for i in 1:length(fnms)
        si = stdm[i]
        print(io, " ", rpad(fnms[i],12))
        @printf(io, " %10f  %10f\n", abs2(si[1]), si[1])
        for j in 2:length(si)
            @printf(io, "             %10f  %10f\n", abs2(si[j]), si[j])
        end
    end
    @printf(io," Number of obs: %d; levels of grouping factor: %d", nobs(m), size(m.u,2))
    println(io); println(io)

    println(io, "Fixed effects parameters:")
    show(io, coeftable(m))
    println(io)
end

Base.std(nm::SimpleNLMM) = scale(nm) * [diag(nm.lambda),1.]


text:
struct Constants{T<:AbstractFloat,Tprog<:AbstractFloat}

    # RUNGE-KUTTA COEFFICIENTS 3rd/4th order including timestep Δt
    RKaΔt::Array{Tprog,1}
    RKbΔt::Array{Tprog,1}

    # BOUNDARY CONDITIONS
    one_minus_α::Tprog      # tangential boundary condition for the ghost-point copy

    # PHYSICAL CONSTANTS
    g::T                    # gravity
    cD::T                   # quadratic bottom friction - incl grid spacing
    rD::T                   # linear bottom friction - incl grid spacing
    γ::T                    # frequency of interface relaxation
    cSmag::T                # Smagorinsky constant
    νB::T                   # biharmonic diffusion coefficient
    rSST::T                 # tracer restoring timescale
    jSST::T                 # tracer consumption timescale
    SSTmin::T               # tracer minimum
    ωFη::Float64            # frequency [1/s] of seasonal surface forcing incl 2π
    ωFx::Float64            # frequency [1/s] of seasonal wind x incl 2π
    ωFy::Float64            # frequency [1/2] of seasonal wind y incl 2π
end

"""Generator function for the mutable struct Constants."""
function Constants{T,Tprog}(P::Parameter,G::Grid) where {T<:AbstractFloat,Tprog<:AbstractFloat}

    # Runge-Kutta 2nd/3rd/4th order coefficients including time step Δt and grid spacing Δ
    # a are the coefficents to sum the rhs on the fly, such that sum=1
    # b are the coefficents for the update that used for a new evaluation of the RHS
    if P.RKo == 2     # Heun's method
        RKaΔt = Tprog.([1/2,1/2]*G.dtint/G.Δ)
        RKbΔt = Tprog.([1]*G.dtint/G.Δ)
    elseif P.RKo == 3     # version 2 / Heun's 3rd order
        RKaΔt = Tprog.([1/4,0.,3/4]*G.dtint/G.Δ)
        RKbΔt = Tprog.([1/3,2/3]*G.dtint/G.Δ)
    elseif P.RKo == 4
        RKaΔt = Tprog.([1/6,1/3,1/3,1/6]*G.dtint/G.Δ)
        RKbΔt = Tprog.([.5,.5,1.]*G.dtint/G.Δ)
    end

    one_minus_α = Tprog(1-P.α)    # for the ghost point copy/tangential boundary conditions
    g = T(P.g)                # gravity - for Bernoulli potential

    # BOTTOM FRICTION COEFFICENTS
    # incl grid spacing Δ for non-dimensional gradients
    cD = T(-G.Δ*P.cD)             # quadratic drag [m]
    rD = T(-G.Δ/(P.τD*24*3600))   # linear drag [m/s]

    # INTERFACE RELAXATION FREQUENCY
    # incl grid spacing Δ for non-dimensional gradients
    γ = T(G.Δ/(P.t_relax*3600*24))    # [m/s]

    # BIHARMONIC DIFFUSION
    cSmag = T(-P.cSmag)   # Smagorinsky coefficient
    νB = T(-P.νB/30000)   # linear scaling based on 540m^s/s at Δ=30km

    # TRACER ADVECTION
    rSST = T(G.dtadvint/(P.τSST*3600*24))    # tracer restoring [1]
    jSST = T(G.dtadvint/(P.jSST*3600*24))    # tracer consumption [1]
    SSTmin = T(P.SSTmin)

    # TIME DEPENDENT FORCING
    ωFη = -2π*P.ωFη/24/365.25/3600
    ωFx = 2π*P.ωFx/24/365.25/3600
    ωFy = 2π*P.ωFy/24/365.25/3600

    return Constants{T,Tprog}(  RKaΔt,RKbΔt,one_minus_α,
                                g,cD,rD,γ,cSmag,νB,rSST,
                                jSST,SSTmin,ωFη,ωFx,ωFy)
end


text:
"""
Provides Genie with encryption and decryption capabilities.
"""
module Encryption

using Genie, App, Nettle

const ENCRYPTION_METHOD = "AES256"

"""
    encrypt{T}(s::T) :: String

Encrypts `s`.
"""
function encrypt{T}(s::T) :: String
  (key32, iv16) = encryption_sauce()
  encryptor = Encryptor(ENCRYPTION_METHOD, key32)

  Nettle.encrypt(encryptor, :CBC, iv16, add_padding_PKCS5(Vector{UInt8}(s), 16)) |> bytes2hex
end

"""
    decrypt(s::String) :: String

Decrypts `s` (a `string` previously encrypted by Genie).
"""
function decrypt(s::String) :: String
  (key32, iv16) = encryption_sauce()
  decryptor = Decryptor(ENCRYPTION_METHOD, key32)
  deciphertext = Nettle.decrypt(decryptor, :CBC, iv16, s |> hex2bytes)

  String(trim_padding_PKCS5(deciphertext))
end

function encryption_sauce() :: Tuple{Vector{UInt8},Vector{UInt8}}
  passwd = App.secret_token()[1:32]
  salt = hex2bytes(App.secret_token()[33:64])

  gen_key32_iv16(Vector{UInt8}(passwd), salt)
end

end


text:
include("vanilla_gradient.jl")

# Guided Backprop

function change_activation(model::Metalhead.InceptionBlock, activation)
    path_1 = Conv(activation, model.path_1.weight, model.path_1.bias, model_path_1.stride, model.path_1.pad, model.path_1.dilation)

    path_2 = (Conv(activation, model.path_2[1].weight, model.path_2[1].bias, model_path_2[1].stride, model.path_2[1].pad, model.path_2[1].dilation),
            Conv(activation, model.path_2[2].weight, model.path_2[2].bias, model_path_2[2].stride, model.path_2[2].pad, model.path_2[2].dilation))

    path_3 = (Conv(activation, model.path_3[1].weight, model.path_3[1].bias, model_path_3[1].stride, model.path_3[1].pad, model.path_3[1].dilation),
            Conv(activation, model.path_3[2].weight, model.path_3[2].bias, model_path_3[2].stride, model.path_3[2].pad, model.path_3[2].dilation))

    path_4 = (model.path_4[1],
            Conv(activation, model.path_4[2].weight, model.path_4[2].bias, model_path_4[2].stride, model.path_4[2].pad, model.path_4[2].dilation))
    Metalhead.InceptionBlock(path_1, path_2, path_3, path_4)
end

# Resnet model is not supported
function change_activation(model::Chain, activation)
  updated_model = []
  for (i, l) in enumerate(model.layers)
    T = typeof(l)
    if T <: Dense
      push!(updated_model, Dense(l.W, l.b, identity))
      push!(updated_model, activation)
    elseif T <: Conv
      push!(updated_model, Conv(identity, l.weight, l.bias, l.stride, l.pad, l.dilation))
      push!(updated_model, activation)
    elseif T <: BatchNorm
      push!(updated_model, BatchNorm(identity, l.β, l.γ, l.μ, l.σ², l.ϵ, l.momentum, l.active))
      push!(updated_model, activation)
    elseif T <: Chain
      push!(updated_model, change_activation(l, activation))
    elseif T <: ConvTranspose
      push!(updated_model, ConvTranspose(identity, l.weight, l.bias, l.stride, l.pad, l.dilation))
      push!(updated_model, activation)
    elseif T <: DepthwiseConv
      push!(updated_model, DepthwiseConv(identity, l.weight, l.bias, l.stride, l.pad, l.dilation))
      push!(updated_model, activation)
    elseif T <: CrossCor
      push!(updated_model, CrossCor(identity, l.weight, l.bias, l.stride, l.pad, l.dilation))
      push!(updated_model, activation)
    elseif T <: Metalhead.Bottleneck # To handle the DenseNet Model
      push!(updated_model, Bottleneck(change_activation(l.layer, activation)))
    elseif T <: Metalhead.InceptionBlock # To handle Inception Layer Model
      push!(updated_model, change_activation(l, activation))
    elseif T <: Metalhead.ResidualBlock
      push!(updated_model, change_activation(l, activation))
    else
      push!(updated_model, l)
    end
  end
  Chain(updated_model...)
end

function guided(x)
    max.(zero(x),x)
end

@adjoint guided(x) = guided(x), c̄ -> (Int.(x .> zero(x)) .* max.(zero(c̄), c̄),)

"""
    viz_guidedbackprop(img, model; top_k=1, target_based=false, target_index=-1)

In this method of Gradient visualization, the backpropagation of RELU activation function are overridden (so called Guided RELU) in such a way that only the non-negative gradients are backpropagated. Here, the guided RELU is applied onto the input gradients which stands as the main difference from the `viz_deconvolution`. Once the activation function is overridden, the gradient is evaluated here using the `viz_backprop` function.

Note: The function currently doesn't support ResNet type models.
See also: [`viz_backprop`](@ref), [`viz_deconvolution`](@ref)
"""
function viz_guidedbackprop(img, model; top_k=1, target_based=false, target_index=-1)
  model1 = change_activation(model, guided)
  print(model1)
  viz_backprop(img, model1; top_k=top_k, target_based=target_based, target_index=target_index)
end

# DeconvNet

function deconv(x)
    max.(zero(x),x)
end

@adjoint deconv(x) = deconv(x), c̄ -> (max.(zero(c̄), c̄),)

"""
    viz_deconvolution(img, model; top_k=1, target_based=false, target_index=-1)

In this method of Gradient visualization, the backpropagation of RELU activation function are overridden (so called Guided RELU) in such a way that only the non-negative gradients are backpropagated. Here, the guided RELU is applied onto the `output` gradients which stands as the main difference from the `viz_guidedbackprop`. Once the activation function is overridden, the gradient is evaluated here using the `viz_backprop` function.

Note: The function currently doesn't support ResNet type models.
See also: [`viz_backprop`](@ref), [`viz_guidedbackprop`](@ref)
"""
function viz_deconvolution(img, model; top_k=1, target_based=false, target_index=-1)
  model = change_activation(model, deconv)
  viz_backprop(img, model; top_k=top_k, target_based=target_based, target_index=target_index)
end


text:
"""
Use line resistance to estimate line length.

In:

* `S_base` [VA], System base VA
* `V_base` [V], Base voltage
* `R_pu` [pu], Line resistance
* `R_cond` [ohms/m], Conductor resistance
* `bundle` [-], Number of conductors per phase

Out: `length` [m], Line length

2015-08-11.
Ported from .m file provided by Jonathon Martin.
Original code by Mads Almassalkhi.
"""
function estimate_length(
    S_base::Float64,
    V_base::Float64,
    R_pu::Float64,
    R_cond::Float64,
    bundle::Int64
    )
    R_base = V_base^2 / S_base
    R_ohms = R_pu * R_base
    line_length = R_ohms / (R_cond / bundle)
    # Ensure length estimate is positive
    default = 1e-4 * 1609.3
    line_length <= 0 && (line_length = default)
    return line_length
end

"""
Estimate line lengths for all branches, given a `network_data` dict. The input
dict is the output of `PowerModels.parse_file`.
"""
function estimate_length(network_data::Dict{String, Any})
    !network_data["per_unit"] && @error "Network data is not in per unit!"

    baseMVA = network_data["baseMVA"]
    S_base = baseMVA * 1e6
    basekv = branch_basekv(network_data)

    branch_data = Vector{Tuple}()
    for (k, v) in network_data["branch"]
        # note: if network_data["per_unit"] == true, then rate_a has been
        # divided by baseMVA.
        push!(branch_data, (k, v["br_r"], v["rate_a"] * baseMVA))
    end

    line_lengths = Dict{String, Float64}()
    for i in 1:length(branch_data)
        id, R_pu, rate_a = branch_data[i]

        V_base = basekv[string(id)] * 1e3 # convert kV to V

        # convert long-term VA rating to current rating
        I_lim = rate_a * 1e6 / V_base

        a = acsr_interpolation(I_lim, V_base)
        line_lengths[string(id)] = estimate_length(S_base, V_base, R_pu, a.R, a.bundle)
    end
    return line_lengths
end


text:
using Oceananigans.Advection
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBoundary
using Oceananigans.Coriolis
using Oceananigans.Operators
using Oceananigans.TurbulenceClosures: ∇_dot_qᶜ, ∂ⱼ_τ₁ⱼ, ∂ⱼ_τ₂ⱼ

@inline half_g_h²(i, j, k, grid, h, g) = @inbounds 1/2 * g * h[i, j, k]^2

@inline x_pressure_gradient(i, j, k, grid, g, h, formulation) = ∂xᶠᶜᶜ(i, j, k, grid, half_g_h², h, g)
@inline y_pressure_gradient(i, j, k, grid, g, h, formulation) = ∂yᶜᶠᶜ(i, j, k, grid, half_g_h², h, g)

@inline x_pressure_gradient(i, j, k, grid, g, h, ::VectorInvariantFormulation) = g * ∂xᶠᶜᶜ(i, j, k, grid, h)
@inline y_pressure_gradient(i, j, k, grid, g, h, ::VectorInvariantFormulation) = g * ∂yᶜᶠᶜ(i, j, k, grid, h)

@inline bathymetry_contribution_x(i, j, k, grid, g, h, hB, formulation) = g * h[i, j, k] * ∂xᶠᶜᶜ(i, j, k, grid, hB)
@inline bathymetry_contribution_y(i, j, k, grid, g, h, hB, formulation) = g * h[i, j, k] * ∂yᶜᶠᶜ(i, j, k, grid, hB)

@inline bathymetry_contribution_x(i, j, k, grid, g, h, hB, ::VectorInvariantFormulation) = g * ∂xᶠᶜᶜ(i, j, k, grid, hB)
@inline bathymetry_contribution_y(i, j, k, grid, g, h, hB, ::VectorInvariantFormulation) = g * ∂yᶜᶠᶜ(i, j, k, grid, hB)

"""
Compute the tendency for the x-directional transport, uh
"""
@inline function uh_solution_tendency(i, j, k, grid,
                                      gravitational_acceleration,
                                      advection,
                                      coriolis,
                                      closure,
                                      bathymetry,
                                      solution,
                                      tracers,
                                      diffusivities,
                                      forcings,
                                      clock,
                                      formulation)

    g = gravitational_acceleration

    return ( - div_mom_u(i, j, k, grid, advection, solution, formulation)
             - x_pressure_gradient(i, j, k, grid, g, solution.h, formulation)
             - x_f_cross_U(i, j, k, grid, coriolis, solution)
             + bathymetry_contribution_x(i, j, k, grid, g, solution.h, bathymetry, formulation)
             - ∂ⱼ_τ₁ⱼ(i, j, k, grid, closure, diffusivities, solution, tracers, clock, nothing)
             + forcings[1](i, j, k, grid, clock, merge(solution, tracers)))
end

"""
Compute the tendency for the y-directional transport, vh.
"""
@inline function vh_solution_tendency(i, j, k, grid,
                                      gravitational_acceleration,
                                      advection,
                                      coriolis,
                                      closure,
                                      bathymetry,
                                      solution,
                                      tracers,
                                      diffusivities,
                                      forcings,
                                      clock,
                                      formulation)

     g = gravitational_acceleration

    return ( - div_mom_v(i, j, k, grid, advection, solution, formulation)
             - y_pressure_gradient(i, j, k, grid, g, solution.h, formulation)
             - y_f_cross_U(i, j, k, grid, coriolis, solution)
             + bathymetry_contribution_y(i, j, k, grid, g, solution.h, bathymetry, formulation)
             - ∂ⱼ_τ₂ⱼ(i, j, k, grid, closure, diffusivities, solution, tracers, clock)
             + forcings[2](i, j, k, grid, clock, merge(solution, tracers)))
end

"""
Compute the tendency for the height, h.
"""
@inline function h_solution_tendency(i, j, k, grid,
                                     gravitational_acceleration,
                                     advection,
                                     coriolis,
                                     closure,
                                     solution,
                                     tracers,
                                     diffusivities,
                                     forcings,
                                     clock,
                                     formulation)

    return ( - div_Uh(i, j, k, grid, advection, solution, formulation)
             + forcings.h(i, j, k, grid, clock, merge(solution, tracers)))
end

@inline function tracer_tendency(i, j, k, grid,
                                 val_tracer_index::Val{tracer_index},
                                 advection,
                                 closure,
                                 solution,
                                 tracers,
                                 diffusivities,
                                 forcing,
                                 clock,
                                 formulation) where tracer_index

    @inbounds c = tracers[tracer_index]

    return ( - div_Uc(i, j, k, grid, advection, solution, c, formulation)
             + c_div_U(i, j, k, grid, solution, c, formulation)
             - ∇_dot_qᶜ(i, j, k, grid, closure, diffusivities, val_tracer_index, solution, tracers, clock, nothing)
             + forcing(i, j, k, grid, clock, merge(solution, tracers))
            )
end


text:
using DynamicPPL
using Distributions
using ForwardDiff
using Tracker
using Zygote

using Random
using Test

dir = splitdir(splitdir(pathof(DynamicPPL))[1])[1]
include(dir*"/test/Turing/Turing.jl")
using .Turing

turnprogress(false)

include("test_util.jl")

@testset "DynamicPPL.jl" begin
    include("utils.jl")
    include("compiler.jl")
    include("varinfo.jl")
    include("sampler.jl")
    include("prob_macro.jl")
    include("independence.jl")
    include("distribution_wrappers.jl")

    include("threadsafe.jl")

    @testset "compat" begin
        include(joinpath("compat", "ad.jl"))
    end
end


text:
@symbol_func function f_BS(cur_reactor::AbstractReactor)
  cur_f = K_BS(cur_reactor)

  cur_f *= cur_reactor.T_bar

  cur_f /= cur_reactor.I_P

  cur_f
end


text:
"""
    assigndist_betabinomial(n, α, β; trunc_lower = -Inf, trunc_upper = Inf)

Assign parameters to a beta binomial distribution with `n` trials and
shape parameters `α` and `β`, optionally truncating the distribution.
"""
function assigndist_betabinomial(n, α, β; trunc_lower = -Inf, trunc_upper = Inf)
    Truncated(BetaBinomial(n, α, β), trunc_lower, trunc_upper)
end

export assigndist_betabinomial


text:
module MangalTestUsers

    using Mangal
    using Test

    @test typeof(users()) <: Vector{MangalUser}

end


text:
# ---
# title: 890. Find and Replace Pattern
# id: problem890
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/find-and-replace-pattern/description/>
# hidden: true
# ---
#
# You have a list of `words` and a `pattern`, and you want to know which words
# in `words` matches the pattern.
#
# A word matches the pattern if there exists a permutation of letters `p` so
# that after replacing every letter `x` in the pattern with `p(x)`, we get the
# desired word.
#
# ( _Recall that a permutation of letters is a bijection from letters to
# letters: every letter maps to another letter, and no two letters map to the
# same letter._ )
#
# Return a list of the words in `words` that match the given pattern.
#
# You may return the answer in any order.
#
#
#
# **Example 1:**
#
#
#
#     Input: words = ["abc","deq","mee","aqq","dkd","ccc"], pattern = "abb"
#     Output: ["mee","aqq"]
#     Explanation: "mee" matches the pattern because there is a permutation {a -> m, b -> e, ...}.
#     "ccc" does not match the pattern because {a -> c, b -> c, ...} is not a permutation,
#     since a and b map to the same letter.
#
#
#
# **Note:**
#
#   * `1 <= words.length <= 50`
#   * `1 <= pattern.length = words[i].length <= 20`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
export _run_vstab, _run_multi_vstab

# mutable struct VDPowerModel <: _PM.AbstractACModel _PM.@pm_fields end

"""
run model for Voltge-Deviation objective with AC Power Flow equations
"""

function _run_vstab(file, model_type::_PM.Type, optimizer; kwargs...)
    return _PM.solve_model(file, model_type, optimizer, _build_vstab; kwargs...)
end

"""
given a JuMP model and a PowerModels network data structure,
builds an Voltge-Deviation formulation of the given data and returns the JuMP model
"""

function _build_vstab(pm::_PM.AbstractPowerModel)

    _PM.variable_bus_voltage(pm)
    _PM.variable_gen_power(pm)
    _PM.variable_branch_power(pm)
    _PM.variable_dcline_power(pm, bounded = false) # TODO: why false?

    objective_vstab(pm)

    _PM.constraint_model_voltage(pm)

    for i in _PM.ids(pm, :ref_buses)
        _PM.constraint_theta_ref(pm, i)
    end

    for i in _PM.ids(pm, :bus)
        _PM.constraint_power_balance(pm, i)
    end

    for (i, branch) in _PM.ref(pm, :branch)
        _PM.constraint_ohms_yt_from(pm, i)
        _PM.constraint_ohms_yt_to(pm, i)

        _PM.constraint_thermal_limit_from(pm, i)
        _PM.constraint_thermal_limit_to(pm, i)
    end

    for i in _PM.ids(pm, :dcline)
        _PM.constraint_dcline_power_losses(pm, i)
    end
end

function objective_vstab(pm::_PM.AbstractPowerModel)
    return JuMP.@objective(pm.model, Min,
        sum((var(pm, :vm, content["element_index"]) - content["value"])^2 for (i, content) in pm.ext[:setpoint_v]))
end

"""
run model for Voltge-Deviation objective with AC Power Flow equations
"""

function _run_multi_vstab(file, model_type::_PM.Type, optimizer; kwargs...)
    return _PM.solve_model(file, model_type, optimizer, _build_multi_vstab; multinetwork=true, kwargs...)
end

"""
given a JuMP model and a PowerModels network data structure,
builds an Voltge-Deviation formulation of the given data and returns the JuMP model
"""

function _build_multi_vstab(pm::_PM.AbstractPowerModel)
    for (n, network) in _PM.nws(pm)
            _PM.variable_bus_voltage(pm, nw=n)
            _PM.variable_gen_power(pm, nw=n)
            _PM.variable_branch_power(pm, nw=n)
            _PM.variable_dcline_power(pm, nw=n)

            _PM.constraint_model_voltage(pm, nw=n)

            for i in ids(pm, :ref_buses, nw=n)
                _PM.constraint_theta_ref(pm, i, nw=n)
            end

            for i in ids(pm, :bus, nw=n)
                _PM.constraint_power_balance(pm, i, nw=n)
            end

            for i in ids(pm, :branch, nw=n)
                _PM.constraint_ohms_yt_from(pm, i, nw=n)
                _PM.constraint_ohms_yt_to(pm, i, nw=n)

                _PM.constraint_voltage_angle_difference(pm, i, nw=n)

                _PM.constraint_thermal_limit_from(pm, i, nw=n)
                _PM.constraint_thermal_limit_to(pm, i, nw=n)
            end

            for i in ids(pm, :dcline, nw=n)
                _PM.constraint_dcline_power_losses(pm, i, nw=n)
            end
        end
        objective_multi_vstab(pm)
end

function objective_multi_vstab(pm::_PM.AbstractPowerModel)
    timestep_ids = [id for id in _PM.nw_ids(pm) if id != 0]
    return JuMP.@objective(pm.model, Min,
        sum(
        sum((var(pm, nw, :vm, content["element_index"]) - content["value"])^2 for (i, content) in pm.ext[:setpoint_v])
        for nw in timestep_ids)
            )
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus <: SwaggerModel

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus(;)
        o = new()
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus = Dict{Symbol,Symbol}()
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus = Dict{Symbol,String}()
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus[property_name]

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus)
    true
end

function validate_property(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus }, name::Symbol, val)
end


text:
# Database URLs
if isempty(ARGS)
    @warn "!!! No database urls provided as command line arguments, proceeding with default urls !!!"
    db_url_building = "sqlite:///C:\\DATA\\Spine\\toolbox\\projects\\spinea4testing\\building_data\\Building_Data.sqlite"
    db_url_combined ="sqlite:///C:\\DATA\\Spine\\toolbox\\projects\\spinea4testing\\combined_data_1\\Combined_Data.sqlite"
else
    if length(ARGS) != 2
        error("!!! Provide only database urls both for the baseline consumption database as well as the target database !!!")
    else
        db_url_building = ARGS[1]
        db_url_combined = ARGS[2]
    end
end

# Load necessary modules
using SpineInterface
using Dates

# Open building database
println("... Loading building model results ...")
using_spinedb(db_url_building)
println("... Done ...")

# Determine relevant indeces
r = first(SpineInterface.report())
dummy_names = [Symbol(75, "FI_Dummy"), Symbol(74, "FI_Dummy")]
units = filter(x->x.name in dummy_names, SpineInterface.unit())
und = filter(x->x.unit in units, SpineInterface.unit__node__direction())
nodes = unique(map(x->x.node, und))
uncd = []
for (u,n,d) in und
    c = first(filter(x->x.node==n, SpineInterface.node__commodity())).commodity
    append!(uncd, [(unit=u, node=n, commodity=c, direction=d)])
end

# Extract the baseline electricity consumption
baseline_electricity_demand = Dict()
for (u,n,c,d) in uncd
    baseline_electricity_demand[(node=n,)] = SpineInterface.flow(report=r, unit=u, node=n, commodity=c, direction=d)
end

# Open combined database
println("... Loading combined database ...")
using_spinedb(db_url_combined)
println("... Done ...")

# Repeat baseline electricity consumption until it's length matches power system data
println("... Repeating baseline electricity consumption to match power system data ...")
for n in nodes
    ind = baseline_electricity_demand[(node=n,)].indexes
    val = baseline_electricity_demand[(node=n,)].values
    ign = baseline_electricity_demand[(node=n,)].ignore_year
    rep = baseline_electricity_demand[(node=n,)].repeat
    while length(ind) < length(SpineInterface.demand(node=n).indexes)
        append!(ind, ind + Hour(length(ind)))
        append!(val, val)
    end
    # Cut out any unnecessary tail
    new_ind = filter(x->x in SpineInterface.demand(node=n).indexes, ind)
    new_val_indexes = findall(x->x in new_ind, ind)
    new_val = val[new_val_indexes]
    # Create the new extended TimeSeries
    baseline_electricity_demand[(node=n,)] = TimeSeries(
        new_ind,
        new_val,
        ign,
        rep
    )
end
println("... Done ...")

# Calculate the new electricity demand
println("... Calculating new electricity demand ...")
new_electricity_demand = Dict{NamedTuple,TimeSeries}()
for (u,n,c,d) in uncd
    ind = SpineInterface.demand(node=n).indexes
    if ind != baseline_electricity_demand[(node=n,)].indexes
        println("Indexes don't match!")
        break
    end
    val = SpineInterface.demand(node=n).values - baseline_electricity_demand[(node=n,)].values
    new_electricity_demand[(node=n,)] = TimeSeries(
        ind,
        val,
        SpineInterface.demand(node=n).ignore_year,
        SpineInterface.demand(node=n).repeat
        )
end
println("... Done ...")

# Overwrite the old demand in the combined database with the new demand
pars = Dict(Symbol("demand") => new_electricity_demand)
println("... Rewriting demand ...")
write_parameters(pars, db_url_combined)
println("... Done ...")


text:
import Base.copy
export PointE, copy, densify, SDCOContext

### Matrix types
const Dense{T} = Array{T, 2} where T<:Number
const DenseSym{T} = Symmetric{T, Dense{T}} where {T<:Real}
const DenseHerm{T} = Hermitian{T, Dense{T}} where {T<:Complex}

const Sparse{T} = SparseArrays.SparseMatrixCSC{T,Int64} where T<:Number
const SparseSym{T} = Symmetric{T, Sparse{T}} where {T<:Real}
const SparseHerm{T} = Hermitian{T, Sparse{T}} where {T<:Complex}

mutable struct PointE{T, U}
    mats::Vector{U}             # Collection of symmetric matrices
    vec::Vector{T}              # R^m array
    dims::Vector{Int}           # Dimension of the sdp cones
end

function PointE(mats::Vector{U}, vec::AbstractArray{T}) where {U<:AbstractMatrix, T<:Number}
    dims = map(x->size(x, 1), mats)

    return PointE(mats, vec, dims)
end

function PointE(dims::AbstractArray{Int}, vecdim::Int, T::DataType, U::DataType)
    mats = Vector{U}(undef, length(dims))

    for (i, dimi) in enumerate(dims)
        if U<:Dense
            mats[i] = zeros(T, dimi, dimi)
        elseif U<: DenseSym
            mats[i] = Symmetric(zeros(T, dimi, dimi))
        elseif U<: DenseHerm
            mats[i] = Hermitian(zeros(T, dimi, dimi))
        elseif U<:Sparse
            mats[i] = spzeros(T, dimi, dimi)
        elseif U<: SparseSym
            mats[i] = Symmetric(spzeros(T, dimi, dimi))
        elseif U<: SparseHerm
            mats[i] = Hermitian(spzeros(T, dimi, dimi))
        else
            @error("PointE(): Unsupported matrix type $U")
        end
    end

    return PointE(mats, zeros(T, vecdim), dims)
end

function copy(x::PointE{T, U}) where {T<:Number, U}
    mats = Vector{U}(undef, length(x.dims))

    for (i, mat) in enumerate(x.mats)
        mats[i] = copy(mat)
    end

    vec = copy(x.vec)
    return PointE(mats, vec)
end

function densify(x::PointE{T, U}) where {T<:Number, U<:Union{AbstractArray, Symmetric}}
    mats = Vector{Dense{T}}([zeros(T, n, n) for n in x.dims])

    for (i, mat) in enumerate(x.mats)
        mats[i] = Matrix(mat)
    end

    return PointE(mats, copy(x.vec))
end

mutable struct PointPrimalDual{T, U}
    x::PointE{T, U}
    y::Vector{T}
    s::PointE{T, U}
end

function PointPrimalDual(dims::AbstractArray{Int}, vecdim::Int, m::Int, T::DataType, U::DataType)
    return PointPrimalDual(PointE(dims, vecdim, T, U), zeros(T, m), PointE(dims, vecdim, T, U))
end

################################################################################
##### Context structure
################################################################################
mutable struct SDCOContext{T, matT}
    c::PointE{T, matT}                 # Objective: linear part
    A::Vector{PointE{T, matT}}         # Ctr : linear operator
    b::Array{T, 1}                  #       rhs

    mr::Int                         # Real valued constraints are indexed from 1 to mr
    mc::Int                         # Clpx valued constraints are indexed from mr+1 to mc

    nsdpvar::Vector{Int}            # Size of the several SDP cones
    nscalvar::Int                   # Number of scalar variables (of type T)
    m::Int                          # Nb of constraints
    nc::Int

    options::OrderedDict
end

function SDCOContext(c::PointE{T, matT}, A::Vector{PointE{T, matT}}, b::Array{T, 1}; options=OrderedDict()) where {T<:Number, matT}
    @assert length(A) == length(b)
    m = length(A)

    nsdpvar = c.dims
    nscalvar = length(c.vec)

    # check all PointE elts are compatible
    for ai in A
        @assert nsdpvar == ai.dims
        @assert nscalvar == length(ai.vec)
    end

    nc = sum(nsdpvar) + nscalvar

    defoptions = get_defaultoption()

    for (opt, val) in options
        if haskey(defoptions, opt)
            defoptions[opt] = val
        else
            @warn "option $opt not supported"
        end
    end

    return SDCOContext(c, A, b, 0, 0, nsdpvar, nscalvar, m, nc, defoptions)
end

function get_defaultoption()
    return OrderedDict{Symbol, Any}(
                    :opt_ε => 1e-10,
                    :opt_maxit => 300,
                    :opt_outlev => 1,
    )
end


text:
using Statistics
using NCDatasets

using ArgParse
using JSON

include("nanop.jl")
include("PCA.jl")

using .PCA

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--data-file"
            help = "Ocean SSTA data file. New variable will be appended."
            arg_type = String
            required = true

        "--output-file"
            help = "Output file."
            arg_type = String
            required = true

        "--domain-file"
            help = "Domain file."
            arg_type = String
            required = true

        "--sparsity"
            help = "Because solving PCA matrix is time consuming. This parameter let user to use coarse data to derive PCA. Sparseity `n` means to skip every `n` other grid point in lon / lat. So sparsity 0 (default) means do not skip any data point. Sparsity 1 means means data density will drop to 1/4 of the original (1/2 * 1/2)."
            arg_type = Int64
            default = 0

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

sparsity = parsed["sparsity"]
inc = sparsity + 1

Dataset(parsed["domain-file"], "r") do ds

    global lon = replace(ds["xc"][:, 1], missing=>NaN)
    global lat = replace(ds["yc"][1, :], missing=>NaN)

    global llon = replace(ds["xc"][:], missing=>NaN)
    global llat = replace(ds["yc"][:], missing=>NaN)

    global mask = replace(ds["mask"][:], missing=>NaN)

    global Nx = length(lon)
    global Ny = length(lat)

    global valid_grids = []

    # Reference: https://www.esrl.noaa.gov/psd/enso/mei/
    # EOF of 30°S-30°N and 100°E-70°W

    global s_j = 1
    global n_j = 1
    global w_i = 1
    global e_i = 1

    for j = 1:Ny
        if - 30.0 <= lat[j]
            global s_j = j
            break
        end
    end
    for j = Ny:-1:s_j
        if 30.0 >= lat[j]
            global n_j = j
            break
        end
    end

    for i = 1:Nx
        if  100.0 <= lon[i]
            global w_i = i
            break
        end
    end
    for i = w_i+1:Nx
        if  290.0 <= lon[i]
            global e_i = i-1
            break
        end
    end

    for i = w_i:inc:e_i, j = s_j:inc:n_j
        if mask[i, j] == 0.0
            push!(valid_grids, (i, j))
        end
    end

end

Dataset(parsed["data-file"], "r") do ds

    global SSTA  = replace(ds["T_ML_MA"][:, :, 1, :], missing=>NaN)
    global months  = ds.dim["time"]

    global EOF_input = zeros(Float64, length(valid_grids), months)
    for i = 1:length(valid_grids)
        EOF_input[i, :] = SSTA[valid_grids[i][1], valid_grids[i][2], :]
    end

end

# El nino 3.4
# 170W - 120W  => [190 - 240]
# 5S - 5N

EN34idx = (abs.(llat) .<= 5) .& (llon .>= 190) .& (llon .<= 240)

EN34 = zeros(Float64, months)
for i = 1:months

    SSTA_slice = view(SSTA, :, :, i)
    EN34[i] = nanmean(SSTA_slice[EN34idx])

end

# Doing PCAs
modes = 2

println("Solving for PCA...")
eigen_vectors = PCA.findPCAs(EOF_input, num=modes)
println("done.")

PCAs = zeros(Float64, Nx, Ny, modes)
PCAs .= NaN

for i = 1:length(valid_grids)
    PCAs[valid_grids[i][1], valid_grids[i][2], :] = eigen_vectors[i, :]
end

ENSO = zeros(Float64, months)

ENSO_mode = view(PCAs, :, :, 1)
for i = 1:length(ENSO)
    ENSO[i] = nansum(ENSO_mode .* SSTA[:, :, i])
end

# Map grid back because of sparsity
if sparsity > 0
    for m = 1:modes, i = w_i:inc:e_i, j=s_j:inc:n_j
        if mask[i, j] == 0.0
            PCAs[i:min(e_i, i+sparsity), j:min(n_j, j+sparsity), m] .= PCAs[i, j, m]
        end
    end
end

Dataset(parsed["output-file"], "c") do ds

    defDim(ds, "Nx", Nx)
    defDim(ds, "Ny", Ny)
    defDim(ds, "time", Inf)
    defDim(ds, "modes", modes)

    for (varname, vardata, vardim, attrib) in [
        ("PCAs",  PCAs, ("Nx", "Ny", "modes",), Dict()),
        ("EN34",  EN34, ("time",), Dict()),
        ("ENSO",  ENSO, ("time",), Dict()),
    ]

        if ! haskey(ds, varname)
            var = defVar(ds, varname, Float64, vardim)
            var.attrib["_FillValue"] = 1e20
        end

        var = ds[varname]

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
"""
getbatch(dataset, i::Int, sz::Int)

Get batch of data from the specified dataset.

 * i - batch index
 * sz - size of the batch

Typical implementation for continuous datasets:

    start = (i-1)*sz + 1
    finish = min(i*sz, length(dataset))
    return start <= finish ? dataset[start:finish] : nothing
"""
function getbatch(a::AbstractArray{T,N}, i::Int, sz::Int; batch_dim=ndims(a)) where {T,N}
    start = (i-1)*sz + 1
    finish = min(i*sz, size(a, batch_dim))
    start > finish && return nothing
    selector = [d == batch_dim ? (start:finish) : (:) for d=1:ndims(a)]
    return @view a[selector...]
end

function getbatch(t::Tuple, i::Int, sz::Int)
    return ((getbatch(a, i, sz) for a in t)...,)
end

mutable struct BatchIter
    data
    sz::Int
end

Base.show(io::IO, it::BatchIter) = print(io, "BatchIter($(typeof(it.data)), sz=$(it.sz))")
batchiter(data; sz::Int=1) = BatchIter(data, sz)

function Base.iterate(it::BatchIter, i::Int=1)
    r = getbatch(it.data, i, it.sz)
    if r == nothing || (r isa Tuple && all(x -> x == nothing, r))
        return nothing
    else
        return r, i+1
    end
end

# TODO: implement `getbatch(data::ImageFolder, i, sz)` with threads

## Supervised learning with input X and output Y

# function partial_fit!(m, X::AbstractArray, Y::AbstractArray, full_loss_fn;
#                       opt=SGD(1e-3), batch_size=100, device=CPU())
#     epoch_loss = 0
#     for (i, (x, y)) in enumerate(batchiter((X, Y), size=batch_size))
#         x = to_device(device, copy(x))
#         y = to_device(device, copy(y))
#         loss, g = grad(full_loss_fn, m, x, y)
#         update!(opt, m, g[2])
#     end
#     return epoch_loss / size(X)[end]
# end

# function fit!(m, X::AbstractArray, Y::AbstractArray, loss_fn;
#               n_epochs=10, batch_size=100, opt=SGD(1e-3), device=CPU(), report_every=1)
#     full_loss_fn = (m, x, y) -> loss_fn(m(x), y)
#     for epoch in 1:n_epochs
#         time = @elapsed epoch_loss =
#             partial_fit!(m, X, Y, full_loss_fn, batch_size=batch_size, device=device)
#         if epoch % report_every == 0
#             println("Epoch $epoch: avg_cost=$(epoch_loss / (size(X,2) / batch_size)), elapsed=$time")
#         end
#     end
#     return m
# end

function fit!(m, X::AbstractArray, Y::AbstractArray, loss_fn;
              n_epochs=10, batch_size=100, opt=SGD(1e-3), device=CPU(), report_every=1)
    f = (m, x, y) -> loss_fn(m(x), y)
    num_batches = size(X)[end] // batch_size
    for epoch in 1:n_epochs
        epoch_loss = 0
        for (i, (x, y)) in enumerate(batchiter((X, Y), sz=batch_size))
            x = device(copy(x))
            y = device(copy(y))
            loss, g = grad(f, m, x, y)
            update!(opt, m, g[2])
            epoch_loss += loss
        end
        if epoch % report_every == 0
            println("Epoch $epoch: avg_loss=$(epoch_loss / num_batches)")
        end
    end
    return m
end

## Unsupervised learning with only X

function partial_fit!(m, X::AbstractArray, loss_fn;
                      opt=SGD(1e-3), batch_size=100, device=CPU())
    epoch_loss = 0
    f = (m, x) -> loss_fn(m(x))
    for (i, x) in enumerate(batchiter(X, sz=batch_size))
        x = device(copy(x))
        loss, g = grad(f, m, x)
        update!(opt, m, g[2])
        epoch_loss += loss
        println("iter $i: loss=$loss")
    end
    return epoch_loss
end

function fit!(m, X::AbstractArray, loss_fn;
              n_epochs=10, batch_size=100, opt=SGD(1e-3), device=CPU(), report_every=1)
    for epoch in 1:n_epochs
        time = @elapsed epoch_loss = partial_fit!(m, X, loss_fn, batch_size=batch_size, device=device)
        if epoch % report_every == 0
            @info("Epoch $epoch: avg_cost=$(epoch_loss / (size(X,2) / batch_size)), elapsed=$time")
        end
    end
    return m
end


text:
@testset "$TEST $G" begin

if !@isdefined(test_rem_edge)
    function test_rem_edge(g, srcv)
        for dstv in collect(neighbors(g, srcv))
            rem_edge!(g, srcv, dstv)
        end
        @test length(neighbors(g,srcv)) == 0
    end
end

g5 = DG(4)
add_edge!(g5,1,2); add_edge!(g5,2,3); add_edge!(g5,1,3); add_edge!(g5,3,4)

g6 = DG(4)
unsafe_add_edge!(g6,1,2); unsafe_add_edge!(g6,2,3);
unsafe_add_edge!(g6,1,3); unsafe_add_edge!(g6,3,4)
@test g6 == g5
rebuild!(g6)
@test g6 == g5

h1 = G(5)
h2 = G(3)
h3 = G()
h4 = DG(7)
h5 = DG()

e1 = E(1,2)
e2 = E(1,3)
e3 = E(1,4)
e4 = E(2,5)
e5 = E(3,5)

g = G(5)
@test add_edge!(g, 1, 2) == (true, edge(g, 1, 2))
for e in [e2,e3,e4,e5]
    @test add_edge!(g, e) == (true, edge(g, src(e), dst(e)))
end

g6 = G(5)
unsafe_add_edge!(g6, 1, 2)
unsafe_add_edge!(g6, e2)
unsafe_add_edge!(g6, e3)
unsafe_add_edge!(g6, e4)
unsafe_add_edge!(g6, e5)
@test g6 == g
rebuild!(g6)
@test g6 == g

h = DG(5)
@test add_edge!(h, 1, 2) == (true, edge(h, 1, 2))
for e in [e2,e3,e4,e5]
    @test add_edge!(h, e) == (true, edge(h, src(e), dst(e)))
end

@test vertices(g) == 1:5
i = 0
for e in edges(g)
    i+=1
end
@test i == 5
@test has_edge(g,3,5)
# @test edges(g) == Set([e1, e2, e3, e4, e5])
# @test Set{E}(edges(g)) == Set([e1, e2, e3, e4, e5])
#
@test degree(g) == [3, 2, 2, 1, 2]
@test in_degree(g) == [3, 2, 2, 1, 2]
@test in_degree(g,1) == 3
@test out_degree(g) == [3, 2, 2, 1, 2]
@test out_degree(g,1) == 3
@test degree(h) == [3, 1, 1, 0, 0]
@test in_degree(h) == [0, 1, 1, 1, 2]
@test in_degree(h,1) == 0
@test out_degree(h) == [3, 1, 1, 0, 0]
@test out_degree(h,1) == 3
@test collect(in_neighbors(h,5)) == adjacency_list(h,:in)[5]  == [2, 3]
@test collect(out_neighbors(h,1)) == adjacency_list(h)[1]  == [2, 3, 4]
@test_throws ErrorException adjacency_list(h, :ciao)
@test_throws MethodError adjacency_list(g, :out)

@test has_edge(g, 1, 2)
@test collect(in_edges(g, 2)) == [e1, reverse(e4)]
@test collect(out_edges(g, 1)) == [e1, e2, e3]
@test collect(out_edges(g, 1)) == collect(edges(g, 1))

@test add_vertex!(g) == nv(g) == 6
@test add_vertices!(g,5) == nv(g) == 11
@test has_vertex(g, 11)
@test ne(g) == 5
@test !is_directed(g)
@test is_directed(h)

@test CompleteGraph(4, G) == CompleteGraph(4, G)
@test CompleteGraph(4, G) != PathGraph(4, G)
@test CompleteDiGraph(4, DG) != PathDiGraph(4, DG)
@test CompleteDiGraph(4, DG) == CompleteDiGraph(4, DG)

@test collect(neighbors(g, 1)) == [2, 3, 4]

@test add_edge!(g, 1, 1) == (true, edge(g, 1, 1))
@test has_self_loops(g)
@test num_self_loops(g) == 1
@test add_edge!(g, 1, 1) == (false, edge(g, 1, 1))
@test rem_edge!(g, 1, 1)
@test !rem_edge!(g, 1, 1)

@test ne(g) == 5
@test rem_edge!(g, 1, 2)
@test ne(g) == 4

@test !rem_edge!(g, 2, 1)
add_edge!(g, 1, 2)
@test ne(g) == 5

@test has_edge(g,2,1)
@test has_edge(g,1,2)
@test rem_edge!(g, 2, 1)
@test add_edge!(h, 1, 1) == (true, edge(g, 1, 1))
@test rem_edge!(h, 1, 1)
@test rem_edge!(h, 1, 2)
@test !rem_edge!(h, 1, 2)
for v in vertices(g)
    test_rem_edge(copy(g),v)
end
for v in vertices(h)
    test_rem_edge(copy(h),v)
end

@test g == copy(g)
@test !(g === copy(g))
g10 = CompleteGraph(5, G)
@test rem_vertex!(g10, 1)
@test g10 == CompleteGraph(4, G)
@test rem_vertex!(g10, 4)
@test g10 == CompleteGraph(3, G)
@test !rem_vertex!(g10, 9)

g = G(2)
add_edge!(g, 1, 2)
@test rem_edge!(g, 1, 2)
@test !rem_edge!(g, 1, 2)
add_edge!(g, 1, 2)
edgs = collect(edges(g))
@test rem_edge!(g, edgs[1])
@test !rem_edge!(g, edgs[1])

g = CompleteGraph(5, G)
for e in edges(g)
    @test typeof(e) == edgetype(g)
end
edgs = collect(edges(g))
@test rem_edge!(g, edgs[1])
@test ne(g) == 9
@test !rem_edge!(g, edgs[1])
@test ne(g) == 9

g = G(10)
add_edge!(g, 1, 2)
add_edge!(g, 1, 3)
ga = G(10)
add_edge!(ga, 1, 2)
add_edge!(ga, 1, 4)
@test !(ga == g)
@test ga != g

g10 = CompleteDiGraph(5, DG)
for e in edges(g10)
    @test typeof(e) == edgetype(g10)
end
for i=2:5
    @test rem_edge!(g10, 1, i)
    @test ne(g10) == 21-i
end
@test in_degree(g10, 1) == 4
@test out_degree(g10, 1) == 0

g = CompleteDiGraph(5, DG)
@test ne(g) == 20
@test nv(g) == 5
clean_vertex!(g, 1)
@test ne(g) == 12
@test nv(g) == 5
@test in_degree(g, 1) == 0
@test out_degree(g, 1) == 0
for i=2:5
    @test in_degree(g, i) == 3
    @test out_degree(g, i) == 3
end

g10 = CompleteDiGraph(5, DG)
@test rem_vertex!(g10, 1)
@test g10 == CompleteDiGraph(4, DG)
rem_vertex!(g10, 4)
@test g10 == CompleteDiGraph(3, DG)
@test !rem_vertex!(g10, 9)
g10 = PathGraph(5, G)
@test rem_vertex!(g10, 5)
@test g10 == PathGraph(4, G)
@test rem_vertex!(g10, 4)
@test g10 == PathGraph(3, G)

g10 = PathDiGraph(5, DG)
@test rem_vertex!(g10, 5)
@test g10 == PathDiGraph(4, DG)
@test rem_vertex!(g10, 4)
@test g10 == PathDiGraph(3, DG)

g10 = PathDiGraph(5, DG)
@test rem_vertex!(g10, 1)
h10 = PathDiGraph(6, DG)
@test rem_vertex!(h10, 1)
@test rem_vertex!(h10, 1)
@test g10 == h10

g10 = CycleGraph(5, G)
@test rem_vertex!(g10, 5)
@test g10 == PathGraph(4, G)

g10 = PathGraph(3, G)
@test rem_vertex!(g10, 2)
@test g10 == G(2)

g10 = PathGraph(4, G)
@test rem_vertex!(g10, 3)
h10 =G(3)
@test add_edge!(h10,1,2) == (true, edge(h10, 1, 2))
@test g10 == h10

g10 = CompleteGraph(5, G)
@test rem_vertex!(g10, 3)
@test g10 == CompleteGraph(4, G)

s = split("$G",'.')[end]
if G <: ANetwork
    @test sprint(show, h1) == s*"(5, 0) with [] graph, [] vertex, [] edge properties"
else
    @test sprint(show, h1) == s*"(5, 0)"
end

g3 = PathGraph(5, G)
@test graph(digraph(g3)) == g3

@test degree(g3, 1) == 1
# @test neighbors(g3, 3) == [2, 4]
@test density(g3) == 0.4

g = G(5)
@test add_edge!(g, 1, 2) == (true, edge(g, 1, 2))

e2 = E(1,3)
e3 = E(1,4)
e4 = E(2,5)
e5 = E(3,5)

for e in [e2,e3,e4,e5]
    @test add_edge!(g, e) == (true, edge(g, src(e), dst(e)))
end

for i in out_neighbors(g, 1)
    @test typeof(i) == V
end

h = DG(10, 20)
for i=1:10
    u = sort(union(in_neighbors(h,i), out_neighbors(h,i)))
    @test sort(collect(all_neighbors(h, i))) == u
end
for i in out_neighbors(h, 1)
    @test typeof(i) == V
end

h = DG(5)
@test add_edge!(h, 1, 2) == (true, edge(h, 1, 2))
for e in [e2,e3,e4,e5]
    @test add_edge!(h, e) == (true, edge(h, src(e), dst(e)))
end

@test collect(out_edges(h, 3)) == collect(edges(h, 3))
@test length(collect(out_edges(h, 3))) == 1
@test length(collect(in_edges(h, 3))) == 1
@test length(collect(all_edges(h, 3))) == 2

for i=1:5
    u = sort(union(in_neighbors(h,i), out_neighbors(h,i)))
    @test sort(collect(all_neighbors(h, i))) == u
end

@test adjacency_list(g)[1] == collect(out_neighbors(g,1)) ==
    collect(in_neighbors(g,1)) == [2,3,4]

e0 = E(2, 3)
s = split("$DG",'.')[end]
# @test sprint(show, h4) == s*"(7, 0)" #TODO
# @test sprint(show, h5) == s*"(0, 0)"
@test has_edge(g, e1)
@test has_edge(h, e1)
@test !has_edge(g, e0)
@test !has_edge(h, e0)

g4 = PathDiGraph(5, DG)
@test degree(g4, 1) == 1
# @test neighbors(g4, 3) == [4]
@test density(g4) == 0.2

adjmx1 = [0 1 0; 1 0 1; 0 1 0] # graph
adjmx2 = [0 1 0; 1 0 1; 1 1 0] # digraph
a1 = G(adjmx1)
a2 = DG(adjmx2)

@test nv(a1) == 3
@test ne(a1) == 2
@test nv(a2) == 3
@test ne(a2) == 5

badadjmx = [ 0 1 0; 1 0 1]
@test_throws ErrorException G(badadjmx)
@test_throws ErrorException G(sparse(badadjmx))
@test_throws ErrorException DG(badadjmx)
@test_throws ErrorException DG(sparse(badadjmx))

@test add_edge!(g, 100, 100) == (false, edge(g, 100, 100))
@test add_edge!(h, 100, 100) == (false, edge(h, 100, 100))

g = G(sparse(adjmx1))
h = DG(sparse(adjmx1))

@test (nv(g), ne(g)) == (3, 2)
@test (nv(h), ne(h)) == (3, 4)
@test graph(h) == g

@test collect(neighbors(WheelDiGraph(10, DG),2)) ==
    collect(out_neighbors(WheelDiGraph(10, DG),2))
@test collect(out_neighbors(WheelDiGraph(10, DG),2)) == [3]

gg = DG(4)
add_edge!(gg, 1, 2)
add_edge!(gg, 1, 3)

g = reverse(gg)
@test ne(g) == ne(gg)
@test E(2,1) in edges(g)
@test !(E(1,2) in edges(g))
@test E(3,1) in edges(g)
@test !(E(1,3) in edges(g))

reverse!(g)
@test g == gg

g = G(4)
add_edge!(g, 1, 2); add_edge!(g, 3, 4);

swap_vertices!(g, 1, 3)
@test ne(g) == 2
@test has_edge(g, 1, 4)
@test has_edge(g, 3, 2)
@test !has_edge(g, 1, 2)
@test !has_edge(g, 3, 4)

add_edge!(g, 1, 2)
@test ne(g) == 3
swap_vertices!(g, 1, 3)
h = G(4)
add_edge!(h, 1, 2); add_edge!(h, 3, 4); add_edge!(h, 3, 2);
@test g == h

g = DG(4)
add_edge!(g, 1, 2); add_edge!(g, 3, 4);

swap_vertices!(g, 1, 3)
@test ne(g) == 2
@test has_edge(g, 1, 4)
@test has_edge(g, 3, 2)
@test !has_edge(g, 1, 2)
@test !has_edge(g, 3, 4)

add_edge!(g, 1, 2)
swap_vertices!(g, 1, 3)
@test !has_edge(g, 1, 4)
@test has_edge(g, 3, 2)
@test has_edge(g, 1, 2)
@test has_edge(g, 3, 4)
@test !has_edge(g, 2, 3)
@test !has_edge(g, 2, 1)
@test !has_edge(g, 4, 3)

g = G(2)
add_edge!(g, 1, 2); add_edge!(g, 1, 1);
swap_vertices!(g, 1, 2)
@test ne(g) == 2
@test has_edge(g, 1, 2)
@test has_edge(g, 2, 2)
@test !has_edge(g, 1, 1)
add_edge!(g, 1, 1)
@test has_edge(g, 1, 1)
@test ne(g) == 3
swap_vertices!(g, 1, 2)
@test ne(g) == 3
@test has_edge(g, 1, 2)
@test has_edge(g, 2, 2)
@test has_edge(g, 1, 1)

g = DG(3)
add_edge!(g, 1, 2); add_edge!(g, 1, 1); add_edge!(g, 3, 1)
swap_vertices!(g, 1, 2)
@test ne(g) == 3
@test !has_edge(g, 1, 2)
@test has_edge(g, 2, 1)
@test has_edge(g, 2, 2)
@test !has_edge(g, 1, 1)
@test has_edge(g, 3, 2)
@test !has_edge(g, 3, 1)
@test !has_edge(g, 1, 3)
@test !has_edge(g, 2, 3)
add_edge!(g, 1, 1)
swap_vertices!(g, 1, 2)
@test ne(g) == 4
@test has_edge(g, 1, 2)
@test !has_edge(g, 2, 1)
@test has_edge(g, 2, 2)
@test has_edge(g, 1, 1)
@test !has_edge(g, 3, 2)
@test has_edge(g, 3, 1)
@test !has_edge(g, 1, 3)
@test !has_edge(g, 2, 3)

g = WheelGraph(10, G)
@test degree(g, 1) == 9
@test degree(g, 10) == 3
swap_vertices!(g, 1, 10)
@test degree(g, 1) == 3
@test degree(g, 10) == 9

g = CompleteGraph(10, G)
vmap = rem_vertices!(g, 6:10)
@test g == CompleteGraph(5, G)
@test typeof(vmap) <: AVertexMap
@test length(vmap) == 0

g = CompleteGraph(10, G)
vmap = rem_vertices!(g, 1:5)
@test g == CompleteGraph(5, G)
@test typeof(vmap) <: AVertexMap
@test length(vmap) ==  5
for v=1:5
    @test vmap[v] == v + 5
end

g = G(10, 20)
g2 = copy(g)
vmap = rem_vertices!(g, 6:10)
for v=10:-1:6
    rem_vertex!(g2, v)
end
@test g == g2
@test typeof(vmap) <: AVertexMap
@test length(vmap) == 0

g = DG(10, 20)
g2 = copy(g)
vmap = rem_vertices!(g, 1:5)
for v=5:-1:1
    rem_vertex!(g2, v)
end
@test g == g2
@test typeof(vmap) <: AVertexMap
@test length(vmap) ==  5
for v=1:5
    @test vmap[v] == v + 5
end

g1 = G(10, 30)
g2 = copy(g1)
@test rem_vertices!(g1, 1, 2, 2) == rem_vertices!(g2, 1, 2)
@test g1 == g2

g = G(10,20)
h = copy(g)
pop_vertex!(g)
rem_vertex!(h, 10)
@test g == h

g = DG(10,20)
h = copy(g)
pop_vertex!(g)
rem_vertex!(h, 10)
@test g == h

A = [1 2 0;
    -1 0 1;
     3  1 2]
g = G(A)
@test nv(g) == 3
@test ne(g) == 5
g = G(A, upper=true)
@test ne(g) == 4
@test !has_edge(g, 1, 3)
@test !has_edge(g, 2, 2)
@test has_edge(g, 1, 1)
@test !has_edge(g, 15, 1)
g = G(A, upper=true, selfedges=false)
@test ne(g) == 2
@test !has_edge(g, 1, 3)
@test !has_edge(g, 2, 2)
@test !has_edge(g, 1, 1)
@test !has_edge(g, 3, 3)
g = G(A, upper=false, selfedges=false)
@test ne(g) == 3
@test has_edge(g, 1, 3)
@test !has_edge(g, 2, 2)
@test !has_edge(g, 1, 1)
@test !has_edge(g, 3, 3)

g = DG(A)
@test nv(g) == 3
@test ne(g) == 7
@test has_edge(g, 1, 1)
@test !has_edge(g, 2, 2)
@test has_edge(g, 3, 3)
@test !has_edge(g, 15, 1)

g = DG(A, selfedges=false)
@test ne(g) == 5
@test !has_edge(g, 1, 1)
@test !has_edge(g, 2, 2)
@test !has_edge(g, 3, 3)

end #testset


text:
module PosteriorSummary

export effective_sample_size,
       hpd_interval

import MCMCDiagnostics.effective_sample_size

function hpd_interval(x::AbstractVector{Float64}; alpha::Float64 = 0.05, sorted::Bool = false)
    n = length(x)
    span = Int(round((1.0 - alpha) * n))
    smallest_rng = Inf
    rng_upper = Inf
    rng_lower = -Inf
    sx::AbstractVector{Float64} = sorted ? x : sort(x)
    lower = 0
    for upper = span:n
        lower = upper - span + 1
        rng = sx[upper] - sx[lower]
        if rng < smallest_rng
            smallest_rng = rng
            rng_lower = sx[lower]
            rng_upper = sx[upper]
        end
    end

    return (lower = rng_lower, upper = rng_upper)
end

end


text:
Random.seed!(91247)
rng1 = MersenneTwister(817)
rng2 = deepcopy(rng1)

@test_throws DomainError PermutedBlock([-0.5, 0.5], 2)
@test_throws DomainError PermutedBlock([-1, 2], 2)
@test_throws DomainError PermutedBlock([1, 2], 4)
@test_throws DomainError PermutedBlock([1, 2, 2], 3)

p0 = [1 / 3, 0, 1 / 3, 1 / 3]
w1 = 1:4
PB0 = PermutedBlock(p0, 3)
PB1 = PermutedBlock(w1, 10)

@test narms(PB0) == 4
@test target(PB0) == p0
@test prob(PB0) ≈ target(PB0)
@test sequence(PB0) == zeros(0)
@test dist(PB0) == zeros(narms(PB0))
@test imbalance(PB0) == 0.0
@test predictability(PB0) == 0.0
@test weights(PB0) == round.(Int, blocksize(PB0) .* target(PB0))

@test target(PB1) ≈ collect(w1 ./ sum(w1))
@test prob(PB1) ≈ target(PB1)

@test begin
    p = [10, 14, 17]
    b = 41
    PermutedBlock(p, b) == PermutedBlock(p ./ sum(p), b)
end

# Check that randomise! does alter state
# Both intialises lead to same sequence
w = 1:4
PB1 = PermutedBlock(w, 20)
PB2 = PermutedBlock(w ./ sum(w), 20)
res1 = randomise!(rng1, PB1, 20)
res2 = randomise!(rng2, PB2, 20)
@test sequence(PB1) == res1
@test sequence(PB2) == res2
@test prob(PB1) ≈ prob(PB2)
@test dist(PB1) == dist(PB2)
@test imbalance(PB0) == 0.0
randomise!(PB1, 19)
@test maximum(prob(PB1)) ≈ 1.0


text:
#################################################################################
#   2D Cavity Flow Problem using NeuralPDE
#   Reference: https://nbviewer.jupyter.org/github/barbagroup/CFDPython/blob/master/lessons/14_Step_11.ipynb
#
#################################################################################
# Andriy: Created by Emmanuel Lujan from Julia Slack
#################################################################################
using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux
using Quadrature, Cubature, Cuba

@parameters t,x,y
@variables u(..),v(..),p(..)
@derivatives Dt'~t
@derivatives Dx'~x
@derivatives Dxx''~x
@derivatives Dy'~y
@derivatives Dyy''~y

# Parameters

rho = 1.0
nu = 0.1
t_max = 0.05 #nit * dt
x_max = 2.0
y_max = 2.0

# Equations, initial and boundary conditions

eqs = [
        ( Dt(u(t,x,y)) + u(t,x,y) * Dx(u(t,x,y)) + v(t,x,y) * Dy(u(t,x,y)) ~
          - 1.0 / rho * Dx(p(t,x,y)) + nu * (Dxx(u(t,x,y)) + Dyy(u(t,x,y))) ),

        ( Dt(v(t,x,y)) + u(t,x,y) * Dx(v(t,x,y)) + v(t,x,y) * Dy(v(t,x,y)) ~
          - 1.0 / rho * Dy(p(t,x,y)) + nu * (Dxx(v(t,x,y)) + Dyy(v(t,x,y))) ),

        ( Dxx(p(t,x,y)) + Dyy(p(t,x,y)) ~ -rho * ( Dx(u(t,x,y)) * Dx(u(t,x,y))
          + 2.0 * Dy(u(t,x,y)) * Dx(v(t,x,y)) + Dy(v(t,x,y)) * Dy(v(t,x,y)) ) )
      ]

bcs = [
        u(0.0,x,y)~0.0,
        v(0.0,x,y)~0.0,
        p(0.0,x,y)~0.0,

        u(t,0.0,y)~0.0,
        u(t,x_max,y)~0.0,
        u(t,x,0.0)~0.0,
        u(t,x,y_max)~1.0,

        v(t,0.0,y)~0.0,
        v(t,x_max,y)~0.0,
        v(t,x,0.0)~0.0,
        v(t,x,y_max)~0.0,

        p(t,x,y_max)~0.0,
        Dy(p(t,x,0.0))~0.0,
        Dx(p(t,0.0,y))~0.0,
        Dx(p(t,x_max,y))~0.0 ]

# Space and time domains
domains = [t ∈ IntervalDomain(0.0,t_max),
           x ∈ IntervalDomain(0.0,x_max),
           y ∈ IntervalDomain(0.0,y_max)]

# Discretization
nx = 41
ny = 41
dx = 2.0 / (nx - 1)
dy = 2.0 / (ny - 1)
dt = 0.001

# Neural network
dim = length(domains)
output = length(eqs)
neurons = 12
chain1 = FastChain( FastDense(dim,neurons,Flux.σ),
                    FastDense(neurons,neurons,Flux.σ),
                    FastDense(neurons,1))
chain2 = FastChain( FastDense(dim,neurons,Flux.σ),
                    FastDense(neurons,neurons,Flux.σ),
                    FastDense(neurons,1))
chain3 = FastChain( FastDense(dim,neurons,Flux.σ),
                    FastDense(neurons,neurons,Flux.σ),
                    FastDense(neurons,1))

strategy = GridTraining(dx=dx)
#strategy = GridTraining(dx=[dt,dx])

discretization = PhysicsInformedNN([chain1,chain2,chain3],strategy=strategy)

pde_system = PDESystem(eqs,bcs,domains,[t,x,y],[u,v,p])
prob = discretize(pde_system,discretization)

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

res = GalacticOptim.solve(prob,Optim.BFGS();cb=cb,maxiters=1000)
#res = GalacticOptim.solve(prob, ADAM(0.001), cb = cb, maxiters=5000)

# Plots

phi = discretization.phi
initθ = discretization.initθ
acum =  [0;accumulate(+, length.(initθ))]
sep = [acum[i]+1 : acum[i+1] for i in 1:length(acum)-1]
minimizers = [res.minimizer[s] for s in sep]
ts,xs,ys = [domain.domain.lower:dx:domain.domain.upper for domain in domains]
u_predict = reshape([ phi[1]([t_max,x,y],minimizers[1])[1] for x in xs for y in ys ], length(xs), length(ys))
v_predict = reshape([ phi[2]([t_max,x,y],minimizers[2])[1] for x in xs for y in ys ], length(xs), length(ys))
p_predict = reshape([ phi[3]([t_max,x,y],minimizers[3])[1] for x in xs for y in ys ], length(xs), length(ys))

using Plots, LaTeXStrings
uv(x, y) = (u_predict[Int(round(x/dx))+1,Int(round(y/dy))+1],
            v_predict[Int(round(x/dx))+1,Int(round(y/dy))+1])
xxs = [x for x in xs for y in ys]
yys = [y for x in xs for y in ys]
plot(xs, ys, p_predict, linetype=:contourf, title = L"$Cavity Flow$")
quiver!(xxs, yys, quiver=uv)
savefig("2d-cavity-flow.svg")

max_stars_repo_name: drozzy/NeuralOceananigans

max_issues_repo_name: drozzy/NeuralOceananigans

max_forks_repo_name: drozzy/NeuralOceananigans


text:
#
# Graph Functionality
#

function _show_conns(io, m, comp_name, which::Symbol)
    datumtype = which == :incoming ? "parameters" : "variables"
    println(io, "   $which $datumtype:")

    conns = get_connections(m, comp_name, which)

    if length(conns) == 0
        println(io, "      none")
    else
        for conn in conns
            if which == :incoming
                println(io, "      - $(conn.src_comp_name).$(conn.dst_par_name)")
            else
                println(io, "      - $(conn.dst_comp_name).$(conn.src_var_name)")
            end
        end
    end
end

function show(io::IO, m::Model)
    println(io, "Model component connections:")

    for (i, comp_name) in enumerate(compkeys(m.md))
        comp_def = compdef(m.md, comp_name)
        println(io, "$i. $(comp_def.comp_id) as :$(comp_def.name)")
        _show_conns(io, m, comp_name, :incoming)
        _show_conns(io, m, comp_name, :outgoing)
    end
end

function get_connections(m::Model, ci::ComponentInstance, which::Symbol)
    return get_connections(m, name(ci.comp), which)
end

function _filter_connections(conns::Vector{InternalParameterConnection}, comp_name::Symbol, which::Symbol)
    if which == :all
        f = obj -> (obj.src_comp_name == comp_name || obj.dst_comp_name == comp_name)
    elseif which == :incoming
        f = obj -> obj.dst_comp_name == comp_name
    elseif which == :outgoing
        f = obj -> obj.src_comp_name == comp_name
    else
        error("Invalid parameter for the 'which' argument; must be 'all' or 'incoming' or 'outgoing'.")
    end

    return collect(Iterators.filter(f, conns))
end

function get_connections(m::Model, comp_name::Symbol, which::Symbol)
    return _filter_connections(internal_param_conns(m.md), comp_name, which)
end

function get_connections(mi::ModelInstance, comp_name::Symbol, which::Symbol)
    return _filter_connections(internal_param_conns(mi.md), comp_name, which)
end


text:
module Searching
export bfs_parents, bfs_tree, dfs_parents, dfs_tree

using ...CSetDataStructures, ..BasicGraphs

"""
    tree(parents)
Convert a parents array into a directed graph.
"""
function tree(parents::AbstractVector{Int})
    n = T(length(parents))
    t = Graph(n)
    for (v, u) in enumerate(parents)
        if u > 0 && u != v
            add_edge!(t, u, v)
        end
    end
    return t
end

"""
    bfs_parents(g, s[; dir=:out])

Perform a breadth-first search of graph `g` starting from vertex `s`.
Return a vector of parent vertices indexed by vertex. If `dir` is specified,
use the corresponding edge direction (`:in` and `:out` are acceptable values).

### Performance
This implementation is designed to perform well on large graphs. There are
implementations which are marginally faster in practice for smaller graphs,
but the performance improvements using this implementation on large graphs
can be significant.
"""
bfs_parents(g::ACSet, s::Int; dir = :out) =
    (dir == :out) ? _bfs_parents(g, s, outneighbors) : _bfs_parents(g, s, inneighbors)

function _bfs_parents(g::ACSet, source, neighborfn::Function)
    n = nv(g)
    visited = falses(n)
    parents = zeros(Int, nv(g))
    cur_level = Int[]
    sizehint!(cur_level, n)
    next_level = Int[]
    sizehint!(next_level, n)
    @inbounds for s in source
        visited[s] = true
        push!(cur_level, s)
        parents[s] = s
    end
    while !isempty(cur_level)
        @inbounds for v in cur_level
            @inbounds @simd for i in neighborfn(g, v)
                if !visited[i]
                    push!(next_level, i)
                    parents[i] = v
                    visited[i] = true
                end
            end
        end
        empty!(cur_level)
        cur_level, next_level = next_level, cur_level
        sort!(cur_level)
    end
    return parents
end

"""
    bfs_tree(g, s[; dir=:out])
Provide a breadth-first traversal of the graph `g` starting with source vertex `s`,
and return a directed acyclic graph of vertices in the order they were discovered.
If `dir` is specified, use the corresponding edge direction (`:in` and `:out` are
acceptable values).
"""
bfs_tree(g::ACSet, s::Integer; dir = :out) = tree(bfs_parents(g, s; dir = dir))

"""
    dfs_parents(g, s[; dir=:out])

Perform a depth-first search of graph `g` starting from vertex `s`.
Return a vector of parent vertices indexed by vertex. If `dir` is specified,
use the corresponding edge direction (`:in` and `:out` are acceptable values).

### Implementation Notes
This version of DFS is iterative.
"""
dfs_parents(g::ACSet, s::Integer; dir=:out) =
    (dir == :out) ? _dfs_parents(g, s, outneighbors) : _dfs_parents(g, s, inneighbors)

function _dfs_parents(g::ACSet, s::Int, neighborfn::Function)
    parents = zeros(Int, nv(g))
    seen = zeros(Bool, nv(g))
    S = [s]
    seen[s] = true
    parents[s] = s
    while !isempty(S)
        v = S[end]
        u = 0
        for n in neighborfn(g, v)
            if !seen[n]
                u = n
                break
            end
        end
        if u == 0
            pop!(S)
        else
            seen[u] = true
            push!(S, u)
            parents[u] = v
        end
    end
    return parents
end

"""
    dfs_tree(g, s)

Return a directed acyclic graph based on
depth-first traversal of the graph `g` starting with source vertex `s`.
"""
dfs_tree(g::AbstractGraph, s::Integer; dir=:out) = tree(dfs_parents(g, s; dir=dir))

end

max_forks_repo_name: epatters/CompCat


text:
include("runtests-prerequisite.jl")

@testset "Test utils.jl - function `string2enum()`" begin
    @enum Fruit apple=10 orange=20 kiwi=30
    @test string2enum(Fruit, "orange") == orange
end

@testset "Test utils.jl - function `int2enum()`" begin
    @enum Fruit apple=10 orange=20 kiwi=30
    @test int2enum(Fruit, 20) == orange
    @test int2enum(Fruit, 30) == kiwi
end

@testset "Test utils.jl - function `dict2vectoroftuples()`" begin
    props = Dict(:a => 1, :b => missing, :c => "de", "d" => 9)
    PostgresORMUtil.dict2vectoroftuples(props)
end

@testset "Test utils.jl - function `vectoroftuples2dict()`" begin
    vectoroftuples = [(:a, 1),(:b, missing),(:c, "de"),("d", 9)]
    PostgresORMUtil.vectoroftuples2dict(vectoroftuples)
end

@testset "Test utils.jl - function `dictstringkeys2symbol()`" begin
    @test dictstringkeys2symbol(Dict("company" => String,
                       "canpost" => Bool)) == Dict(:company => String,
                                          :canpost => Bool)
end

@testset "Test utils.jl - function `dict2namedtuple()`" begin
    @test dict2namedtuple(Dict(:company => String,
                       :canpost => Bool)) == (company = String, canpost = Bool)
end

@testset "Test utils.jl - function `namedtuple2dict()`" begin
    @test namedtuple2dict((company = String,
                           canpost = Bool)) == Dict(:company => String,
                                              :canpost => Bool)
end

@testset "Test utils.jl - function `tovector()`" begin
    @test tovector(("1","2")) == ["1","2"]
    @test tovector(("1")) == ["1"]
    @test_throws ArgumentError tovector((a = "1", b= "2"))
    @test tovector(:toto) == [:toto]
    @test tovector("toto";elementstype = Symbol) == [:toto]
    @test tovector((:toto,:tata);elementstype = Symbol) == [:toto,:tata]
    @test tovector(:toto;elementstype = Symbol) == [:toto]
end

@testset "Test utils.jl - function `get_nonmissing_typeof_uniontype()`" begin
    @test get_nonmissing_typeof_uniontype(Union{String,Missing}) == String
    @test get_nonmissing_typeof_uniontype(Union{Missing,String,Nothing}) == String
    @test get_nonmissing_typeof_uniontype(Union{Bool,Missing,Nothing}) == Bool
    @test get_nonmissing_typeof_uniontype(Union{Vector{String},Missing,Nothing}) == Vector{String}
end

@testset "Test utils.jl - function `dataframe2vector_of_namedtuples()`" begin

    df = DataFrame(SepalLength = [5.1,3],
                    SepalWidth = [3.5,4],
                    Species = ["setosa","versicolor"])

    df = @from i in df begin
         @select {i.SepalLength,i.SepalWidth,i.Species}
         @collect DataFrame
    end

    result = dataframe2vector_of_namedtuples(df)
    @test result[1] == (SepalLength = 5.1, SepalWidth = 3.5, Species = "setosa")

end

@testset "Test utils.jl - function `diff_dict()`" begin

    old_dict = Dict("attr1" => 4,"attr2" => 2, "attr4" => 8, "attr5" => missing, "attr6" => missing)
    new_dict = Dict("attr1" => 5,"attr3" => 6, "attr4" => missing, "attr5" => 7, "attr6" => missing)

    (ismissing(old_dict["attr6"])  == ismissing(new_dict["attr6"]))

    diff_result  = diff_dict(old_dict, new_dict)

    @test diff_result["attr1"] == (old = 4, new = 5)
    @test diff_result["attr2"].old == 2
    @test ismissing(diff_result["attr2"].new)
    @test ismissing(diff_result["attr3"].old)
    @test diff_result["attr3"].new == 6

end

@testset "Test utils.jl - function `string2zoneddatetime`" begin
    PostgresORMUtil.string2zoneddatetime("2019-09-03T11:00:00.000Z")
end

@testset "Test utils.jl - function `postgresql_string_array_2_string_vector`" begin
    str = "{\"Deleted. 3 Scenes (a,2,e)\",\"Behind the Scenes\",\"\"}"
    res = PostgresORMUtil.postgresql_string_array_2_string_vector(str)
    @test res == ["Deleted. 3 Scenes (a,2,e)","Behind the Scenes",""]

    str = "{aa,\"Deleted. 3 Scenes (a,2,e)\",bb,\"Behind the Scenes\",\"\",cc}"
    res = PostgresORMUtil.postgresql_string_array_2_string_vector(str)
    @test res == ["aa","Deleted. 3 Scenes (a,2,e)","bb","Behind the Scenes","","cc"]

    str = "{aa,bb,cc,dd}"
    res = PostgresORMUtil.postgresql_string_array_2_string_vector(str)
    @test res == ["aa","bb","cc","dd"]

    # This one fails
    str = "{\"bla bla\",bb,cc}"
    res = PostgresORMUtil.postgresql_string_array_2_string_vector(str)
    @test res == ["bla bla","bb","cc"]

end

@testset "Test utils.jl - function `getpropertiesvalues`" begin

    mutable struct MyNewStruct <: PostgresORM.IEntity

      actor_id::Union{Missing,Int32}
      first_name::Union{Missing,String}

      MyNewStruct(args::NamedTuple) = MyNewStruct(;args...)
      MyNewStruct(;
            actor_id = missing,
            first_name = missing,
     ) = (
                x = new(missing, missing);
                x.actor_id = actor_id;
                x.first_name = first_name;
                return x
                )
    end

    myNewStruct = MyNewStruct(actor_id = 3,first_name="toto")

    @test PostgresORMUtil.getpropertiesvalues(myNewStruct,[:actor_id,:first_name]) == [3,"toto"]
end

@testset "Test utils.jl - function `setpropertiesvalues`" begin

    mutable struct MyNewStruct <: PostgresORM.IEntity

      actor_id::Union{Missing,Int32}
      first_name::Union{Missing,String}

      MyNewStruct(args::NamedTuple) = MyNewStruct(;args...)
      MyNewStruct(;
            actor_id = missing,
            first_name = missing,
     ) = (
                x = new(missing, missing);
                x.actor_id = actor_id;
                x.first_name = first_name;
                return x
                )
    end

    myNewStruct = MyNewStruct(actor_id = 3,first_name="toto")
    PostgresORMUtil.setpropertiesvalues!(myNewStruct,[:actor_id,:first_name],[4,"tata"])
    @test myNewStruct.actor_id == 4
    @test myNewStruct.first_name == "tata"
end

@testset "Test utils.jl - function `remove_spaces_and_split`" begin
    @test PostgresORMUtil.remove_spaces_and_split("rfr ") == ["rfr"]
    PostgresORMUtil.remove_spaces_and_split("  rfr, ded") == ["rfr","ded"]
end


text:
using StataDTAFiles, Test, Dates, StrFs
using StataDTAFiles: LSF, verifytag, read_header, read_map, read_variable_types,
    read_variable_names, read_sortlist, read_formats, TIMESTAMPFMT, readrow
using Parameters: @unpack
import Tables

testdata = joinpath(@__DIR__, "data", "testdata.dta")

≅(a, b) = a == b
≅(a::AbstractVector, b::AbstractVector) = all(a .≅ b)
≅(a::Tuple, b::Tuple) = all(a .≅ b)
≅(a::NamedTuple, b::NamedTuple) = all(map(≅, a,  b))
≅(::Missing, ::Missing) = true

@testset "reading tags" begin
    @test verifytag(IOBuffer("<atag>"), "atag") == nothing
    @test verifytag(IOBuffer("</btag>"), "btag", true) == nothing
    @test_throws ErrorException verifytag(IOBuffer("noopening>"), "noopening")
    @test_throws EOFError verifytag(IOBuffer("<eof"), "eof")
    @test_throws ErrorException verifytag(IOBuffer("<a>"), "b")
end

@testset "reading header" begin
    dta = open(DTAFile, testdata)
    @test dta.header.release == 118
    @test dta.header.variables == 3
    @test dta.header.observations == 10
    @test dta.header.label == ""
    @test dta.map.eof == filesize(testdata)
    @test Base.IteratorEltype(typeof(dta)) ≡ Base.HasEltype()
    @test eltype(dta) == NamedTuple{(:a, :b, :c), Tuple{Union{Missing, Float32},
                                                        Union{Missing, Int16},
                                                        StrF{2}}}
    @test Base.IteratorSize(typeof(dta)) ≡ Base.HasLength()
    @test length(dta) == 10
    @test dta.sortlist == []
    @test dta.formats == ["%9.0g", "%9.0g", "%9s"]
    @test @inferred(collect(dta)) ≅ [(a = i > 7 ? missing : Float32(i),
                                     b = i ≤ 2 ? missing : Int16(i),
                                     c = StrF{2}(string(i))) for i in 1:10]
    close(dta)
end

@testset "reading header" begin
    str = open(repr, DTAFile, testdata)
    r_header = raw"^Stata DTA file 118, 3 vars in 10 rows, .*\n\s+not sorted\n"
    r_vars = raw"\s+a::Union\{Missing,\s*Float32\}.*\n\s+b::Union\{Missing,\s*Int16\}.*\n\s+c::StrF\{2\}.*$"
    @test occursin(Regex(r_header * r_vars), str)
end

@testset "type stability" begin
    # Test that readrow is inferred. Tests use internals, which are not part of the API.
    dta = open(DTAFile, testdata)
    @unpack boio = dta
    seek(boio, dta.map.data)
    verifytag(boio, "data")
    for _ in 1:length(dta)
        @test @inferred(readrow(boio, eltype(dta))) isa NamedTuple
    end
end

@testset "timestamp parsing" begin
    @test DateTime("04 Jul 2032 04:23", TIMESTAMPFMT) ==
        DateTime("04 Jul 2032 04:23", TIMESTAMPFMT) ==
        DateTime("2032-07-04T04:23:00")
end

@testset "elapsed dates" begin
    @test elapsed_days(0) == Date(1960, 1, 1)
    @test elapsed_days(1) == Date(1960, 1, 2)
    @test elapsed_days(-1) == Date(1959, 12, 31)
    @test elapsed_days(365) == Date(1960, 12, 31)
end

@testset "tables interface" begin
    dta = open(DTAFile, testdata)
    @test Tables.istable(typeof(dta)) == true
    @test Tables.rowaccess(typeof(dta)) == true
    @test Tables.rows(dta) ≡ dta
    @test Tables.schema(dta) ≡ Tables.Schema((:a, :b, :c),
                                             (Union{Missing, Float32},
                                              Union{Missing, Int16},
                                              StrF{2}))
    @test Tables.columntable(dta) ≅
        (a = Union{Missing, Float32}[1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, missing, missing, missing],
         b = Union{Missing, Int16}[missing, missing, 3, 4, 5, 6, 7, 8, 9, 10],
         c = StrF{2}["1", "2", "3", "4", "5", "6", "7", "8", "9", "10"])
end


text:
using Documenter
using PowerSystems
using PowerSimulations
using Literate
using PowerGraphics
using DataStructures
const PG = PowerGraphics

folders = Dict(
#  "Operations" => readdir("docs/src/Operations"),
#  "Simulations" => readdir("docs/src/Simulations"),
# "PowerGraphics" => readdir("docs/src/PowerGraphics"),
)

for (name, folder) in folders
    for file in folder
        outputdir = joinpath(pwd(), "docs/src/howto")
        inputfile = joinpath(pwd(), "docs/src/$name/$file")
        Literate.markdown(inputfile, outputdir)
    end
end
if isfile("docs/src/howto/.DS_Store.md")
    rm("docs/src/howto/.DS_Store.md")
end

makedocs(
    sitename = "PowerSimulations.jl",
    format = Documenter.HTML(
        mathengine = Documenter.MathJax(),
        prettyurls = get(ENV, "CI", nothing) == "true",
    ),
    modules = [PowerSimulations],
    strict = true,
    authors = "Jose Daniel Lara, Clayton Barrows and Dheepak Krishnamurthy",
    pages = Any[
        "Introduction" => "index.md",
        #"Quick Start Guide" => "qs_guide.md",
        "Logging" => "man/logging.md",
        "Operation Model" => "man/op_problem.md",
        # "How To" => Any[
        #     "Set Up Plots" => "howto/3.0_set_up_plots.md",
        #     "Make Stack Plots" => "howto/3.1_make_stack_plots.md",
        #     "Make Bar Plots" => "howto/3.2_make_bar_plots.md",
        #     "Make Fuel Plots" => "howto/3.3_make_fuel_plots.md",
        #     "Make Forecast Plots" => "howto/3.4_make_forecast_plots.md",
        #     "Plot Fewer Variables" => "howto/3.5_plot_fewer_variables.md",
        #     "Plot Multiple Results" => "howto/3.6_plot_multiple_results.md",
        # ],
        "Simulation Recorder" => "man/simulation_recorder.md",
        "Model References" => Any["Hydro Models" => "ref/hydro.md"],
        "API" => Any["PowerSimulations" => "api/PowerSimulations.md"],
    ],
)

deploydocs(
    repo = "github.com/NREL-SIIP/PowerSimulations.jl.git",
    target = "build",
    branch = "gh-pages",
    devbranch = "master",
    devurl = "dev",
    versions = ["stable" => "v^", "v#.#"],
)


text:
using HTTP, JSON

const HOST = "127.0.0.1"
const PORT = 8080
const ROUTER = HTTP.Router()

function show_factorial(req::HTTP.Request)
  headers = [                                       # 7
    "Access-Control-Allow-Origin" => "*",
    "Access-Control-Allow-Methods" => "POST, OPTIONS"
  ]
  if HTTP.method(req) == "OPTIONS"                  # 8
    return HTTP.Response(200, headers)
  end
  body = parse(Int64, String(HTTP.body(req)))       # 4
  fact = factorial(big(body))                       # 5
  HTTP.Response(200, headers; body = string(fact))  # 6
end

function render(req::HTTP.Request)
    HTTP.Response(200, read(raw"index.html"))
end

HTTP.@register(ROUTER, "POST", "/factorial", show_factorial)  # 3
HTTP.@register(ROUTER, "GET", "/", render)        # 2

HTTP.serve(ROUTER, HOST, PORT)                    # 1


text:
function laszip_error(laszip_obj::Ptr{Cvoid})
    errstr = Ref(Cstring(C_NULL))
    laszip_get_error(laszip_obj, errstr)
    if errstr[] != C_NULL
        error(unsafe_string(errstr[]))
    end
    nothing
end

macro check(obj, ex)
    return :( $(esc(ex)) == 0 ? nothing : laszip_error($(esc(obj))) )
end


text:
using Plots, Printf
pyplot()
viz = true

@views function heat_1D()
	# physics
	lx  = 10.0
	λ   = 1.0
	ρCp = 1.0
	nt  = 200
	# numerics
	nx  = 127
	∂x  = lx/nx
	xc  = LinRange(∂x/2, lx-∂x/2, nx)
	🔥  = zeros(nx  )
	qx  = zeros(nx+1)
	🔥  = exp.(.-(xc.-lx./2.0).^2)
	∂t  = ∂x^2/ρCp/λ/2.1
	# action
	t0  = Base.time()
	for it = 1:nt
		qx[2:end-1] .= .-λ.*diff(🔥)./∂x
		🔥          .= 🔥 .- ∂t./ρCp.*diff(qx)./∂x
	end
	@printf("T_eff = %1.2e GB/s \n", (2/1e9*nx*sizeof(lx))/((Base.time()-t0)/nt))
	if viz display(plot(xc, 🔥, legend=false, xlabel="lx", ylabel="heat", title="diffusion")) end
	return
end

@time heat_1D()


text:
using Documenter
using FundamentalsNumericalComputation

Documenter.Writers.HTMLWriter.HTML(sidebar_sitename=false)
makedocs(
    sitename = "FNC Functions",
    format = Documenter.HTML(),
    modules = [FundamentalsNumericalComputation],
    pages = [
        "index.md",
        "Functions" => "functions.md",
    ]
)

deploydocs(
    devbranch = "main",
    repo = "github.com/fncbook/FundamentalsNumericalComputation.jl.git",
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#