text:
## data.jl Julia code to handle matrix-type data on disc
using Statistics
using Distributed

## report the kind of Data structure from the type instance
kind(d::Data{T,S}) where {T,S<:AbstractString} = :file
kind(d::Data{T,Matrix{T}}) where {T} = :matrix
Base.eltype(d::Data{T}) where {T} = T

## constructor for a plain matrix.  rowvectors: data points x represented as rowvectors
Data(x::Matrix{T}) where {T} = Data(Matrix{T}[x])

## constructor for a vector of files
## Data([strings], type, loadfunction)
function Data(files::Vector{S}, datatype::DataType, load::Function) where {S<:AbstractString}
    Data(files, datatype, @compat Dict(:load => load))
end

## default load function
function _load(file::AbstractString)
    load(file, "data")
end

## default size function
function _size(file::AbstractString)
    jldopen(file) do fd
        size(fd["data"])
    end
end

## courtesy compatible save for a matrix
function JLD.save(file::AbstractString, x::Matrix)
    save(file,"data", x)
end

## Data([strings], type; load=loadfunction, size=sizefunction)
function Data(files::Vector{S}, datatype::DataType; kwargs...) where {S<:AbstractString}
    all([isa((k,v), (Symbol,Function)) for (k,v) in kwargs]) || error("Wrong type of argument", args)
    d = Dict{Symbol,Function}([kwargs...])
    if !haskey(d, :load)
        d[:load] = _load
        d[:size] = _size
    end
    Data(files, datatype, d)
end

## constructor for a plain file.
Data(file::AbstractString, datatype::DataType, load::Function) = Data([file], datatype, load)
Data(file::AbstractString, datatype::DataType; kwargs...) = Data([file], datatype; kwargs...)

## is this really a shortcut?
API(d::Data, f::Symbol) = d.API[f]

function Base.getindex(x::Data, i::Int)
    if kind(x) == :matrix
        x.list[i]
    elseif kind(x) == :file
        x.API[:load](x.list[i])
    else
        error("Unknown kind")
    end
end

## A range as index (including [1:1]) returns a Data object, not the data.
function Base.getindex(x::Data{T,VT}, r::UnitRange{Int}) where {T,VT}
    Data{T, VT}(x.list[r], x.API)
end

## define an iterator for Data
Base.length(x::Data) = length(x.list)
function Base.iterate(x::Data, state=1)
    count = state

    if(count > length(x.list))
        return nothing
    else
        return (x[count], count+1)
    end
end

## This function is like pmap(), but executes each element of Data on a predestined
## worker, so that file caching at the local machine is beneficial.
## This is _not_ dynamic scheduling, like pmap().
function dmap(f::Function, x::Data)
    if kind(x) == :file
        nₓ = length(x)
        w = workers()
        nw = length(w)
        worker(i) = w[1 .+ ((i-1) % nw)]
        results = cell(nₓ)
        getnext(i) = x.list[i]
        load = x.API[:load]
        @sync begin
            for i = 1:nₓ
                @async begin
                    next = getnext(i)
                    results[i] = remotecall_fetch(worker(i), s->f(load(s)), next)
                end
            end
        end
        results
    elseif kind(x) == :matrix
        pmap(f, x)
    else
        error("Unknown kind")
    end
end

## this is like mapreduce(), but works on a data iterator
## for every worker, the data is reduced immediately, so that
## we only keep a list of (possibly large results) of the size
## of the array (and not the size of the data list)
function dmapreduce(f::Function, op::Function, x::Data)
    nₓ = length(x)
    nw = nworkers()
    results = Array{Any}(undef, nw) ## will contain pointers for parallel return value.
    valid = Any[false for i=1:nw] # can't use bitarray as parallel return value, must be pointers
    id=0
    nextid() = (id += 1)
    @sync begin
        for (wi,wid) in enumerate(workers())
            @async begin
                while true
                    i = nextid()
                    if i > nₓ
                        break
                    end
                    if kind(x) == :matrix
                        r = remotecall_fetch(f, wid, x[i])
                    else
                        r = remotecall_fetch(s->f(x.API[:load](s)), wid, x.list[i])
                    end
                    if valid[wi]
                        results[wi] = op(results[wi], r)
                    else
                        results[wi] = r
                        valid[wi] = true
                    end
                end
            end
        end
    end
    reduce(op, results[findall(valid)])
end

## stats: compute nth order stats for array (this belongs in stats.jl)
function stats(x::Matrix{T}, order::Int=2; kind=:diag, dim=1) where {T<:AbstractFloat}
    if dim==1
        n, d = size(x)
    else
        d, n = size(x)
    end
    if kind == :diag
        if order == 2
            return n, vec(sum(x, dims=dim)), vec(sum(abs2, x, dims=dim))
        elseif order == 1
            return n, vec(sum(x, dims=dim))
        else
            sx = zeros(T, order, d)
            for j=1:d
                for i=1:n
                    if dim==1
                        xi = xp = x[i,j]
                    else
                        xi = xp = x[j,i]
                    end
                    sx[1,j] += xp
                    for o=2:order
                        xp *= xi
                        sx[o,j] += xp
                    end
                end
            end
            return tuple([n, map(i->vec(sx[i,:]), 1:order)...]...)
        end
    elseif kind == :full
        order == 2 || error("Can only do covar starts for order=2")
        ## lazy implementation
        sx = vec(sum(x, dims=dim))
        sxx = x' * x
        return n, sx, sxx
    else
        error("Unknown kind")
    end
end

## Helper functions for stats tuples:
## This relies on sum(::Tuple), which sums over the elements of the tuple.
import Base: +
function +(a::Tuple, b::Tuple)
    length(a) == length(b) || error("Tuples must be of same length in addition")
    tuple(map(sum, zip(a,b))...)
end
Base.zero(t::Tuple) = map(zero, t)

## this function calls dmap as an option for parallelism
function stats(d::Data, order::Int=2; kind=:diag, dim=1)
    s = dmap(x->stats(x, order, kind=kind, dim=dim), d)
    if dim==1
        return reduce(+, s)
    else
        ## this is admittedly hairy: vertically concatenate each element of stats
        n = s[1][1]
        st = map(i->reduce(vcat, [x[i] for x in s]), 1+(1:order))
        return tuple(n, st...)
    end
end

## helper function to get summary statistics in traditional shape
function retranspose(x::Array, dim::Int)
    if dim==1
        return x'
    else
        return x
    end
end

## sum, mean, var
function Base.sum(d::Data)
    s = zero(eltype(d))
    for x in d
        s += sum(x)
    end
    return s
end

Base.sum(d::Data, dim::Int) = retranspose(stats(d,1, dim=dim)[2], dim)

function Statistics.mean(d::Data)
    n, sx = stats(d, 1)
    sum(sx) / (n*length(sx))
end

 function Statistics.mean(d::Data, dim::Int)
     n, sx = stats(d, 1, dim=dim)
     return retranspose(sx ./ n, dim)
end

function Statistics.var(d::Data)
    n, sx, sxx = stats(d, 2)
    n *= length(sx)
    ssx = sum(sx)                       # keep type stability...
    ssxx = sum(sxx)
    μ = ssx / n
    return (ssxx - n*μ^2) / (n - 1)
end

function Statistics.var(d::Data, dim::Int)
    n, sx, sxx = stats(d, 2, dim=dim)
    μ = sx ./ n
    return retranspose((sxx - n*μ.^2) ./ (n-1), dim)
end

function Statistics.cov(d::Data)
    n, sx, sxx = stats(d, 2, kind=:full)
    μ = sx ./ n
    (sxx - n*μ*μ') ./ (n-1)
end

## this is potentially very slow because it reads all file just to find out the size
function Base.size(d::Data)
    if kind(d) == :file && :size in keys(d.API)
        s = map(d.API[:size], d.list) # don't run parallel for fileserver performance
    else
        s = dmap(size, d)       # loads all data in workers, this may thrash the fileserver
    end
    nrow, ncol = s[1]
    ok = true
    for i in 2:length(s)
        ok &= s[i][2] == ncol
        nrow += s[i][1]
    end
    if !ok
        error("Inconsistent number of columns in data")
    end
    nrow, ncol
end

function Base.size(d::Data, dim::Int)
    if dim==2
        size(d[1],dim)
    else
        size(d)[dim]
    end
end

Base.collect(d::Data) = vcat([x for x in d]...)

function Base.convert(::Type{Data{Td}}, d::Data{Ts}) where {Td,Ts}
    Td == Ts && return d
    if kind(d) == :file
        api = copy(d.API)
        _load = api[:load] # local copy
        api[:load] = x -> Array{Td}(_load(x))
        return Data(d.list, Td, api)
    elseif kind(d) == :matrix
        Data(Matrix{Td}[Matrix{Td}(x) for x in d.list])
    else
        error("Unknown kind")
    end
end


text:
#___  ____       ___ _____   _   _       _ _
#|  \/  | |     |_  /  __ \ | | | |     (_| |
#| .  . | |       | | /  \/ | | | |_ __  _| |_ ___
#| |\/| | |       | | |     | | | | '_ \| | __/ _ \
#| |  | | |___/\__/ | \__/\ | |_| | | | | | || (_) |
#_______\_____\____/ \____/  _____|_|___|_____\_____ _____ _____
#| ___ \        (_)         | | \ \ / / / __  |  _  / __  |  _  |
#| |_/ _ __ ___  _  ___  ___| |_ \ V /  `' / /| |/' `' / /| |/' |
#|  __| '__/ _ \| |/ _ \/ __| __|/   \    / / |  /| | / / |  /| |
#| |  | | | (_) | |  __| (__| |_/ /^\ \ ./ /__\ |_/ ./ /__\ |_/ /
#\_|  |_|  \___/| |\___|\___|\__\/   \/ \_____/\___/\_____/\___/
#              _/ |
#             |__/
#
# This code is part of the proposal of the team "MLJC UniTo" - University of Turin
# for "ProjectX 2020" Climate Change for AI.
# The code is licensed under MIT 3.0
# Please read readme or comments for credits and further information.

# Compiler: Julia 1.5

# Short description of this file: Level set implementation for Isom Creek simulation
#with a final boundary condition

using NeuralPDE
using Quadrature, Cubature, Cuba
using Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux
using Plots
using PyPlot
using DelimitedFiles

print("Precompiling Done")

##  DECLARATIONS
@parameters t x y θ
@variables u(..)
@derivatives Dt'~t
@derivatives Dx'~x
@derivatives Dy'~y

## DOMAINS AND OPERATORS

#Real world data to model mesh scaling

Lx = 10000   #Domain x width in meters
Ly = 10000   #Domain y width in meters
Lt = 3600    #Total simulation time in seconds

xRes = 30    #Resolution of x axis in meters, mesh cell width
yRes = 30    #Resolution of y axis in meters, mesh cell width
tRes = 60    #Time resultion, seconds

# Discretization
xwidth      = 10.0      #ft
ywidth      = 10.0
tmax        = 10.0      #min
xScale      = 1.0
yScale      = 1.0
xMeshNum    = max(round(Lx/xRes),500)
yMeshNum    = max(round(Ly/yRes),500)
tMeshNum    = max(round(Lt/tRes),20)
dx  = xwidth/xMeshNum
dy  = ywidth/yMeshNum
dt  = tmax/tMeshNum

#Domain, fire position and shape

shape      = ["zeroIsVertex","zeroIsCenter"]
x0         = [5.0]      #Fire ingnition coordinates
y0         = [5.0]
xSpread    = [1.0]      #Fire shape factors
ySpread    = [0.7]
tIgnition  = [10.0]      #Fire's time of ignition
amplitude  = [2.5]      #Fire's initial spread (radius if circle)

domainShape = shape[1]

if domainShape == shape[1]
    domains = [t ∈ IntervalDomain(0.0,tmax),
           x ∈ IntervalDomain(0.0,xwidth),
           y ∈ IntervalDomain(0.0,ywidth)]
           xs = 0.0 : dx : xwidth
           ys = 0.0 : dy : ywidth
elseif domainShape == shape[2]
    domains = [t ∈ IntervalDomain(0.0,tmax),
           x ∈ IntervalDomain(-xwidth*0.5,xwidth*0.5),
           y ∈ IntervalDomain(-ywidth*0.5,ywidth*0.5)]
           xs = -xwidth*0.5: dx : xwidth*0.5
           ys = -ywidth*0.5: dy : ywidth*0.5
end

#Terrain interpolation parameters

p0 = 316.977
p1 = -31.7087
p2 = 42.2348
p3 = -12.7828
p4 = 1.26132
p5 = -0.0386497
p6 = 38.5245
p7 = 9.65925
p8 = -4.35081
p9 = 0.554218
p10 = -0.0245952
p11 = -21.5973
p12 = 1.79071
p13 = 2.49292
p14 = -0.231523

#Terrain gradient's components
Dxz = 0.01*(p1 + 2*p2*x + 3*p3*x*x + 4*p4*x*x*x + 5*p5*x*x*x*x + p11*y + 2*p12*y*x + p13*y*y + 2*p14*y*y*x)
Dyz = 0.01*(p6 + 2*p7*y + 3*p8*y*y + 4*p9*y*y*y + 5*p10*y*y*y*y + p11*x + p12*x*x + 2*p13*x*y + 2*p14*x*x*y)

#Domain discretization
x_s     = [x for x in xs][1:end-1]
y_s     = [y for y in ys][1:end-1]
hgt(x,y) = 0.01*(p0 + p1*x + p2*x*x + p3*x*x*x + p4*x*x*x*x + p5*x*x*x*x*x + p6*y + p7*y*y + p8*y*y*y + p9*y*y*y*y + p10*y*y*y*y*y +
           p11*x*y + p12*x*x*y + p13*x*y*y + p14*x*x*y*y)
z_s     = [hgt(x,y) for x in x_s for y in y_s]
z_s     = reshape(z_s,(500,500))

Plots.plot(x_s,y_s,z_s, st=:surface)    #terrain plot

Uwind = [0.0, 0.0]  #wind vector

gn   = (Dx(u(t,x,y,θ))^2 + Dy(u(t,x,y,θ))^2)^0.5 #gradient's norm
∇u   = [Dx(u(t,x,y,θ)), Dy(u(t,x,y,θ))]
∇z   = [Dxz,Dyz]
n    = ∇u/gn              #normal versor
normalized = ((Uwind[1]*n[1] + Uwind[2]*n[2])^2)^0.5 #inner product between wind and normal vector

## FUEL PARAMETERS

#Fuel parameters given by namelist.fire from WRF-Sfire
windrf = [0.36, 0.36, 0.44,  0.55,  0.42,  0.44,  0.44, 0.36, 0.36, 0.36,  0.36,  0.43,  0.46, 1e-7]
fgi    = [0.166, 0.897, 1.076, 2.468, 0.785, 1.345, 1.092, 1.121, 0.780, 2.694, 2.582, 7.749, 13.024, 1.e-7]
fueldepthm = [0.305, 0.305, 0.762, 1.829, 0.61,  0.762, 0.762, 0.061, 0.061, 0.305, 0.305, 0.701, 0.914, 0.305]
savr = [3500., 2784., 1500., 1739., 1683., 1564., 1562., 1889., 2484., 1764., 1182., 1145., 1159., 3500.]
fuelmce = [0.12, 0.15, 0.25, 0.20, 0.20, 0.25, 0.40, 0.30, 0.25, 0.25, 0.15, 0.20, 0.25, 0.12]
fueldens = [32.,32.,32.,32.,32.,32.,32.,32.,32.,32.,32.,32.,32.,32.] #from namelist.fire: "! 32 if solid, 19 if rotten"
st = [0.0555, 0.0555, 0.0555, 0.0555, 0.0555, 0.0555, 0.0555, 0.0555, 0.0555, 0.0555, 0.0555, 0.0555, 0.0555, 0.0555]
se = [0.010, 0.010, 0.010, 0.010, 0.010, 0.010, 0.010, 0.010, 0.010, 0.010, 0.010, 0.010, 0.010, 0.010]
cmbcnst  = 17.433e+06
fuelmc_g = 0.08  #or 0.09
weight = [7.,  7.,  7., 180., 100., 100., 100., 900., 900., 900., 900., 900., 900., 7.]

FuelNumber = 9
fuel_scale = 900

a     = windrf[FuelNumber]          #those numbers should be dependent on position
zf    =
z0    =
w     = weight[FuelNumber]/60               #seconds to minutes
wl    = fgi[FuelNumber]*0.204816144         #kg*m^-2 to lb*ft^-2
δm    = fueldepthm[FuelNumber]*3.28084      #m to ft
sigma = savr[FuelNumber]/3.28084            #m^-1 to ft^-1
Mx    = fuelmce[FuelNumber]
ρP    = fueldens[FuelNumber]*0.0624279606   #kg*m^-3 to lb*ft^-3
ST    = st[FuelNumber]
SE    = se[FuelNumber]
h     = cmbcnst*2.6839192e-5    #J*m^-3 to BTU*ft^-3
Mf    = fuelmc_g

## FIRE SPREAD RATE EQUATIONS

river_1 = 1 - 1/(1 + exp(-(x - r0)*1000))=#
tanϕ = sum(∇z.*n)
βop  = 3.348*sigma^(-0.8189)        #from Rothermel, eq (37)
U    = normalized                   #wind correction factor
w0   = wl/(1 + Mf)
ρb   = w0/δm                        #different from paper for units reasons
β    = ρb/ρP
ξ    = exp((0.792 + 0.618*sigma^0.5)*(β+0.1))/(192 + 0.25965*sigma)
ηs   = 0.174*(SE^(-0.19))
ηM   = 1 - 2.59*Mf/Mx + 5.11*(Mf/Mx)^2 - 3.52*(Mf/Mx)^3
wn   = w0/(1 + ST)
Γmax = (sigma^(1.5))/(495 + 0.594*sigma^(1.5))
A    = 1/(4.77*sigma^(0.1) - 7.27)
Γ    = Γmax*(β/βop)^A*exp(A*(1 - β/βop))
ϵ    = exp(-138/sigma)
Qig  = 250*β + 1116*Mf
C    = 7.47*exp(-0.133*sigma^0.55)
Ua   = a*U
E    = 0.715*exp(-0.000359*sigma)
IR   = Γ*wn*h*ηM*ηs
R0   = IR*ξ/(ρb*ϵ*Qig)              #spread rate without wind
ϕw   = C*max(Ua^β, (β/βop)^E)       #wind factor
ϕS   = 5.275*β^(-0.3)*tanϕ^2        #slope factor
S    = fuel_scale*R0*(1 + ϕw + ϕS)         #fire spread rate

eq = Dt(u(t,x,y,θ)) + S*gn ~ 0      #LEVEL SET EQUATION

initialCondition = (((xScale*(x-x0[1]))^2)*xSpread[1] + ((yScale*(y-y0[1]))^2)*ySpread[1])^0.5 - amplitude[1]   #Distance from ignition

#Multiple ignition points
#=
if length(x0) > 2
    for b = 2:length(x0)
        initialCondition = min(initialCondition, (((xScale*(x-x0[b]))^2)*xSpread[b] + ((yScale*(y-y0[b]))^2)*ySpread[b])^0.5 - amplitude[b])
    end
end
=#

bcs = [u(tIgnition[1],x,y,θ) ~ initialCondition]  #from literature

## NEURAL NETWORK
n = 16   #neuron number
maxIters = 3000     #number of iterations

chain = FastChain(FastDense(3,n,Flux.σ),FastDense(n,n,Flux.σ),FastDense(n,1))   #Neural network from Flux library

q_strategy = NeuralPDE.QuadratureTraining(algorithm =CubaCuhre(),reltol=1e-8,abstol=1e-8,maxiters=100)  #Training strategy

discretization = NeuralPDE.PhysicsInformedNN([dt,dx,dy],chain,strategy = q_strategy)

indvars = [t,x,y]   #phisically independent variables
depvars = [u]       #dependent (target) variable

dim = length(domains)

losses = []
cb = function (p,l)     #loss function handling
    println("Current loss is: $l")
    append!(losses, l)
    return false
end

pde_system = PDESystem(eq, bcs, domains, indvars, depvars)
prob = discretize(pde_system, discretization)

a_1 = time_ns()

res = GalacticOptim.solve(prob, GalacticOptim.ADAM(0.08), cb = cb, maxiters=maxIters) #allow_f_increase = false,

b_1 = time_ns()
print(string("Training time = ",(b_1-a_1)/10^9))
#initθ = res.minimizer

pars = open(readdlm,"/Julia_implementation/LevelSetEq/Isom Creek/params_level_set_Isom_Creek.txt")  #to import parameters from previous training (change also line 249 accordingly)
pars

phi = discretization.phi

##VISUALIZATION

extrapolate  = false     #allows extrapolating values outside the original domain
printBCSComp = true     #prints initial condition comparison and training loss plot

tStepFactor = 2 #Used to tune the time scale, if =tMeshNum/tmax the time step is the physical unit time
FPS = 5     #GIF frame per second

if extrapolate
    timeFactor  = 2 #used to extrapolate the prediction outside the domain
    xAxisFactor = 1.25 #IF IsZeroCenter THE RESULTING DOMAIN WILL BE (xAxisFactor * yAxisFactor times)^2 TIMES LARGER !!!
    yAxisFactor = 1.25
else
    timeFactor  = 1 #used to extrapolate the prediction outside the domain
    xAxisFactor = 1 #IF IsZeroCenter THE RESULTING DOMAIN WILL BE (xAxisFactor * yAxisFactor times)^2 TIMES LARGER !!!
    yAxisFactor = 1
end

if domainShape == shape[1]
    xs = 0.0 : dx : xwidth*xAxisFactor
    ys = 0.0 : dy : ywidth*yAxisFactor
elseif domainShape == shape[2]
    xs = -xwidth*0.5*xAxisFactor : dx : xwidth*0.5*xAxisFactor
    ys = -ywidth*0.5*yAxisFactor : dy : ywidth*0.5*yAxisFactor
end
ts = 1 : dt*tStepFactor : tmax*timeFactor

u_predict = [reshape([first(phi([t,x,y],pars)) for x in xs for y in ys], (length(xs),length(ys))) for t in ts]  #matrix of model's prediction
outfile = "u_predict.txt"
writedlm(outfile, u_predict)

maxlim = maximum(maximum(u_predict[t]) for t = 1:length(ts))
minlim = minimum(minimum(u_predict[t]) for t = 1:length(ts))

result = @animate for time = 1:length(ts)   #Animation of the level set function
    Plots.plot(xs, ys, u_predict[time],st=:surface,camera=(30,30), zlim=(minlim,maxlim), clim=(minlim,maxlim),
                title = string("ψ: max = ",round(maxlim, digits = 3)," min = ", round(minlim, digits = 3),"\\n t = ",
                round((time - 1)/tMeshNum*tStepFactor*tmax, digits = 3)))
end
gif(result, "isom_creek_surface_forensic.gif", fps = FPS)

result_level = @animate for time = 1:length(ts)     #Animation of the level set contour at z=0
    Plots.contour(xs, ys, u_predict[time::Int], levels = [0], title = string("Fireline \\n t = ",
    round((time - 1)/tMeshNum*tStepFactor*tmax, digits = 3)), legend = false, size = (600,600))
end
gif(result_level, "isom_creek_contour_forensic.gif", fps = FPS)

if printBCSComp
    zbcs(x,y) = (((xScale*(x-x0[1]))^2)*xSpread[1] + ((yScale*(y-y0[1]))^2)*ySpread[1])^0.5 - amplitude[1]

    z_s = reshape([zbcs(x,y) for x in xs for y in ys], (length(xs),length(ys)))
    target = reshape(z_s, (length(xs),length(ys)))
    diff = (u_predict[end] - target).^2
    MSE = sum(diff)/(length(xs)*length(ys))     #Mean square difference

    bcsPlot = Plots.plot(xs,ys,z_s, st=:surface,  title = "Initial Condition")    #camera=(30,30)
    bcsPredict = Plots.plot(xs, ys, u_predict[end],st=:surface, zlim=(minlim,maxlim), clim=(minlim,maxlim),
        title = string("ψ: max = ",round(maxlim, digits = 3)," min = ", round(minlim, digits = 3),"\\n t = ",0))    #initial condition prediction

    bcsDiff = Plots.plot(xs,ys,diff, st=:surface,  title = string("MSE = ", MSE))   #difference between prediction and target

    bcsFirelinePredict = Plots.contour(xs, ys, u_predict[end], levels = [0], title = string(" Fireline \\n t = ", 0)) #initial fireline prediction

    bcsFireline = Plots.contour(xs, ys, z_s, levels = [0], title = "BCS fireline ignition")     #target initial fireline

    trainingPlot = Plots.plot(1:(maxIters + 1), losses, yaxis=:log, title = string("Training time = 270 s",
        "\\n Iterations: ", maxIters, "   NN: 3>16>1"), ylabel = "log(loss)", legend = false) #loss plot

    bcsComparisonPlots = Plots.plot(bcsPlot, bcsPredict, bcsDiff, bcsFireline,bcsFirelinePredict, trainingPlot, size = (1500,600))
    Plots.savefig("isom_creek_bcs_comparison.pdf")
    bcsComparisonPlots
end

##IMPORT THE WRF OUTPUT

tensor = readdlm("/WRF_tensor_new/tensor_isom.txt")
tensor = reshape(tensor, (1440,1440,73))
tensor = permutedims(tensor, [2,1,3])

#Plots.contour(tensor[:,:,30], levels = 0.1:0.5)

# axis_scale_to_mesh_cell_ID()
xs =  0 : dx*144 : 1440
ys =  0 : dx*144 : 1440

#=
gif_tensor = @animate for time = 1:31
    Plots.plot(tensor[:,:,time], legend = false, levels = [0.01], size = (600,600))
end
gif(gif_tensor, "WRF_out_one_fire.gif", fps = FPS)
=#

##PLOTS FOR COMPARISON BETWEEN WRF AND PINNs

for i = 3:20
    timel = i
    tensor_p = Plots.contour(tensor[:,:,Int(floor((8 + timel)))], levels = [0.01], tick = true, grid = true, size = (400,400), colorbar=false, color="red", label=["WRF output"],
            xlabel = "Space domain [m]", legend=true, title=string("Isom Creek Fire at t = ", ((8+timel)*15), " min"))
    pred_p   = Plots.contour!(xs, ys, u_predict[Int(floor(timel))], levels = [0], tick = true, grid = true, size = (400,400), color="blue", label=["PINNs output"], legend=true)
    Plots.savefig(string("/Visualization/Levelset isom creek/", i ,".png"))
end

#h_tensor = heatmap(tensor_p)
#h_pred = heatmap(u_predict[1],levels = [0.6])

##SAVE PARAMETERS TO REPRODUCE THE RESULTS WITHOUT TRAINING AGAIN
param = initθ

outfile = "/Julia_implementation/LevelSetEq/Isom Creek/params_level_set_Isom_Creek.txt"
open(outfile, "w") do f
  for i in param
    println(f, i)
  end
end

ext: jl


text:
using GoldenSequences
using Base.Iterators
using LinearAlgebra
using Statistics
using Test
using Random
using Base.MathConstants
min_dist(x) = minimum(norm(x[i] .- x[j]) for (i,j) in product(eachindex(x), eachindex(x)) if i != j)

@testset "traits" begin
    @test eltype(GoldenSequence(3)) === Tuple{Float64, Float64, Float64}
    for g in [GoldenSequence(3), GoldenIntSequence(UInt, 3), GoldenCartesianSequence((13,14))]
        @test eltype(g) == eltype(collect(take(g, 3)))
    end
end

Random.seed!(1)
for d in 1:6
    @testset "dim d=$d" begin
        n = 5000
        xrand = [ntuple(d) do i rand() end for i in 1:n]
        xgolden = collect(take(GoldenSequence(d), n))
        @test all(isa.(xgolden, Ref(eltype(GoldenSequence(d)))))
        @test eltype(xgolden) == eltype(GoldenSequence(d))
        xgolden_ = reshape(reinterpret(Float64, xgolden), (d, n))
        xrand_ = reshape(reinterpret(Float64, xrand), (d, n))

        # Var(rand()) = 1/12

        @test all(xgolden_[1:d, 1] .== xgolden[1])
        @test norm(mean(xgolden_, dims=2) .- 0.5) < 1/sqrt(n)
        @test norm(cov(xgolden_') - I/12) < 1/sqrt(n)
        @test norm(mean(xgolden_, dims=2) .- 0.5) < norm(mean(xrand_, dims=2) .- 0.5)

        @test min_dist(xgolden) > min_dist(xrand)
        #@test norm(cov(xgolden_') - I/12) <  norm(cov(xrand_') - I/12)
    end
end

@testset "Numbers" begin
    @test collect(take(GoldenSequence(0.0), 10))[2] ≈ 1/golden
end
function f()
    g = GoldenSequence((0.,0.,0.))
    x, s = iterate(g)
    xsum = x[1] + x[2] + x[3]
    for i in 1:1000-1
        x, s = iterate(g, s)
        xsum += x[1] + x[2] + x[3]
    end
    xsum
end

@testset "Allocation" begin
    @test f() ≈ 1498.3407019608956
    @test  (@allocated f()) < 200
end
nothing


text:
module StarAlgebras

using SparseArrays
import LinearAlgebra

export StarAlgebra, AlgebraElement
export aug, basis, coeffs, star, supp

include("bases.jl")

include("mstructures.jl")
include("mtables.jl")

include("types.jl")
include("algebra_elts.jl")
include("arithmetic.jl")
include("show.jl")

end


text:
using IPPCore
import IPPDSP
using Base.Test

# write your own tests here

for T = IPPDSP.conv_types
    a = convert(Vector{eval(T)}, [ x = 1:10 ] )
    b = convert(Vector{eval(T)}, [ x = 11:21 ] )
    reference = conv( a, b )
    IPPResult = IPPDSP.conv( a, b )
    @test length( reference ) == length( IPPResult )
    for i = 1:length( reference )
        @test isapprox( reference[i], IPPResult[i] )
    end
end

for T = IPPDSP.xcorr_types
    a = convert(Vector{eval(T)}, [ x = 1:10 ] )
    b = convert(Vector{eval(T)}, [ x = 1:10 ] )
    reference = xcorr( a, b )
    IPPResult = IPPDSP.xcorr( a, b, lowlag = -9 )
    @test length( reference ) == length( IPPResult )
    for i = 1:length( reference )
        @test isapprox( reference[i], IPPResult[i] )
    end
end


text:
function solve(case, pitches)
    streak = 1
    dir = nothing
    last = pitches[1]
    broken = 0

    for elem in pitches[2:end]
        if elem > last
            dir == :D && (streak = 1)
            dir = :U
        elseif elem < last
            dir == :U && (streak = 1)
            dir = :D
        end

        if elem != last
            if streak >= 4
                streak = 0
                broken += 1
            end

            streak += 1
        end

        last = elem
    end

    println("Case #$(case): $(broken)")
end

function main()
    T = parse(Int, readline(stdin))
    for case in 1:T
        K = parse(Int, readline(stdin))
        pitches = parse.(Int, stdin |> readline |> split)
        solve(case, pitches)
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
module LatticeTools

include("imports.jl")

include("Basic/coord.jl")
include("Basic/hypercube.jl")  # replacement for hypercube
include("Basic/unitcell.jl")
include("Basic/lattice.jl")

include("SymmetryOperation/abstractoperation.jl")
# include("SymmetryOperation/identityoperation.jl")
include("SymmetryOperation/translationoperation.jl")
include("SymmetryOperation/pointoperation.jl")
include("SymmetryOperation/spaceoperation.jl")
# include("SymmetryOperation/directproductoperation.jl")
# include("SymmetryOperation/productoperation.jl")

include("Symmetry/abstractsymmetry.jl")
include("Symmetry/translationsymmetry.jl")
include("Symmetry/pointsymmetry.jl")
include("Symmetry/compatibility.jl")
include("Symmetry/littlesymmetry.jl")

include("SymmetryEmbedding/sitepermutation.jl")
include("SymmetryEmbedding/sitemap.jl")
include("SymmetryEmbedding/localunitary.jl")
include("SymmetryEmbedding/symmetryembedding.jl")

include("Irrep/irrep.jl")

export PointSymmetryDatabase
include("Symmetry/pointsymmetrydatabase.jl")
export IrrepDatabase
include("Irrep/irrepdatabase.jl")

include("Symmorphic/symmorphicsymmetry.jl")
include("Symmorphic/symmorphicsymmetryembedding.jl")
include("Symmorphic/symmorphicirrep.jl")

include("momentumpath.jl")

end # module


text:
# Autogenerated wrapper script for LLVM_full_jll for aarch64-apple-darwin
export clang, libclang, libclang_cpp, libllvm, liblto, llc, llvm_config, llvm_mca, mlir, opt

using Zlib_jll
JLLWrappers.@generate_wrapper_header("LLVM_full")
JLLWrappers.@declare_executable_product(clang)
JLLWrappers.@declare_library_product(libclang, "@rpath/libclang.dylib")
JLLWrappers.@declare_library_product(libclang_cpp, "@rpath/libclang-cpp.dylib")
JLLWrappers.@declare_library_product(libllvm, "@rpath/libLLVM.dylib")
JLLWrappers.@declare_library_product(liblto, "@rpath/libLTO.dylib")
JLLWrappers.@declare_executable_product(llc)
JLLWrappers.@declare_executable_product(llvm_config)
JLLWrappers.@declare_executable_product(llvm_mca)
JLLWrappers.@declare_library_product(mlir, "@rpath/libMLIR.dylib")
JLLWrappers.@declare_executable_product(opt)
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll)
    JLLWrappers.@init_executable_product(
        clang,
        "tools/clang",
    )

    JLLWrappers.@init_library_product(
        libclang,
        "lib/libclang.dylib",
        nothing,
    )

    JLLWrappers.@init_library_product(
        libclang_cpp,
        "lib/libclang-cpp.dylib",
        nothing,
    )

    JLLWrappers.@init_library_product(
        libllvm,
        "lib/libLLVM.dylib",
        nothing,
    )

    JLLWrappers.@init_library_product(
        liblto,
        "lib/libLTO.dylib",
        nothing,
    )

    JLLWrappers.@init_executable_product(
        llc,
        "tools/llc",
    )

    JLLWrappers.@init_executable_product(
        llvm_config,
        "tools/llvm-config",
    )

    JLLWrappers.@init_executable_product(
        llvm_mca,
        "tools/llvm-mca",
    )

    JLLWrappers.@init_library_product(
        mlir,
        "lib/libMLIR.dylib",
        nothing,
    )

    JLLWrappers.@init_executable_product(
        opt,
        "tools/opt",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
###################################################
## network/solvers/pathchoice.jl
## Model that incorporates path choice
###################################################

"""
    pathChoice :
    optimize travel times to minimize error from data with given paths
    Model from Zhan et al.
"""
function pathChoice(s::IterativeState; args...)
    g = s.data.network.graph
    paths = s.paths
    tripData = s.trips
    roads = s.data.network.roads

    beta1 = 0.275 / 60.
    beta2 = 1.563 / 1609.344

    MAX_SPEED = 30 * 1609.344 / 3600
    MIN_SPEED = 1 * 1609.344 / 3600

    lambda = 1e2
    DROP = 0.1
    BOOST = 10.
    MAX_LAMBDA = 1e16
    MIN_LAMBDA = 1e-16

    oldObjective = 0
    objective = 1e16

    # current times of edges
    t = s.timings.times
    theta = 1.
    # length of each path
    d = zeros(length(tripData), s.pathsPerTrip)
    for i in eachindex(tripData), j = 1:s.pathsPerTrip
        @inbounds d[i,j] = sum(paths[i][j][edge] * roads[src(edge), dst(edge)].distance
                               for edge=keys(paths[i][j]))
    end
    # time of each path
    g = zeros(length(tripData), s.pathsPerTrip)
    # probability of each path
    p = zeros(length(tripData), s.pathsPerTrip)
    # normalizing factor
    S = zeros(length(tripData))
    # residuals
    r = zeros(length(tripData))
    # put all the path/probability/objective calculations in a method
    function calculate!()
        for i in eachindex(tripData), j = 1:s.pathsPerTrip
            @inbounds g[i,j] = sum(paths[i][j][edge] * t[src(edge), dst(edge)]
                                   for edge=keys(paths[i][j]))
        end
        p = exp(theta .* (-beta1 .* g - beta2 .* d))
        # normalizing factor
        S = sum(p, 2)
        @simd for i = eachindex(tripData)
            @inbounds r[i] = tripData[i].time - dot(p[i,:], g[i,:])/S[i]
        end
    end
    # perform these calculations
    calculate!()
    # derivative of normalizing factor
    dSdt = zeros(length(tripData), length(roads) + 1)
    # derivative matrix
    J = zeros(length(tripData), length(roads) + 1)

    while abs(oldObjective - objective) > 1e-6 * oldObjective
        oldObjective = objective
        # derivative of normalizing factor
        for i in eachindex(tripData)
            for (k, (src, dst)) in enumerate(keys(roads))
                @inbounds dSdt[i,k] = -beta1 * sum(p[i,m] * get(paths[i][m], Edge(src, dst), 0)
                                         for m = 1:s.pathsPerTrip)
            end
            dSdt[i,end] = - dot(p[i,:], beta1 * g[i,:] + beta2 * d[i,:])
        end
        # derivative matrix
        for i = eachindex(tripData)
            for (k, (src, dst)) in enumerate(keys(roads))
                @inbounds J[i,k] = 1/S[i] ^ 2 * sum(
                                          [(p[i,m] * S[i] * get(paths[i][m], Edge(src, dst), 0)
                                            * (1 - beta1 * theta * g[i,m])) for m=1:s.pathsPerTrip]
                                          - g[i,:] .* p[i,:] * dSdt[i,k]
                                          )
            end
            J[i,end] = - dot(p[i,:] .* g[i,:]/S[i],
                             beta1 * g[i,:] + beta2 * d[i,:] + dSdt[i,end]/S[i])
        end
        diagTerms = vec(sumabs2(J,1))
        @simd for i in eachindex(diagTerms)
            @inbounds if diagTerms[i] < 1e-6
                diagTerms[i] = 1e-6
            end
        end
        println("*** Solving LM ***")
        while true
            # solve Levenberg-Marquardt system
            updateDir = (J'*J + lambda * diagm(diagTerms)) \ (J'*r)
            # save old parameters
            oldTheta = theta
            oldTimes = copy(t)
            # update parameters
            for (k, (src, dst)) in enumerate(keys(roads))
                t[src, dst] += updateDir[k]
                if t[src, dst] < roads[src, dst].distance/MAX_SPEED
                    t[src, dst] = roads[src, dst].distance/MAX_SPEED
                elseif t[src,dst] > roads[src, dst].distance/MIN_SPEED
                    t[src, dst] = roads[src, dst].distance/MIN_SPEED
                end
            end
            theta += updateDir[end]
            # see what things looklike now
            calculate!()
            objective = sumabs2(r)
            # check if we made progress
            println("λ = ", lambda)
            println("Old objective: ", oldObjective)
            println("New objective: ", objective)
            if objective <= oldObjective
                if lambda < MAX_LAMBDA && abs(objective - oldObjective) < 1e-3 * oldObjective
                    println("Reselecting lambda")
                    lambda = min(lambda * BOOST, MAX_LAMBDA) # need steeper descent
                    # reset parameters
                    t = oldTimes
                    theta = oldTheta
                    calculate!()
                else
                    if abs(objective - oldObjective) > 1e-2 * oldObjective
                        # can decrease lambda: fast convergence
                        lambda = max(lambda * DROP, MIN_LAMBDA)
                    end
                    break
                end
            else
                println("Reselecting lambda")
                if lambda < MAX_LAMBDA
                    lambda = lambda * BOOST # need steeper descent
                else
                    break
                end
                # reset parameters
                t = oldTimes
                theta = oldTheta
                calculate!()
            end
        end
    end
    println("θ = $(theta)")
    return t
end


text:
function gausshermite( n::Integer )
    x,w = unweightedgausshermite(n)
    w .*= exp.(-x.^2)
    x, w
end
function unweightedgausshermite( n::Integer )
    # GAUSSHERMITE(n) COMPUTE THE GAUSS-HERMITE NODES AND WEIGHTS IN O(n) time.
    if n < 0
        x = (Float64[],Float64[])
        return x
    elseif n == 0
        x = (Float64[],Float64[])
        return x
    elseif n == 1
        x = ([0.0],[sqrt(pi)])
        return x
    elseif n <= 20
       # GW algorithm
       x = hermpts_gw( n )
    elseif n <= 200
       # REC algorithm
       x = hermpts_rec( n )
    else
       # ASY algorithm
       x = hermpts_asy( n )
    end

    if mod(n,2) == 1                              # fold out
        w = [flipdim(x[2][:],1); x[2][2:end]]
        x = [-flipdim(x[1],1) ; x[1][2:end]]
    else
        w = [flipdim(x[2][:],1); x[2][:]]
        x = [-flipdim(x[1],1) ; x[1]]
    end
    w .*= sqrt(π)/sum(exp.(-x.^2).*w)
    (x, w)
end

function hermpts_asy( n::Integer )
    # Compute Hermite nodes and weights using asymptotic formula

    x0 = HermiteInitialGuesses( n ) # get initial guesses
    t0 = x0./sqrt(2n+1)
    theta0 = acos.(t0)               # convert to theta-variable
    val = x0;
    for k = 1:20
        val = hermpoly_asy_airy(n, theta0);
        dt = -val[1]./(sqrt(2).*sqrt(2n+1).*val[2].*sin.(theta0))
        theta0 .-= dt;                        # Newton update
        if norm(dt,Inf) < sqrt(eps(Float64))/10
           break
        end
    end
    t0 = cos.(theta0)
    x = sqrt(2n+1)*t0                          #back to x-variable
    w = x.*val[1] .+ sqrt(2).*val[2]
    w .= 1 ./ w.^2;            # quadrature weights

    (x, w)
end

function hermpts_rec( n::Integer )
    # Compute Hermite nodes and weights using recurrence relation.

    x0 = HermiteInitialGuesses( n )
    x0 .*= sqrt(2)
    val = x0
    for _ = 1:10
        val = hermpoly_rec.(n, x0)
        dx = first.(val)./last.(val)
        dx[ isnan.( dx ) ] .= 0
        x0 .= x0 .- dx
        if norm(dx, Inf)<sqrt(eps(Float64))
            break
        end
    end
    x0 ./= sqrt(2)
    w = 1 ./ last.(val).^2           # quadrature weights

    x = (x0, w)
end

function hermpoly_rec( n::Integer, x0)
    # HERMPOLY_rec evaluation of scaled Hermite poly using recurrence
    n < 0 && throw(ArgumentError("n = $n must be positive"))
    # evaluate:
    w = exp(-x0^2 / (4*n))
    wc = 0 # 0 times we've applied wc
    Hold = one(x0)
    # n == 0 && return (Hold, 0)
    H = x0
    for k = 1:n-1
        Hold, H = H, (x0*H/sqrt(k+1) - Hold/sqrt(1+1/k))
        while abs(H) ≥ 100 && wc < n # regularise
            H *= w
            Hold *= w
            wc += 1
        end
        k += 1
    end
    for _ = wc+1:n
        H *= w
        Hold *= w
    end

    # return (value, derivative):
    val = (H, -x0*H + sqrt(n)*Hold)
end

function hermpoly_rec( r::Base.OneTo, x0)
    isempty(r) && return [1.0]
    n = maximum(r)
    # HERMPOLY_rec evaluation of scaled Hermite poly using recurrence
    n < 0 && throw(ArgumentError("n = $n must be positive"))
    n == 0 && return [exp(-x0^2 / 4)]
    p = max(1,floor(Int,x0^2/100))
    w = exp(-x0^2 / (4*p))
    wc = 0 # 0 times we've applied wc
    ret = Vector{Float64}()
    Hold = one(x0)
    push!(ret, Hold)
    H = x0
    push!(ret, H)
    for k = 1:n-1
        Hold, H = H, (x0*H/sqrt(k+1) - Hold/sqrt(1+1/k))
        while abs(H) ≥ 100 && wc < p # regularise
            ret .*= w
            H *= w
            Hold *= w
            wc += 1
        end
        push!(ret, H)
        k += 1
    end
    ret .*= w^(p-wc)

    ret
end

hermpoly_rec( r::AbstractRange, x0) = hermpoly_rec(Base.OneTo(maximum(r)), x0)[r.+1]

function hermpoly_asy_airy(n::Integer, theta)
    # HERMPOLY_ASY evaluation hermite poly using Airy asymptotic formula in
    # theta-space.

    musq = 2n+1;
    cosT = cos.(theta)
    sinT = sin.(theta)
    sin2T = 2 .* cosT.*sinT
    eta = 0.5 .* theta .- 0.25 .* sin2T
    chi = -(3*eta/2).^(2/3)
    phi = (-chi./sinT.^2).^(1/4)
    C = 2*sqrt(pi)*musq^(1/6)*phi
    Airy0 = real.(airyai.(musq.^(2/3).*chi))
    Airy1 = real.(airyaiprime.(musq.^(2/3).*chi))

    # Terms in (12.10.43):
    a0 = 1; b0 = 1
    a1 = 15/144; b1 = -7/5*a1
    a2 = 5*7*9*11/2/144^2; b2 = -13/11*a2
    a3 = 7*9*11*13*15*17/6/144^3
    b3 = -19/17*a3

    # u polynomials in (12.10.9)
    u0 = 1; u1 = (cosT.^3-6*cosT)/24
    u2 = @. (-9*cosT^4 + 249*cosT^2 + 145)/1152
    u3 = @. (-4042*cosT^9+18189*cosT^7-28287*cosT^5-151995*cosT^3-259290*cosT)/414720

    #first term
    A0 = 1
    val = A0*Airy0

    #second term
    B0 = @. -(a0*phi^6 * u1+a1*u0)/chi^2
    val .+=  B0.*Airy1./musq.^(4/3)

    # third term
    A1 = @. (b0*phi^12 * u2 + b1*phi^6 * u1 + b2*u0)/chi^3
    val .+= A1.*Airy0/musq.^2

    # fourth term
    B1 = @. -(phi^18 * u3 + a1*phi^12 * u2 + a2*phi^6 * u1 + a3*u0)/chi^5
    val .+= B1.*Airy1./musq.^(4/3+2)

    val .= C.*val

    ## Derivative

    eta = .5*theta - .25*sin2T
    chi = -(3*eta/2).^(2/3)
    phi = (-chi./sinT.^2).^(1/4)
    C = sqrt(2*pi)*musq^(1/3)./phi

    # v polynomials in (12.10.10)
    v0 = 1;
    v1 = @. (cosT^3+6*cosT)/24
    v2 = @. (15*cosT^4-327*cosT^2-143)/1152
    v3 = @. (259290*cosT + 238425*cosT^3 - 36387*cosT^5 + 18189*cosT^7 - 4042*cosT^9)/414720

    # first term
    C0 = -(b0*phi.^6 .* v1 .+ b1.*v0)./chi
    dval = C0.*Airy0/musq.^(2/3)

    # second term
    D0 =  a0*v0
    dval = dval + D0*Airy1

    # third term
    C1 = @. -(phi^18 * v3 + b1*phi^12 * v2 + b2*phi^6 * v1 + b3*v0)/chi^4
    dval = dval + C1.*Airy0/musq.^(2/3+2)

    #fourth term
    D1 = @. (a0*phi^12 * v2 + a1*phi^6 * v1 + a2*v0)/chi^3
    dval = dval + D1.*Airy1/musq.^2

    dval = C.*dval

    val = (val, dval)
end

let T(t) = @. t^(2/3)*(1+5/48*t^(-2)-5/36*t^(-4)+(77125/82944)*t^(-6) -108056875/6967296*t^(-8)+162375596875/334430208*t^(-10))
    global function HermiteInitialGuesses( n::Integer )
        #HERMITEINTITIALGUESSES(N), Initial guesses for Hermite zeros.
        #
        # [1] L. Gatteschi, Asymptotics and bounds for the zeros of Laguerre
        # polynomials: a survey, J. Comput. Appl. Math., 144 (2002), pp. 7-27.
        #
        # [2] F. G. Tricomi, Sugli zeri delle funzioni di cui si conosce una
        # rappresentazione asintotica, Ann. Mat. Pura Appl. 26 (1947), pp. 283-300.

        # Error if n < 20 because initial guesses are based on asymptotic expansions:
        @assert n>=20

        # Gatteschi formula involving airy roots [1].
        # These initial guess are good near x = sqrt(n+1/2);
        if mod(n,2) == 1
            m = (n-1)>>1
            bess = (1:m)*pi
            a = .5
        else
            m = n>>1
            bess = ((0:m-1) .+ 0.5)*pi
            a = -.5
        end
        nu = 4*m + 2*a + 2

        airyrts = -T(3/8*pi*(4*(1:m) .- 1))

        airyrts_exact = [-2.338107410459762           # Exact Airy roots.
            -4.087949444130970
            -5.520559828095555
            -6.786708090071765
            -7.944133587120863
            -9.022650853340979
            -10.040174341558084
            -11.008524303733260
            -11.936015563236262
            -12.828776752865757]
        airyrts[1:10] = airyrts_exact  # correct first 10.

        x_init = sqrt.(abs.(nu .+ (2^(2/3)).*airyrts.*nu^(1/3) .+ (1/5*2^(4/3)).*airyrts.^2 .* nu^(-1/3) .+
            (11/35-a^2-12/175).*airyrts.^3 ./ nu .+ ((16/1575).*airyrts.+(92/7875).*airyrts.^4).*2^(2/3).*nu^(-5/3) .-
            ((15152/3031875).*airyrts.^5 .+ (1088/121275).*airyrts.^2).*2^(1/3).*nu^(-7/3)))
        x_init_airy = real( flipdim(x_init,1) )

        # Tricomi initial guesses. Equation (2.1) in [1]. Originally in [2].
        # These initial guesses are good near x = 0 . Note: zeros of besselj(+/-.5,x)
        # are integer and half-integer multiples of pi.
        # x_init_bess =  bess/sqrt(nu).*sqrt((1+ (bess.^2+2*(a^2-1))/3/nu^2) );
        Tnk0 = fill(pi/2,m)
        nu = (4*m+2*a+2)
        rhs = ((4*m+3) .- 4*(1:m))/nu*pi

        for k = 1:7
            val = Tnk0 .- sin.(Tnk0) .- rhs
            dval = 1 .- cos.(Tnk0)
            dTnk0 = val./dval
            Tnk0 = Tnk0 .- dTnk0
        end

        tnk = cos.(Tnk0./2).^2
        x_init_sin = @. sqrt(nu*tnk - (5 / (4 * (1-tnk)^2) - 1 / (1 - tnk)-1 + 3*a^2)/3 / nu)

        # Patch together
        p = 0.4985+eps(Float64)
        x_init = [x_init_sin[1:convert(Int,floor(p*n))] ;
        x_init_airy[convert(Int,ceil(p*n)):end]]

        if mod(n, 2) == 1
            x_init = [0 ; x_init]
            x_init = x_init[1:m+1]
        else
            x_init = x_init[1:m]
        end

        return x_init
    end
end

function hermpts_gw( n::Integer )
    # Golub--Welsch algorithm. Used here for n<=20.

    beta = sqrt.(0.5 .* (1:n-1))              # 3-term recurrence coeffs
    T = SymTridiagonal(zeros(n), beta)  # Jacobi matrix
    (D, V) = eigen(T)                      # Eigenvalue decomposition
    indx = sortperm(D)                  # Hermite points
    x = D[indx]
    w = sqrt(pi)*V[1,indx].^2            # weights

    # Enforce symmetry:
    ii = floor(Int, n/2)+1:n
    x = x[ii]
    w = w[ii]
    return (x,exp.(x.^2).*w)
end


text:
### parsing.jl
#
# Copyright (C) 2016, 2017 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: uncertainty, error propagation, physics
#
# This file is a part of Measurements.jl.
#
# License is MIT "Expat".
#
### Commentary:
#
# This file defines a method for `measurement` which enables parsing of a string
# into a `Measurement` object.
#
### Code:

# Regxep matching the exponent part: "e-4"
const rxp_expn = "e[+-]?[0-9]+"
# Regexp matching the complete number: "5.6e-7"  or "48"
const rxp_full_numb = "[0-9.]+" * "(?:" * rxp_expn * ")?"
# Regexp matching the complete number, including the sign: "-12.3e4"
const rxp_full_numb_sign = "[+-]?" * rxp_full_numb

# Regexp matching a measurement, with error expressed within parentheses:
# "123(45)e6".  Captures:
#  1: the whole nominal value
#  2: the decimal part of the nominal value (optional)
#  3: the uncertainty on the last digits
#  4: a global exponent factor (optional)
const rxp_error_with_parentheses =
    Regex("^[ \\t]*([+-]?[0-9]+(\\.[0-9]+)?)\\(([0-9]+)\\)(" *
          rxp_expn * ")?[ \\t]*\$")

# Regexp matching a global exponent: "(...)e-1".  Captures:
#  1: the rest of the string
#  2: the global exponent
const rxp_global_exponent = Regex("^[ \\t]*\\((.*)\\)(" *
                                  rxp_expn * ")[ \\t]*\$")

# Regexp matching a measurement, with error expressed within plus-minus sign:
# "12.3e-4 ± 5.6e-7".  Captures:
#  1: the nominal value
#  2: the uncertainty
#  3: global exponent factor (optional)
const rxp_error_with_pm =
    Regex("^[ \\t]*(" * rxp_full_numb_sign * ")[ \\t]*(?:\\+/?-|±)[ \\t]*(" *
          rxp_full_numb * ")[ \\t]*\$")

# Regexp matching number without uncertainty "123.4e5"
const rxp_no_error = Regex("^[ \\t]*(" * rxp_full_numb_sign * ")[ \\t]*\$")

"""
    measurement(string) -> Measurement{Float64}

Parse the string and return a `Measurement{Float64}` object.

Examples of valid strings and the resulting `Measurement{Float64}` are:

```jldoctest
julia> measurement("-123.4(56)")
-123.4 ± 5.6

julia> measurement("+1234(56)e-1")
123.4 ± 5.6000000000000005

julia> measurement("12.34e-1 +- 0.56e1")
1.234 ± 5.6

julia> measurement("(-1.234 ± 0.056)e2")
-123.4 ± 5.6000000000000005

julia> measurement("1234e-1 +/- 5.6e0")
123.4 ± 5.6

julia> measurement("-1234e-1")
-123.4 ± 0.0
```
"""
measurement(str::AbstractString) = parse(Measurement{Float64}, str)

function Base.parse(::Type{Measurement{T}}, str::S) where {T<:AbstractFloat, S<:AbstractString}
    m = match(rxp_error_with_parentheses, str)
    if m !== nothing # "123(45)e6"
        val_str::S, val_dec, err_str::S, expn = m.captures
    else
        m = match(rxp_global_exponent, str)
        if m === nothing # There is no global exponent
            expn = nothing
        else # "(...)e5"
            str, expn = m.captures
        end
        m = match(rxp_error_with_pm, str)
        if m !== nothing # "12.3e-4 ± 5.6e-7"
            val_str, err_str, val_dec = m.captures..., nothing
        else
            m = match(rxp_no_error, str) # Match "123.4e5"
            if m !== nothing === expn # Exclude the case "(123.4e5)e6"
                # measurement(val) returns val ± 0.  For consistency,
                # measurement("val") should give the same result.
                val_str, err_str, val_dec, expn =
                    m.captures[1], "0", nothing, nothing
            else
                throw(ArgumentError("cannot parse $(repr(str)) as Measurement{$T}"))
            end
        end
    end
    val::T = parse(T, val_str)
    err::T = parse(T, err_str)
    if val_dec !== nothing # The nominal value has a decimal part
        # Divide the uncertainty by 10^(number_of_decimal_digits)
        err /= 10 ^ (length(val_dec) - 1)
    end
    if expn !== nothing # There is a global exponent factor
        # Parse to a number
        fact = parse(T, "1" * expn)
        val *= fact
        err *= fact
    end
    return measurement(val, err)
end


text:
"""
Calculates the From - To complex power flow (Flow injected at the bus) of branch of type
TapTransformer
"""
function flow_val(b::TapTransformer)
    !get_available(b) && return 0.0
    Y_t = get_series_admittance(b)
    c = 1 / get_tap(b)
    arc = get_arc(b)
    V_from = arc.from.magnitude * (cos(arc.from.angle) + sin(arc.from.angle) * 1im)
    V_to = arc.to.magnitude * (cos(arc.to.angle) + sin(arc.to.angle) * 1im)
    I = (V_from * Y_t * c^2) - (V_to * Y_t * c)
    flow = V_from * conj(I)
    return flow
end

"""
Calculates the From - To complex power flow (Flow injected at the bus) of branch of type
Line
"""
function flow_val(b::ACBranch)
    !get_available(b) && return 0.0
    Y_t = get_series_admittance(b)
    arc = PowerSystems.get_arc(b)
    V_from = arc.from.magnitude * (cos(arc.from.angle) + sin(arc.from.angle) * 1im)
    V_to = arc.to.magnitude * (cos(arc.to.angle) + sin(arc.to.angle) * 1im)
    I = V_from * (Y_t + (1im * PowerSystems.get_b(b).from)) - V_to * Y_t
    flow = V_from * conj(I)
    return flow
end

"""
Calculates the From - To complex power flow (Flow injected at the bus) of branch of type
Transformer2W
"""
function flow_val(b::Transformer2W)
    !get_available(b) && return 0.0
    Y_t = get_series_admittance(b)
    arc = get_arc(b)
    V_from = arc.from.magnitude * (cos(arc.from.angle) + sin(arc.from.angle) * 1im)
    V_to = arc.to.magnitude * (cos(arc.to.angle) + sin(arc.to.angle) * 1im)
    I = V_from * (Y_t + (1im * get_primary_shunt(b))) - V_to * Y_t
    flow = V_from * conj(I)
    return flow
end

function flow_val(b::PhaseShiftingTransformer)
    error("Systems with PhaseShiftingTransformer not supported yet")
    return
end

"""
Calculates the From - To complex power flow using external data of voltages of branch of type
TapTransformer
"""
function flow_func(b::TapTransformer, V_from::Complex{Float64}, V_to::Complex{Float64})
    !get_available(b) && return (0.0, 0.0)
    Y_t = get_series_admittance(b)
    c = 1 / get_tap(b)
    I = (V_from * Y_t * c^2) - (V_to * Y_t * c)
    flow = V_from * conj(I)
    return real(flow), imag(flow)
end

"""
Calculates the From - To complex power flow using external data of voltages of branch of type
Line
"""
function flow_func(b::ACBranch, V_from::Complex{Float64}, V_to::Complex{Float64})
    !get_available(b) && return (0.0, 0.0)
    Y_t = get_series_admittance(b)
    I = V_from * (Y_t + (1im * get_b(b).from)) - V_to * Y_t
    flow = V_from * conj(I)
    return real(flow), imag(flow)
end

"""
Calculates the From - To complex power flow using external data of voltages of branch of type
Transformer2W
"""
function flow_func(b::Transformer2W, V_from::Complex{Float64}, V_to::Complex{Float64})
    !get_available(b) && return (0.0, 0.0)
    Y_t = get_series_admittance(b)
    I = V_from * (Y_t + (1im * get_primary_shunt(b))) - V_to * Y_t
    flow = V_from * conj(I)
    return real(flow), imag(flow)
end

function flow_func(
    b::PhaseShiftingTransformer,
    V_from::Complex{Float64},
    V_to::Complex{Float64},
)
    error("Systems with PhaseShiftingTransformer not supported yet")
    return
end

"""
Updates the flow on the branches
"""
function _update_branch_flow!(sys::System)
    for b in get_components(ACBranch, sys)
        S_flow = get_available(b) ? flow_val(b) : 0.0 + 0.0im
        set_active_power_flow!(b, real(S_flow))
        set_reactive_power_flow!(b, imag(S_flow))
    end
end

"""
Obtain total load on bus b
"""
function _get_load_data(sys::System, b::Bus)
    active_power = 0.0
    reactive_power = 0.0
    for l in get_components(ElectricLoad, sys)
        !get_available(l) && continue
        if (l.bus == b)
            if isa(l, FixedAdmittance)
                # Assume v rated = 1.0
                active_power += real(get_Y(l))
                reactive_power += imag(get_Y(l))
            else
                active_power += get_active_power(l)
                reactive_power += get_reactive_power(l)
            end
        end
    end
    return active_power, reactive_power
end

"""
Updates system voltages and powers with power flow results
"""
function _write_pf_sol!(sys::System, nl_result)
    function _is_available_source(x, bus::Bus)
        return get_available(x) && x.bus == bus && !isa(x, ElectricLoad)
    end

    result = round.(nl_result.zero; digits = 7)
    buses = enumerate(sort!(collect(get_components(Bus, sys)), by = x -> get_number(x)))

    for (ix, bus) in buses
        if bus.bustype == BusTypes.REF
            P_gen = result[2 * ix - 1]
            Q_gen = result[2 * ix]
            devices =
                get_components(StaticInjection, sys, x -> _is_available_source(x, bus))
            sum_basepower = sum(get_base_power.(devices))
            for d in devices
                part_factor = get_base_power(d) / sum_basepower
                set_active_power!(d, P_gen * part_factor)
                set_reactive_power!(d, Q_gen * part_factor)
            end
        elseif bus.bustype == BusTypes.PV
            Q_gen = result[2 * ix - 1]
            θ = result[2 * ix]
            devices =
                get_components(StaticInjection, sys, x -> _is_available_source(x, bus))
            sum_basepower = sum(get_base_power.(devices))
            for d in devices
                set_reactive_power!(d, Q_gen * get_base_power(d) / sum_basepower)
            end
            bus.angle = θ
        elseif bus.bustype == BusTypes.PQ
            Vm = result[2 * ix - 1]
            θ = result[2 * ix]
            set_magnitude!(bus, Vm)
            set_angle!(bus, θ)
        end
    end

    _update_branch_flow!(sys)
    return
end

"""
Return power flow results in dictionary of dataframes.
"""
function _write_results(sys::System, nl_result)
    @info("Voltages are exported in pu. Powers are exported in MW/MVAr.")
    result = round.(nl_result.zero; digits = 7)
    buses = sort!(collect(get_components(Bus, sys)), by = x -> get_number(x))
    N_BUS = length(buses)
    bus_map = Dict(buses .=> 1:N_BUS)
    sys_basepower = get_base_power(sys)
    sources = get_components(StaticInjection, sys, d -> !isa(d, ElectricLoad))
    Vm_vect = fill(0.0, N_BUS)
    θ_vect = fill(0.0, N_BUS)
    P_gen_vect = fill(0.0, N_BUS)
    Q_gen_vect = fill(0.0, N_BUS)
    P_load_vect = fill(0.0, N_BUS)
    Q_load_vect = fill(0.0, N_BUS)

    for (ix, bus) in enumerate(buses)
        P_load_vect[ix], Q_load_vect[ix] = _get_load_data(sys, bus) .* sys_basepower
        if bus.bustype == BusTypes.REF
            Vm_vect[ix] = get_magnitude(bus)
            θ_vect[ix] = get_angle(bus)
            P_gen_vect[ix] = result[2 * ix - 1] * sys_basepower
            Q_gen_vect[ix] = result[2 * ix] * sys_basepower
        elseif bus.bustype == BusTypes.PV
            Vm_vect[ix] = get_magnitude(bus)
            θ_vect[ix] = result[2 * ix]
            for gen in sources
                !get_available(gen) && continue
                if gen.bus == bus
                    P_gen_vect[ix] += get_active_power(gen) * sys_basepower
                end
            end
            Q_gen_vect[ix] = result[2 * ix - 1] * sys_basepower
        elseif bus.bustype == BusTypes.PQ
            Vm_vect[ix] = result[2 * ix - 1]
            θ_vect[ix] = result[2 * ix]
            for gen in sources
                !get_available(gen) && continue
                if gen.bus == bus
                    P_gen_vect[ix] += get_active_power(gen) * sys_basepower
                    Q_gen_vect[ix] += get_reactive_power(gen) * sys_basepower
                end
            end
        end
    end

    branches = get_components(ACBranch, sys)
    N_BRANCH = length(branches)
    P_from_to_vect = fill(0.0, N_BRANCH)
    Q_from_to_vect = fill(0.0, N_BRANCH)
    P_to_from_vect = fill(0.0, N_BRANCH)
    Q_to_from_vect = fill(0.0, N_BRANCH)
    for (ix, b) in enumerate(branches)
        !get_available(b) && continue
        bus_f_ix = bus_map[get_arc(b).from]
        bus_t_ix = bus_map[get_arc(b).to]
        V_from = Vm_vect[bus_f_ix] * (cos(θ_vect[bus_f_ix]) + sin(θ_vect[bus_f_ix]) * 1im)
        V_to = Vm_vect[bus_t_ix] * (cos(θ_vect[bus_t_ix]) + sin(θ_vect[bus_t_ix]) * 1im)
        P_from_to_vect[ix], Q_from_to_vect[ix] = flow_func(b, V_from, V_to) .* sys_basepower
        P_to_from_vect[ix], Q_to_from_vect[ix] = flow_func(b, V_to, V_from) .* sys_basepower
    end

    bus_df = DataFrames.DataFrame(
        bus_number = get_number.(buses),
        Vm = Vm_vect,
        θ = θ_vect,
        P_gen = P_gen_vect,
        P_load = P_load_vect,
        P_net = P_gen_vect - P_load_vect,
        Q_gen = Q_gen_vect,
        Q_load = Q_load_vect,
        Q_net = Q_gen_vect - Q_load_vect,
    )

    branch_df = DataFrames.DataFrame(
        line_name = get_name.(branches),
        bus_from = get_number.(get_from.(get_arc.(branches))),
        bus_to = get_number.(get_to.(get_arc.(branches))),
        P_from_to = P_from_to_vect,
        Q_from_to = Q_from_to_vect,
        P_to_from = P_to_from_vect,
        Q_to_from = Q_to_from_vect,
        P_losses = P_from_to_vect + P_to_from_vect,
        Q_losses = Q_from_to_vect + Q_to_from_vect,
    )
    DataFrames.sort!(branch_df, [:bus_from, :bus_to])

    return Dict("bus_results" => bus_df, "flow_results" => branch_df)
end

"""
Solves a the power flow into the system and writes the solution into the relevant structs.
Updates generators active and reactive power setpoints and branches active and reactive
power flows (calculated in the From - To direction) (see
[`flow_val`](@ref))

Supports solving using Finite Differences Method (instead of using analytic Jacobian)
by setting finite_diff = true.
Supports passing NLsolve kwargs in the args. By default shows the solver trace.

Arguments available for `nlsolve`:

* `method` : See NLSolve.jl documentation for available solvers
* `xtol`: norm difference in `x` between two successive iterates under which
  convergence is declared. Default: `0.0`.
* `ftol`: infinite norm of residuals under which convergence is declared.
  Default: `1e-8`.
* `iterations`: maximum number of iterations. Default: `1_000`.
* `store_trace`: should a trace of the optimization algorithm's state be
  stored? Default: `false`.
* `show_trace`: should a trace of the optimization algorithm's state be shown
  on `STDOUT`? Default: `false`.
* `extended_trace`: should additifonal algorithm internals be added to the state
  trace? Default: `false`.

## Examples
```julia
solve_powerflow!(sys)
# Passing NLsolve arguments
solve_powerflow!(sys, method = :newton)
# Using Finite Differences
solve_powerflow!(sys, finite_diff = true)

```

"""
function solve_powerflow!(system::System; finite_diff = false, kwargs...)
    #Save per-unit flag
    settings_unit_cache = deepcopy(system.units_settings.unit_system)
    #Work in System per unit
    set_units_base_system!(system, "SYSTEM_BASE")
    res = _solve_powerflow(system, finite_diff; kwargs...)
    if res.f_converged
        _write_pf_sol!(system, res)
        @info("PowerFlow solve converged, the results have been stored in the system")
        #Restore original per unit base
        set_units_base_system!(system, settings_unit_cache)
        return res.f_converged
    end
    @error("The powerflow solver returned convergence = $(res.f_converged)")
    set_units_base_system!(system, settings_unit_cache)
    return res.f_converged
end

"""
Similar to solve_powerflow!(sys) but does not update the system struct with results.
Returns the results in a dictionary of dataframes.

## Examples
```julia
res = solve_powerflow(sys)
# Passing NLsolve arguments
res = solve_powerflow(sys, method = :newton)
# Using Finite Differences
res = solve_powerflow(sys, finite_diff = true)

```

"""
function solve_powerflow(system::System; finite_diff = false, kwargs...)
    #Save per-unit flag
    settings_unit_cache = deepcopy(system.units_settings.unit_system)
    #Work in System per unit
    set_units_base_system!(system, "SYSTEM_BASE")
    res = _solve_powerflow(system, finite_diff; kwargs...)
    if res.f_converged
        @info("PowerFlow solve converged, the results are exported in DataFrames")
        df_results = _write_results(system, res)
        #Restore original per unit base
        set_units_base_system!(system, settings_unit_cache)
        return df_results
    end
    @error("The powerflow solver returned convergence = $(res.f_converged)")
    set_units_base_system!(system, settings_unit_cache)
    return res.f_converged
end

function _solve_powerflow(system::System, finite_diff::Bool; kwargs...)
    ybus_kw = [:check_connectivity, :connectivity_method]
    ybus_kwargs = (k for k in kwargs if first(k) ∈ ybus_kw)
    kwargs = (k for k in kwargs if first(k) ∉ ybus_kw)

    buses = sort!(collect(get_components(Bus, system)), by = x -> get_number(x))
    N_BUS = length(buses)

    # assumes the ordering in YBus is the same as in the buses.
    Yb = Ybus(system; ybus_kwargs...).data
    a = collect(1:N_BUS)
    I, J, V = SparseArrays.findnz(Yb)
    neighbors = [Set{Int}([i]) for i in 1:N_BUS]
    for nz in eachindex(V)
        push!(neighbors[I[nz]], J[nz])
        push!(neighbors[J[nz]], I[nz])
    end
    x0 = zeros(N_BUS * 2)

    #Create Jacobian structure
    J0_I = Int[]
    J0_J = Int[]
    J0_V = Float64[]

    for ix_f in a
        F_ix_f_r = 2 * ix_f - 1
        F_ix_f_i = 2 * ix_f

        for ix_t in neighbors[ix_f]
            X_ix_t_fst = 2 * ix_t - 1
            X_ix_t_snd = 2 * ix_t
            b = get_bustype(buses[ix_t])
            #Set to 0.0 only on connected buses
            if b == BusTypes.REF
                if ix_f == ix_t
                    #Active PF w/r Local Active Power
                    push!(J0_I, F_ix_f_r)
                    push!(J0_J, X_ix_t_fst)
                    push!(J0_V, 0.0)
                    #Rective PF w/r Local Reactive Power
                    push!(J0_I, F_ix_f_i)
                    push!(J0_J, X_ix_t_snd)
                    push!(J0_V, 0.0)
                end
            elseif b == BusTypes.PV
                #Active PF w/r Angle
                push!(J0_I, F_ix_f_r)
                push!(J0_J, X_ix_t_snd)
                push!(J0_V, 0.0)
                #Reactive PF w/r Angle
                push!(J0_I, F_ix_f_i)
                push!(J0_J, X_ix_t_snd)
                push!(J0_V, 0.0)
                if ix_f == ix_t
                    #Reactive PF w/r Local Reactive Power
                    push!(J0_I, F_ix_f_i)
                    push!(J0_J, X_ix_t_fst)
                    push!(J0_V, 0.0)
                end
            elseif b == BusTypes.PQ
                #Active PF w/r VoltageMag
                push!(J0_I, F_ix_f_r)
                push!(J0_J, X_ix_t_fst)
                push!(J0_V, 0.0)
                #Reactive PF w/r VoltageMag
                push!(J0_I, F_ix_f_i)
                push!(J0_J, X_ix_t_fst)
                push!(J0_V, 0.0)
                #Active PF w/r Angle
                push!(J0_I, F_ix_f_r)
                push!(J0_J, X_ix_t_snd)
                push!(J0_V, 0.0)
                #Reactive PF w/r Angle
                push!(J0_I, F_ix_f_i)
                push!(J0_J, X_ix_t_snd)
                push!(J0_V, 0.0)
            end
        end
    end
    J0 = SparseArrays.sparse(J0_I, J0_J, J0_V)

    # Use vectors to cache data for closure
    # These should be read only
    P_GEN_BUS = fill(0.0, N_BUS)
    Q_GEN_BUS = fill(0.0, N_BUS)
    P_LOAD_BUS = fill(0.0, N_BUS)
    Q_LOAD_BUS = fill(0.0, N_BUS)
    P_net = fill(0.0, N_BUS)
    Q_net = fill(0.0, N_BUS)
    Vm = fill(0.0, N_BUS)
    θ = fill(0.0, N_BUS)

    state_variable_count = 1
    sources = get_components(StaticInjection, system, d -> !isa(d, ElectricLoad))

    for (ix, b) in enumerate(buses)
        bus_angle = get_angle(b)::Float64
        Vm[ix] = bus_voltage_magnitude = get_magnitude(b)::Float64
        P_GEN_BUS[ix] = 0.0
        Q_GEN_BUS[ix] = 0.0
        get_ext(b)["neighbors"] = neighbors[ix]
        for gen in sources
            !get_available(gen) && continue
            if gen.bus == b
                P_GEN_BUS[ix] += get_active_power(gen)
                Q_GEN_BUS[ix] += get_reactive_power(gen)
            end
        end

        P_LOAD_BUS[ix], Q_LOAD_BUS[ix] = _get_load_data(system, b)
        if b.bustype == BusTypes.REF
            injection_components = get_components(
                StaticInjection,
                system,
                d -> get_available(d) && get_bus(d) == b,
            )
            isempty(injection_components) && throw(
                IS.ConflictingInputsError(
                    "The slack bus does not have any injection component. Power Flow can not proceed",
                ),
            )
            θ[ix] = get_angle(b)::Float64
            x0[state_variable_count] = P_GEN_BUS[ix]
            x0[state_variable_count + 1] = Q_GEN_BUS[ix]
            state_variable_count += 2
        elseif b.bustype == BusTypes.PV
            x0[state_variable_count] = Q_GEN_BUS[ix]
            x0[state_variable_count + 1] = bus_angle
            state_variable_count += 2
        elseif b.bustype == BusTypes.PQ
            x0[state_variable_count] = bus_voltage_magnitude
            x0[state_variable_count + 1] = bus_angle
            state_variable_count += 2
        else
            throw(ArgumentError("Bustype not recognized"))
        end
    end

    @assert state_variable_count - 1 == N_BUS * 2
    bus_types = get_bustype.(buses)
    function pf!(F::Vector{Float64}, X::Vector{Float64})
        for (ix, b) in enumerate(bus_types)
            if b == BusTypes.REF
                # When bustype == REFERENCE Bus, state variables are Active and Reactive Power Generated
                P_net[ix] = X[2 * ix - 1] - P_LOAD_BUS[ix]
                Q_net[ix] = X[2 * ix] - Q_LOAD_BUS[ix]
            elseif b == BusTypes.PV
                # When bustype == PV Bus, state variables are Reactive Power Generated and Voltage Angle
                P_net[ix] = P_GEN_BUS[ix] - P_LOAD_BUS[ix]
                Q_net[ix] = X[2 * ix - 1] - Q_LOAD_BUS[ix]
                θ[ix] = X[2 * ix]
            elseif b == BusTypes.PQ
                # When bustype == PQ Bus, state variables are Voltage Magnitude and Voltage Angle
                P_net[ix] = P_GEN_BUS[ix] - P_LOAD_BUS[ix]
                Q_net[ix] = Q_GEN_BUS[ix] - Q_LOAD_BUS[ix]
                Vm[ix] = X[2 * ix - 1]
                θ[ix] = X[2 * ix]
            end
        end

        # F is active and reactive power balance equations at all buses
        for ix_f in a
            S_re = -P_net[ix_f]
            S_im = -Q_net[ix_f]
            for ix_t in neighbors[ix_f]
                gb = real(Yb[ix_f, ix_t])
                bb = imag(Yb[ix_f, ix_t])
                if ix_f == ix_t
                    S_re += Vm[ix_f] * Vm[ix_t] * gb
                    S_im += -Vm[ix_f] * Vm[ix_t] * bb
                else
                    S_re +=
                        Vm[ix_f] *
                        Vm[ix_t] *
                        (gb * cos(θ[ix_f] - θ[ix_t]) + bb * sin(θ[ix_f] - θ[ix_t]))
                    S_im +=
                        Vm[ix_f] *
                        Vm[ix_t] *
                        (gb * sin(θ[ix_f] - θ[ix_t]) - bb * cos(θ[ix_f] - θ[ix_t]))
                end
            end
            F[2 * ix_f - 1] = S_re
            F[2 * ix_f] = S_im
        end
    end

    function jsp!(J::SparseArrays.SparseMatrixCSC{Float64, Int}, X::Vector{Float64})
        for ix_f in a
            F_ix_f_r = 2 * ix_f - 1
            F_ix_f_i = 2 * ix_f

            for ix_t in neighbors[ix_f]
                X_ix_t_fst = 2 * ix_t - 1
                X_ix_t_snd = 2 * ix_t
                b = bus_types[ix_t]

                if b == BusTypes.REF
                    # State variables are Active and Reactive Power Generated
                    # F[2*i-1] := p[i] = p_flow[i] + p_load[i] - x[2*i-1]
                    # F[2*i] := q[i] = q_flow[i] + q_load[i] - x[2*i]
                    # x does not appears in p_flow and q_flow
                    if ix_f == ix_t
                        J[F_ix_f_r, X_ix_t_fst] = -1.0
                        J[F_ix_f_i, X_ix_t_snd] = -1.0
                    end
                elseif b == BusTypes.PV
                    # State variables are Reactive Power Generated and Voltage Angle
                    # F[2*i-1] := p[i] = p_flow[i] + p_load[i] - p_gen[i]
                    # F[2*i] := q[i] = q_flow[i] + q_load[i] - x[2*i]
                    # x[2*i] (associated with q_gen) does not appear in q_flow
                    # x[2*i] (associated with q_gen) does not appear in the active power balance
                    if ix_f == ix_t
                        #Jac: Reactive PF against local active power
                        J[F_ix_f_i, X_ix_t_fst] = -1.0
                        #Jac: Active PF against same Angle: θ[ix_f] =  θ[ix_t]
                        J[F_ix_f_r, X_ix_t_snd] =
                            Vm[ix_f] * sum(
                                Vm[k] * (
                                    real(Yb[ix_f, k]) * -sin(θ[ix_f] - θ[k]) +
                                    imag(Yb[ix_f, k]) * cos(θ[ix_f] - θ[k])
                                ) for k in neighbors[ix_f] if k != ix_f
                            )
                        #Jac: Reactive PF against same Angle: θ[ix_f] = θ[ix_t]
                        J[F_ix_f_i, X_ix_t_snd] =
                            Vm[ix_f] * sum(
                                Vm[k] * (
                                    real(Yb[ix_f, k]) * cos(θ[ix_f] - θ[k]) -
                                    imag(Yb[ix_f, k]) * -sin(θ[ix_f] - θ[k])
                                ) for k in neighbors[ix_f] if k != ix_f
                            )
                    else
                        g_ij = real(Yb[ix_f, ix_t])
                        b_ij = imag(Yb[ix_f, ix_t])
                        #Jac: Active PF against other angles θ[ix_t]
                        J[F_ix_f_r, X_ix_t_snd] =
                            Vm[ix_f] *
                            Vm[ix_t] *
                            (g_ij * sin(θ[ix_f] - θ[ix_t]) + b_ij * -cos(θ[ix_f] - θ[ix_t]))
                        #Jac: Reactive PF against other angles θ[ix_t]
                        J[F_ix_f_i, X_ix_t_snd] =
                            Vm[ix_f] *
                            Vm[ix_t] *
                            (g_ij * -cos(θ[ix_f] - θ[ix_t]) - b_ij * sin(θ[ix_f] - θ[ix_t]))
                    end
                elseif b == BusTypes.PQ
                    # State variables are Voltage Magnitude and Voltage Angle
                    # Everything appears in everything
                    if ix_f == ix_t
                        #Jac: Active PF against same voltage magnitude Vm[ix_f]
                        J[F_ix_f_r, X_ix_t_fst] =
                            2 * real(Yb[ix_f, ix_t]) * Vm[ix_f] + sum(
                                Vm[k] * (
                                    real(Yb[ix_f, k]) * cos(θ[ix_f] - θ[k]) +
                                    imag(Yb[ix_f, k]) * sin(θ[ix_f] - θ[k])
                                ) for k in neighbors[ix_f] if k != ix_f
                            )
                        #Jac: Active PF against same angle θ[ix_f]
                        J[F_ix_f_r, X_ix_t_snd] =
                            Vm[ix_f] * sum(
                                Vm[k] * (
                                    real(Yb[ix_f, k]) * -sin(θ[ix_f] - θ[k]) +
                                    imag(Yb[ix_f, k]) * cos(θ[ix_f] - θ[k])
                                ) for k in neighbors[ix_f] if k != ix_f
                            )

                        #Jac: Reactive PF against same voltage magnitude Vm[ix_f]
                        J[F_ix_f_i, X_ix_t_fst] =
                            -2 * imag(Yb[ix_f, ix_t]) * Vm[ix_f] + sum(
                                Vm[k] * (
                                    real(Yb[ix_f, k]) * sin(θ[ix_f] - θ[k]) -
                                    imag(Yb[ix_f, k]) * cos(θ[ix_f] - θ[k])
                                ) for k in neighbors[ix_f] if k != ix_f
                            )
                        #Jac: Reactive PF against same angle θ[ix_f]
                        J[F_ix_f_i, X_ix_t_snd] =
                            Vm[ix_f] * sum(
                                Vm[k] * (
                                    real(Yb[ix_f, k]) * cos(θ[ix_f] - θ[k]) -
                                    imag(Yb[ix_f, k]) * -sin(θ[ix_f] - θ[k])
                                ) for k in neighbors[ix_f] if k != ix_f
                            )
                    else
                        g_ij = real(Yb[ix_f, ix_t])
                        b_ij = imag(Yb[ix_f, ix_t])
                        #Jac: Active PF w/r to different voltage magnitude Vm[ix_t]
                        J[F_ix_f_r, X_ix_t_fst] =
                            Vm[ix_f] *
                            (g_ij * cos(θ[ix_f] - θ[ix_t]) + b_ij * sin(θ[ix_f] - θ[ix_t]))
                        #Jac: Active PF w/r to different angle θ[ix_t]
                        J[F_ix_f_r, X_ix_t_snd] =
                            Vm[ix_f] *
                            Vm[ix_t] *
                            (g_ij * sin(θ[ix_f] - θ[ix_t]) + b_ij * -cos(θ[ix_f] - θ[ix_t]))

                        #Jac: Reactive PF w/r to different voltage magnitude Vm[ix_t]
                        J[F_ix_f_i, X_ix_t_fst] =
                            Vm[ix_f] *
                            (g_ij * sin(θ[ix_f] - θ[ix_t]) - b_ij * cos(θ[ix_f] - θ[ix_t]))
                        #Jac: Reactive PF w/r to different angle θ[ix_t]
                        J[F_ix_f_i, X_ix_t_snd] =
                            Vm[ix_f] *
                            Vm[ix_t] *
                            (g_ij * -cos(θ[ix_f] - θ[ix_t]) - b_ij * sin(θ[ix_f] - θ[ix_t]))
                    end
                else
                    error("Undefined Conditional")
                end
            end
        end
    end

    res = similar(x0)
    pf!(res, x0)
    if sum(res) > 10 * (N_BUS * 2)
        _, ix = findmax(res)
        bus_no = get_number(buses[ix])
        @warn "Initial guess provided results in a large initial residual. Largest residual at bus $bus_no"
    end

    if finite_diff
        res = NLsolve.nlsolve(pf!, x0; kwargs...)
    else
        F0 = similar(x0)
        df = NLsolve.OnceDifferentiable(pf!, jsp!, x0, F0, J0)
        res = NLsolve.nlsolve(df, x0; kwargs...)
    end

    return res
end


text:
using Test
using Distributions
using StatsBase
using Random

using EnsembleKalmanProcesses.ParameterDistributionStorage

@testset "ParameterDistributionStorage" begin
    @testset "ParameterDistributionType" begin
        # Tests for the ParameterDistributionType
        d = Parameterized(Gamma(2.0, 0.8))
        @test d.distribution == Gamma(2.0, 0.8)

        d = Samples([1 2 3; 4 5 6])
        @test d.distribution_samples == [1 2 3; 4 5 6]
        d = Samples([1 4; 2 5; 3 6]; params_are_columns = false)
        @test d.distribution_samples == [1 2 3; 4 5 6]
        d = Samples([1, 2, 3, 4, 5, 6])
        @test d.distribution_samples == [1 2 3 4 5 6]
        d = Samples([1, 2, 3, 4, 5, 6]; params_are_columns = false)
        @test d.distribution_samples == reshape([1, 2, 3, 4, 5, 6], :, 1)

    end
    @testset "ConstraintType" begin
        # Tests for the ConstraintType
        # The predefined transforms
        c1 = bounded_below(0.2)
        @test isapprox(c1.constrained_to_unconstrained(1.0) - (log(1.0 - 0.2)), 0)
        @test isapprox(c1.unconstrained_to_constrained(0.0) - (exp(0.0) + 0.2), 0)

        c2 = bounded_above(0.2)
        @test isapprox(c2.constrained_to_unconstrained(-1.0) - (log(0.2 - -1.0)), 0)
        @test isapprox(c2.unconstrained_to_constrained(10.0) - (0.2 - exp(10.0)), 0)

        c3 = bounded(-0.1, 0.2)
        @test isapprox(c3.constrained_to_unconstrained(0.0) - (log((0.0 - -0.1) / (0.2 - 0.0))), 0)
        @test isapprox(c3.unconstrained_to_constrained(1.0) - ((0.2 * exp(1.0) + -0.1) / (exp(1.0) + 1)), 0)
        @test_throws DomainError bounded(0.2, -0.1)

        #an example with user defined invertible transforms
        c_to_u = (x -> 3 * x + 14)
        u_to_c = (x -> (x - 14) / 3)

        c4 = Constraint(c_to_u, u_to_c)
        @test isapprox(c4.constrained_to_unconstrained(5.0) - c_to_u(5.0), 0)
        @test isapprox(c4.unconstrained_to_constrained(5.0) - u_to_c(5.0), 0)

    end

    @testset "ParameterDistribution" begin

        d = Parameterized(Normal(0, 1))
        c = no_constraint()
        name = "unconstrained_normal"
        u = ParameterDistribution(d, c, name)
        @test u.distributions == [d]
        @test u.constraints == [c]
        @test u.names == [name]

        # Tests for the ParameterDistribution
        d = Parameterized(MvNormal(4, 0.1))
        c = [no_constraint(), bounded_below(-1.0), bounded_above(0.4), bounded(-0.1, 0.2)]

        name = "constrained_mvnormal"
        u = ParameterDistribution(d, c, name)
        @test u.distributions == [d]
        @test u.constraints == c
        @test u.names == [name]
        @test_throws DimensionMismatch ParameterDistribution(d, c[1:3], name)
        @test_throws DimensionMismatch ParameterDistribution(d, c, [name, "extra_name"])

        # Tests for the ParameterDistribution
        d1 = Parameterized(MvNormal(4, 0.1))
        c1 = [no_constraint(), bounded_below(-1.0), bounded_above(0.4), bounded(-0.1, 0.2)]
        name1 = "constrained_mvnormal"
        u1 = ParameterDistribution(d1, c1, name1)

        d2 = Samples([1.0 3.0 5.0 7.0; 9.0 11.0 13.0 15.0])
        c2 = [bounded(10, 15), no_constraint()]
        name2 = "constrained_sampled"
        u2 = ParameterDistribution(d2, c2, name2)

        u = ParameterDistribution([d1, d2], [c1, c2], [name1, name2])
        @test u.distributions == [d1, d2]
        @test u.constraints == cat([c1, c2]..., dims = 1)
        @test u.names == [name1, name2]
    end

    @testset "getter functions" begin
        # setup for the tests:
        d1 = Parameterized(MvNormal(4, 0.1))
        c1 = [no_constraint(), bounded_below(-1.0), bounded_above(0.4), bounded(-0.1, 0.2)]
        name1 = "constrained_mvnormal"
        u1 = ParameterDistribution(d1, c1, name1)

        d2 = Samples([1 2 3 4])
        c2 = [bounded(10, 15)]
        name2 = "constrained_sampled"
        u2 = ParameterDistribution(d2, c2, name2)

        u = ParameterDistribution([d1, d2], [c1, c2], [name1, name2])

        # Test for get_dimension(s)
        @test get_total_dimension(u1) == 4
        @test get_total_dimension(u2) == 1
        @test get_total_dimension(u) == 5
        @test get_dimensions(u1) == [4]
        @test get_dimensions(u2) == [1]
        @test get_dimensions(u) == [4, 1]

        # Tests for get_name
        @test get_name(u1) == [name1]
        @test get_name(u) == [name1, name2]

        # Tests for get_n_samples
        @test typeof(get_n_samples(u)[name1]) <: String
        @test get_n_samples(u)[name2] == 4

        # Tests for get_distribution
        @test get_distribution(d1) == MvNormal(4, 0.1)
        @test get_distribution(u1)[name1] == MvNormal(4, 0.1)
        @test typeof(get_distribution(d2)) == Array{Int64, 2}
        @test typeof(get_distribution(u2)[name2]) == Array{Int64, 2}

        d = get_distribution(u)
        @test d[name1] == MvNormal(4, 0.1)
        @test typeof(d[name2]) == Array{Int64, 2}

        # Test for get_all_constraints
        @test get_all_constraints(u) == cat([c1, c2]..., dims = 1)
    end

    @testset "statistics functions" begin

        # setup for the tests:
        d1 = Parameterized(MvNormal(4, 0.1))
        c1 = [no_constraint(), bounded_below(-1.0), bounded_above(0.4), bounded(-0.1, 0.2)]
        name1 = "constrained_mvnormal"
        u1 = ParameterDistribution(d1, c1, name1)

        d2 = Samples([1 2 3 4])
        c2 = [bounded(10, 15)]
        name2 = "constrained_sampled"
        u2 = ParameterDistribution(d2, c2, name2)

        d3 = Parameterized(Beta(2, 2))
        c3 = [no_constraint()]
        name3 = "unconstrained_beta"
        u3 = ParameterDistribution(d3, c3, name3)

        u = ParameterDistribution([d1, d2], [c1, c2], [name1, name2])

        d4 = Samples([1 2 3 4 5 6 7 8; 8 7 6 5 4 3 2 1])
        c4 = [no_constraint(), no_constraint()]
        name4 = "constrained_MVsampled"
        v = ParameterDistribution([d1, d2, d3, d4], [c1, c2, c3, c4], [name1, name2, name3, name4])

        # Tests for sample distribution
        seed = 2020
        Random.seed!(seed)
        s1 = rand(MvNormal(4, 0.1), 1)
        Random.seed!(seed)
        @test sample_distribution(u1) == s1

        Random.seed!(seed)
        s1 = rand(MvNormal(4, 0.1), 3)
        Random.seed!(seed)
        @test sample_distribution(u1, 3) == s1

        Random.seed!(seed)
        idx = StatsBase.sample(collect(1:size(d2.distribution_samples)[2]), 1)
        s2 = d2.distribution_samples[:, idx]
        Random.seed!(seed)
        @test sample_distribution(u2) == s2

        Random.seed!(seed)
        idx = StatsBase.sample(collect(1:size(d2.distribution_samples)[2]), 3)
        s2 = d2.distribution_samples[:, idx]
        Random.seed!(seed)
        @test sample_distribution(u2, 3) == s2

        Random.seed!(seed)
        s1 = sample_distribution(u1, 3)
        s2 = sample_distribution(u2, 3)
        Random.seed!(seed)
        s = sample_distribution(u, 3)
        @test s == cat([s1, s2]..., dims = 1)

        #Test for get_logpdf
        @test_throws ErrorException get_logpdf(u, zeros(get_total_dimension(u)))
        x_in_bd = [0.5]
        Random.seed!(seed)
        lpdf3 = logpdf.(Beta(2, 2), x_in_bd)[1] #throws deprecated warning without "."
        Random.seed!(seed)
        @test isapprox(get_logpdf(u3, x_in_bd) - lpdf3, 0.0; atol = 1e-6)
        @test_throws DimensionMismatch get_logpdf(u3, [0.5, 0.5])

        #Test for get_cov, get_var
        block_cov = cat([get_cov(d1), get_var(d2), get_var(d3), get_cov(d4)]..., dims = (1, 2))
        @test isapprox(get_cov(v) - block_cov, zeros(get_total_dimension(v), get_total_dimension(v)); atol = 1e-6)
        #Test for get_mean
        means = cat([get_mean(d1), get_mean(d2), get_mean(d3), get_mean(d4)]..., dims = 1)
        @test isapprox(get_mean(v) - means, zeros(get_total_dimension(v)); atol = 1e-6)

    end

    @testset "transform functions" begin
        #setup for the tests
        tol = 1e-8
        d1 = Parameterized(MvNormal(4, 0.1))
        c1 = [no_constraint(), bounded_below(-1.0), bounded_above(0.4), bounded(-0.1, 0.2)]
        name1 = "constrained_mvnormal"
        u1 = ParameterDistribution(d1, c1, name1)

        d2 = Samples([1.0 3.0 5.0 7.0; 9.0 11.0 13.0 15.0])
        c2 = [bounded(10, 15), no_constraint()]
        name2 = "constrained_sampled"
        u2 = ParameterDistribution(d2, c2, name2)

        u = ParameterDistribution([d1, d2], [c1, c2], [name1, name2])

        x_unbd = rand(MvNormal(6, 3), 1000)  #6 x 1000
        # Tests for transforms
        x_real_constrained1 = mapslices(x -> transform_unconstrained_to_constrained(u1, x), x_unbd[1:4, :]; dims = 1)
        @test isapprox(
            x_unbd[1:4, :] -
            mapslices(x -> transform_constrained_to_unconstrained(u1, x), x_real_constrained1; dims = 1),
            zeros(size(x_unbd[1:4, :]));
            atol = tol,
        )

        x_real_constrained2 = mapslices(x -> transform_unconstrained_to_constrained(u2, x), x_unbd[5:6, :]; dims = 1)
        @test isapprox(
            x_unbd[5:6, :] -
            mapslices(x -> transform_constrained_to_unconstrained(u2, x), x_real_constrained2; dims = 1),
            zeros(size(x_unbd[5:6, :]));
            atol = tol,
        )

        x_real = mapslices(x -> transform_unconstrained_to_constrained(u, x), x_unbd; dims = 1)
        x_unbd_tmp = mapslices(x -> transform_constrained_to_unconstrained(u, x), x_real; dims = 1)
        @test isapprox(x_unbd - x_unbd_tmp, zeros(size(x_unbd)); atol = tol)

        # Tests transforms for other input structures
        @test isapprox(transform_unconstrained_to_constrained(u1, x_unbd[1:4, :]), x_real_constrained1; atol = tol)
        @test isapprox(
            x_unbd[1:4, :] - transform_constrained_to_unconstrained(u1, x_real_constrained1),
            zeros(size(x_unbd[1:4, :]));
            atol = tol,
        )
        @test isapprox(transform_unconstrained_to_constrained(u2, x_unbd[5:6, :]), x_real_constrained2; atol = tol)
        @test isapprox(
            x_unbd[5:6, :] - transform_constrained_to_unconstrained(u2, x_real_constrained2),
            zeros(size(x_unbd[5:6, :]));
            atol = tol,
        )
        @test isapprox(
            transform_unconstrained_to_constrained(u2, [x_unbd[5:6, :], x_unbd[5:6, :]]),
            [x_real_constrained2, x_real_constrained2];
            atol = tol,
        )
    end
end


text:
export iauNum00a
"""
Form the matrix of nutation for a given date, IAU 2000A model.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   date1,date2  double          TT as a 2-part Julian Date (Note 1)

Returned:
   rmatn        double[3][3]    nutation matrix

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. The matrix operates in the sense V(true) = rmatn * V(mean), where
   the p-vector V(true) is with respect to the true equatorial triad
   of date and the p-vector V(mean) is with respect to the mean
   equatorial triad of date.

   3. A faster, but slightly less accurate result (about 1 mas), can be
   obtained by using instead the iauNum00b function.

Called:
   iauPn00a     bias/precession/nutation, IAU 2000A

Reference:

   Explanatory Supplement to the Astronomical Almanac,
   P. Kenneth Seidelmann (ed), University Science Books (1992),
   Section 3.222-3 (p114).

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauNum00a(double date1, double date2, double rmatn[3][3])
function iauNum00a(date1::Real, date2::Real)
   # Allocate return values
   rmatn = zeros(Float64, 3, 3)

   ccall((:iauNum00a, libsofa_c), Cvoid,
        (Cdouble, Cdouble, Ptr{Cdouble},),
        convert(Float64, date1), convert(Float64, date2),
        rmatn)

   return SMatrix{3,3}(rmatn')
end